Nuclear Instruments and Methods in Physics Research B 194 (2002) 123–131 www.elsevier.com/locate/nimb
A Monte Carlo track structure code for low energy protons S. Endo
a,*
, E. Yoshida a, H. Nikjoo b, S. Uehara c, M. Hoshi d, M. Ishikawa d, K. Shizuma a
a
d
Department of Quantum Energy Applications, Graduate School of Engineering, Hiroshima University, 1-4-1, Kagamiyama, Higashi–Hiroshima 739-8527, Japan b Radiation & Genome Stability Unit, Medical Research Council, Harwell, Didcot, Oxfordshire OX11 0RD, UK c School of Health Sciences, Kyushu University, Maidashi, Higashi-ku, Fukuoka 812-8592, Japan Research Institute for Radiation Biology and Medicine, Hiroshima University, 1-2-3, Kasumi, Minami-ku, Hiroshima 734-8553, Japan Received 12 March 2001; received in revised form 9 October 2001
Abstract A code is described for simulation of protons (100 eV to 10 MeV) track structure in water vapor. The code simulates molecular interaction by interaction for the transport of primary ions and secondary electrons in the form of ionizations and excitations. When a low velocity ion collides with the atoms or molecules of a target, the ion may also capture or lose electrons. The probabilities for these processes are described by the quantity cross-section. Although proton track simulation at energies above Bragg peak (>0.3 MeV) has been achieved to a high degree of precision, simulations at energies near or below the Bragg peak have only been attempted recently because of the lack of relevant cross-section data. As the hydrogen atom has a different ionization cross-section from that of a proton, charge exchange processes need to be considered in order to calculate stopping power for low energy protons. In this paper, we have used state-ofthe-art Monte Carlo track simulation techniques, in conjunction with the published experimental and established theoretical data, to develop a model for the extension of the proton track simulation to the low energy region. Data are presented on charge-state-fraction, proton stopping power, range and averaged energy producing an ion pair (Wvalues) in a mixture of hydrogen (H2 ) and Oxygen (O2 =2) gas. The results are compared with the published experimental data. Ó 2002 Elsevier Science B.V. All rights reserved. Keywords: Track simulation; Proton track; Proton stopping power; W value
1. Introduction For many years at Hiroshima University, we have studied the radiobiological effects of atomic-
*
Corresponding author. Tel.: +81-824-24-7612; fax: +81824-22-7192. E-mail address:
[email protected] (S. Endo).
bomb neutrons using a neutron generator (HIRRAC) [1]. Recent developments in track-structure studies and biophysical modeling [2–5] have made it possible to make an attempt to simulate tracks of low energy neutrons at a molecular level using Monte Carlo track structure simulation methods. In general, neutron collisions with atomic nuclei lead to recoil atomic ions and nuclear reaction products in the form of secondary charged particles.
0168-583X/02/$ - see front matter Ó 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 5 8 3 X ( 0 2 ) 0 0 4 9 7 - 4
124
S. Endo et al. / Nucl. Instr. and Meth. in Phys. Res. B 194 (2002) 123–131
The energies of these particles vary according to the energy of the incident neutrons. In hydrogen containing media such as water and tissue the most important interaction is elastic scattering with hydrogen nuclei (protons), which accounts for more than 90% of energy transfer. Protons are, therefore, considered to be the most important recoil particles for estimating neutron induced radiation effect. The recoil protons and other ions set in motion in the cell, in turn, interact with the biomolecule leading to DNA damage and subsequent biological lesions. Most of these recoil protons are low energy particles below 1 MeV energy. Therefore, the motivation for this work arises for a need to simulate the tracks of low energy protons emitted in neutron interactions in tissue as these simulations are needed in biophysical modeling. Although proton track simulation at energies above the Bragg peak has been achieved to a high degree of precision [2] simulations at energies near or below the Bragg peak have not been attempted because of the lack of relevant cross-section data. As the hydrogen atom has a different ionization cross-section from that of the proton, charge exchange processes (CEP) need to be considered in order to calculate the electronic stopping power for low energy protons. In this paper, we have used state-of-the-art Monte Carlo track simulation techniques, in conjunction with the published experimental and established theoretical data, to develop a model for the extension of the proton track simulation in the low energy region. Data are presented on charge-state-fraction, proton electronic stopping power, range and averaged energy producing an ion pair (W-values) in a mixture of hydrogen (H2 ) and Oxygen (O2 =2) gas at the density of 334:3 1020 molecules=cm3 . Results compare well with the published experimental data.
cross-sections were calculated using the Rutherford-formula taking into account the screening parameters given by Moliere [8] (T P 80 keV). Below 80 keV, the cross-section formula proposed by Poter and Jump [9] was used. The excitation cross-sections down to 10 eV were calculated using the empirical formula given by Berger and Wang [10]. Below 10 eV, the excitation cross-sections were used with a function given by Zaider et al. [7] which takes into account the vibrational crosssections. Energy transfer for the excitation event is obtained by averaged values as kexc of which kinetic and potential events are selected by ejected electron energy and electron binding energy, according to fractions of partial ionization crosssections [6]. The fractions of partial ionization cross-sections and the binding energies (Bi ) for water molecule-subshells 1b1 , 3a1 , 1b2 , 2a1 and for the oxygen O1s or 1a1 state are taken from the published data in [6] (Bi ¼ 12:62, 14.75, 18.51, 32.4 and 539.7 eV for the five sub-shells, respectively). Calculation of secondary electron energies was carried out using a formula given by Green and Sawada [11]. A response function taken from the experimental data of Opal et al. [12] was used for the angular distribution of ejected electrons. 2.2. Proton code In this work, the CEP is taken into account in track simulation in the low energy region. In the absence of experimental data, electron capture and electron stripping cross-sections were calculated for hydrogen and oxygen by assuming water vapor composition as a mixture of H2 and O2 /2. Total cross-sections for proton (rT p ) and hydrogen (rT H ) are given by rT p ¼ rT p þ rion p þ rexc p þ relas p þ r10 ;
ð1aÞ
2. Materials and methods
rT H ¼ rT H þ rion H þ rexc H þ relas H þ r01 ;
ð1bÞ
2.1. Electron code
where rion , rexc , relas , r10 and r01 are the crosssections of ionization, excitation, elastic, capture by proton and stripping by hydrogen atom, respectively. In this code, we have ignored the higher order of charge-exchange processes such as p !
An electron transport code was developed according to formulations published by Uehara et al. [6] and Zaider et al. [7]. The elastic scattering
S. Endo et al. / Nucl. Instr. and Meth. in Phys. Res. B 194 (2002) 123–131
H , H ! H2þ , etc. For a particle of energy E and a charge state C, the path length (x) to the next interaction is sampled by x¼
1 logðnÞ; rT ðE; CÞn
where the fitting parameters K and C are listed in Table 1. The calculated cross-sections are shown in Fig. 1. The experimental data for ionization, excitation, electron capture and electron stripping [16–19] are also plotted in the same figure. In the absence of published experimental data, the elastic cross-section for incident proton and hydrogen atoms were assumed to be a universal potential and an effective charge depending on the incident particle velocity [20]. The magnitude of this elastic cross-sections is close to a result by Uehara et al. using the impact parameter methods [21]. Also, to obtain the total cross-section for incident hydrogen, it is assumed that a cross-section ratio of rðH þ moleculeÞ=rðp þ moleculeÞ scales to rðH þ H2 Þ=rðp þ H2 Þ, where the elastic crosssections for p þ H2 and H þ H2 reactions are calculated by Krstic and Schultz [22]. The ratio was fitted by the function rH =rp ¼ 1 þ 0:0224 logðEÞ þ 0:01285 logðEÞ2 where E is the proton energy. The elastic cross-sections for protons and hydrogen atoms are shown in Fig. 2. For the excitation cross-sections, we adopted the formulation by Miller and Green [14]. Table 2 presents the 12 parameters used for the calculations of the 12 excited states of water molecules as used in our calculations. Ejected electron spectra produced by protons and hydrogen atoms are computed using the Rudd’s formula [16]. To save computation time, angular distributions of ejected electrons were kinematically calculated by
ð2Þ
where n is a random number uniformly distributed in the range (0,1) and n ¼ 334:3 1020 molecules=cm3 is the atomic density of water vapor. The cross-sections for ionization by protons and hydrogen atoms, and the electron capture by protons are assumed to be the cross-sections for a mixture of H2 and O2 =2 as water vapor [13], and electron stripping of hydrogen atoms in water vapor [14]. The cross-sections rion p , rion H and r01 , are obtained from X
rðT Þ ¼ r0
m
ðZaÞ ðT IÞ ; J ðXþmÞ þ T ðXþmÞ
ð3Þ
where r0 is 1 1016 cm2 , I is the ionization threshold, X and m are initial and asymptotic powers, and J is an adjustable parameter [13]. The value of the parameters given by Green and McNeal [13] and Miller and Green [14] are listed in Table 1. The difference between the ionization cross-sections for water vapor and a mixture of H2 and O2 =2 is estimated to be less than 16% [15]. The cross-section for electron capture by protons (r10 ) is presented as follows: X
rðT Þ ¼ r0
m
ðZaÞ ðT IÞ
X
J ðXþmÞ þ T ðXþmÞ þ ðZaÞ T m ðE=CÞ
K
;
125
M þ me cos h ¼ 2M
ð4Þ
rffiffiffiffiffiffiffiffiffi MQ ; me E
ð5Þ
Table 1 Ionization parameters by protons and hydrogen atoms, electron capture by protons, and stripping by hydrogen atoms [13,14] Parameters Proton ionization Hydrogen ionization Electron capture by proton Stripping of hydrogen
H2 O2 H2 O2 H2 O2 H2 O
m
J (keV)
a (keV)
X
K
C
1.02 0.82 1.04 0.43 2.0 2.0 0.943
47.62 69.08 23.68 73.11 1.215 0.057 27.7
171.5 105.6 47.45 73.88 4084.0 1038.0 79.3
0.75 0.75 0.75 0.75 0.271 0.258 0.652
– – – – 4.8 3.5 –
– – – – 75.8 125.0 –
126
S. Endo et al. / Nucl. Instr. and Meth. in Phys. Res. B 194 (2002) 123–131
Fig. 1. Ionization (p for proton, H for hydrogen atom), excitation, capture (capt.) and stripping (strip.) cross-sections. Comparison with this code and experimental data [16–19].
Fig. 2. Elastic cross-sections of proton and hydrogen atom in matter. Plus symbol is the calculated result assuming impact parameter by Uehara [21].
where me is the electron rest mass, M and E are the mass and kinetic energy of the incident particle,
and the energy transfer Q is defined in term of the ejected electron energy e and binding energy B,
S. Endo et al. / Nucl. Instr. and Meth. in Phys. Res. B 194 (2002) 123–131
127
Table 2 Excitations induced by protons [13] Transition
W (eV)
a (keV)
J (keV)
X
m
m2 m1 m3 m1 þ m3 Dissociative Dissociative Diffuse bands H Lyman a H Lyman b OH A þ B [3] C þ B [3]
0.198 0.453 0.455 0.899 7.4 9.67 13.32 17.0 19.0 10.0 10.0 10.0
0.0422 0.0910 0.0936 0.116 0.912 1.47 7.28 9.31 1.17 0.324 0.249 0.555
0.627 1.43 1.48 2.85 14.4 18.7 26.6 34.1 34.1 7.09 19.6 21.6
6.0 6.0 6.0 6.0 0.75 0.8 0.75 0.75 0.75 1.0 0.75 0.75
2.0 2.0 2.0 2.0 1.0 1.0 1.0 1.0 1.0 5.35 1.0 1.0
Q ¼ e þ B:
ð6Þ
For the width of the emission angle, angular distributions of ejected electrons were calculated by the function 2
dr CðT Þ / ; dXdT fðh h0 Þ2 þ CðT Þ2 g
ð7Þ
CðT Þ ¼ aT b ;
ð8Þ
where a ¼ 3:702, b ¼ 0:4969 and h0 was determined by Eq. (5). The energy transfer for electron capture (Qc ) is given by Q c ¼ BW BH ;
ð9Þ
3.1. Electron radial distribution A part of the electron transport in the present code is very similar to the well-established code kurbuc [6]. We indicate an example of a radial distribution of energy deposition at the center of spheres with spherical shells at 1 nm intervals from the origin for 1 keV electron track comparison with the code kurbuc [6] as shown in Fig. 3. Both this code and kurbuc give similar results for the distribution. 3.2. Charge state fraction The charge state fraction is calculated from the number of protons (Np ) and hydrogen atoms (NH )
where BW and BH are electron binding energies of the water molecule (Bi ¼ 12:62, 14.75, 18.51, 32.4 and 539.7 eV for the 5 sub-shells, respectively) and the hydrogen atom (13.6 eV), respectively. The energy transfer for the stripping process (Qs ) by ejected electron (energy denoted by Te ) is given by Q s ¼ BH þ T e :
ð10Þ
3. Results and discussion In order to check the accuracy and consistency of the code, a number of physical quantities such as charge state fraction, electronic stopping power, proton range and W-values for protons are calculated and compared with the experimental data.
Fig. 3. Radial distribution from electron transport code compared to result for kurbuc [6] and this code.
128
S. Endo et al. / Nucl. Instr. and Meth. in Phys. Res. B 194 (2002) 123–131
Fig. 4. Charge state fraction for H2 gas and mixed gas of H2 –O2 =2 as a function of proton energy. Plots by open and closed circles are experimental charge state fractions discussed by Fainstein et al. [23]. H0 and Hþ mean hydrogen atom and proton, respectively. Data are well reproduced by this calculation for the H2 gas.
calculated from the relationships F ðHþ Þ ¼ Np = ðNp þ NH Þ and F ðH0 Þ ¼ 1 F ðHþ Þ, ignoring the higher order exchange processes. The charge state fractions F ðHþ Þ and F ðH0 Þ as a function of proton energy in hydrogen gas and a mixed gas are shown in Fig. 4. The fraction of protons decreases very rapidly at proton energies below 100 keV, reducing to nearly 10% at 10 keV. In the case of the charge state fraction in hydrogen gas, the calculated result is consistent with the experimental data in hydrogen gas [23] as plotted in Fig. 4. The present calculations are in reasonable agreement with the experimental data.
3.3. Stopping power and range for protons The calculated electronic stopping power and range are shown in Figs. 5 and 6 and in comparison with the ICRU data [24]. The calculations are presented for electronic and nuclear stopping powers showing good agreement with the ICRU data. There are some discrepancies around the 1 keV and 10 MeV regions. This is due to differences between the magnitudes of cross-sections of rðH2 OÞ and rðH2 Þ þ rðO2 Þ=2. Similar tendencies
for proton ranges are observed between the ICRU data and this calculation. 3.4. W value The W-value, which is defined as the averaged energy to produce an ion-pair, is approximately constant at about 31 eV for protons down to 100 keV as shown in Fig. 7. Below 100 keV, the Wvalue increases rapidly and is greater than 100 eV at proton energies below 1 keV. As experimental W-values for water vapor are not available in this energy range, comparisons have been made with those for tissue equivalent (TE) gas [25–29]. Our calculations agree well with the experimental data for proton energies down to 100 keV. At lower energies, calculated W-values deviated by about 10% from the experimental values. 3.5. Radial distribution of protons Radial distributions of absorbed energy were calculated by 1000 full slowing down tracks of 10, 50 keV, 0.1, 0.3 and 1 MeV protons at the center axis of cylindrical shells. The distributions of absorbed energy are plotted in Fig. 8 and compared
S. Endo et al. / Nucl. Instr. and Meth. in Phys. Res. B 194 (2002) 123–131
129
Fig. 5. Comparison between calculated and data [24] for protons both for the electronic- and nuclear-stopping power.
Fig. 6. Comparison of proton range between the ICRU [24] and this code.
with previous experimental data for 1 MeV protons by Wingate and Baum [30]. The calculated result fits well with the data. The maximum range for 0.1, 0.3 and 1.0 MeV protons correspond to the maximum energies of ejected electrons (0.022, 0.058 and 0.22 lm, respectively). At the incident
proton energies lower than 0.1 MeV, the Rudd formula takes into account for the high energy electron tail around about 200 eV to reproduced experimental data. Therefore, the maximum ranges in the low proton energy region (
