- Email: [email protected]
Monte Carlo simulation of cosmogenic nuclide production
Nuclear Instruments and Methods North-Holland, Amsterdam
in Physics
MONTE CARLO SIMULATION J. MASARIK,
P. EMRICH
Research
OF COSMOGENIC
*, P. POVINEC
Comenius University, Faculty of Mathematics
483
B17 (1986) 483-489
NUCLIDE
and S. TOUR
and Physics, Department
PRODUCTION
**
of Nuclear Physics, 842 15 Bmtislava,
CSSR
A method for calculation of absolute production rates of cosmogenic nuclides, based on the Monte Carlo simulation of hadron cascade process in matter is presented. The experimental depth profiles of 22Na, 26A1 and 53Mn in lunar samples are compared with results of computations. A disagreement between theoretical and experimental results for 26A1 and 53Mn has been found.
1. Introduction The study of the interactions of cosmic ray particles with objects of extraterrestrial origin in meteorites, lunar samples and cosmic dust have supplied important information on time and space variations in cosmic ray intensity and spectra, as well as on the exposure histories of these objects. Information on the history of cosmic rays is obtained from the analysis of the content of cosmogenic stable and radioactive nuclides in extraterrestrial objects. For the evaluation of concentrations of cosmogenic nuclides in these objects it is necessary to know or to estimate their production rates. Knowing the flux, the composition and the energy spectra of cosmic ray particles of galactic and solar origin, the composition of target nuclei and cross sections of nuclear reactions, it is possible to calculate production rates of cosmogenic nuclides. By comparing the measured activities of radionuclides of different half-lives with independently determined production rates it is possible to estimate the average cosmic ray intensity over different periods of time. Various methods for calculation of cosmogenic nuclide production rates have been developed [l-4]. In this article we describe our Monte Carlo (MC) approach and we compare theoretical and experimental production rates of 22Na, 26A1 and 53Mn in lunar samples.
2. Physical processes and procedures used in calculations In our calculations incoming primary particles in the energy range from one to tens of GeV are considered. These particles undergo inelastic collisions with target
* Present address: Ministry of Education SSR, Bratislava, CSSR. ** Present address: Joint Institute for Nuclear Research, Dubna,
USSR.
0168-583X/86/$03.50 0 Elsevier Science Publishers (North-Holland Physics Publishing Division)
B.V.
nuclei. The characteristic feature of collisions at these energies is the origin of new hadrons whose average number is approximately proportional to the logarithm of the primary energy. They have enough energy to undergo further inelastic collisions and to produce the next generation of secondary particles. After production of several generations of particles the cascade process comes to an end when the energies are too small for further particle production, The energy losses due to elastic scattering are negligible in comparison with the ionization and nuclei excitation energy losses. High energy collisions in which secondary particles are produced and also collisions at lower energies transform the target nuclei mainly through spallation and evaporation processes to radioactive or stable nuclei. Amongst strongly interacting hadrons it is sufficient to consider only protons, neutrons and pions, because only these appear in significant numbers. Electron-photon showers are initiated in great numbers by decay of neutral Q’ mesons. In most materials the radiation length is much shorter than the absorption length for strongly interacting particles and hence the electromagnetic cascades play only a minor role in the development of the hadronic cascade. They are, however, very important for energy deposition. The basic problem of hadron transport in solid matter is the solution of the transport equation. The transport equation of hadron cascade in solid matter represents a type of Boltzman kinetic equation [5]. This equation can be solved analytically or stochastically. In this paper we only deal with the stochastical approach. In our problem we consider two basic types of particles. Firstly hadrons, which form the hadron cascade, and secondly the target nuclei, which are steady elements with well known distribution functions. In the arbitrary unit of phase space we consider all processes yielding the generation or absorption of particles. In the case of the hadron cascade it is enough to consider the following processes [6]: a) inelastic nuclear interactions, IV. COSMOGENIC NUCLIDES
2. Musarik ef a/. / Simulation for cosmogenic m&de production
484
b) elastic scattering, c) ionization and excitation of target atoms, d) decay of particles. The solution of the problem of hadron cascade in matter using a Monte Carlo (MC) method is in the final analysis reduced to the calculation of following integrals (1)
Rij=(X,
Wl(k’(xo,
y=(s,T,r),
x’=(u),
E, 0’)
and
Qi(x)=Qi(E,
w)S(r,--r)
is the primary source of i-type particles with energy E and the direction w, r, is the position vector of the surface at which the primary particle enters the target, cj(E, T, o, 7) is the number of i-type particles with energy E and direction w produced in inelastic collisions of j-type particles with energy T and the direction T with target nuclei. The calculation of the functional (1) using the MC method is based on the estimation of the average random variable 2 (an estimator) defined on the region SE of all possible realizations of Markovian process [7,8] which represents the process of the hadron cascade. Therefore we can write
Xwjk)(xo, x1,
E R&x)
i
W:k)(xO, . ..xk)h(xj).
~~52~ = 0
Xk).
(4)
The choice of ulk) is limited only by the variance condition (2). A comparison of the explicit expressions for E(Z) and (h, $1 leads to the following choice of the w$“’ 181
. ..dx.=O,
(5)
for i = 1, 2, 3, . . . . The estimator can be represented for example by the function Z, defined on the region 9, which we obtain from (3) when we introduce wj’) = 1 (for i, k = 0, 1, 2...). It is seen from eq. (4) that the weight function W>k) connects the real transport process represented by the transport eq. (1) with the chosen Markovian process. Calculation of the functional (2) is organized as follows. On the basis of the distributions p*(x), p(x, y) and of chosen Markovian process we generate random trajectories a: = (x$/j.. . xij)) and then we calculate the quantity (61 for each ej. From the quantity Zi we go to an estimation of the value E(Z) = (h, 111) by generation of the quantity
This quantity converges for large N to E(Z) because 2, = Z(ej) and the trajectories aj are generated in agreement with the measure of the probability P.
3. The structure of the program The method of calculating production rates of cosmogenic nuclides is based on the following expression
(3)
i=O
where Rn,k(cr)=l
. . .
(2)
where (Yis the random realization of Markovian process (a trajectory) and P is the probability distribution defined on s2. The introduction of Markovian process and a system (L?, P) enables to express the basic idea of a stochastic calculation contained in eq. (2). The estimator Z, which is a real function defined on 52, must represent a true estimate of (h, J/). In general we can express Z as a function of the trajectory ar as
k=O
$o;;;.*‘Kk-lsk , pk-l, k
$,gI%.
(h.~)=E(Z)=SZ(a)dP(a),
Z(o)=
._.Xk)==
’ ‘.
xdxi+l E, o),
Xl,
Wik)gi+ ,=F+r (R,_,Wlk’Ri, i-t . . ’ ‘k-1, k
Y)=/Gi(x, Y)F,i(E,T, 0, r)T
where x=(r,
choose in the form
LYE&
is the characteristic function on &, and Wjk)(x,, xk) is a function that it is advantageous to Xl,“’
X
IVNd3rqk(r, E’, E),
(7)
where Q, is the production rate in a chosen region of space V,, N(E) is the energetic spectrum of primary particles, u$( E’) is the cross section for the production of x-type nuclide by k-type incident particle with the energy E’, o$“‘(E’) is the inelastic collision cross sec-
J. Masarik et al. / Simulation for cosmogenic nuclide production
tion of k-type particle with energy E’, $(r, E’, E) is the inelastic collision density for the k-type particle with the energy E’ at the position r. The calculation is performed by program units CASCA, LOWENDAT, COMDAT, NUPROD and SUMENW. The basic problem, the solution of the transport equation using the MC method is realized in the program CASCA. The output of this program is the quantity G+(E)
d3*dE&(r,
=/,
E’, E),
‘I
(8)
where F,j = (E,, E,,,). V,, which gives the number of inelastic collisions of k-type particle from the i th-energetic region in the space volume 5. LOWENDAT gives basic information about the target material, input particles and it also calculates the data required for the low energy particle simulation. In the program unit NUPROD the depth dependence and the energy dependence of the production rates of cosmogenic nuclides with consideration of contribution from the evaporation processes are calculated as p/(E)
= C #ij,e;(Ei)/ep’(E,). k.
(9)
i
The program COMDAT calculates the formation probabilities of cosmogenic nuclides as a function of the type and the energy of cascade particles and the energy of the incoming particles. The program SUMENW calculates the production rates of nuclides for a given input spectrum of primary cosmic rays. Communication between program units is realized throught the data libraries.
4. Simulation of the cascade process In this section we briefly describe the important part of the simulation process used in the program CASCA. We have shown that the calculation of the functional $J,~~ (eq. (8)) using the MC method is based on the estimation of the average quantity estimator which is defined on the region s2 of all possible trajectories of chosen Markovian process. Owing to this fact it is necessary to solve the following two problems: a) to construct the estimator for the inelastic collision number as a function defined on the space region 0. b) to construct a common Markovian process. We have used within our calculations two expressions for the considered estimators. Firstly, the estimator is defined as Z+(a)
= N,
where N=l, 2, 3-.. is the number sions of cascade particles of k-type cascade and represented by the tree The second type of estimator used in
(10) of inelastic collibelonging to the in the region Fij. our calculation is
defined
485
as
(11) m=l
where M is the number of cascade particles belonging to the cascade, represented by the tree 01, the range of which crosses the region r,j; w,$,‘)is the probality that the M th particle of type I, whose range crosses the region rii, undergoes an inelastic interaction in this region. From the definitions of the estimators follows that in the case of the estimator Z a particle contributes to the collision density only at the end of its path. In the case of the estimator Y the investigated particle contributes continuously to the collision density along its path. The simulation of Markovian process is started by a choice of primary particle characteristics. The particle history is followed taking into account its multiple Coulomb scattering up to the moment of its inelastic interaction, its stopping or leaving the investigated target region. In the case of inelastic collisions of investigated particles the type of target nucleus is determined and only one of the secondary particles is chosen for the next investigation. The succesive choice of always a higher generation of secondary particles is performed up to the case when the generated particle is characterized by one of the following conditions: a) the momentum of generated particles is less than per = 0.3 GeV/c, b) the particle momentum is in the interaction place less than the threshold momentum for production of the next generation of secondary particles, c) the particle leaves the target region. Except for case c) the investigated particle is followed within the low energy simulation. When all possible secondary particles are realized according to the amount of energy reserved for them, the program goes on to investigate the next generation of secondary particles. The sampling of the cascade process is stopped when all particles generated in interactions of the primary particle are realized. As we do not have full information about inelastic collisions we have used for their modelling procedures which are realized by satisfying the following conditions: a) in all collisions the law of energy conservation is fulfilled. b) average values of multiplicities, calculated for the large number of collisions, must be in agreement with the prescribed values. c) all secondary particles are generated according to the following equations [9]
IV. COSMOGENIC NUCLIDES
486
J. Masarik et al. / Simulation for cosmogenic nuclide production
cases for which experimental values of total cross sections were not available, we calculated them using the following inte~olation [lo]
- proton and neutron production in nucleon-nucleus collision dZN dp,, dp,
(16)
=- ~(l~~~,,+~p:,)p,[exp(-a,p:) where
+expt-a6p,>.a5],
_ pion
production in nucleon-nucleus
d=N
sl
dp,, dp, +a,
~1?l
(
al exp”-“2/E2)Pll)p,
exp(-osp,>I)/(
pi +pl
collision
+mZ)“’
03)
collision
‘32/a, A,/A,)
u&J = [email protected]
d2N
) ’
ai\i = up*/1 2.
In the case of known chemical composition of the target effective values of cross sections have been used
dp,, dp, XP,
led
0;’ E c
-+p?l.)
+ a5 exp(-~6pl]
I
04)
PII < 0,
Pll> 0,
- pion production
d2N =[
a5
in pion-nucleus exp( -4zpi)
collision +al]oip,
Xexd-wJ x (pi t-p:
Wn”kn)
n
E0
dpri dp,
log(
where u1 and o2 are cross sections on nuclei with mass numbers A, and A,, respectively (A, I; A s AZ). Nucleon and pion cross sections have been calculated using the relation
exp( -a~p2,) 1
- nucleon production in pion-nucleus
log( B =
where k is the type of the incident particle {proton, neutron, charged pion), w, is the relative statistical weight of n-type of target nuclei which gives the probability of its occurrence, ekn is the inelastic interaction cross section for a k-type particle with an n-type nucleus. The elastic scattering cross sections have been calculated using the formula Q(A)
+m2,y2,PII'
= 6.38A1”4[mb].
The differential angular distribution of secondary particles is well described by the expression
0,
due’ = exp( - 15.7~~8~A’/~), da
Xexd--a4p) x(p:,+p?
+mZ,y,
PII < 0, 09
where EC is the total collision energy in PCMS, pII, pI are the lon~tu~n~ and perpendicular components of the secondary particle momentum in PCMS, M, is the mass of pioa, a, . . ’ a9 are the parameters obtained from experimental data fitting. When the number of generated input particles or the number of realized collisions is equal to the required value, the inelastic collision densities as a function of primary particle momentum and space position are normalized to one input particle. 5. Auxiliary computations For modelling of the cascade process it is necessary to know the total inelastic cross sections for the interaction of protons, neutrons and charged pions with the target nuclei. In the region above 10 GeV/c we have supposed that the total cross section is constant. In
where A is the target nuclei mass number, p is the particle momentum in GeV/c, 8 is the scattering angle in radians. Elastic scattering is considered only for the perpendicular incoming particles. In other cases (secondary particles, an isotropic flux of primary particles) the production angles or input angles are much greater than the elastic scattering angles and therefore elastic scattering is neglected. The simulation of the inelastic interaction of particles with target nuclei is an important part of the simulation of the cascade development. To be able to generate the secondary particles initiated in the inelastic interaction it is necessary to know the momentum spectra of the secondary particles and their production angles as a function of the type and momentum of the incident particle and of the type of the target nucleus. In our calculations we have used experimentally obtained functions which give the secondary particle distribution as a function of the longitudial p,, and perpendicular p I component of their momenta in the hadron-nucleon centre of mass system. The multiplici-
J. Masarik et al. / Simulation for cosmogenic nuclide production ties of each type of secondary particles originating in hadron-nucleon collision have been calculated using the following equation
d’“: K,, = k/E.-
,
‘dpdQ
dp dQ,
>9
where Kij are the multiplicities without the correction on the nucleus excitation, E,, is the nucleus excitation energy. The choice of a random angle for the primary particle direction which has undergone multiple Coulomb scattering is realized using the Gaussian distribution with quadratic angle deviaton given by the Rossi formula [ll] (82)~~
=
=E
t
+ (E-E,)(B-E,) 3 - E,
= A/100
6. Results and discussion The MC method has been applied to production rate calculations of 22Na 26A1 and 53Mn in lunar samples. 22Na was chosen for testing and calibration reasons while 26A1 and 53Mn for investigation of the cosmic ray (CR) behaviour in the past. The calculations have been carried out for targets having the form of an infinite layer of material. The thickness of the layer was chosen to be 9000 kgme2. The effective chemical composition was taken to be in accordance with that of real lunar soils from Apollo 15 [12-141. For the purpose of investigation of production rate depth profiles the target was divided into 25 sublayers. An isotropic flux of CR protons is assumed to impinge on the target. The CR spectrum used in the calculations was taken from ref. [15]. The form of the spectrum takes into account the influence of solar activity on the galactic CR intensity through the modulation parameter $I. T(T+2E,,)(T+++m)-’
($)( $)“i[rad],
‘(T’
where pc is the momemtum of the considered particle in MeV, I,, is the radiation length of the considered material, 1 is the range between the input point and the position, in which we consider the angle of particle direction. Multiple Coulomb scattering has been considered for the secondary particles for the same reasons as for elastic scattering. All particles which take part in the hadron cascade lose their energy in these processes through: a) the ionization and excitation of atoms (ionization losses), b) the production of cascade particles (protons, neutrons, charged pions), c) the production of neutral pions, d) the excitation of nuclei. The ionization energy losses were calculated using the Bethe-Bloch formula with a correction for the effects of material density and of close binding. The excitation energy in GeV is approximated by the following experimental expression E ex
where E is the incident energy in the laboratory system, E, is the threshold energy of a particle able to develop a cascade (E, z 0.05 GeV), B is the maximum of both values A/100 and E,.
(17)
where E, is the kinetic energy of an incident particle in the laboratory system, Ej is the kinetic (in the case of pion - total) energy of j-type secondary particles in the laboratory system, d2N/dpdQ is the distribution function of j-type generated particle number as the function of secondary particle momentum p, D is the solid angle. As a part of the energy of the incoming particle is used for the target nucleus excitation, it is necessary to make corrections in the K,, multiplicities, eq. (17). We have used the correction K,J = KFJ (1 - E.JE,
487
E,cE<3GeV E>3GeV,
‘)=A
’
(T++)(T+2E,,++)
where T (MeV) is the proton kinetic energy, E, its rest energy, $I the modulation parameter (MeV), A = 9.9 X
L
0
I
100
200
300
I I 400 gcm2
Fig. 1. A comparison of calculated and measured depth dependence of production rate of ‘*Na in Apollo 15 samples. The experimental data have been taken from ref. (131. (a: + = 650, b: $I = 450, c: +J= 200 MeV). IV. COSMOGENIC NUCLIDES
J. Masarik et at. / ~irn~~ati~n for cosmogenic m&de pr~~ct~on
488
Bij
kg” Ffl r
4.0
3.0
2.0
1.0 1
0
I
100
I
Sbll
300
I
400 gcni2
Fig. 2. A comparison of calculated and measured depth dependence of production rate of 26A1in Apollo 15 samples. The experimental data have been taken from ref. [13,14]. (a: + = 100, b: #I = 200, c: rp= 450, d: 9 =i 650 MeV).
y = 2.65, for the 108; 172= 780 exp( - 2.5 X 10e4T); unmodulated GCR flux rp = 100 Me?‘, for the solar minimum GCR flux # = 300 MeV, for the solar maximum GCR flu 9 = 900 MeV and for the average GCR flux Cp= 650 MeV. In the calculation, the spectrum was corrected for the non-proton contribution, which can be estimated to be as high as 35-40%. Fig. 1 shows a comparison of calculated and measured depth dependences of the production rate of *‘Na in Apollo 15 samples. The experimental data have been taken from ref. [13]. The best compromise between the theoretical and experimental data was found for the modulation parameter + = 650 MeV. This represents an average galactic cosmic ray flux a few years before the sample collection. in the case of 26A1, which integrates the galactic cosmic ray flux over 1 Myr, fig. 2 shows a disagreement between the theoretical and experimental data. Although the theoretical curve follows the shape of the experimental dependence, at small depths there is an overproduction of 26A1.The best compromise between the theoretical and experimental data was found for the modulation parameter @= 200 MeV. Even larger disagreement has been found for the production rate of 53Mn which integrates the galactic cosmic ray flux over a few million years. Fig. 3 shows that an overproduction of s3Mn in comparison with the calculated production rate due to the average galactic cosmic ray intensity + = 650 MeV has been found. Although there are large uncertainties in the measured rates the and calculated values of 53Mn production observed
effect seems to overlap these errors. The shape
I
0
100
a
200
1
300
4oo .?
Fig. 3. A comparison of calculated and measured 53Mn production rate in Apollo 15 samples. The experimental data have been taken from ref. [14]. (a: + = 100, b: + = 450, c: QI= 650 MeV).
of the theoretical curve would imply that 53Mn was mostly produced by an unmodulated galactic cosmic ray flux #I = 100 MeV. Similar conclusions were obtained in our previous work where different methods of calculations were used [16,17].
References PI R.C. Reedy and J.R. Arnold, J. Geophys. Rev. 77 (1972) 537.
PI AX. Lavrukhina and G.K. Ustinova, Astron. J. 44 (1966) 1081. 131T.W. Armstrong and R.G. Alsmiller, Proc. Apollo 12 Lunar Sci. Conf., Geochim. Cosmochim. Acta, suppl 2 (1971) 1729. 141 ES. Light, M. Merker, R.B. Mendel and S.A. Korf, Acta Phys., suppl. 2 (1970) 745. PI N.P. Kalashnikov, VS. Remizovich and M.I. Riazanov, Stolknovenya bystrykh zaryazhenykh chastic v tverdykh telakh (Atomizdat, Moskva, 1980). PI R.G. Alsmiller, Nucl. Sci. Eng. 27 (1967) 158. [71 J. Spanier and E.M. Celbard, Monte Carlo Principles and Neutron Transport Problems (Addison-Wesley, Massachusetts, 1969). PI SM. Jermakov, Monte Carlo Metody i Smeshnye Problemy (Nauka, Moskva, 1975). 191 S.Tok&r, Ph.D. Thesis, Comer&s University, Bratislava (1983). 1101 J. Ranft and J.T. Routi, Comput. Phys. Commun. 7 (1972) 101. 1111 S. Rossi, Kosmicheskiye luchi (Atomizdat, Moskva, 1975). WI L.A. Rancitelli, J.S. Fruchter, D.M. Felix, R.W. Perkins
J. Masarik
489
et al. / Simulation for cosmogenic nuclide production
and N.A. Wogman, Proc. 6th Lunar. Sci. Conf. (NASA, Houston, 1975) p. 1876. [13] K. Nishiizumi, J. Klein, R. Middleton and J.R. Arnold, Earth and Planet. Sci. Lett. 70 (1984) 164. [14] M. Imamura, R.C. Finkel and M. Wahlen. Earth and Planet. Sci. Lett. 20 (1973) 107.
[15] G. Bonino. G. Castagnoli and M. Dardo. Cosmic Ray Conf. (IUPAP, Kyoto, 1979) vol. [16] P. Povinec, Izv. Akad. Nauk SSSR, ser. fiz. 41 [17] A.K. Lavrukhina, P. Povinec, G.K. Ustinova, skie issledovanya 22 (1984) 110.
IV. COSMOGENIC
Proc. 16th 12, p. 155. (1977) 383. Kosmiche-
NUCLIDES
in Physics
MONTE CARLO SIMULATION J. MASARIK,
P. EMRICH
Research
OF COSMOGENIC
*, P. POVINEC
Comenius University, Faculty of Mathematics
483
B17 (1986) 483-489
NUCLIDE
and S. TOUR
and Physics, Department
PRODUCTION
**
of Nuclear Physics, 842 15 Bmtislava,
CSSR
A method for calculation of absolute production rates of cosmogenic nuclides, based on the Monte Carlo simulation of hadron cascade process in matter is presented. The experimental depth profiles of 22Na, 26A1 and 53Mn in lunar samples are compared with results of computations. A disagreement between theoretical and experimental results for 26A1 and 53Mn has been found.
1. Introduction The study of the interactions of cosmic ray particles with objects of extraterrestrial origin in meteorites, lunar samples and cosmic dust have supplied important information on time and space variations in cosmic ray intensity and spectra, as well as on the exposure histories of these objects. Information on the history of cosmic rays is obtained from the analysis of the content of cosmogenic stable and radioactive nuclides in extraterrestrial objects. For the evaluation of concentrations of cosmogenic nuclides in these objects it is necessary to know or to estimate their production rates. Knowing the flux, the composition and the energy spectra of cosmic ray particles of galactic and solar origin, the composition of target nuclei and cross sections of nuclear reactions, it is possible to calculate production rates of cosmogenic nuclides. By comparing the measured activities of radionuclides of different half-lives with independently determined production rates it is possible to estimate the average cosmic ray intensity over different periods of time. Various methods for calculation of cosmogenic nuclide production rates have been developed [l-4]. In this article we describe our Monte Carlo (MC) approach and we compare theoretical and experimental production rates of 22Na, 26A1 and 53Mn in lunar samples.
2. Physical processes and procedures used in calculations In our calculations incoming primary particles in the energy range from one to tens of GeV are considered. These particles undergo inelastic collisions with target
* Present address: Ministry of Education SSR, Bratislava, CSSR. ** Present address: Joint Institute for Nuclear Research, Dubna,
USSR.
0168-583X/86/$03.50 0 Elsevier Science Publishers (North-Holland Physics Publishing Division)
B.V.
nuclei. The characteristic feature of collisions at these energies is the origin of new hadrons whose average number is approximately proportional to the logarithm of the primary energy. They have enough energy to undergo further inelastic collisions and to produce the next generation of secondary particles. After production of several generations of particles the cascade process comes to an end when the energies are too small for further particle production, The energy losses due to elastic scattering are negligible in comparison with the ionization and nuclei excitation energy losses. High energy collisions in which secondary particles are produced and also collisions at lower energies transform the target nuclei mainly through spallation and evaporation processes to radioactive or stable nuclei. Amongst strongly interacting hadrons it is sufficient to consider only protons, neutrons and pions, because only these appear in significant numbers. Electron-photon showers are initiated in great numbers by decay of neutral Q’ mesons. In most materials the radiation length is much shorter than the absorption length for strongly interacting particles and hence the electromagnetic cascades play only a minor role in the development of the hadronic cascade. They are, however, very important for energy deposition. The basic problem of hadron transport in solid matter is the solution of the transport equation. The transport equation of hadron cascade in solid matter represents a type of Boltzman kinetic equation [5]. This equation can be solved analytically or stochastically. In this paper we only deal with the stochastical approach. In our problem we consider two basic types of particles. Firstly hadrons, which form the hadron cascade, and secondly the target nuclei, which are steady elements with well known distribution functions. In the arbitrary unit of phase space we consider all processes yielding the generation or absorption of particles. In the case of the hadron cascade it is enough to consider the following processes [6]: a) inelastic nuclear interactions, IV. COSMOGENIC NUCLIDES
2. Musarik ef a/. / Simulation for cosmogenic m&de production
484
b) elastic scattering, c) ionization and excitation of target atoms, d) decay of particles. The solution of the problem of hadron cascade in matter using a Monte Carlo (MC) method is in the final analysis reduced to the calculation of following integrals (1)
Rij=(X,
Wl(k’(xo,
y=(s,T,r),
x’=(u),
E, 0’)
and
Qi(x)=Qi(E,
w)S(r,--r)
is the primary source of i-type particles with energy E and the direction w, r, is the position vector of the surface at which the primary particle enters the target, cj(E, T, o, 7) is the number of i-type particles with energy E and direction w produced in inelastic collisions of j-type particles with energy T and the direction T with target nuclei. The calculation of the functional (1) using the MC method is based on the estimation of the average random variable 2 (an estimator) defined on the region SE of all possible realizations of Markovian process [7,8] which represents the process of the hadron cascade. Therefore we can write
Xwjk)(xo, x1,
E R&x)
i
W:k)(xO, . ..xk)h(xj).
~~52~ = 0
Xk).
(4)
The choice of ulk) is limited only by the variance condition (2). A comparison of the explicit expressions for E(Z) and (h, $1 leads to the following choice of the w$“’ 181
. ..dx.=O,
(5)
for i = 1, 2, 3, . . . . The estimator can be represented for example by the function Z, defined on the region 9, which we obtain from (3) when we introduce wj’) = 1 (for i, k = 0, 1, 2...). It is seen from eq. (4) that the weight function W>k) connects the real transport process represented by the transport eq. (1) with the chosen Markovian process. Calculation of the functional (2) is organized as follows. On the basis of the distributions p*(x), p(x, y) and of chosen Markovian process we generate random trajectories a: = (x$/j.. . xij)) and then we calculate the quantity (61 for each ej. From the quantity Zi we go to an estimation of the value E(Z) = (h, 111) by generation of the quantity
This quantity converges for large N to E(Z) because 2, = Z(ej) and the trajectories aj are generated in agreement with the measure of the probability P.
3. The structure of the program The method of calculating production rates of cosmogenic nuclides is based on the following expression
(3)
i=O
where Rn,k(cr)=l
. . .
(2)
where (Yis the random realization of Markovian process (a trajectory) and P is the probability distribution defined on s2. The introduction of Markovian process and a system (L?, P) enables to express the basic idea of a stochastic calculation contained in eq. (2). The estimator Z, which is a real function defined on 52, must represent a true estimate of (h, J/). In general we can express Z as a function of the trajectory ar as
k=O
$o;;;.*‘Kk-lsk , pk-l, k
$,gI%.
(h.~)=E(Z)=SZ(a)dP(a),
Z(o)=
._.Xk)==
’ ‘.
xdxi+l E, o),
Xl,
Wik)gi+ ,=F+r (R,_,Wlk’Ri, i-t . . ’ ‘k-1, k
Y)=/Gi(x, Y)F,i(E,T, 0, r)T
where x=(r,
choose in the form
LYE&
is the characteristic function on &, and Wjk)(x,, xk) is a function that it is advantageous to Xl,“’
X
IVNd3rqk(r, E’, E),
(7)
where Q, is the production rate in a chosen region of space V,, N(E) is the energetic spectrum of primary particles, u$( E’) is the cross section for the production of x-type nuclide by k-type incident particle with the energy E’, o$“‘(E’) is the inelastic collision cross sec-
J. Masarik et al. / Simulation for cosmogenic nuclide production
tion of k-type particle with energy E’, $(r, E’, E) is the inelastic collision density for the k-type particle with the energy E’ at the position r. The calculation is performed by program units CASCA, LOWENDAT, COMDAT, NUPROD and SUMENW. The basic problem, the solution of the transport equation using the MC method is realized in the program CASCA. The output of this program is the quantity G+(E)
d3*dE&(r,
=/,
E’, E),
‘I
(8)
where F,j = (E,, E,,,). V,, which gives the number of inelastic collisions of k-type particle from the i th-energetic region in the space volume 5. LOWENDAT gives basic information about the target material, input particles and it also calculates the data required for the low energy particle simulation. In the program unit NUPROD the depth dependence and the energy dependence of the production rates of cosmogenic nuclides with consideration of contribution from the evaporation processes are calculated as p/(E)
= C #ij,e;(Ei)/ep’(E,). k.
(9)
i
The program COMDAT calculates the formation probabilities of cosmogenic nuclides as a function of the type and the energy of cascade particles and the energy of the incoming particles. The program SUMENW calculates the production rates of nuclides for a given input spectrum of primary cosmic rays. Communication between program units is realized throught the data libraries.
4. Simulation of the cascade process In this section we briefly describe the important part of the simulation process used in the program CASCA. We have shown that the calculation of the functional $J,~~ (eq. (8)) using the MC method is based on the estimation of the average quantity estimator which is defined on the region s2 of all possible trajectories of chosen Markovian process. Owing to this fact it is necessary to solve the following two problems: a) to construct the estimator for the inelastic collision number as a function defined on the space region 0. b) to construct a common Markovian process. We have used within our calculations two expressions for the considered estimators. Firstly, the estimator is defined as Z+(a)
= N,
where N=l, 2, 3-.. is the number sions of cascade particles of k-type cascade and represented by the tree The second type of estimator used in
(10) of inelastic collibelonging to the in the region Fij. our calculation is
defined
485
as
(11) m=l
where M is the number of cascade particles belonging to the cascade, represented by the tree 01, the range of which crosses the region r,j; w,$,‘)is the probality that the M th particle of type I, whose range crosses the region rii, undergoes an inelastic interaction in this region. From the definitions of the estimators follows that in the case of the estimator Z a particle contributes to the collision density only at the end of its path. In the case of the estimator Y the investigated particle contributes continuously to the collision density along its path. The simulation of Markovian process is started by a choice of primary particle characteristics. The particle history is followed taking into account its multiple Coulomb scattering up to the moment of its inelastic interaction, its stopping or leaving the investigated target region. In the case of inelastic collisions of investigated particles the type of target nucleus is determined and only one of the secondary particles is chosen for the next investigation. The succesive choice of always a higher generation of secondary particles is performed up to the case when the generated particle is characterized by one of the following conditions: a) the momentum of generated particles is less than per = 0.3 GeV/c, b) the particle momentum is in the interaction place less than the threshold momentum for production of the next generation of secondary particles, c) the particle leaves the target region. Except for case c) the investigated particle is followed within the low energy simulation. When all possible secondary particles are realized according to the amount of energy reserved for them, the program goes on to investigate the next generation of secondary particles. The sampling of the cascade process is stopped when all particles generated in interactions of the primary particle are realized. As we do not have full information about inelastic collisions we have used for their modelling procedures which are realized by satisfying the following conditions: a) in all collisions the law of energy conservation is fulfilled. b) average values of multiplicities, calculated for the large number of collisions, must be in agreement with the prescribed values. c) all secondary particles are generated according to the following equations [9]
IV. COSMOGENIC NUCLIDES
486
J. Masarik et al. / Simulation for cosmogenic nuclide production
cases for which experimental values of total cross sections were not available, we calculated them using the following inte~olation [lo]
- proton and neutron production in nucleon-nucleus collision dZN dp,, dp,
(16)
=- ~(l~~~,,+~p:,)p,[exp(-a,p:) where
+expt-a6p,>.a5],
_ pion
production in nucleon-nucleus
d=N
sl
dp,, dp, +a,
~1?l
(
al exp”-“2/E2)Pll)p,
exp(-osp,>I)/(
pi +pl
collision
+mZ)“’
03)
collision
‘32/a, A,/A,)
u&J = [email protected]
d2N
) ’
ai\i = up*/1 2.
In the case of known chemical composition of the target effective values of cross sections have been used
dp,, dp, XP,
led
0;’ E c
-+p?l.)
+ a5 exp(-~6pl]
I
04)
PII < 0,
Pll> 0,
- pion production
d2N =[
a5
in pion-nucleus exp( -4zpi)
collision +al]oip,
Xexd-wJ x (pi t-p:
Wn”kn)
n
E0
dpri dp,
log(
where u1 and o2 are cross sections on nuclei with mass numbers A, and A,, respectively (A, I; A s AZ). Nucleon and pion cross sections have been calculated using the relation
exp( -a~p2,) 1
- nucleon production in pion-nucleus
log( B =
where k is the type of the incident particle {proton, neutron, charged pion), w, is the relative statistical weight of n-type of target nuclei which gives the probability of its occurrence, ekn is the inelastic interaction cross section for a k-type particle with an n-type nucleus. The elastic scattering cross sections have been calculated using the formula Q(A)
+m2,y2,PII'
= 6.38A1”4[mb].
The differential angular distribution of secondary particles is well described by the expression
0,
due’ = exp( - 15.7~~8~A’/~), da
Xexd--a4p) x(p:,+p?
+mZ,y,
PII < 0, 09
where EC is the total collision energy in PCMS, pII, pI are the lon~tu~n~ and perpendicular components of the secondary particle momentum in PCMS, M, is the mass of pioa, a, . . ’ a9 are the parameters obtained from experimental data fitting. When the number of generated input particles or the number of realized collisions is equal to the required value, the inelastic collision densities as a function of primary particle momentum and space position are normalized to one input particle. 5. Auxiliary computations For modelling of the cascade process it is necessary to know the total inelastic cross sections for the interaction of protons, neutrons and charged pions with the target nuclei. In the region above 10 GeV/c we have supposed that the total cross section is constant. In
where A is the target nuclei mass number, p is the particle momentum in GeV/c, 8 is the scattering angle in radians. Elastic scattering is considered only for the perpendicular incoming particles. In other cases (secondary particles, an isotropic flux of primary particles) the production angles or input angles are much greater than the elastic scattering angles and therefore elastic scattering is neglected. The simulation of the inelastic interaction of particles with target nuclei is an important part of the simulation of the cascade development. To be able to generate the secondary particles initiated in the inelastic interaction it is necessary to know the momentum spectra of the secondary particles and their production angles as a function of the type and momentum of the incident particle and of the type of the target nucleus. In our calculations we have used experimentally obtained functions which give the secondary particle distribution as a function of the longitudial p,, and perpendicular p I component of their momenta in the hadron-nucleon centre of mass system. The multiplici-
J. Masarik et al. / Simulation for cosmogenic nuclide production ties of each type of secondary particles originating in hadron-nucleon collision have been calculated using the following equation
d’“: K,, = k/E.-
,
‘dpdQ
dp dQ,
>9
where Kij are the multiplicities without the correction on the nucleus excitation, E,, is the nucleus excitation energy. The choice of a random angle for the primary particle direction which has undergone multiple Coulomb scattering is realized using the Gaussian distribution with quadratic angle deviaton given by the Rossi formula [ll] (82)~~
=
=E
t
+ (E-E,)(B-E,) 3 - E,
= A/100
6. Results and discussion The MC method has been applied to production rate calculations of 22Na 26A1 and 53Mn in lunar samples. 22Na was chosen for testing and calibration reasons while 26A1 and 53Mn for investigation of the cosmic ray (CR) behaviour in the past. The calculations have been carried out for targets having the form of an infinite layer of material. The thickness of the layer was chosen to be 9000 kgme2. The effective chemical composition was taken to be in accordance with that of real lunar soils from Apollo 15 [12-141. For the purpose of investigation of production rate depth profiles the target was divided into 25 sublayers. An isotropic flux of CR protons is assumed to impinge on the target. The CR spectrum used in the calculations was taken from ref. [15]. The form of the spectrum takes into account the influence of solar activity on the galactic CR intensity through the modulation parameter $I. T(T+2E,,)(T+++m)-’
($)( $)“i[rad],
‘(T’
where pc is the momemtum of the considered particle in MeV, I,, is the radiation length of the considered material, 1 is the range between the input point and the position, in which we consider the angle of particle direction. Multiple Coulomb scattering has been considered for the secondary particles for the same reasons as for elastic scattering. All particles which take part in the hadron cascade lose their energy in these processes through: a) the ionization and excitation of atoms (ionization losses), b) the production of cascade particles (protons, neutrons, charged pions), c) the production of neutral pions, d) the excitation of nuclei. The ionization energy losses were calculated using the Bethe-Bloch formula with a correction for the effects of material density and of close binding. The excitation energy in GeV is approximated by the following experimental expression E ex
where E is the incident energy in the laboratory system, E, is the threshold energy of a particle able to develop a cascade (E, z 0.05 GeV), B is the maximum of both values A/100 and E,.
(17)
where E, is the kinetic energy of an incident particle in the laboratory system, Ej is the kinetic (in the case of pion - total) energy of j-type secondary particles in the laboratory system, d2N/dpdQ is the distribution function of j-type generated particle number as the function of secondary particle momentum p, D is the solid angle. As a part of the energy of the incoming particle is used for the target nucleus excitation, it is necessary to make corrections in the K,, multiplicities, eq. (17). We have used the correction K,J = KFJ (1 - E.JE,
487
E,cE<3GeV E>3GeV,
‘)=A
’
(T++)(T+2E,,++)
where T (MeV) is the proton kinetic energy, E, its rest energy, $I the modulation parameter (MeV), A = 9.9 X
L
0
I
100
200
300
I I 400 gcm2
Fig. 1. A comparison of calculated and measured depth dependence of production rate of ‘*Na in Apollo 15 samples. The experimental data have been taken from ref. (131. (a: + = 650, b: $I = 450, c: +J= 200 MeV). IV. COSMOGENIC NUCLIDES
J. Masarik et at. / ~irn~~ati~n for cosmogenic m&de pr~~ct~on
488
Bij
kg” Ffl r
4.0
3.0
2.0
1.0 1
0
I
100
I
Sbll
300
I
400 gcni2
Fig. 2. A comparison of calculated and measured depth dependence of production rate of 26A1in Apollo 15 samples. The experimental data have been taken from ref. [13,14]. (a: + = 100, b: #I = 200, c: rp= 450, d: 9 =i 650 MeV).
y = 2.65, for the 108; 172= 780 exp( - 2.5 X 10e4T); unmodulated GCR flux rp = 100 Me?‘, for the solar minimum GCR flux # = 300 MeV, for the solar maximum GCR flu 9 = 900 MeV and for the average GCR flux Cp= 650 MeV. In the calculation, the spectrum was corrected for the non-proton contribution, which can be estimated to be as high as 35-40%. Fig. 1 shows a comparison of calculated and measured depth dependences of the production rate of *‘Na in Apollo 15 samples. The experimental data have been taken from ref. [13]. The best compromise between the theoretical and experimental data was found for the modulation parameter + = 650 MeV. This represents an average galactic cosmic ray flux a few years before the sample collection. in the case of 26A1, which integrates the galactic cosmic ray flux over 1 Myr, fig. 2 shows a disagreement between the theoretical and experimental data. Although the theoretical curve follows the shape of the experimental dependence, at small depths there is an overproduction of 26A1.The best compromise between the theoretical and experimental data was found for the modulation parameter @= 200 MeV. Even larger disagreement has been found for the production rate of 53Mn which integrates the galactic cosmic ray flux over a few million years. Fig. 3 shows that an overproduction of s3Mn in comparison with the calculated production rate due to the average galactic cosmic ray intensity + = 650 MeV has been found. Although there are large uncertainties in the measured rates the and calculated values of 53Mn production observed
effect seems to overlap these errors. The shape
I
0
100
a
200
1
300
4oo .?
Fig. 3. A comparison of calculated and measured 53Mn production rate in Apollo 15 samples. The experimental data have been taken from ref. [14]. (a: + = 100, b: + = 450, c: QI= 650 MeV).
of the theoretical curve would imply that 53Mn was mostly produced by an unmodulated galactic cosmic ray flux #I = 100 MeV. Similar conclusions were obtained in our previous work where different methods of calculations were used [16,17].
References PI R.C. Reedy and J.R. Arnold, J. Geophys. Rev. 77 (1972) 537.
PI AX. Lavrukhina and G.K. Ustinova, Astron. J. 44 (1966) 1081. 131T.W. Armstrong and R.G. Alsmiller, Proc. Apollo 12 Lunar Sci. Conf., Geochim. Cosmochim. Acta, suppl 2 (1971) 1729. 141 ES. Light, M. Merker, R.B. Mendel and S.A. Korf, Acta Phys., suppl. 2 (1970) 745. PI N.P. Kalashnikov, VS. Remizovich and M.I. Riazanov, Stolknovenya bystrykh zaryazhenykh chastic v tverdykh telakh (Atomizdat, Moskva, 1980). PI R.G. Alsmiller, Nucl. Sci. Eng. 27 (1967) 158. [71 J. Spanier and E.M. Celbard, Monte Carlo Principles and Neutron Transport Problems (Addison-Wesley, Massachusetts, 1969). PI SM. Jermakov, Monte Carlo Metody i Smeshnye Problemy (Nauka, Moskva, 1975). 191 S.Tok&r, Ph.D. Thesis, Comer&s University, Bratislava (1983). 1101 J. Ranft and J.T. Routi, Comput. Phys. Commun. 7 (1972) 101. 1111 S. Rossi, Kosmicheskiye luchi (Atomizdat, Moskva, 1975). WI L.A. Rancitelli, J.S. Fruchter, D.M. Felix, R.W. Perkins
J. Masarik
489
et al. / Simulation for cosmogenic nuclide production
and N.A. Wogman, Proc. 6th Lunar. Sci. Conf. (NASA, Houston, 1975) p. 1876. [13] K. Nishiizumi, J. Klein, R. Middleton and J.R. Arnold, Earth and Planet. Sci. Lett. 70 (1984) 164. [14] M. Imamura, R.C. Finkel and M. Wahlen. Earth and Planet. Sci. Lett. 20 (1973) 107.
[15] G. Bonino. G. Castagnoli and M. Dardo. Cosmic Ray Conf. (IUPAP, Kyoto, 1979) vol. [16] P. Povinec, Izv. Akad. Nauk SSSR, ser. fiz. 41 [17] A.K. Lavrukhina, P. Povinec, G.K. Ustinova, skie issledovanya 22 (1984) 110.
IV. COSMOGENIC
Proc. 16th 12, p. 155. (1977) 383. Kosmiche-
NUCLIDES