RADIATION RESEARCH 170, 244–259 (2008)

0033-7587/08 $15.00
䉷 2008 by Radiation Research Society.
All rights of reproduction in any form reserved.

Development of PARMA: PHITS-based Analytical Radiation Model in
the Atmosphere
Tatsuhiko Sato,a,1 Hiroshi Yasuda,b Koji Niita,c Akira Endoa and Lembit SihverdaJapan Atomic Energy Agency (JAEA); b National Institute of Radiological Sciences (NIRS); c Research Organization for Information Science andTechnology (RIST); and d Chalmers University of Technology

Several calculation codes, e.g. EPCARD (2), CARI-6 (3)and PCAIRE (4), have been developed to estimate the aircrew dose. They can calculate route doses, the dose enroute between two airports, by specifying the flight conditions, such as the departure and destination airports. The
accuracy of the calculations was well verified by experimental data. However, the calculated dose is intrinsicallyderived from a non-physical quantity: the fluence-to-dose
conversion coefficient, which depends significantly on its
calculation model properties such as the dose type (the effective dose or the ambient dose equivalent), the radiation
weighting factor, the characterization of the human model,
and the irradiation geometry. Hence calculation results reflect the radiation protection policy adopted in the code. It
is therefore worthwhile to develop a new route-dose calculation code that explicitly determines the atmosphericcosmic-ray spectra on flight routes and combines these datawith user-specified fluence-to-dose conversion coefficients.
Establishment of a reliable model for calculating the cosmic-ray spectra under any global conditions is the key issue
in the development. Estimation of the spectra is also important for research in astrophysics and elementary particle
physics.
A number of studies have been carried out to build themodel. The most of the earlier work was devoted to theconstruction of semi-empirical or analytical models. For instance, O’Brien et al. developed a deterministic code LUIN
(5) based on an analytical two-component solution of the
Boltzmann transport equation, which is employed in theroute-dose calculation code CARI-6 (3). However, mostthis research has shifted to the development of Monte Carlocode that can be used in the simulation of atmosphericpropagation of cosmic rays, owing to the dramatic improvement of computer performance in the last decade. Several
simulation codes such as CORSIKA (6), COSMOS (7),PLANETOCOSMICS (8) and FLUKA (9, 10) were developed or used for this purpose. Some of their simulation
results, e.g. ref. (11), proved their excellent ability to reproduce measured neutron spectra at flight altitudes, which
are the most important quantity to be reproduced in routedose calculation. However, those Monte Carlo codes are
rarely incorporated directly into a route-dose calculation

Sato, T., Yasuda, H., Niita, K., Endo, A. and Sihver, L.
Development of PARMA: PHITS-based Analytical Radiation
Model in the Atmosphere. Radiat. Res. 170, 244–259 (2008).
Estimation of cosmic-ray spectra in the atmosphere has
been essential for the evaluation of aviation doses. We therefore calculated these spectra by performing Monte Carlo simulation of cosmic-ray propagation in the atmosphere using the
PHITS code. The accuracy of the simulation was well verified
by experimental data taken under various conditions, evennear sea level. Based on a comprehensive analysis of the simulation results, we proposed an analytical model for estimating the cosmic-ray spectra of neutrons, protons, helium ions,
muons, electrons, positrons and photons applicable to any location in the atmosphere at altitudes below 20 km. Our model,
named PARMA, enables us to calculate the cosmic radiation
doses rapidly with a precision equivalent to that of the Monte
Carlo simulation, which requires much more computational
time. With these properties, PARMA is capable of improving
the accuracy and efficiency of the cosmic-ray exposure dose
estimations not only for aircrews but also for the public on
the ground. 䉷 2008 by Radiation Research Society

INTRODUCTION

At high altitude, aircrews are exposed to a high level ofcosmic radiations. Protection for crew members against
these terrestrial cosmic rays has been widely discussed
since the publication of ICRP publication 60 (1), in which
this aircrew exposure is recognized as an occupational hazard. As a result of this discussion, many countries haveissued regulations (or recommendations) for an annual dose
limitation for aircrews. Development of an aviation dose
calculation code is indispensable for following the regulations, since the doses depend on the altitude, the geomagnetic location and the solar activity (referred to here as
global conditions) along the flight routes in a complicatedmanner, and it is impractical to measure the doses for all
flight conditions.
Address for correspondence: Research Group for Radiation Protection, Division of Environment and Radiation Sciences, Nuclear Scienceand Engineering Directorate, Japan Atomic Energy Agency, Tokai, Naka,
Ibaraki, 319-1195, Japan; e-mail: sato.tatsuhiko@jaea.go.jp.
1244

DEVELOPMENT OF PARMA

code, since it is extremely time-consuming to perform
Monte Carlo simulation of the cosmic-ray propagation foreach route-dose calculation even using the latest computer.
For example, it takes approximately half a day to calculate
terrestrial cosmic-ray spectra at a certain location by Monte
Carlo simulation using our parallel computer with 24 CPUs,
and route-dose calculation directly based on the Monte Carlo simulation is expected to be much more time-consuming
due to the variety of operational flight conditions. Assumption or parameterization is thus required to allow the Monte
Carlo-obtained spectra to be used in route-dose calculation.
With these problems in mind, we calculated the terrestrial
cosmic-ray spectra by performing the Monte Carlo simulation of cosmic-ray propagation in the atmosphere by theParticle and Heavy Ion Transport code System PHITS (12)coupled with the nuclear data library JENDL-High-Energy
File (JENDL/HE) (13, 14). Based on a comprehensive analysis of the simulation results, we proposed an analyticalmodel for estimating the atmospheric cosmic-ray spectra
for neutrons, protons, helium ions, muons, electrons, positrons and photons applicable to any global conditions ataltitudes below 20 km. The model was designated PARMA,or PHITS-based Analytical Radiation Model in the Atmosphere. The details of the simulation procedure together
with the calculated cosmic-ray neutron spectra were published in our previous paper (15). This paper therefore focuses on describing the results for other particles, including
the derivation and verification of PARMA for these particles.
MONTE CARLO SIMULATION
Simulation Procedure
The simulation procedure for the atmospheric propagation of cosmic
rays is basically the same as that described in our previous paper (15)except for the source-term determination. In the simulation, the atmosphere was divided into 28 concentric spherical shells, and its maximum
altitude was assumed to be 86 km. The densities and temperatures ofeach shell were determined by referring to the U.S. Standard Atmosphere
1976. Note that argon was replaced by the atom with the same mass
number, calcium, in our simulation, since JENDL/HE did not include thedata for argon. The Earth was represented as a sphere with a radius of
6378.14 km, and its composition was presumed to be the same as thatof the air at sea level to obtain the atmospheric cosmic-ray spectra under
the ideal condition, i.e. without the disturbance of the local geometryeffect. The particles reaching 1000 g/cm2 below the ground level werediscarded in the simulation to reduce the computational time. Note thatthe composition of the soil significantly influences the neutron spectra at
the ground level (16, 17), and we analyzed the dependence of the spectraon the composition by changing the water density in the ground (15).
However, the effects of the composition on the other particle spectra are
expected to be much smaller compared to the neutron case, since there
are few albedo particles other than neutrons.
In the simulation, cosmic rays were incident from the top of the atmosphere assumed in the virtual Earth system, i.e. from the altitude of
86 km. The galactic cosmic-ray (GCR) protons and heavy ions with energies and charges up to 200 GeV/nucleon and 28 (nickel), respectively,
were considered as the source particles. The GCR spectra around theEarth can be estimated from their local interstellar (LIS) spectra, considering the modulation due to the solar wind magnetic field, so-called solar

245

FIG. 1. Calculated GCR proton and helium-ion spectra above theEarth’s atmosphere in comparison with the experimental data obtained by
the BESS spectrometer during the last solar period. The values in theparentheses indicate the force field potential calculated from several neutron monitor count rates at each experimental date.
modulation. In the determination of the source particle spectra in our
simulation, we employed the LIS spectra calculated by the Nymmik model (18) coupled with modified parameters. The solar modulation was considered based on the force field formalism (19, 20), using the force fieldpotential that is occasionally called the heliocentric potential (5, 21).
Figure 1 shows the calculated GCR proton and helium-ion spectra
above the Earth’s atmosphere compared to the corresponding experimental data obtained by the BESS spectrometer (22) during the last solar
period. The numerical values of the force field potential at each experimental date were calculated from the count rates of several neutron monitors located all over the world (23); their relationship will be presented
in a future paper. Note that our calculation procedure for estimating thenumerical value of the force field potential is different from that for the
heliocentric potential (5), although the results are close to each other. We
therefore adopted the name ‘‘force field potential’’ instead of ‘‘heliocentric potential’’ in this paper to prevent confusion between the two quantities. It is found from the graph that the lower energy fluxes decrease
with an increase in the modulation potential, and the calculation can
reproduce the experimental tendency very well. We therefore concluded
that the calculated GCR spectra are precise enough to be used in making
the source determination in the atmospheric propagation simulation ofcosmic rays.
The Monte Carlo simulations were carried out for five force field potentials—400, 600, 900, 1200 and 1800 MV—and 18 geomagnetic fields
with the vertical cut-off rigidities from 0 to 17 GV. The azimuth andzenith dependences of the cut-off rigidity were considered by assumingthat the geomagnetism can be simply expressed by a dipole magnet, as
described in ref. (17).
The atmospheric propagation of the incident cosmic rays and theirassociated cascades was simulated by the PHITS code, which can deal
with the transport of all kinds of hadrons and heavy ions with energies
up to 200 GeV/nucleon. PHITS can also treat the production, transportand decay of photons, electrons, positrons, pions, muons, kaons and various resonant states. In the simulation, two models, JENDL/HE and INC,
the widely used model of the intranuclear cascade (24), were alternatively
employed for simulating nuclear reactions induced by neutrons and protons below 3 GeV. An advantage of JENDL/HE compared to INC is thatit can precisely calculate the yields of high-energy secondary particles
knocked out from light ions such as nitrogen and oxygen, which are thedominant components of the atmosphere. Owing to this property, thesimulation using JENDL/HE can reproduce the experimental data of cosmic-ray neutron spectra very well even near sea level, as shown in our246

SATO ET AL.

FIG. 2. Calculated cosmic-ray spectra at a typical flight altitude (A) and near sea level (B). The values of d andrc are the atmospheric depth and the cut-off rigidity, respectively, while smin and smax indicate the solar minimum andmaximum conditions with the force field potential 400 and 1200 MV, respectively. These notations are also used in
the other figures.

previous paper (15). However, the current version of JENDL/HE written
in PHITS-readable format did not include the pion-production channels,
and consequently the spectra of pions and the particles associated with
their decay—muons, electrons, positrons and photons—could not be determined precisely by the simulation. We therefore decided to adopt theresults obtained by the simulation employing JENDL/HE for neutron,
proton and helium-ion spectra and the simulation employing INC for
muon, electron, positron and photon spectra in the analysis throughout
this paper.
Results of the Simulation
As examples of the simulation results, the cosmic-ray spectra at a typical flight altitude and near sea level are shown in Fig. 2. The statisticaluncertainties in the values obtained by the simulation are generally small,less than approximately 5% and 20% for the flight altitude and sea-level
data, respectively, except for helium-ion spectra. The corresponding spectra calculated by the analytical model proposed in the next section,
PARMA, are also plotted in the figure. The discussions about the comparison between the results obtained by the Monte Carlo simulation and
PARMA are given in the next section.
In general, the fluxes of the terrestrial cosmic rays at flight altitude areapproximately 100 times larger than those near sea level except for thoseof muons. This tendency is consistent with the well-known fact that dose
rates at conventional flight altitudes are approximately 100 times higher
than those at sea level. The muon flux does not increase with rising
altitude very much in comparison to other particles because of theirstrong penetrability in the atmosphere.
The energy ranges of the electron and positron spectra obtained by thesimulation have an upper limit of 1 GeV, since the current version of
PHITS cannot handle the transport of electrons and positrons above thatenergy. When such particles are created in the PHITS calculation, theirconsequent transports are simulated under the assumption that the energyof the produced particle is equal to 1 GeV, and its importance is weighted
by the ratio of the real energy to 1 GeV.
Comparison with Experimental Data
To examine the usefulness of the calculated cosmic-ray spectra in dose
estimation, it is ideal to compare the angular-integrated spectra obtainedby Monte Carlo simulation with the corresponding experimental data, as
shown in our previous paper dealing with neutrons. In the cases of otherparticles, however, there are few experimental data that can be used forthis purpose, since most of these measurements have been performed for

elucidating angular-differential spectra, this information being of prime
importance in astrophysics and elementary particle physics. We therefore
compared simulated and experimental angular-differential spectra to verify the accuracy of the Monte Carlo simulation.
The results of the comparisons are shown in Fig. 3. Panel A shows
the angular-integrated neutron spectra in the unit of cm⫺2 s⫺1 lethargy⫺1,
which had been reported in our previous paper (15) in a different unit,
in comparison with the measurement (25). Panels B to F, respectively,
show the Monte Carlo-obtained proton, muon, photon and electron spectra for the vertical direction in the unit of cm⫺2 s⫺1 sr⫺1 MeV⫺1, plotted
together with the corresponding measured data (26–34) under a variety
of experimental conditions. Excellent agreement between the calculated
and measured spectra was observed except for the electron data, indicating the adequacy of the assumptions adopted in the Monte Carlo simulation. We therefore concluded that the cosmic-ray spectra estimated byour simulation are precise enough to be subjected to systematic analysis
for developing an analytical model of these spectra. The discrepancy inthe electron spectra is presumed to be attributable to the assumption madein the treatment of high-energy electron transport adopted in PHITS, as
described before. The slight disagreement observed in the photon spectra
at sea level is probably due to the effect of photons emitted from radioisotopes, which was not considered in our simulation.

ANALYTICAL MODEL: PARMA

General Description of PARMA
The analytical model proposed in this section enables us
to estimate the cosmic-ray proton, helium-ion, electron,
positron, photon and muon spectra with energies from 1
MeV to 200 GeV. The unit of the obtained spectra is cm⫺2s⫺1 MeV⫺1, supplying the force field potential in MV, thevertical cut-off rigidity in GV, and the kinetic energy in
MeV except for the case of helium ions, where the kinetic
energy is given in MeV/nucleon to estimate its spectra incm⫺2 s⫺1 (MeV/nucleon)⫺1. Although the full two-dimensional distribution of the cut-off rigidity was considered in
the PHITS simulation, we simply adopted the vertical cutoff rigidity as an index for accessing the simulation results.
The associating errors are expected to be small except forthe magnetic equator region, as discussed in ref. (5). This

DEVELOPMENT OF PARMA

247

FIG. 3. Comparison between the calculated and experimental cosmic-ray spectra in the atmosphere. Panel A
shows the angular-integrated spectra in the unit of cm⫺2 s⫺1 lethargy⫺1, while panels B–F show the spectra for thevertical direction in the unit of cm⫺2 s⫺1 sr⫺1 MeV⫺1.

model coupled with that for neutrons proposed in our previous paper (15) was given the name PARMA.
In the development of PARMA, the Monte Carlo-obtained spectra for the force field potentials 400 and 1200
MV were regarded as the data for the solar minimum andmaximum conditions, respectively, although the highest
force field potential adopted in our Monte Carlo simulation
was 1800 MV. The data for the other force field potentials
were used only for the determination of the solar-modulation dependence of the secondary particle fluxes, as dis-cussed later in this section. The Monte Carlo-obtained spectra at the altitudes above 20 km (⬃59 g/cm2) were notconsidered in the derivation of PARMA for the following
two reasons: (1) the equilibrium between the numbers ofincoming and outgoing particles, which is a necessary condition for calculating lower-energy particle fluxes by ourmodel, is not established at the higher altitudes, and (2)commercial flights never exceed an altitude of 20 km. Thus
the applicable altitude range of PARMA is limited to 20
km.

248

SATO ET AL.

FIG. 4. Comparison between the cosmic-ray proton and helium-ion spectra obtained by our Monte Carlo simulation and PARMA. The left and right panels show the atmospheric depth and the cut-off rigidity dependence of thespectra, respectively.

Consideration of Proton and Helium-Ion Spectra
Figure 4 shows the cosmic-ray proton and helium-ion
spectra obtained by the Monte Carlo simulation, together
with the corresponding data calculated by PARMA. The
left and right panels show the atmospheric depth and the
TABLE 1
Numerical Values of the Parameters a used in Eqs.
(1) to (3) and (11) to (13) for Estimating the
Atmospheric Proton and Helium-Ion Spectra
Parameter⫺1a1 (cm g MeV/nucleon)a2
a3 (cm2 g⫺1)
a4 (cm2 g⫺1)
a5
a6
a7 (s⫺1m⫺2 sr⫺1 GV⫺1)
a8 (s⫺1 m⫺2 sr⫺1 GV⫺1)a9 (MeV/nucleon)a10 (MeV/nucleon)
a11
a12
a13
a14 (cm2 g⫺1 MeV/nucleon)
a15 (MeV/nucleon)
a16 (MeV/nucleon)
2

Protons

Helium ions

2.12
0.445
0.01010.03962.924
2.708
1.27 ⫻ 104
4.83 ⫻ 103
3.28 ⫻ 104
7.44 ⫻ 103
3.46
1.68
1.37
2.07
108
2.30 ⫻ 103

17.6
0.438
0.0121
0.04341.841
2.646
2.36 ⫻ 103
432
6.06 ⫻ 103
2.41 ⫻ 103
3.33
11.7
0.967
3.20
15.0
853

vertical cut-off rigidity dependences of the spectra, respectively. It is found from the data in the figure that the spectra
can be divided into higher- and lower-energy components,
although they are not clearly distinguished in the case of
protons. The two components consist predominantly of the
primary cosmic rays and their secondary particles producedin the atmosphere, respectively. The switching energy between the two components depends on the cut-off rigidity,
as discussed later in this section. We first established themathematical functions to estimate the primary and secondary particle spectra separately and then constructed analytical models for predicting the whole spectrum by combining them.
Primary Particle Spectra for Protons and Helium Ions
Considering the energy loss and nuclear interactions in
the atmosphere, the primary proton and helium-ion spectra␾pri can be assumed to be expressed by␾ pri (s, d, E)⫽ ␾ TOA (s, E ⫹ a1 d)⫻ [a2 exp(⫺a3 d) ⫹ (1 ⫺ a2 )exp(⫺a4 d)].

(1)where s, d and E denote the force field potential, the atmospheric depth and the kinetic energy per nucleon, re-DEVELOPMENT OF PARMA

spectively, and ␾TOA is the spectra at the top of the atmosphere. The effect of the magnetosphere was not consideredin this equation and will be taken into account in the equation for combining the primary spectra with the secondary
one, i.e. in Eq. (11) proposed later in this section.
The parameter a1 is related to the particle stopping powerof the atmosphere, and consequently E ⫹ a1d indicates the
kinetic energy at the top of the atmosphere. The two exponential decay terms represent the decrease in the primary-particle fluxes due to nuclear interactions, and the parameters a2 to a4 are related to primary particle macroscopic
cross sections in the atmosphere. Hence the parameters ai
should depend on the particle energy in a complicated manner, but we assumed that they were constant for each particle and determined their numerical values from the leastsquares fitting of their high-energy spectra (above several
GeV, depending of the cut-off rigidity) as obtained by theMonte Carlo simulation, assuming that there is no secondary particle in such a high-energy region. The results of thefitting are summarized in Table 1.
In the least-squares fitting, ␾TOA was obtained from the
LIS spectra calculated by the Nymmik model (18) coupledwith modified parameters, considering the solar modulation
based on the force field formalism, as described in the previous section. It can be written as
␾ TOA (s, E) ⫽[ ]2

C(ELIS )[␤(ELIS )] a5 R(E),
[R(ELIS )] a6
R(ELIS )(2)

where ELIS is the kinetic energy at LIS, ␤ is the speed ofthe particle relative to light, and R is the rigidity of theparticle in GV, which can be obtained from the equation R⫽ 0.001 ⫻ 兹(AE)2 ⫹ 2AEmE /Z, where A, Z and Em are
the mass and charge number and the rest mass of the particle, respectively. The former part of the right side of Eq.
(2), C␤a5 R⫺a6, expresses the LIS spectra, while the latterdoes the effect of the solar modulation based on the force
field formalism. Note that E and Em should be supplied in
MeV/nucleon and MeV, respectively. Based on the force
field formalism, ELIS can be simply determined by E ⫹ sZ/
A, where s is the force field potential expressed in MV.
In the original Nymmik model, the parameter C is regarded as constant. However, Eq. (2) coupled with a constant C parameter cannot reproduce the proton and heliumion spectra at the top of the atmosphere measured with the
BESS spectrometer during the variety of the solar-modulation conditions over wide energy ranges very much. We
therefore presumed that the parameter C depends on theparticle energy and their relationship can be expressed by
the sigmoid function:
C(E) ⫽ a7 ⫹a8.
1 ⫹ exp[(E ⫺ a9 )/a10 ](3)

where a7 to a10 are constant parameters. The numerical values of a5, a8, a9 and a10 were determined from the leastsquares fitting of the BESS data. In the least-squares fitting,

249

FIG. 5. Normalized proton and helium-ion spectra by their fluxes atthe energy of 1 MeV/nucleon.

the parameters a6 and a7 were fixed to force ␾TOA calculated
by Eq. (2) asymptotic to that obtained from the Honda model (35) in the high-energy region. It should be noted that
this calculation method was introduced only for reproducing the BESS data and is not fully based on a theoreticalmodel. The results of the fitting are also summarized inTable 1. Note that Eqs. (2) and (3) coupled with the parameters listed in Table 1 give the particle spectra in the unitof s⫺1 m⫺2 sr⫺1 GV⫺1. Thus they should be multiplied by(4␲ ⫺ ⍀E) ⫻ (A/Z)/␤ ⫻ 10⫺7, where ⍀E is the solid angleof the Earth from a point at the top of the atmosphere in
our simulation, to convert that unit into the one used in theleast-squares fitting of the terrestrial proton and helium-ion
spectra—cm⫺2 s⫺1 (MeV/nucleon)⫺1. The numerical valueof ⍀E is equal to 1.675␲.
Secondary Particle Spectra for Protons and Helium Ions
To estimate the influence of the global conditions on therelative shapes of the secondary spectra, we normalized thespectra to their fluxes at the energy of 1 MeV/nucleon, at
which the contributions of the primary particles to the fluxes are almost negligible. By analyzing the normalized spectral shape, we found that the proton data depend only on
the atmospheric depth and are almost independent of the
cut-off rigidity and the force field potential and that the datafor helium ions are independent of all global conditions.
Figure 5 shows examples of the normalized spectra. The
proton data shown in the figure were obtained by averaging
over all Monte Carlo results at all the altitudes, while thedata for helium ions were the mean spectra of all global
conditions. It is found from the data in the figure that thenormalized spectra in lower- and higher-energy regions canbe expressed simply by power functions of the particle energy.
Based on these considerations, we introduced the following
function to estimate the secondary particle spectra ␾sec:

250

SATO ET AL.

TABLE 2
Numerical Values of the Coefficients of Third Order of the Polynomial Function—theParameters c in Eq. (5)—Obtained from the Least-Squares Fittings
ParameterProtons

Helium ionsElectrons

Positrons

Photons

b1
b2
b3
b4
b1
b2
b3
b4
b1
b2
b3
b4
b1
b2
b3
b4
b1
b2
b3
b4

ci1

ci2

ci3

ci4

1.26
0.438
1.81 ⫻ 10⫺4
1.71
1.00
0.881
1.80 ⫻ 10⫺4
4.77
6.44
⫺0.8940.002311.13
2.29⫺0.4550.00387
1.28
15.8
⫺1.25
0.0121
0.825

0.00323⫺5.58 ⫻ 10⫺4
⫺5.18 ⫻ 10⫺7
7.16 ⫻ 10⫺44.81 ⫻ 10⫺67.84 ⫻ 10⫺77.59 ⫻ 10⫺10⫺9.32 ⫻ 10⫺7

2.28 ⫻ 10⫺9⫺3.87 ⫻ 10⫺10⫺3.82 ⫻ 10⫺13
5.27 ⫻ 10⫺100.0266⫺3.58 ⫻ 10⫺4
8.07 ⫻ 10⫺6
6.64 ⫻ 10⫺40.0124⫺7.14 ⫻ 10⫺4⫺5.82 ⫻ 10⫺68.82 ⫻ 10⫺40.00963
8.44 ⫻ 10⫺42.47 ⫻ 10⫺50.00134⫺5.97 ⫻ 10⫺52.37 ⫻ 10⫺7⫺2.00 ⫻ 10⫺8
⫺1.04 ⫻ 10⫺6
⫺2.94 ⫻ 10⫺5
1.17 ⫻ 10⫺66.76 ⫻ 10⫺9⫺1.31 ⫻ 10⫺6
⫺2.22 ⫻ 10⫺5
⫺2.14 ⫻ 10⫺6
⫺5.33 ⫻ 10⫺8
⫺2.18 ⫻ 10⫺6

4.39 ⫻ 10⫺8
5.39 ⫻ 10⫺11
1.89 ⫻ 10⫺11
3.13 ⫻ 10⫺10
2.08 ⫻ 10⫺8
⫺5.13 ⫻ 10⫺102.93 ⫻ 10⫺12
4.22 ⫻ 10⫺10
2.53 ⫻ 10⫺8
1.36 ⫻ 10⫺92.79 ⫻ 10⫺11
1.15 ⫻ 10⫺9Note. The parameters b for helium ions are assumed to be independent of the atmospheric depth, i.e. ci2 ⫽ ci3 ⫽ci4 ⫽ 0.

␾ sec (s, rc , d, E) ⫽ ⌽N (s, rc , d)b1 (d)E b2(d),1 ⫹ b3 (d)E b4(d)(4)

where ⌽N is the flux used for the normalization of the spectra, i.e. flux at 1 MeV, rc is the cut-off rigidity, and bi arefree fitting parameters. For b3 K 1, ␾sec is proportional tothe power functions Eb2 and Eb2⫺b4 in lower- and higherenergy regions, respectively. The parameters bi for reproducing the proton spectra are assumed to be dependent ond, while those for the helium-ion spectra are assumed to beconstant. If ␾sec can be fitted perfectly by the equation, thenb1 should be equal to 1. However, b1 was also regarded asa fitting parameter in the case of protons, since the normalized proton spectra in the very low-energy region could
not be expressed by a simple power function. The numerical values of bi were determined from the least-squares
fitting of the normalized spectra. The curves resulting fromthe fitting curves are also plotted in Fig. 5, and they clearly
demonstrate the ability of the equation to reproduce thesimulation results.
For expressing the dependences of bi for protons on theatmospheric depth d, we simply introduced the third-orderpolynomial function:bi (d) ⫽ ci1 ⫹ ci2 d ⫹ ci3 d 2 ⫹ ci4 d 3,(5)

where the parameters cij are the fitting constants. Table 2
summarizes the numerical values of the parameters cij obtained from the least-squares fitting of the bi data together
with those of the constant bi parameters for reproducing thesecondary helium-ion spectra.
The flux used for the normalization ⌽N can be obtainedin the same manner as the low-energy atmospheric neutron
fluxes, whose calculation procedure was described in detailin our previous paper (15). Since equilibrium between thenumbers of incoming and outgoing particles is almost established at low energies, ⌽N for the solar minimum condition, ⌽Nmin, can be expressed by⌽Nmin (rc , d)
⫽ g1min (rc ){exp[⫺g2min (rc )d]⫺ g3min (rc )exp[⫺g4min (rc )d]},(6)

where g1 to g4 are parameters depending on the cut-off rigidity rc. The normalization flux for the solar maximum
condition, ⌽Nmax, can be also estimated from Eq. (6) byreplacing the subscript ‘‘min’’ with ‘‘max’’ in the equation.
For expressing the dependence of gi on rc, we introduced
the sigmoid function coupled with the linear term:gi (rc ) ⫽ hi1 ⫹ hi 2 rc ⫹hi 3,
1 ⫹ exp[(rc ⫺ hi 4 )/hi 5 ](7)where hi1 to hi5 are free parameters. The numerical valuesof the parameters hij for the solar minimum and maximum
conditions were determined from the least-squares fitting of
the Monte Carlo-obtained ⌽N data for the force field potentials 400 and 1200 MV, respectively. The results of thisleast-squares fitting are summarized in Table 3. It should
be noted that the cosmic-ray spectra below the altitude of
20 km are almost independent of the cut-off rigidity for rc⬍ 1 GV, since incident cosmic rays with rigidity below 1
GV together with their secondary particles generally cannot

251

DEVELOPMENT OF PARMATABLE 3
Numerical Values of the Parameters h used in Eq. (7) for Expressing the Cut-off Rigidity Dependenceof the g Data
ParameterProtons

Helium ionsElectrons

Positrons

Photons

g1min (cm⫺2 s⫺1)
g1max (cm⫺2 s⫺1)
g2min (cm2 g⫺1)
g2max (cm2 g⫺1)
g3min
g3max
g4min (cm2 g⫺1)
g4max (cm2 g⫺1)
g5
g6(cm2 g⫺1)
g1min(cm⫺2 s⫺1)
g1max (cm⫺2 s⫺1)
g2min (cm2 g⫺1)
g2max (cm2 g⫺1)
g3min
g3max
g4min (cm2 g⫺1)
g4max (cm2 g⫺1)
g5
g6(cm2 g⫺1)
g1min (cm⫺2 s⫺1)
g1max (cm⫺2 s⫺1)
g2min (cm2 g⫺1)
g2max (cm2 g⫺1)
g3min
g3max
g4min (cm2 g⫺1)
g4max (cm2 g⫺1)
g5
g6(cm2 g⫺1)
g1min (cm⫺2 s⫺1)
g1max (cm⫺2 s⫺1)
g2min (cm2 g⫺1)
g2max (cm2 g⫺1)
g3min
g3max
g4min (cm2 g⫺1)
g4max (cm2 g⫺1)
g5
g6 (cm2 g⫺1)
g1min (cm⫺2 s⫺1)
g1max (cm⫺2 s⫺1)
g2min (cm2 g⫺1)
g2max (cm2 g⫺1)
g3min
g3max
g4min (cm2 g⫺1)
g4max (cm2 g⫺1)
g5
g6 (cm2 g⫺1)

hi1

hi2 (GV⫺1)hi3hi4 (GV)

hi5 (GV)

0.002440.00255
0.00779
0.00768
0.963
0.974
0.007810.00735
0.191
5.71 ⫻ 10⫺4⫺2.00 ⫻ 10⫺5⫺5.21 ⫻ 10⫺50.00566
0.005240.925
0.918
0.008930.008310.212
5.57 ⫻ 10⫺40.04660.01250.006900.006481.04
1.06
0.01610.01660.464⫺8.21 ⫻ 10⫺50.02840.007020.00655
0.00633
1.07
1.07
0.01700.01670.463⫺7.12 ⫻ 10⫺5
0.528
0.150
0.00532
0.00523
1.07
1.06
0.01390.01400.464⫺5.75 ⫻ 10⫺5⫺6.03 ⫻ 10⫺5
⫺7.18 ⫻ 10⫺5
⫺9.58 ⫻ 10⫺6
⫺2.41 ⫻ 10⫺60.00160
0.00106
9.71 ⫻ 10⫺11
2.56 ⫻ 10⫺50.07036.13 ⫻ 10⫺6
1.79 ⫻ 10⫺6
3.81 ⫻ 10⫺6
7.51 ⫻ 10⫺59.97 ⫻ 10⫺50.00260
0.00332⫺5.75 ⫻ 10⫺5
⫺3.19 ⫻ 10⫺50.0769⫺1.81 ⫻ 10⫺5⫺0.00111
4.30 ⫻ 10⫺4
⫺8.92 ⫻ 10⫺6
9.61 ⫻ 10⫺60.005540.00416⫺1.64 ⫻ 10⫺10⫺3.07 ⫻ 10⫺50.0255⫺1.07 ⫻ 10⫺5
⫺7.92 ⫻ 10⫺4
2.24 ⫻ 10⫺46.30 ⫻ 10⫺6
1.63 ⫻ 10⫺50.005360.00532⫺6.99 ⫻ 10⫺6
⫺2.52 ⫻ 10⫺80.0255⫺1.13 ⫻ 10⫺5⫺0.0146
0.00277
4.99 ⫻ 10⫺55.20 ⫻ 10⫺50.00291
0.00300⫺1.71 ⫻ 10⫺5⫺2.47 ⫻ 10⫺50.0255⫺1.19 ⫻ 10⫺50.002200.00146
6.22 ⫻ 10⫺4
6.64 ⫻ 10⫺4⫺0.0712⫺0.0214
8.24 ⫻ 10⫺40.00125⫺0.6455.47 ⫻ 10⫺49.01 ⫻ 10⫺5
9.88 ⫻ 10⫺50.00275
0.00309⫺0.828⫺0.133
7.47 ⫻ 10⫺49.70 ⫻ 10⫺4⫺0.6204.56 ⫻ 10⫺40.02340.04106.86 ⫻ 10⫺48.85 ⫻ 10⫺4⫺0.0808⫺0.05130.002530.00104⫺0.3300.001030.00963
0.02230.001190.00109⫺0.240⫺0.06780.00255
0.00213⫺0.3290.001020.208
0.414
0.00271
0.00249⫺0.145⫺0.03180.00135
7.20 ⫻ 10⫺4⫺0.3299.58 ⫻ 10⫺4

6.68
6.92
7.78
7.75
2.23
3.01
8.51
8.19
2.03
1.11
6.87
8.26
8.20
8.21
⫺0.637
2.26
5.53
7.85
2.47
0.943
6.80
9.16
7.31
8.93
4.49
5.70
4.61
6.63
3.79
1.19
6.43
8.92
6.88
8.61
1.92
5.37
4.84
6.48
3.78
1.04
5.86
8.70
7.90
9.01
0.00622
4.84
5.38
7.27
3.79
1.41

0.932
0.994
1.85
1.94
0.788
0.918
2.31
2.94
1.30
0.837
4.82
3.81
4.96
4.67
1.90
1.41
1.34
3.74
1.43
1.25
2.69
4.35
3.51
4.22
1.56
2.11
3.42
1.89
1.33
4.86
2.16
4.29
4.54
4.47
2.87
1.47
3.08
3.62
1.33
4.85
2.58
4.69
7.21
7.25
3.42
1.75
2.16
1.70
1.33
4.67

reach these lower altitudes. Hence we regarded gi for rc ⬍1 GV as the constant that obtained from Eq. (7) by substituting 1.0 for rc; i.e., the minimum of the applicable cutoff rigidity in Eq. (7) is 1 GV. This limitation is also appliedin Eq. (17) defined later in this paper.
To determine ⌽N for arbitrary solar conditions, we assumed that the dependence of ⌽N on the force field potential s can be written as⌽N (s, rc, d) ⫽ f1(rc, d) ⫹ f2(rc, d)s f3(rc,d),(8)where fi are the parameters depending on rc and d. The
power index of s, f3(rc, d), generally increases with the riseof the atmospheric depth. We assumed that their relation
can be represented by the linear function
f3(rc, d) ⫽ g5(rc) ⫹ g6(rc)d,(9)where g5 and g6 are parameters depending on rc as ex-

252

SATO ET AL.

FIG. 6. Proton fluxes at 1 MeV/nucleon, ⌽N, as a function of (panel A) atmospheric depth and (panel B) force
field potential.

pressed by Eq. (7). The numerical values of the parameters
h5j and h6j were obtained from the least-squares fitting of
the Monte Carlo data of ⌽N for all force field potentials—
400, 600, 900, 1200 and 1800 MV. The results of this leastsquares fitting are also provided in Table 3. The parameters
f1(rc, d) and f2(rc, d) were determined by solving the simultaneous equations
f 3(r c ,d)⌽Nmin (rc , d) ⫽ f 1 (rc , d) ⫹ f 2 (rc , d)smin
andf 3(r c ,d)⌽Nmax (rc , d) ⫽ f 1 (rc , d) ⫹ f 2 (rc , d)smax,

influencing the smoothness of the spectrum switching. The
functions of the hyperbolic tangent were introduced only
for expressing the gradual switching of the two spectra, andthe form of the equation has little physical meaning. In
general, the switching energy for the primary spectra Es1 isequal to that for the secondary spectra Es2 and can be determined from rc and d by the equationEs(rc, d) ⫽ a13[Ec(rc) ⫺ a14d],(10)

using f3, and ⌽Nmin and ⌽Nmax calculated by Eqs. (9) and(6), respectively. Substituting the obtained fi parameters
into Eq. (8), ⌽N for arbitrary solar conditions can be estimated.
As examples of the calculation results, the dependencesof ⌽N for proton on the atmospheric depth and the force
field potential are shown in Fig. 6 in comparison with thecorresponding Monte Carlo-obtained data. It is evident
from the figure that the analytical and Monte Carlo results
agree closely with each other, indicating the validity of theequations for expressing the changes in ⌽N. Substituting bi
and ⌽N calculated by Eqs. (5) and (8), respectively, into
Eq. (4), we can analytically estimate the secondary protonand helium-ion spectra.
Combining Primary and Secondary Spectra of Protons
and Helium Ions
Considering the influence of the cut-off rigidity, the terrestrial cosmic-ray proton and helium-ion spectra ␾ can be
estimated from their primary and secondary spectra by␾(s, rc , d, E)⫽ ␾ pri (s, d, E)[tanh{a11 [E/Es1 (rc , d) ⫺ 1]} ⫹ 1]/2
⫹ ␾ sec (s, rc , d, E)[tanh{a12 [1 ⫺ E/Es2 (rc , d)]} ⫹ 1]/2,(11)
where Esi is the switching energy between the primary andsecondary spectra, and a11 and a12 are constant parameters(12)

where a13 and a14 are constant parameters, and Ec corresponds to the cut-off energy of the particle at the top of
the atmosphere. This can be obtained by Ec ⫽[兹(1000 ⫻ rc Z) 2 ⫹ E2m ⫺ Em]/A, where rc is given in GV
to obtain Ec in MeV/nucleon. However, it is obvious thatmost of the lower-energy particles are produced in the atmosphere, i.e. are secondary particles, even for very lower
cut-off rigidity cases, and hence we introduced minimum
values of Es for the primary and secondary particles in theequationsEs1 (rc , d) ⫽ max[a15 , Es (rc , d)]
Es2 (rc , d) ⫽ max[a16 , Es (rc , d)],(13)

respectively, where a15 and a16 represent the minimum energy. The numerical values of the parameters a11 to a16 were
determined from the least-squares fitting of the Monte Carlo-obtained spectra, using the primary and secondary spectra calculated using Eqs. (1) and (4), respectively. The results of the fitting are also summarized in Table 1.
The PARMA results for proton and helium-ion spectra
shown in Figs. 2 and 4 were obtained by substituting all
the parameters and equations given in this section into Eq.
(11). The figures clearly indicate the ability of PARMA toreproduce the Monte Carlo simulation results except for the
proton spectra below 10 MeV and the helium-ion spectra
at intermediate energies, around 100 MeV/nucleon. The
disagreement in the helium-ion spectra is attributed to thefact that helium ions with such an intermediate energy weregenerally produced by the pre-equilibrium or nuclear fragmentation processes, but we assumed that all secondary he-DEVELOPMENT OF PARMA

253

FIG. 7. Comparison between the cosmic-ray electron, positron and photon spectra obtained by our Monte Carlosimulation and PARMA. The left and right panels show the atmospheric depth and the cut-off rigidity dependenceof the spectra, respectively.

lium ions were created by the evaporation process in thederivation of PARMA. However, the doses from protons
below 10 MeV and helium ions between 20 MeV/nucleonand 1 GeV/nucleon are less than 1% of their total values
in most cases. Hence the effects of these discrepancies on
the dose estimation can be considered to be negligible.
Electron, Positron and Photon Spectra
Figure 7 shows the cosmic-ray electron, positron andphoton spectra obtained by the Monte Carlo simulation,
together with the corresponding data calculated by
PARMA. In the simulation, all electrons, positrons and photons were generated in the Earth system, i.e. were second-ary particles, since only protons and heavy ions were considered to be source particles incident to the atmosphere.
We therefore assumed that the electron, positron and photon
spectra can be estimated in a manner similar to the secondary proton spectra described before. The difference between the calculation procedures for secondary proton
spectra and the others is that we normalized the electron,
positron and photon spectra to their fluxes at 10 MeV instead of those at 1 MeV. This is because the photon fluxes
below a few MeV are influenced by the cut-off rigidity andthus are inadequate for use in the normalization of the spectra. The parameters cij and hij used in Eqs. (5) and (7),respectively, for estimating the atmospheric electron, posi-

254

SATO ET AL.

FIG. 8. Comparison between the cosmic-ray muon spectra obtained by our Monte Carlo simulation and PARMA.
The left and right panels show the atmospheric depth and the cut-off rigidity dependence of the spectra, respectively.

tron and photon spectra are also listed in Tables 2 and 3,
where their numerical values were determined from theleast-squares fittings of the corresponding Monte Carlodata.
The PARMA results for electron, positron and photon
spectra shown in Figs. 2 and 7 were obtained by substituting the results of Eqs. (5) and (8) into Eq. (4), using parameters listed in Tables 2 and 3. It is evident from thegraphs that PARMA can reproduce the corresponding Monte Carlo data for energies below 1 GeV. Above 1 GeV, on
the other hand, the analytical fluxes obtained for photons
as well as electrons and positrons are generally larger thanthe Monte Carlo data. Thus their doses estimated by
PARMA are approximately 10% higher than those obtainedby the Monte Carlo simulation. However, this disagreement
is not directly connected to a decrease in the reliability of
PARMA with respect to the aircrew dose estimation, since
their spectra in the high-energy region obtained by theMonte Carlo simulation are supposed to be underestimatedfor two reasons: (1) High-energy electrons and positrons
were not transported, and (2) incident cosmic rays withenergy above 200 GeV/nucleon, which can produce a large
number of high-energy photons, electrons and positrons,
were not considered in our simulation.
Muon Spectra
Figure 8 shows the dependence of muon spectra on theatmospheric depth and the cut-off rigidity. Although all of

the cosmic-ray muons are secondary particles, as are theelectrons, positrons and photons, the muon spectra cannot
be estimated by the same procedure, i.e. normalizing theirspectra to their flux at a certain energy. This is because theatmospheric depth dependence of the muon fluxes even atlower energies cannot be expressed by Eq. (6) because thestrong penetrability of muons prevents equilibrium between
the numbers of incoming and outgoing particles. We therefore constructed the model for predicting atmospheric muon
spectra based on that for estimating underground muon
spectra proposed by Lipari et al. (36).
According to their model, underground muon spectra dN/dE can be expressed bydN⫽ ␣1 exp[(1 ⫺ ␣3 )␣2 X]
dE

冦

⫻ E⫹

⫺␣3

冧 , (14)

␣4 X[1 ⫺ exp(⫺␣2 X)]
␣2 X

where X denotes the depth of rock and ␣i are constants.
Note that the notations of some parameters are changed
from the original equation to avoid the duplicate use of anotation in this paper. For the boundary condition, they
assumed that the muon spectrum at the ground level canbe expressed as ␣1 E⫺␣3, where ␣3 ⫽ 3.7, close to the spectrum expected for muons generated by pion and kaon decays in the atmosphere. The parameters ␣2 and ␣4 are re-255

DEVELOPMENT OF PARMA

lated to the radiation and ionization energy loss of highenergy muons, respectively, and their stopping power canbe written as ␣2E ⫹ ␣4.
This equation, however, cannot be applied directly to theestimation of the atmospheric muon spectra, since not only
deceleration but also production and decay of muons must
be considered in their calculation. To take the production
into account, ␣1 should be replaced by functions that depend on the atmospheric depth in a complicated manner.
Furthermore, the assumption that ␣4 is a constant is inadequate for our purpose, since it depends significantly on ␤
at lower energies and increases slightly with a logarithmic
increase of muon energy at higher energies, as expectedfrom the Bethe-Bloch formula. On the other hand, ␣2X is
expected to be very small for atmospheric muons, since thetrack length of a muon in the atmosphere is at most 103
g/cm2, which is much shorter than 1/␣2. Consequently, theapproximations of exp[(1 ⫺ ␣3) ␣2 X] ⬃ 1 and [1 ⫺exp(⫺␣2 X)] ⬃ ␣2 X are established.
From these considerations, we proposed a best-fit curveto the atmospheric muon fluxes for the solar minimum condition, ␾␮min, as generated by the equation
␾ ␮min (rc , d, E)[]t (r , d) ⫹ t2min (rc , d)log 10 (E)⫽ ⌽␮ (d) E ⫹ 1min c
␤ t3min(rc ,d)

⫺␣3

,
(15)

where ⌽␮ is an index for the high-energy muon fluxes at
the atmospheric depth d, and ti are parameters related tothe mean ionization energy loss of muons during the transport through the atmosphere. Note that ⌽␮ depends only on
the atmospheric depth, since the high-energy muon fluxes
are independent of the solar modulation and the cut-off rigidity. In the derivation of this equation, we assumed thatthe mean energy loss due to the ionization can be expressedby a linear function of the common logarithm of E divided
by a power function of ␤, since it depends on the global
condition in a complicated manner and cannot be calculated
directly by theoretical formulas such as the Bethe-Blochformula. The muon fluxes for the solar maximum condition,
␾␮max, can be also estimated from Eq. (15) by replacing thesubscript ‘‘min’’ with ‘‘max’’ in the equation.
For a best-fit curve to ⌽␮, we introduced a function withthe form similar to Eq. (6) as written by⌽␮(d) ⫽ u1[exp(⫺u2 d) ⫺ u3 exp(⫺u4 d)] ⫹ u5,(16)

where ui are free parameters. For that to ti, we simply adopted the fourth order of the polynomial function:ti (rc , d) ⫽ vi1 (rc ) ⫹ vi 2 (rc )d ⫹ vi 3 (rc )d 2
⫹ vi 4 (rc )d 3 ⫹ vi 5 (rc )d 4 ,(17)where vij are parameters related to the solar modulation andthe cut-off rigidity. For expressing the dependence of vij on

rc, we employed a function with the same form as Eq. (7)as written by
vi j (rc ) ⫽ wi j1 ⫹ wi j2 rc ⫹wi j3,
1 ⫹ exp[(rc ⫺ wi j4 )/wi j 5 ](18)

where wijk are free parameters. The parameters vij for rc ⬍1 GV are regarded as the constant values of those at rc ⫽1 GV, as described before. The numerical values of theparameters ui together with wijk for the solar minimum andmaximum conditions were determined from the leastsquares fitting of the Monte Carlo-obtained ␾␮ data for theforce field potentials 400 and 1200 MV, respectively. In thefitting, the numerical value of ␣3 is fixed at 3.7, following
the original Lipari model. The atmospheric muon spectra
for arbitrary solar conditions can be estimated from ␾␮minand ␾␮max calculated by Eq. (15), assuming that ␾␮ depends
on the force field potential in the same manner of ⌽N as
written by Eq. (8). For muons, the numerical values of theparameters h5j and h6j in Eq. (7) for determining g5 and g6
in Eq. (9) were obtained from the least-squares fitting of
the Monte Carlo-obtained energy-integrated ␾␮ for all force
field potentials—400, 600, 900, 1200 and 1800 MV. The
results of these least-squares fittings are summarized in Table 4.
The PARMA results for muon spectra shown in Figs. 2
and 8 were calculated by substituting the evaluated ui inthis section into Eq. (15). It is evident from the graphs that
the PARMA results closely agree with the Monte Carlo data
except for energies below 10 MeV and above 20 GeV. The
discrepancy in the low-energy region is not important in
the dose estimation, since doses from such low-energy muons are negligibly small—much less than 1% of their total
values. That in the high-energy region is attributed to thefact that such high-energy muons are generally producedby nuclear reactions caused by cosmic rays with energies
above 200 GeV, which are not considered in the Monte
Carlo simulation. According to theory, the high-energy
muon spectra can be simply expressed by a power functionof the muon energy, and hence the PARMA results are
more reliable than the corresponding Monte Carlo data in
the high-energy region. This tendency is also verified by
the comparison of the Monte Carlo-obtained and experimental muon spectra shown in Fig. 3C and D, where theMonte Carlo simulation underestimates the experimental
data in the high-energy region.
Comparison with Monte Carlo Simulation in Terms of
Dose Estimation
To verify the agreement between PARMA and the Monte
Carlo simulation in dose estimation, the ratios of the doses
calculated by PARMA to those by the Monte Carlo simulation were evaluated for 1620 global conditions: five force
field potentials from 400 to 1800 MV, 18 geomagnetic
fields with vertical cut-off rigidities from 0 to 17 GV, and18 altitude ranges from sea level to 20 km. In the dose

256

SATO ET AL.

TABLE 4
Numerical Values of the Parameters ui, wijk and hij used in Eqs. (16), (18) and (7), respectively, to Estimate theAtmospheric Muon Spectra␮⫹Parameter⌽␮Parameter

u1
6.26 ⫻ 109
wij1

u2
0.00343
wij2 (GV⫺1)

u3
1.01
wij3

u4
0.00418
wij4 (GV)

u5
3.75 ⫻ 108
wij5 (GV)

v11minv11maxv12minv12maxv13minv13maxv14minv14maxv15minv15maxv21minv21maxv22minv22maxv23minv23maxv24minv24maxv25minv25maxv31minv31maxv32minv32maxv33minv33maxv34minv34maxv35minv35max2.05 ⫻ 103
2.39 ⫻ 103
⫺5.67
⫺5.62
0.01170.0117⫺2.31 ⫻ 10⫺6
⫺2.24 ⫻ 10⫺6
1.74 ⫻ 10⫺91.75 ⫻ 10⫺9
84.8
87.3
3.41
3.41
⫺0.00332⫺0.00331⫺2.68 ⫻ 10⫺6
⫺2.68 ⫻ 10⫺6
2.33 ⫻ 10⫺92.32 ⫻ 10⫺9
0.760
0.923
0.00206
8.44 ⫻ 10⫺4⫺5.96 ⫻ 10⫺6
⫺3.91 ⫻ 10⫺6
6.46 ⫻ 10⫺9
1.99 ⫻ 10⫺9
⫺3.21 ⫻ 10⫺12⫺1.78 ⫻ 10⫺12126118⫺0.655⫺0.6580.001570.00157⫺7.60 ⫻ 10⫺7
⫺7.56 ⫻ 10⫺7
⫺2.22 ⫻ 10⫺10⫺2.26 ⫻ 10⫺10
⫺5.77
⫺5.90
0.07870.0785⫺1.49 ⫻ 10⫺4
⫺1.49 ⫻ 10⫺4
⫺8.88 ⫻ 10⫺8⫺8.81 ⫻ 10⫺8
1.49 ⫻ 10⫺10
1.49 ⫻ 10⫺10⫺0.0180⫺0.02966.17 ⫻ 10⫺5
1.34 ⫻ 10⫺4
⫺1.48 ⫻ 10⫺7
⫺2.88 ⫻ 10⫺7
⫺9.28 ⫻ 10⫺12
3.57 ⫻ 10⫺10
5.46 ⫻ 10⫺14⫺3.17 ⫻ 10⫺14⫺1.01 ⫻ 103
⫺949
3.59
3.28
⫺0.0125⫺0.0124
2.48 ⫻ 10⫺5
2.48 ⫻ 10⫺5
⫺1.69 ⫻ 10⫺8
⫺1.69 ⫻ 10⫺8370377⫺0.520⫺0.5230.001850.00185⫺2.71 ⫻ 10⫺6
⫺2.71 ⫻ 10⫺6
1.20 ⫻ 10⫺9
1.20 ⫻ 10⫺9⫺0.273⫺0.4280.001780.00181⫺4.13 ⫻ 10⫺6
⫺2.49 ⫻ 10⫺6
1.74 ⫻ 10⫺92.95 ⫻ 10⫺99.21 ⫻ 10⫺13
4.79 ⫻ 10⫺13

6.18
7.04
1.31
1.06
3.26
3.31
4.94
4.89
5.12
5.18
4.81
4.59
6.87
6.84
7.02
7.02
7.04
7.04
7.04
7.05
11.3
9.66
7.55
9.26
7.53
9.74
23.6
10.4
7.54
4.21

3.47
3.84
3.22
3.34
3.65
3.58
3.80
3.80
4.39
4.40
3.36
3.39
1.09
1.09
0.607
0.611
0.4685
0.472
0.364
0.367
5.39
4.00
3.93
2.44
4.39
1.49
1.67
1.94
2.66
0.747

Parameterg5
g6(cm2g⫺1)

hi1
0.506
1.39 ⫻ 10⫺4hi2 (GV⫺1)
0.01306.95 ⫻ 10⫺6hi3
⫺0.394
7.47 ⫻ 10⫺4hi4 (GV)
4.12
3.72

hi5 (GV)
1.33
1.97

calculation, we adopted the fluence to effective dose conversion coefficients for the isotropic irradiation geometry
calculated by PHITS (37, 38) coupled with the updated
radiation weighting factor defined in ICRP publication 103(39). It should be noted that the Monte Carlo simulation
took more than 1 month using a parallel computer with 24
CPUs, while the dose calculation by PARMA took only 10s using a conventional PC. This difference clearly indicatesthe efficiency of PARMA when it is adopted in the routedose calculation.
Figure 9 shows the distributions of the ratios of the doses
calculated by PARMA to those obtained from the Monte
Carlo simulation for each particle type as well as that for
the total dose. It should be noted that the total doses obtained by the Monte Carlo simulation include the doses
from particles whose spectrum model was not developed inthis paper, i.e. pions and heavy ions with Z ⱖ 3, but theircontributions are generally negligible—less than 1% of thetotal. It is evident from this graph that the total doses cal-culated by PARMA and the Monte Carlo simulation agreewith each other within 5% for more than 99% of the global
conditions, indicating the equivalence of the PARMA andthe Monte Carlo simulation with respect to the accuracy ofdose estimation. The distribution of this ratio for helium
ions is broad in comparison to those for other particles, but
most of the scattered data are for lower altitudes, where thecontribution of helium ions to the total dose is negligible.
The ratios for electrons, positrons and photons are generally
greater than 1 by approximately 10%, since PARMA gives
higher values for these fluxes above 1 GeV than to theMonte Carlo simulation, as described before.
CONCLUSIONS

Monte Carlo simulations were performed for estimatingcosmic-ray spectra in the atmosphere using the PHITS
code. Excellent agreement was observed between the calculated and measured spectra for various conditions, except

257

DEVELOPMENT OF PARMATABLE 4
Extended␮⫺u1
5.82 ⫻ 109
wij1

u2
0.00362
wij2 (GV⫺1)

u3
1.02
wij3

u4
0.00451
wij4 (GV)

u5
3.20 ⫻ 108
wij5 (GV)

2.09 ⫻ 103
2.42 ⫻ 103
⫺5.66
⫺5.61
0.01180.0118⫺2.59 ⫻ 10⫺6
⫺2.52 ⫻ 10⫺6
1.87 ⫻ 10⫺9
1.86 ⫻ 10⫺9
85.9
87.0
3.42
3.42
⫺3.33 ⫻ 10⫺3
⫺3.32 ⫻ 10⫺3
⫺2.69 ⫻ 10⫺6
⫺2.68 ⫻ 10⫺6
2.34 ⫻ 10⫺92.33 ⫻ 10⫺9
0.787
0.814
0.002146.65 ⫻ 10⫺4⫺6.05 ⫻ 10⫺6⫺3.80 ⫻ 10⫺66.68 ⫻ 10⫺92.75 ⫻ 10⫺9
⫺3.10 ⫻ 10⫺12⫺1.81 ⫻ 10⫺12121
112⫺0.650⫺0.6510.00158
0.00158⫺7.99 ⫻ 10⫺7⫺7.93 ⫻ 10⫺7
⫺1.98 ⫻ 10⫺10⫺2.01 ⫻ 10⫺10
⫺5.86
⫺5.88
0.07900.0787⫺1.49 ⫻ 10⫺4⫺1.50 ⫻ 10⫺4
⫺9.00 ⫻ 10⫺8⫺8.93 ⫻ 10⫺8
1.50 ⫻ 10⫺10
1.50 ⫻ 10⫺10⫺0.0180⫺0.0248
4.99 ⫻ 10⫺5
1.35 ⫻ 10⫺4
⫺1.36 ⫻ 10⫺7
⫺2.92 ⫻ 10⫺7
1.09 ⫻ 10⫺12
3.35 ⫻ 10⫺103.80 ⫻ 10⫺14⫺4.17 ⫻ 10⫺14

⫺929
⫺895
3.58
3.31
⫺0.0125⫺0.0125
2.54 ⫻ 10⫺5
2.53 ⫻ 10⫺5⫺1.71 ⫻ 10⫺8
⫺1.70 ⫻ 10⫺8369
372
⫺0.529⫺0.5320.001860.00186⫺2.71 ⫻ 10⫺6
⫺2.71 ⫻ 10⫺6
1.19 ⫻ 10⫺9
1.20 ⫻ 10⫺9⫺0.304⫺0.3110.001430.00184⫺3.94 ⫻ 10⫺6
⫺2.58 ⫻ 10⫺6
1.58 ⫻ 10⫺92.31 ⫻ 10⫺9
7.46 ⫻ 10⫺13
4.63 ⫻ 10⫺13

6.86
7.45
0.89
0.78
3.44
3.48
4.95
4.92
5.12
5.14
4.82
4.68
6.88
6.86
7.04
7.04
7.05
7.05
7.05
7.04
14.5
10.6
8.10
9.29
7.83
9.67
22.7
10.3
7.85
4.54

3.29
3.55
3.73
3.76
3.65
3.59
3.72
3.74
4.24
4.27
3.30
3.30
1.06
1.09
0.602
0.601
0.464
0.465
0.356
0.362
5.61
3.61
3.46
2.39
4.34
1.38
1.99
1.37
2.00
0.479

hi1
0.565
8.80 ⫻ 10⫺5hi2 (GV⫺1)
0.0121⫺3.89 ⫻ 10⫺6hi3⫺0.357
4.91 ⫻ 10⫺4hi4 (GV)
4.73
4.51

hi5 (GV)
1.46
1.72

for the electron data. Further study is needed to improve
the accuracy of the PHITS simulation in regard to electron
and positron transport, since the cosmic-ray electron spectra
obtained by the simulation do not agree well with someexperimental data due to the uncertain calculation technique for dealing with the transport of high-energy electrons and positrons in the code. Based on a comprehensive
analysis of the simulation results, we proposed an analyticalmodel for estimating the atmospheric cosmic-ray spectra
for neutrons, protons, helium ions, muons, electrons, positrons and photons that is applicable to any global conditions
at altitudes below 20 km. The PARMA model enables us
to calculate the cosmic-radiation doses instantaneously with
precision equivalent to that of the Monte Carlo simulationthat requires much computational time, although some discrepancies are observed between their calculated spectra forcertain particle types and energies. One shortcoming of
PARMA is that it adopts the vertical cut-off rigidity insteadof its full two-dimensional distribution in the considerationof the magnetosphere on the cosmic-ray spectra, but theresulting errors are expected to be small except for the magnetic equator region, as discussed in ref. (5). We therefore
conclude that PARMA can improve the accuracy and efficiency of cosmic-ray exposure dose estimations not only
for aircrews but also for the public on the ground.
For the practical use of PARMA, software based on themodel that we named EXPACS was developed for calculating atmospheric cosmic-ray spectra. It has been released
to the public online (40). In the near future, PARMA will
be incorporated into the Japanese Internet System for Calculation of Aviation Route Doses (JISCARD) (41) and used
for adhering to the dose limit (⬍5 mSv/year) recommended
for aircrews of Japanese airline companies. The accuracyof the updated JISCARD is currently being evaluated bycomparing its calculated route doses with the corresponding
experimental data under various flight conditions, and theresults will be presented in our forthcoming paper. The simulation technique established by this work is also capable

258

SATO ET AL.

FIG. 9. Distributions of the ratios of doses calculated by PARMA to
those obtained by our Monte Carlo simulation.

of contributing to the estimation of the atmospheric cosmicray spectra under solar-geomagnetic storm conditions, anarea requiring additional study. Application of this work to
the estimation of cosmic-ray spectra on the Martian surface
will be of great interest in future NASA human space explorations.
ACKNOWLEDGMENTS
We would like to thank Dr. Y. Watanabe and Dr. T. Fukahori for theirsupport in incorporating JENDL/HE into PHITS, Dr. M. Honda and Dr.
K. O’Brien for their advice on constructing the model for estimating thecosmic-ray spectra at the top of the atmosphere, and Dr. T. Nakamura for
his valuable comments on this work. We also wish to thank the technicalstaff of the CCSE office of JAEA for their help in performing the Monte
Carlo simulation. This work was partially supported by ‘‘Ground-based
Research Program for Space Utilization’’ promoted by the Japan Space
Forum.
Received: May 18, 2007; accepted: March 24, 2008

REFERENCES
1. ICRP, 1990 Recommendations of the International Commission on
Radiological Protection. Publication 60, Annals of the ICRP, Vol.
21, No. 1–3, Pergamon Press, Oxford, 1991.
2. H. Schraube, V. Mares, S. Roesler and W. Heinrich, Experimental
verification and calculation of aviation route doses. Radiat. Prot. Dosimetry 86, 309–315 (1999).
3. W. Friedberg, K. Copeland, F. E. Duke, K. O’Brien and E. B. Darden,
Guidelines and technical information provided by the U.S. Federal
Aviation Administration to promote radiation safety for air carriercrew members. Radiat. Prot. Dosimetry 86, 323–327 (1999).

4. B. J. Lewis, L. G. I. Bennett, A. R. Green, A. Butler, M. Desormeaux,F. Kitching, M. J. McCall, B. Ellaschuk and M. Pierre, Aircrew dosimetry using the predictive code for aircrew radiation exposure(PCAIRE). Radiat. Prot. Dosimetry 116, 320–326 (2005).
5. K. O’Brien, W. Friedberg, H. H. Sauer and D. F. Smart, Atmosphericcosmic rays and solar energetic particles at aircraft altitudes. Environ.
Int. 22 (Suppl. 1), S9–S44 (1996).
6. D. Heck, J. Knapp, J. N. Capdevielle, G. Schatz and T. Thouw, CORSIKA: A Monte Carlo Code to Simulate Extensive Air Showers. Report FZKA 6019, Forschungszentrum Karlsruhe, 1998. http://
www.ik.fzk.de/corsika/.
7. COSMOS: A versatile Monte Carlo simulation code for propagationof cosmic rays in the atmosphere and near earth environment.
http://cosmos.n.kanagawa-u.ac.jp/cosmosHome/index.html.
8. PLANETOCOSMICS GEANT4. http://cosray.unibe.ch/⬃laurent/
planetocosmics/.
9. A. Fasso, A. Ferrari and P. R. Sala, Electron-photon transport in FLUKA: status. In Proceedings of the Monte Carlo 2000 Conference,
Lisbon, Oct. 23–26 (A. Kling, F. Barao, M. Nakagawa, L. Tavora
and P. Vaz, Eds.), pp. 159–164. Springer-Verlag, Berlin, 2001.
10. A. Fasso, A. Ferrari, J. Ranft and P. R. Sala, FLUKA: status andprospective for hadronic applications. In Proceedings of the Monte
Carlo 2000 Conference, Lisbon, Oct. 23–26 (A. Kling, F. Barao, M.
Nakagawa, L. Tavora, P. Vaz, Eds.), pp. 955–960. Springer-Verlag,
Berlin, 2001.
11. J. M. Clem, G. De Angelis, P. Goldhagen and J. W. Wilson, New
calculations of the atmospheric cosmic radiation field—results forneutron spectra. Radiat. Prot. Dosimetry 110, 423–428 (2004).
12. H. Iwase, K. Niita and T. Nakamura, Development of a generalpurpose particle and heavy ion transport Monte Carlo code. J. Nucl.
Sci. Technol. 39, 1142–1151 (2002).
13. T. Fukahori, Y. Watanabe, N. Yoshizawa, F. Maekawa, S. Meigo, C.
Konno, N. Yamano, A. Yu. Konobeyev and S. Chiba, JENDL HighEnergy File. J. Nucl. Sci. Technol. (Suppl. 2), 25–30 (2002).
14. Y. Watanabe, T. Fukahori, K. Kosako, N. Shigyo, T. Murata, N. Yamano, T. Hino, K. Maki, H. Nakashima, N. Odano and S. Chiba,
Nuclear data evaluations for JENDL high-energy file. In Proceedingsof International Conference on Nuclear Data for Science and Technology, pp. 326–331. AIP CP769, American Institute of Physics,
Melville, NY, 2005.
15. T. Sato and K. Niita, Analytical functions to predict cosmic-ray neutron spectra in the atmosphere. Radiat. Res. 166, 544–555 (2006).
16. K. O’Brien, H. A. Sandmeier, G. E. Hansen and J. E. Campbell,
Cosmic-ray-induced neutron background sources and fluxes for geometries of air over water, ground, iron, and aluminum. J. Geophys.
Res. 83, 114–120 (1978).
17. S. Roesler, W. Heinrich and H. Schraube, Calculation of radiation
fields in the atmosphere and comparison to experimental data. Radiat.
Res. 149, 87–97 (1998).
18. R. A. Nymmik, M. I. Panasyuk, T. I. Pervaja and A. A. Suslov, A
model of galactic cosmic ray fluxes. Nucl. Tracks Radiat. Meas. 20,
427–429 (1992).
19. L. J. Gleeson and W. I. Axford, Modulation of galactic cosmic rays.
Astrophys. J. 154, 1011–1026 (1968).
20. R. A. Caballero-Lopez and H. Moraal, Limitations of the force field
equation to describe cosmic ray modulation. J. Geophys. Res. 109,
A01101 (2004).
21. K. O’Brien, E. Felsberger and P. Kindl, Application of the heliocentric potential to aircraft dosimetry. Radiat. Prot. Dosimetry 116, 336–
342 (2005).
22. Y. Shikaze, S. Haino, K. Abe, H. Fuke, T. Hams, K. C. Kim, Y.
Makida, S. Matsuda, J. W. Mitchell and K. Yoshimura, Measurementsof 0.2–20 GeV/n cosmic-ray proton and helium spectra from 1997
through 2002 with the BESS spectrometer. Astropart. Phys. 28, 154–
167 (2007).
23. Index to /ftp/pub/WDCCR. http://www.env.sci.ibaraki.ac.jp/ftp/pub/
WDCCR.

DEVELOPMENT OF PARMA

24. H. W. Bertini, Low-energy intranuclear cascade calculation. Phys.
Rev. 131, 1801–1821 (1963).
25. P. Goldhagen, J. M. Clem and J. W. Wilson, The energy spectrum of
cosmic-ray induced neutrons measured on an airplane over a wide
range of altitude and latitude. Radiat. Prot. Dosimetry 110, 387–392(2004).
26. N. M. Kocharian, G. S. Saakian and Z. A. Kirakosian, Energy spectraand nuclear interactions of cosmic-ray particles. Sov. Phys. JETP 35,
933–942 (1959).
27. I. S. Diggory, J. R. Hook, I. A. Jenkins and K. E. Turver, The momentum spectra of nuclear active particles in the cosmic radiation atsea level I. Experimental data. J. Phys. A Math. Nucl. Gen. 7, 741–
764 (1974).
28. H. B. Barber, T. Bowen, D. A. Delise, E. W. Jenkins, J. J. Jones, R. M.
Kalbach and A. E. Pifer, Predictions and measurements of mass spectra of the charged nucleonic component of cosmic rays at mountain
altitude. Phys. Rev. 22, 2667–2687 (1980).
29. T. Sanuki, M. Fujikawa, H. Matsunaga, K. Abe, K. Anraku, H. Fuke,
S. Haino, M. Imori, K. Izumi and K. Yoshimura, Measurement ofcosmic-ray proton and antiproton spectra at mountain altitude. Phys.
Lett. 577, 10–17 (2003).
30. J. Kremer, M. Boezio, M. L. Ambriola, G. Barbiellini, S. Bartalucci,
R. Bellotti, D. Bergstrom, U. Bravar, F. Cafagna and N. Zamp, Measurements of ground-level muons at two geomagnetic locations.
Phys. Rev. Lett. 83, 4241–4244 (1999).
31. T. Sanuki, M. Fujikawa, K. Abe, K. Anraku, Y. Asaoka, H. Fuke, S.
Haino, M. Imori, K. Izumi and K. Yoshimura, Measurement of atmospheric muon spectra at mountain altitude. Phys. Lett. B 541, 234–
242 (2002).
32. S. Haino, T. Sanuki, K. Abe, K. Anraku, Y. Asaoka, H. Fuke, M.
Imori, A. Itasaki, T. Maeno and K. Yoshimura, Measurements of259

primary and atmospheric cosmic-ray spectra with the BESS-TeV
spectrometer. Phys. Lett. B 594, 35–46 (2004).
33. S. Cecchini, T. Chiarusi, A. Pagliarin and G. Giovannini, On thealtitude dependence of ␥-rays spectra in the Earth’s atmosphere. In
Proceedings of 29th International Cosmic Ray Conference, Vol. 2,
pp. 349–352. Pune, India, 2005.
34. G. J. Fulks and P. Meyer, Cosmic ray electrons in the atmosphere. J.
Geophys. 40, 751–759 (1974).
35. M. Honda, T. Kajita, K. Kasahara, S. Midorikawa, T. Sanuki, J. Nishimura and A. Okada, Reducing uncertainty in atmospheric neutrino
flux prediction. In Proceedings of 30th International Cosmic Ray
Conference, 0274. Merida, Mexico, 2007.
36. P. Lipari and T. Stanev, Propagation of multi-TeV muons. Phys. Rev.
D 44, 3543–3554 (1991).
37. T. Sato, S. Tsuda, Y. Sakamoto, Y. Yamaguchi and K. Niita, Conversion coefficients from fluence to effective dose for heavy ions withenergies up to 3 GeV/A. Radiat. Prot. Dosimetry 106, 137–144(2003).
38. T. Sato, S. Tsuda, Y. Sakamoto, Y. Yamaguchi and K. Niita, Analysisof dose-LET distribution in the human body irradiated by high energy
hadrons. Radiat. Prot. Dosimetry 106, 145–153 (2003).
39. ICRP, The 2007 Recommendations of the International Commission
on Radiological Protection. ICRP Publication 103, Annals of the
ICRP, Vol. 37, No. 2–4, Elsevier, Oxford, 2007.
40. EXPACS, Excel-based Program for Calulating Atmospheric Cosmicray
Spectrum.
http://www.jaea.go.jp/04/nsed/ers/radiation/rpro/
EXPACS/expacs-eng.html.
41. JISCARD, Japanese Internet System for Calculation of Aviation
Route Doses. http://www.nirs.go.jp/research/jiscard/information/index.shtml. [in Japanese]