RADIATION RESEARCH 166, 544–555 (2006)

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

Analytical Functions to Predict Cosmic-Ray Neutron Spectrain the Atmosphere
Tatsuhiko Satoa,1 and Koji NiitabaJapan Atomic Energy Agency (JAEA) and b Research Organization for Information Science and Technology (RIST), Ibaraki, Japan

ical two-component solution of the Boltzmann transport
equation. The code is capable of estimating the altitudedependences of aircrew doses and integrated neutron fluxes
precisely and is adopted in the route dose calculation code
CARI-6 (3). Several authors (4–6) adopted a Monte Carloparticle transport code, FLUKA (7, 8), for the calculationof the cosmic-ray propagation in the atmosphere. Their calculation methods were successful in reproducing the neutron spectra measured at high altitudes, and a similar model
was employed in the European Program Package for the
Calculation of Aviation Route Doses, EPCARD (9). However, the cosmic-ray neutron spectra depend not only on theatmospheric depth, cut-off rigidity and solar modulation(referred to here as global conditions) but also on the structure of the aircraft (10) and the water density around thepoint of interest (4) (referred to here as local geometries)
in an intricate manner, and none of the existing models canreproduce the measured neutron spectra at any location andtime with satisfactory accuracy.
With these situations in mind, we have calculated thecosmic-ray neutron spectra by performing a Monte Carloparticle transport simulation in the atmosphere based on theParticle and Heavy Ion Transport code System, PHITS (11),using the latest version of the nuclear data library JENDLHigh-Energy File (JENDL/HE) (12, 13). In the simulation,
cosmic rays with charges up to 28 (Ni) were incident on
the Earth system represented by the concentric sphericalshells, and the neutron spectra in each shell, i.e. a certain
altitude range, were obtained. Based on a comprehensive
analysis of the simulation results, we proposed analytical
functions to predict the cosmic-ray neutron spectra at any
global condition at altitudes below 20 km, considering the
local geometry effect.
This paper describes the details of the simulation procedure, together with the comparison of results between thecalculated and experimental neutron spectra. The derivationand verification of the analytical function are also discussed
in detail below.

Sato, T. and Niita, K. Analytical Functions to Predict Cosmic-Ray Neutron Spectra in the Atmosphere. Radiat. Res. 166,
544–555 (2006).
Estimation of cosmic-ray neutron spectra in the atmosphere
has been an essential issue in the evaluation of the aircrew
doses and the soft-error rates of semiconductor devices. We
therefore performed Monte Carlo simulations for estimating
neutron spectra using the PHITS code in adopting the nuclear
data library JENDL-High-Energy file. Excellent agreements
were observed between the calculated and measured spectra
for a wide altitude range even at the ground level. Based on
a comprehensive analysis of the simulation results, we propose
analytical functions that can predict the cosmic-ray neutron
spectra for any location in the atmosphere at altitudes below
20 km, considering the influences of local geometries such as
ground and aircraft on the spectra. The accuracy of the analytical functions was well verified by various experimental
data. 䉷 2006 by Radiation Research Society

INTRODUCTION

In the last decade, protection for aircrews against terrestrial cosmic rays has been one of the most extensively discussed issues in dosimetry. Furthermore, increasing attention has been paid to the soft errors of semiconductor devices induced by the cosmic rays even at the ground level,
since the miniaturization of the devices causes a large decrease in their critical charges. These radiation effects are
triggered predominantly by neutrons produced by nuclear
reactions between the cosmic rays and atmospheric components. Therefore, estimation of cosmic-ray neutron spectra in the atmosphere is an essential issue in the evaluationof the aircrew doses and the soft-error rates (SERs).
A number of studies have been devoted to the estimationof neutron spectra. For instance, Ziegler (1) proposed anempirical formula for predicting the high-energy neutron
spectra to estimate SERs at ground level. O’Brien et al. (2)developed a deterministic code LUIN based on an analyt-

MONTE CARLO SIMULATION
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.
1

Incident Cosmic Rays
In the simulation, cosmic rays were incident on the Earthsystem at the altitude of 86 km up to which the US-Stan544

PREDICTION OF COSMIC-RAY NEUTRON SPECTRA

dard-Atmosphere-1976 has atmospheric data. Protons, ␣particles and heavy ions with charges up to 28 (Ni) wereconsidered as the source particles, although the contributions of heavy ions to the cosmic-ray neutron spectra are
generally small.
The spectra of the incident particles were estimated by
the CREME96 code (14). In the code, the effect of thegeomagnetic field on the spectra was considered by specifying the McIlwain L, where the vertical cut-off rigidity rc
in units of GV is related to this parameter by the simple
formularc ⫽ 14.5/L2.

(1)

It should be noted that the output fluxes of CREME96 are
the averaged values for all directions after accounting forthe solid Earth’s shadowing effect. Hence the values should
be converted into downward fluxes by multiplying the factor of 4␲/(4␲ ⫺ ⍀E), where ⍀E is the solid angle of the
Earth from a point at the altitude of 86 km, to be employedin our simulation.
Atmospheric Propagation
The atmospheric propagation of the incident cosmic rays
and their associated cascades was simulated by the PHITS
code, which can deal with the transports of all kinds ofhadrons and heavy ions with energies up to 200 GeV/nucleon. PHITS is based on the high-energy hadron transport
code NMTC/JAM (15) and incorporates the MCNP4C code(16) for low-energy neutron transports and the JAERI
Quantum Molecular Dynamics (JQMD) model (17) for
simulating nucleus-nucleus interactions. PHITS can also
treat the transport of photons, electrons and muons, but they
were not considered in our current simulation, since ourprimary interest was in the evaluation of the neutron spectrum, in which they play a negligible role. Note that suchelectromagnetic radiations have a significant contribution
to the total aircrew dose, and a simulation to estimate theirspectra in aircraft is to be performed in the future.
For nuclear reactions induced by neutrons and protons
below 3 GeV, the latest version of JENDL/HE was adopted
in our simulation instead of the widely used model of theintranuclear cascade (INC) (18). This is because the INC
model tends to overestimate the yields of high-energy secondary particles knocked out from light nuclei such as nitrogen and oxygen, which are the dominant components of
the atmosphere. This tendency causes the overprediction ofneutron fluences in deep-penetration calculations.
The atmosphere was divided into 28 concentric sphericalshells. The densities and temperatures of each shell were
determined referring to the US-Standard-Atmosphere-1976.
The atmosphere is assumed to be composed of 75.4% nitrogen, 23.3% oxygen and 1.3% argon by mass above thealtitude of 2 km, and additionally, 0.06% hydrogen by mass
below this altitude due to the existence of water vapor. Notethat argon was replaced by the atom with the same mass

545

number—calcium—in our simulation, since JENDL/HE
does not yet include the data for argon. The Earth wasrepresented as a sphere with the radius of 6378.14 km, andits composition was assumed to be 59.2% oxygen, 28.0%
silicon, 10.6% aluminum and 2.2% hydrogen by mass. This
constitution corresponds to 60% SiO2, 20% Al2O3 and 20%
H2O by mass. The particles arriving at 1000 g/cm2 below
the ground level were discarded in the simulation for reducing the computational time, since there are few albedo
neutrons from so deep underground to the atmosphere.
Comparison with Experiments
Figure 1 shows the comparisons of the calculated neutron
spectra with the corresponding experimental data obtained
by Goldhagen et al. (19) and Nakamura et al. (20). Notethat all the neutron spectra are expressed in terms of perunit lethargy throughout this paper. The statistical errors in
the values obtained by the simulation are generally small—less than approximately 5% and 20% for the high-altitude
and ground-level data, respectively, except for very highand low energies. The results at the altitudes below 2 m(⬃1035 g/cm2) are depicted in the figure at the ground levelto consider the effect of the Earth’s albedo neutrons precisely. The spectra predicted by the analytical functions
proposed in the next section are also plotted in the figures.
Two peaks around 1 MeV and 100 MeV can be observedin every spectrum. The former is attributed to neutrons
emitted by the evaporation process, while the latter is attributed to those produced by the pre-equilibrium and intranuclear cascade processes; these are called the evaporation and the pre-equilibrium peak, respectively, in this paper. The peaks at the thermal energy can be found only inthe spectra at the ground level, since they are composed
predominantly of the Earth’s albedo neutrons.
It is evident from the figure that our simulation can reproduce the experimental data for a variety of global conditions very well. We therefore concluded that the cosmicray neutron spectra estimated by our simulation are precise
enough to be used in systematic analysis for developing
analytical functions. It should be noted that slight overestimations by the simulation can be observed in the evaporation peaks at high altitudes. This discrepancy is due toour ignorance of the shielding effect of aircraft in the spectra; the details of the effect will be discussed in the next
section.
ANALYTICAL FUNCTION

Global Condition Dependence of Cosmic-Ray Neutron
Spectrum
To analyze the dependences of the cosmic-ray neutron
spectra on the global conditions, the simulations of the atmospheric propagation in semi-infinite atmosphere were
performed to obtain reference spectra without the disturbance of the local geometry effects. The semi-infinite at-

546

SATO AND NIITA

FIG. 1. Calculated and experimental neutron spectra in the atmosphere for various global conditions. Panels A and B show the comparison withthe data measured by Goldhagen et al. (19) and Nakamura et al. (20), respectively. The values of d and rc are the atmospheric depth and the cut-off
rigidity, respectively, while smin and smax indicate the solar minimum and maximum, respectively. These denotations are also used in the other figures.
Equations (10) and (18) were adopted in the calculation of the analytical values for the ground level and high altitudes, respectively.

mosphere consists of the dry air only, even at subterranean
regions down to 1000 g/cm2 below sea level. The simulations were carried out for the solar minimum and solar
maximum with modulation potentials of approximately 465

and 1700 MV, respectively, and 15 geomagnetic fields withthe vertical cut-off rigidities from 0.1 to 14 GV.
Figure 2 provides graphical presentations of the global
condition dependences of the neutron spectra in the semi-

FIG. 2. Neutron flux ratios normalized at (panel A) sea level and (panel B) the solar minimum with the cut-off rigidity 1 GV, obtained by the Monte
Carlo simulation for the semi-infinite atmosphere. Panel A shows the atmospheric-depth dependence of the neutron spectra, fixing the cut-off rigidity
and the solar condition, while panel B depicts spectral dependence on the cut-off rigidity and the solar condition at a conventional commercial flight
altitude. The gray lines in panel A indicate the spectra for the altitudes over the Pfotzer maximum.

547

PREDICTION OF COSMIC-RAY NEUTRON SPECTRA

trons fluxes is established, the differential equation of ⌽Lwith respect to the depth can be expressed by⳵⌽L (s, rc , d)
⫽ Np (s, rc , d) ⫺ ⌺A⌽L (s, rc , d),⳵d(3)

where NP and ⌺A⌽L, respectively, are the number of produced and absorbed neutrons with energies below 15 MeV
in the unit volume, with ⌺A denoting the macroscopic absorption cross-section in the atmosphere. The parameter NP
is significantly related to the high-energy particle fluxes,
and we assumed that the dependence of NP on the depth
can be represented by an exponential decay curve. Basedon the assumption, ⌽L can be expressed by⌽L (s, rc , d) ⫽ a1 (s, rc ){exp[⫺a2 (s, rc )d]⫺ a3 (s, rc )exp[⫺a4 (s, rc )d]},FIG. 3. Integrated fluxes of neutrons with energies below 15 MeV ⌽L
under the solar minimum and with several cut-off rigidities.

infinite atmosphere. It is found from the figures that theshapes of the neutron spectra below 15 MeV are almost
independent of the global conditions, except for those at
the altitudes over Pfotzer maximum (⬃18 km) as denotedby gray lines in Fig. 2A. This is because equilibrium between production and deceleration of neutrons is established in the lower energy regions. We therefore express
the neutron spectra in the semi-infinite atmosphere ␾Inf as
the product of the basic spectrum ␾B and the integrated
neutron flux below 15 MeV denoted as ⌽L as follows:␾Inf(s, rc, d, E) ⫽ ␾B(s, rc, d, E)⌽L(s, rc, d),(2)where s, rc, d and E denote the solar modulation potential,
vertical cut-off rigidity, atmospheric depth and neutron energy, respectively; we analyze the dependences of ⌽L and␾B on the parameters separately. The advantage of introducing ␾B is that the spectrum below 15 MeV is almost
independent of the global conditions for altitudes under 20km, and it can be expressed by a relatively simple equation,
as described later in this section. Since commercial flights
never exceed the altitude of 20 km, we will discard theneutron spectra above the altitude in the following analysis
for the purpose of simplification and propose analytical
functions for predicting the neutron spectra at altitudes under 20 km.
Low-Energy Neutron Flux
The low-energy neutron fluxes ⌽L obtained by the simulations under the solar minimum with several cut-off rigidities are shown in Fig. 3 as a function of the atmosphericdepth. As expected from previous studies, ⌽L reaches its
maximum at around 100 g/cm2 and then decreases with theincrease in the atmospheric depth. Assuming that the equilibrium between incoming and outgoing low-energy neu-

(4)

where a1 to a4 are the parameters depending on the modulation potential and cut-off rigidity. In the derivation of
the equation, we also assumed that no neutron exists at the
top of the atmosphere, i.e. ⌽L(s, rc, 0) ⫽ 0. If the all assumptions were completely true, then a3 and a4 should beequal to 1 and ⌺A, respectively. However, the validity ofthese assumptions is not always established, and hence a3
and a4 are also regarded as the global condition-dependent
parameters in our model. The values of the parameters afor each condition were determined by least-squares fittingof the ⌽L data.
Figure 4 plots the parameters obtained as a function of
the cut-off rigidity and solar modulation. It is seen fromthe figure that the differentials of the parameters with respect to the cut-off rigidity generally become the maximum
at rc around 4 GV. The solar modulation affects the parameters only at lower cut-off rigidities except for a1, since theincident particle spectra calculated by CREME96 are almost independent of the solar condition at higher cut-off
rigidities. We therefore introduced a sigmoid function to
generate the best-fitting curve to the parameters:ai (s, rc ) ⫽ bi1 (s) ⫹bi2 (s),1 ⫹ exp[(rc ⫺ bi3 )/bi 4 ](5)where bi1 to bi 4 are free parameters, and only b11 and bi 2
are related to the solar modulation; i.e., b21, b31, b41, bi 3 andbi 4 are the constant. It should also be noted that the sigmoid
function was also adopted in the work of Lewis et al. (21)for expressing the dependence of aircrew doses on the cutoff rigidity. Table 1 summarizes the numerical values of theparameters b obtained from the least-squares fitting of thea data shown in Fig. 4. We assumed that the dependencesof b11 and bi 2 on the solar modulation s can be expressed
as a linear relationship:b(s) ⫽b(smin )(smax ⫺ s) ⫹ b(smax )(s ⫺ smin ),smax ⫺ smin

(6)

where smin and smax denote the minimum and maximum ofs, 465 and 1700 MV, respectively. The validity of the lin-

548

SATO AND NIITATABLE 1
Numerical Values of the Parameters b used in Eq.
(5) Obtained from the Least-Squares Fitting of thea Data Shown in Fig. 4
Unit of ai⫺2a1 (cm⫺1s )a2 (cm2 g⫺1)

bi1

bi2

bi3 (GV)

bi4 (GV)

13.9(12.9)
0.00706

25.5(15.7)
6.73 ⫻ 10⫺4(5.70 ⫻ 10⫺4)
0.292(⫺0.210)0.00582(0.00441)

5.62

1.79

5.99

1.94

0.99

2.24

2.24

2.66

a3

0.975

a4 (cm2 g⫺1)0.00840

Notes. The values in parentheses are the data for the solar maximum,
while the others are those for the solar minimum or the solar-independent
parameters. The units of bi1 and bi2 are the same as those of ai.

¯ B that can be obtained by averaging ␾B for all thetrum ␾
cases—18 altitude ranges below 20 km, two solar conditions, and 15 geomagnetic fields. The graphical presentationof ␾B is given in Fig. 5. To reproduce ␾B by an analytical
function, we propose the following equation composed ofthree terms that represent the evaporation and pre-equilibrium peaks, and the pedestal part, respectively:冢 冣 exp 冢 c 冣 ⫹ c exp 冦 2[log (c )]E⫹ c log 冢 冣冦1 ⫹ tanh c log 冢 冣 冧[ cE ]c⫻ 冦1 ⫺ tanh c log 冢 冣 冧 ,[ cE ]␾¯ B (E) ⫽ c1

E
c2
7

⫺E

c3

冧

⫺[log 10 (E) ⫺ log 10 (c5 )] 2

4

2

10

10

9

8

6

2

10

10

11

10

(7)

12

FIG. 4. The dependences of parameters a used in Eq. (4) on the cutoff rigidity. The solid and open circles are the values obtained by theleast-squares fitting of the Monte Carlo simulation results for the solar
minimum and maximum, respectively, while the solid lines denote thevalues calculated by the sigmoid function of Eq. (5).

where c1 to c12 are constant parameters. The parameter c3
is fixed at 0.5 for the conventional evaporation spectrum,
but it is regarded as a free parameter in our model since
the peak has been disturbed by the deceleration and absorption before the equilibrium is accomplished. Note that

ear-relationship assumption could not be verified, since weperformed the Monte Carlo simulation only for the solar
minimum and maximum conditions. However, the differences in ⌽L between the two conditions are at most 50%,
and the error associated with introducing the assumption
can be considered small.
The calculated a parameters using Eq. (5) are also depicted in Fig. 4. The agreements between the data are quitesatisfactory for all the cases. The low-energy neutron flux
⌽L can be estimated analytically by substituting the parameters a obtained from Eq. (5) into Eq. (4); the results arealso shown in Fig. 3. It is evident from this graph that theanalytical function can reproduce the Monte Carlo resultsvery well.
Mean Basic Spectrum
Although the basic spectrum ␾B depends on the global
conditions, it is worthwhile to analyze the mean basic spec-

FIG. 5. The mean basic spectrum ␾¯ B. The gray line indicates the averaged spectrum for the Monte Carlo results of ␾B for all the global
conditions except for altitudes over 20 km, while the dashed line denotes
the values calculated by Eq. (7).

PREDICTION OF COSMIC-RAY NEUTRON SPECTRA

549

TABLE 2
Numerical Values of the Parameters c used in Eq.
(7) Obtained from the Least-Squares Fitting of the␾¯ B Data Shown in Fig. 5
Unit of ci⫺1c1 (lethargy )c2 (MeV)c3c4 (lethargy⫺1)ac5 (MeV)c6 (MeV)c7 (lethargy⫺1)c8 (MeV)c9c10 (MeV)
c11
c12 (MeV)a
a

ci
0.229
2.31
0.721
0.0516126
2.17
0.00108
3.33 ⫻ 10⫺12
1.62
9.59 ⫻ 10⫺8
1.48
299

Dependent on the global conditions assumed for estimating ␾B.

the second and third terms in the equation were introduced
only for reproducing the simulation result, and little physical meaning exists in their forms.
The numerical values of the parameters c obtained fromthe least-squares fittings of the ␾B data shown in Fig. 5 arelisted in Table 2. The ␾B calculated from Eq. (7) is depicted
in Fig. 5. The figure clearly indicates that the analytical
function can reproduce the simulation results very well.
Global Condition Dependence of Basic Spectrum
As mentioned before, the basic spectrum ␾B depends onthe global conditions at energies above 15 MeV. It can be
deduced from Fig. 2 that the global conditions affect ␾Bwith regard to the height of the pre-equilibrium peak andthe high-energy edge of the spectra. Figure 2B also showsthat the influence of the solar condition on ␾B is limited to
the lower cut-off rigidity cases. We therefore surmised thatc4 and c12, which are respectively related to the magnitudeof the pre-equilibrium peak and the high-energy edge ofthe pedestal part, are dependent on both atmospheric depth
and cut-off rigidity, and the others are the constant for all
the conditions.
Figure 6 plots the values of c4 and c12 obtained from theleast-squares fitting of ␾B after fixing the other c parameters
for several cut-off rigidities as a function of the atmosphericdepth. In general, c4 rises with an increase in the atmospheric depth and a decrease in the cut-off rigidity, whereas
the relation is reversed for c12. For the purpose of estimating
the depth dependences of c4 and c12, we introduced the following functions:c4(rc, d) ⫽ a5(rc) ⫹a6d
and1 ⫹ a7 exp(a8d)FIG. 6. The parameters c4 and c12 used in Eq. (7) for representing the
global-condition dependence of the neutron spectra in the semi-infinite
atmosphere.

of the parameters can be obtained by the least-squares fitting of the parameters c4 and c12. The dependences of a5,
a9, a10 and a11 on the cut-off rigidity are illustrated in Fig.
7. For best-fit curves to reproduce the cut-off rigidity dependences of the parameters, we also employed the sigmoid
function defined in Eq. (5). Table 3 summarizes the numerical values of the parameters b obtained from the leastsquares fitting of the a data shown in Fig. 7, together withthe constant values of a6, a7, a8 and a12. The calculated
results from the analytical functions adopting these parameters are also shown in Figs. 6 and 7. It is clear from thediagrams that the analytical functions can reproduce thesimulation result fairly well.
The basic spectrum ␾B for each condition can be estimated using Eq. (7), using the c4 and c12 values obtained
from Eqs. (8) and (9), respectively, instead of the data listed
in Table 2. The neutron spectra in the semi-infinite atmosphere ␾Inf can be predicted by multiplying ␾B by the lowenergy neutron flux ⌽L obtained from Eq. (4), as indicatedby Eq. (2).
Comparison with Simulation Results(8)c12(rc, d) ⫽ a9(rc){exp[⫺a10(rc)d] ⫹ a11(rc)exp(⫺a12d)},(9)where a5 to a12 are free parameters, and only a5, a9, a10 anda11 are related to the cut-off rigidity. The numerical values

Figure 8 provides graphical presentations of the comparisons between the neutron spectra in the semi-infinite atmosphere obtained by the simulation and by the analytical
function. It is evident from the graphs that the results estimated by the analytical function agree with those by thesimulation very well for a variety of global conditions.

550

SATO AND NIITA

influenced by the local structures such as the ground conditions, buildings and aircraft. The water density around thepoint of interest is one of the key quantities in the estimation of spectra at the ground level (4), while the configuration of aircraft is of great importance at higher altitudes(10). We therefore analyzed the dependences of the neutron
spectra on the water density and the aircraft structure as
representative of the local geometry effects.
Estimation of Neutron Spectra at Ground Level

FIG. 7. The dependences of the parameters a used in Eqs. (8) and (9)on the cut-off rigidity.

Equivalent agreements can be observed in the comparisonsin other cases, and this enables us to conclude that theanalytical function is capable of reproducing the simulation
results of the neutron spectra in the semi-infinite atmosphere for any global condition for altitudes under 20 km.
Local Geometry EffectIn practice, the cosmic-ray neutron spectra are different
from those in the semi-infinite atmosphere, since they are

To analyze the water-density dependence of the neutron
spectra at the ground level, simulations were performed thatwere similar to those described in the Monte Carlo Simulation section, but with the weight fraction of water, w, inthe subterranean region changed from 0 (dry ground) to 1
(pure water). In the simulation, the composition of the restof the ground is assumed to be 75% SiO2 and 25% Al2O3
by mass; e.g., the ground consists of 60% SiO2, 20% Al2O3
and 20% H2O by mass for w ⫽ 0.2, which corresponds tothe condition described in the Monte Carlo Simulation section. The neutron spectra for a particular global condition—
altitudes below 2 m from the ground level, the solar maximum, and 10 GV cut-off rigidity—were estimated for eachground condition.
Figures 9 and 10 show the absolute and relative valuesof the neutron spectra compared to the corresponding spectrum in the semi-infinite atmosphere. The figures clearly
indicate the significant dependence of the neutron spectra
on w especially for lower energies. The low-energy neutron
fluxes become larger with the existence of dry ground,since it can be regarded as a reflector of the neutrons. On
the other hand, the fluxes decrease with an increase of w,
since water works as an absorber of the neutrons. The
ground also causes the thermal neutron peak, which cannot
be observed in the neutron spectra in the semi-infinite atmosphere.
For the purpose of reproducing the neutron spectra at the
ground level ␾G by an analytical function, we propose thefollowing equation:
␾G (s, rc , d, E, w)⫽ ⌽L (s, rc , d)[␾B (s, rc , d, E) f G (E, w) ⫹ ␾T (E, w)],

TABLE 3
Numerical Values of the Parameters a and b used in Eq. (8) and (9) Obtained from theLeast-Squares Fittings
Unit of ai

Functionai or bi1

a5 (lethargy⫺1)
a6 (cm2 g⫺1 lethargy⫺1)a7
a8 (cm2 g⫺1)
a9 (MeV)a10 (cm2 g⫺1)a11
a12 (cm2 g⫺1)

Eq. (5)constant
constant
constant
Eq. (5)Eq. (5)Eq. (5)constant⫺0.00701
1.71 ⫻ 10⫺40.530
0.00136642
0.001121.26
0.0133Note. The units of bi1 and bi2 are the same as those of ai.

bi2

bi3 (GV)

bi4 (GV)

0.0258

10.9

2.38

⫺189
1.81 ⫻ 10⫺4⫺0.958

2.32
8.84
3.18

0.897
0.587
1.47

(10)

PREDICTION OF COSMIC-RAY NEUTRON SPECTRA

551

FIG. 8. Calculated neutron spectra in the semi-infinite atmosphere ␾Inf. Panel A shows the altitude dependence of the neutron spectra, fixing the
cut-off rigidity and solar condition, while panel B depicts spectra dependence on the cut-off rigidity and the solar condition at a conventional commercial
flight altitude.

where fG(E, w) denotes the functions representing the disturbance of the spectrum due to the local geometry effect,and ␾T(E, w) expresses the thermal neutron peak. For neutron energies above 1 eV, ␾T becomes negligible and fG(E,w) corresponds to the value plotted in Fig. 10, i.e., the ratioof neutron flux with the local geometry effect to neutron
flux in the semi-infinite atmosphere ␾G/␾Inf. The reason fortaking fG and ␾T to be independent of the global conditions
is that the local geometry effect influences these values
predominantly in the low-energy region, where the shapeof the basic spectrum ␾B is independent of the conditions.
We introduced a best-fit curve to the common logarithmof fG(E, w) as generated by the equationlog 10 [ f G (E, w)] ⫽ g1 ⫹ g2 log 10

[ ][

冦

[ ]冧]

E
E
1 ⫺ tanh g4 log 10
g3 (w)
g5 (w)

where g1 to g5, are free parameters determined by the leastsquares fitting of the ratios of the neutron fluxes shown in
Fig. 10. Note that this equation was introduced only forreproducing the simulation results, and little physical meaning was hidden in the form. The numerical values of theparameter g are given in Table 4 except for those of g3 andg5. In the fitting, the flux ratios for neutron energies below
1 eV were ignored, since the contributions of thermal neutrons are separately expressed by ␾T(E, w). The ratios forenergies above 500 MeV were also excluded from the fitting because of their poor statistics.
The parameters g3 and g5, are related to w, and this dependence is presented graphically in Fig. 11. For the curves, (11)

FIG. 9. Calculated neutron spectra at the ground level in comparison
with the corresponding spectrum in the semi-infinite atmosphere.

FIG. 10. The ratios of the calculated neutron fluxes at the ground levelto the corresponding fluxes in the semi-infinite atmosphere.

552

SATO AND NIITATABLE 4
Numerical Values of the Parameters g used in Eqs.
(11) and (17) Obtained from the Least-Squares
Fittings of the Neutron Flux Ratios shown in Figs.
10 and 12, Respectively
Unit of gi

gi

gi1a

gi2b

g1
g2
g3 (MeV)g4
g5 (MeV)⫺0.0235⫺0.0129
Eq. (12)0.969
Eq. (13)⫺0.0379⫺0.0179
6.03 ⫻ 10⫺8
1.36
5.52

⫺0.121⫺0.0220
5.42 ⫻ 10⫺10
1.40
6.45

a
b

Location of a pilot.
Location of a cabin passenger.

best fitting g3 and g5 respectively, we introduced the functionh32
w ⫹ h33(12)g5(w) ⫽ h51 ⫹ h52w ⫹ h53w2,(13)

log 10 [g3 (w)] ⫽ h31
and the third-order polynomial

respectively, where the parameters h are the fitting constants. Table 5 summarizes the numerical values of the parameters h obtained from the least-squares fitting of theg(w) data shown in Fig. 11, together with those of the constant g parameters. The calculated results of the analytical
functions adopting these parameters are also depicted in
Figs. 10 and 11. The figures clearly indicate the suitabilityof the analytical functions for reproducing the simulation
results.
The spectrum of the thermal neutron peak can be expressed by冢 冣 冢 冣

E
⫺E
exp
,
ET
ET
2

␾T (E, w) ⫽ g6 (w)(14)

where g6 denotes the magnitude of the peak and ET represents the thermal energy. The values of g6 can be determined by the equation

冕1eV

g6 (w) ⫽0[␾G (s, rc , d, E, w)/⌽L (s, rc , d) ⫺ ␾B (s, rc , d, E) f G (E, w)] dE/E

冕

1eV

(E/ET ) 2 exp(⫺E/ET ) dE/E

FIG. 11. The dependences of the parameters g used in Eqs. (11) and(14) on the weight fraction of water in the ground.

as shown in Fig. 10. For a best-fit curve to g6, we introduce
the function
g6 (w) ⫽h6l ⫹ h62 exp(⫺h63 w),
1 ⫹ h64 exp(⫺h65 w)(16)

where the parameters h6l are fitting constants determined by
the least-squares fitting of the g6 data. The numerical valuesof h6l are also given in Table 5, while the graphical fitting
result is depicted in Fig. 11.
The neutron spectra at the ground level ␾G can be estimated from Eq. (10) by substituting the calculated valuesof the analytical functions for ⌽L, ␾B, fG and ␾T. The results
are presented graphically in Fig. 9. The graph indicates thatthe analytical functions can express the water-density dependence of the neutron spectra satisfactorily. It should be

0

(15)

to reproduce the total flux of thermal neutrons using theanalytical functions. The calculated results setting ET ⫽0.025 eV are also shown in Fig. 11. Note that the Monte
Carlo simulation results were employed for the ␾G data inthe calculation, while the values obtained from the analytical functions were adopted for ⌽L, ␾B and fG. It is clear
from this graph that g6 generally becomes larger with adecrease of w, probably due to the lower number of thermal
neutrons absorbed by water, but the relationship is reversed
at the lowest values of w because of a rapid increase of fG

TABLE 5
Numerical Values of the Parameters h used in Eqs.
(12), (13) and (15) Obtained from the LeastSquares Fitting of the g(w) Data Shown in Fig. 11
Unit of gi

hi1

hi2

hi3

hi4

hi5

g3 (MeV)g5 (MeV)g6 (lethargy⫺1)

⫺25.2
0.348
0.118

2.73
3.35
0.144

0.0715
⫺1.57
3.87

0.653

42.8

Notes. The units of h51, h52, h53, h61 and h62 are the same as those of
the corresponding g parameters. The value of g3 in unit of MeV can be
estimated from Eq. (12) with the given h3l values.

PREDICTION OF COSMIC-RAY NEUTRON SPECTRA

553

where l indicates the location of interest, i.e. of a pilot ora passenger, g1l to g5l are the fitting constants, and mA andmA340 denote the mass of the aircraft being considered andAirbus-340, respectively. Although the detailed dependenceof the spectrum on the aircraft size has not been well investigated, we simply assumed that the common logarithmof fAl was proportional to mA since the influence of theaircraft structure on the spectrum generally could be expected to become larger with the mass. The values of theparameters g were determined by the least-squares fittingof the data shown in Fig. 12 and by setting mA ⫽ mA340.
The obtained values were also listed in Table 4, while thecalculated results using the parameters are depicted in Fig.
12. It is evident from the graph that the analytical function
is capable of reproducing the simulation results well.
The neutron spectra in aircraft ␾Al can be simply estimated by the equation
FIG. 12. The ratios of the calculated neutron fluxes inside Airbus-340
at the positions of a pilot and a cabin passenger to those in the freeatmosphere. The numerical values denoted by the symbols were calculated from the data shown in figure 4 (A) of ref. (10).

noted that not only the weight fraction of water in groundbut also the structures of surroundings such as concrete walland flooring are relevant to the neutron spectra. Furthermore, ground generally consists of a variety of materials
with several layers, and hence the assumption of groundwith a uniform composition is over simplified for somecases. Therefore, the weight fraction of water used in thefunctions should be considered not as a purely physicalquantity but an adjustable index for representing the existence of various materials, especially for hydrogen, aroundthe point of interest.
Estimation of Neutron Spectra in Aircraft
The influence of aircraft on the neutron spectra was wellanalyzed by Ferrari et al. (10). In their work, the particle
transport simulations in the mathematical model of Airbus340 were performed with the FLUKA code, and the neutron
spectra at the several locations in the aircraft such as the
cockpit and the cabin were estimated. Figure 12 shows theratios of their calculated neutron spectra inside the aircraftat the positions of a pilot and a cabin passenger to those in
the free atmosphere, these data were taken from Figure 4(A) in their paper (10).
A comparison of Figs. 10 and 12 indicates that the aircraft disturbs the neutron spectra in a manner similar to theground except for the thermal energies. For the function fAl
that represents the disturbance of the spectrum due to aircraft, we introduced an equation similar to that describingfG:
log 10 [ f Al (E, mA )]
⫽

冢

[ ]冦

[

冢 冣]冧冣 , (17)MAEEg ⫹ g2l log 10
1 ⫺ tanh g4l log 10
mA340 1l
g3l
g5l

␾Al (s, rc , d, E, mA ) ⫽ ⌽L (s, rc , d)␾B (s, rc , d, E) f Al (E, mA ).

(18)

This equation has no term representing the thermal neutron
peak, which was considered separately in the estimation of
␾G as expressed by Eq. (10). This is justified by the factthat aircraft generally are not large enough to form a thermal neutron peak, as indicated from experimental data (19,20). It should be noted that not only the location and themass of aircraft but also detailed configuration characteristics, such as the number of passengers and the amount offuel, affect the neutron spectra, and further studies are desired to develop a more sophisticated model for describing
the influence of the aircraft structure on the spectra.
Comparison with Experimental Data
The neutron spectra calculated by Eqs. (10) and (18) for
the ground level and high altitudes, respectively, are shownin Fig. 1 for the corresponding global conditions to theexperimental data. The weight fraction of water was set at
0.2 for the calculation at the ground level. At high altitudes,
the calculations were performed for the pilot location, andthe masses of the aircrafts were assumed to be a quarter of
mA340, apparently since the airplanes used in their measurements—ER-2 and DC-8—were smaller than Airbus-340.
It is evident from the figures that the neutron spectra of
the experiments and the calculations are in excellent agreement. Furthermore, the analytical functions are substantially superior to the Monte Carlo simulation in reproducing
experimental data at lower energies, since they can consider
the local geometry effect on the spectra even at high altitudes.
The accuracy of the analytical functions was also verified
by comparing their calculation results with the experimental
data by Kowatari et al. (22), as shown in Fig. 13. These
researchers measured the neutron spectra at the ground level at different altitudes while climbing up Mt. Fuji in Japan.
The weight fraction of water w was also set at 0.2 for the
analytical calculation. It can be seen clearly from the graph

554

SATO AND NIITA

portance in aircrew dosimetry, can also be done by applying the Monte Carlo simulation technique established bythis work. The simulation technique is also capable of contributing to the estimation of the particle spectra in the atmosphere at solar-geomagnetic storm conditions, an area
requiring additional study. Application of this work to theestimation of neutron 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 and Dr. A. Endo, Dr. M.
Kowatari, Dr. T. Nunomiya, Dr. H. Yasuda and Dr. T. Nakamura for theiradvice on this work. We also wish to thank the technical staff of the
CCSE office of JAEA for their help in performing the Monte Carlo simulation.
FIG. 13. Neutron spectra calculated with Eq. (10) compared to themeasured data in the Mt. Fuji area at the altitudes of 20, 1020 and 2400m, where the corresponding atmospheric depths are approximately 1030,
915 and 770 g/cm2, respectively.

that the analytical function can reproduce the experimental
data very well. Note that the corresponding Monte Carlodata have not been generated, since it is very time-consuming to perform the simulations with different altitudes of
the ground level.
From these considerations, we concluded that our analytical functions can predict cosmic-ray neutron spectrawith enough accuracy for any global condition except for
altitudes over 20 km, taking the local geometry effects into
account.
CONCLUSION

Monte Carlo simulations were performed for estimating
the cosmic-ray neutron spectra in the atmosphere by using
the PHITS code adopting JENDL/HE. Excellent agreement
was observed between the calculated and measured spectra
even at the ground level. Based on the comprehensive analysis of the simulation results, we proposed analytical functions that can predict the cosmic-ray neutron spectra forany time and location at the altitudes under 20 km, considering the local geometry effects. The accuracy of the analytical functions was clearly verified by comparing theirresults with various experimental data.
In the future, the systematic analysis of doses estimatedby our analytical functions is needed for their validation inpractical use. Further studies are desired to improve themodel for describing the local geometry effect on the neutron spectra, since that proposed in this paper incorporates
several simplifications. More experimental data for neutron
spectra at locations with detailed information on the surrounding structures are required for the purpose, especiallyfor those inside aircraft. Estimation of the photon andcharged-particle spectra at high altitudes, also of great im-

Received: February 28, 2006; accepted: May 11, 2006

REFERENCES
1. J. F. Ziegler, Terrestrial cosmic ray intensities. IBM J. Res. Dev. 42,
117–139 (1998).
2. 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).
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. Dosim. 86, 323–327 (1999).
4. 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).
5. A. Ferrari, M. Pelliccioni and T. Rancati, Calculation of the radiation
environment caused by galactic cosmic rays for determining air crew
exposure. Radiat. Prot. Dosim. 93, 101–114 (2001).
6. J. M. Clem, G. De Angelis, P. Goldhagen and J. W. Wilson, New
calculations of the atmospheric cosmic radiation field–results for neutron spectra. Radiat. Prot. Dosim. 110, 423–428 (2004).
7. 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.
8. 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 and P. Vaz, Eds.), pp. 955–960. Springer-Verlag, Berlin, 2001.
9. H. Schraube, V. Mares, S. Roesler and W. Heinrich, Experimental
verification and calculation of aviation route doses. Radiat. Prot. Dosim. 86, 309–315 (1999).
10. A. Ferrari, M. Pelliccioni and R. Villari, Evaluation of the influenceof aircraft shielding on the aircrew exposure through an aircraft mathematical model. Radiat. Prot. Dosim. 108, 91–105 (2004).
11. 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).
12. 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).
13. Y. Watanabe, T. Fukahori, K. Kosako, N. Shigyo, T. Murata, N. Yamano, T. Hino, K. Maki, H. Nakashima and S. Chiba, Nuclear data
evaluations for JENDL high-energy file. In Proceedings of Interna-

PREDICTION OF COSMIC-RAY NEUTRON SPECTRA

tional Conference on Nuclear Data for Science and Technology, pp.
326–331. American Institute of Physics, New York, 2005.
14. A. J. Tylka, J. H. Adams, Jr., P. R. Boberg, B. Brownstein, W. F.
Dietrich, E. O. Flueckiger, E. L. Petersen, M. A. Shea, D. F. Smartand E. C. Smith, CREME96: A revision of the cosmic ray effects on
micro-electronics code. IEEE Trans. Nucl. Sci. 44, 2150–2160(1997).
15. K. Niita, S. Meigo, H. Takada and Y. Ikeda, High Energy Particle
Transport Code NMTC/JAM. JAERI-Data/Code 2001-007, Japan
Atomic Energy Research Institute, Ibaraki, 2001.
16. J. F. Briesmeister, MCNP—A General Monte Carlo n-Particle Transport Code. LA-12625-M, Los Alamos National Laboratory, Los Alamos, NM, 1997.
17. K. Niita, S. Chiba, T. Maruyama, H. Takada, T. Fukahori, Y. Nakahara and A. Iwamoto, Analysis of the (N,xN⬘) reactions by quantum
molecular dynamics plus statistical decay model. Phys. Rev. C52,
2620–2635 (1995).

555

18. H. W. Bertini, Low-energy intranuclear cascade calculation. Phys.
Rev. 131, 1801–1821 (1963).
19. P. Goldhagen, J. M. Clem and J. W. Wilson, The energy spectrum of
cosmic-ray induced neutrons measured on an airplane over a wide
rage of altitude and latitude. Radiat. Prot. Dosim. 110, 387–392(2004).
20. T. Nakamura, T. Nunomiya, S. Abe, K. Terunuma and H. Suzuki,
Sequential measurements of cosmic-ray neutron spectrum and dose
rate at sea level in Sendai, Japan. J. Nucl. Sci. Technol. 42, 843–853(2005).
21. B. J. Lewis, L. G. I. Bennett, A. R. Green, M. J. McCall, B. Ellaschuk, A. Butler and M. Pierre, Galactic and solar radiation exposure
to aircrew during a solar cycle. Radiat. Prot. Dosim. 102, 207–227(2002).
22. M. Kowatari, K. Nagaoka, S. Satoh, Y. Ohta, J. Abukawa, S. Tachimori and T. Nakamura, Evaluation of the altitude variation of thecosmic-ray induced environmental neutrons in the Mt. Fuji area. J.
Nucl. Sci. Technol. 42, 495–502 (2005).