Improvement of photonuclear reaction model below 140 MeV in the PHITS code

Abstract

The photonuclear reaction model in the particle and heavy ion transport code system (PHITS) code is improved for incident photon energies below 140 MeV. Japanese Evaluated Nuclear Data Library (JENDL) Photonuclear Data File 2004 (JENDL/PD-2004) is adopted to determine the total reaction cross section. The statistical decay model after excitation of the nucleus in PHITS is improved by modifying the decay widths for light nuclei under the isospin selection rule to reasonably reproduce the proton and neutron emission in the de-excitation process. The quasideuteron disintegration process is newly introduced into PHITS to handle the photonuclear reaction up to 140 MeV of incident photon energy. The accuracy of the improvements was verified by comparison with the experimental literature data. The improved PHITS can contribute to various practical applications such as neutron dose estimation in X-ray therapy.

Keywords:

• PHITS
• photon
• photonuclear reaction
• JENDL/PD-2004
• giant resonance
• statistical decay model
• evaporation model
• GEM
• quasideuteron disintegration
• JQMD

1. Introduction

The electron linear accelerator (linac) is widely employed in radiation therapy. In the linac facility, high-energy photons produced by electron bremsstrahlung can induce the photonuclear reaction, which generates secondary particles such as protons and neutrons. Thus, it is necessary in terms of radiation protection and shielding to estimate the dose rate from the secondary particles that originate from the photonuclear reaction.

Photonuclear reactions also occur in high-intensity laser facilities and synchrotron radiation facilities. In laser facilities, when the high-intensity laser hits a material, it generates high electric fields that can produce high-energy electrons inside the material [1]. These electrons can also produce photons that are capable of inducing the photonuclear reaction. In synchrotron facilities, high-energy photons are used for studies on element and structural analyses. In these facilities, protons and neutrons are produced through the photonuclear reaction. For the shielding design of these facilities, a particle transport simulation code that can handle photonuclear reactions is needed.

The particle and heavy ion transport code system (PHITS) [2] is one of the Monte-Carlo calculation codes to be utilized for this purpose [3]. PHITS was able to simulate a photonuclear reaction only for the giant resonance region using the generalized evaporation model (GEM) [4]. However, GEM cannot distinguish among the nuclear excitations caused by photons and other particles. The de-excitation of nuclei excited by photons is different from that induced by protons or neutrons at low excitation energies owing to the difference in the isospins at the excited nuclei. This distinction causes different particles to be output from the evaporation process. In addition, the photonuclear reaction cross section adopted in PHITS was derived from a simple empirical equation that poorly reproduces the evaluated data for some nuclei, especially for those with a small mass.

In this study, we improve the photonuclear reaction model in PHITS with respect to the following three points: (1) the total photonuclear reaction cross sections are replaced by the evaluated photonuclear reaction data contained in JENDL/PD-2004 [5]; (2) the evaporation model for the giant resonance of some light nuclei is modified considering the branching ratios calculated using the isospin selection rule; and (3) the quasideuteron disintegration process, which is the dominant reaction mechanism arising from 20 MeV of incident photon energy, is implemented in the Japan Atomic Energy Research Institute (JAERI) quantum molecular dynamics (JQMD) model [6] of PHITS. These improvements are implemented in PHITS 2.64 and later versions.

2. Improvement of the PHITS code

2.1. Photonuclear reaction cross section/JENDL Photonuclear Data File 2004

PHITS had employed the Lorentzian function with three parameters for representing the total photonuclear reaction cross section [7]: (1) $$\sigma \left(E\right)={\sigma }_0\frac{E^2{\Gamma }^2}{{\left(E_0^2-E^2\right)}^2+E^2{\Gamma }^2},$$(1) where E stands for the incident photon energy, namely, the excitation energy. σ₀, Γ, and E₀ are the resonance cross section, resonance width, and resonance energy, respectively. It is known that these parameters used in the former PHITS cannot reproduce the evaluated data for some nuclei, especially the lighter ones with a mass below 20. Therefore, we decided to adopt the JENDL Photonuclear Data File 2004 (JENDL/PD-2004) [5], which contains data for 68 nuclei selected for practical applications such as shielding calculations, for determining the total photonuclear cross section.

For other nuclei, we continue to use the parameterized Lorentzian function so as to obtain their photonuclear reaction cross section. Their parameters were re-evaluated to reproduce the 68 nuclei's data available in JENDL/PD-2004. In the re-evaluation, we simply assumed that Γ and E₀ are given by functions of their mass, as expressed by Equations (2(2) $$\begin{array}{ccc}\hfill E_0\left(A\right)& =& E_1+E_2{A}^{E_3},\hfill \end{array}$$(2) ) and (3(3) $$\begin{array}{ccc}\hfill \Gamma \left(A\right)& =& {\Gamma }_1+{\Gamma }_2{A}^{{\Gamma }_3}.\hfill \end{array}$$(3) ), where A is the mass number, and E_i and Γ_i (i = 1 − 3) are the fitting parameters. (2) $$\begin{array}{ccc}\hfill E_0\left(A\right)& =& E_1+E_2{A}^{E_3},\hfill \end{array}$$(2) (3) $$\begin{array}{ccc}\hfill \Gamma \left(A\right)& =& {\Gamma }_1+{\Gamma }_2{A}^{{\Gamma }_3}.\hfill \end{array}$$(3) On the other hand, we considered the dependence of the shell energies of nuclei on σ₀ because the cross sections of the giant resonance can be affected by the spherical state of each nucleus [8]. Therefore, we represent σ₀ as the function of not only A but also Z, namely, σ₀(Z, A), where Z stands for the atomic number. As in [9], we express σ₀ with σ_i (i = 1 − 3): (4) $${\sigma }_0(Z,A)={\sigma }_1{Z}^{{\sigma }_2}exp\left\{{\sigma }_3{E}_{\mathrm{sh}}(Z,A)\right\},$$(4) where E_sh stands for the shell energy of each nucleus. In the fitting, we employed the revised 2005 version of the Koura, Tachibana, Uno, Yamada mass formula [10] to obtain the value of E_sh(Z, A). The numerical values of the fitting parameters in Equations (2(2) $$\begin{array}{ccc}\hfill E_0\left(A\right)& =& E_1+E_2{A}^{E_3},\hfill \end{array}$$(2) )–(4(4) $${\sigma }_0(Z,A)={\sigma }_1{Z}^{{\sigma }_2}exp\left\{{\sigma }_3{E}_{\mathrm{sh}}(Z,A)\right\},$$(4) ) were determined by the least-square fitting of E₀(A), Γ(A), and σ₀(Z, A) given in JENDL/PD-2004. The quality of their fitting is shown in Figures 1–3. Figure 4 shows all the values of σ₀(Z, A) that Equation (4(4) $${\sigma }_0(Z,A)={\sigma }_1{Z}^{{\sigma }_2}exp\left\{{\sigma }_3{E}_{\mathrm{sh}}(Z,A)\right\},$$(4) ) can give.

Figure 1. The least-square fitting of E₀(A). E₀ obtained from JENDL/PD-2004 are represented by circles. The dashed line shows Equation (2(2) $$\begin{array}{ccc}\hfill E_0\left(A\right)& =& E_1+E_2{A}^{E_3},\hfill \end{array}$$(2) ) after least-square fitting.

Figure 2. The least-square fitting of Γ(A). Γ obtained from JENDL/PD-2004 are represented by circles. The dashed line shows Equation (3(3) $$\begin{array}{ccc}\hfill \Gamma \left(A\right)& =& {\Gamma }_1+{\Gamma }_2{A}^{{\Gamma }_3}.\hfill \end{array}$$(3) ) after least-square fitting.

Figure 3. The least-square fitting of σ₀(Z, A). σ₀ obtained from JENDL/PD-2004 are shown by circles. The fitted values using Equation (4(4) $${\sigma }_0(Z,A)={\sigma }_1{Z}^{{\sigma }_2}exp\left\{{\sigma }_3{E}_{\mathrm{sh}}(Z,A)\right\},$$(4) ) are represented by triangles.

Figure 4. All the σ₀(Z, A) values that Equation (4(4) $${\sigma }_0(Z,A)={\sigma }_1{Z}^{{\sigma }_2}exp\left\{{\sigma }_3{E}_{\mathrm{sh}}(Z,A)\right\},$$(4) ) can give.

All the obtained parameters are shown in Table 1.

Table 1. The obtained parameters for Equations (2(2) $$\begin{array}{ccc}\hfill E_0\left(A\right)& =& E_1+E_2{A}^{E_3},\hfill \end{array}$$(2) )–(4(4) $${\sigma }_0(Z,A)={\sigma }_1{Z}^{{\sigma }_2}exp\left\{{\sigma }_3{E}_{\mathrm{sh}}(Z,A)\right\},$$(4) ).

	Ei	Γi	σi
i = 1	−45.0	−1400	0.294
2	81.0	1420	1.75
3	−0.0608	−0.00184	0.00105

2.2. Light nuclei's branching ratios in the evaporation model

When a photon induces the giant-dipole resonance on a nucleus, the so-called isovector giant-dipole resonance, the excited nucleus takes the state of its isospin T = 1. The residual nucleus, after emitting the alpha particle with isospin T = 0, has to take the state of T = 1 in accordance with the isospin conservation law. However, the residual nucleus rarely can take this state T = 1 because its level density for T = 1 is too small to take, even in its ground state. Meanwhile, for emitting protons or neutrons whose isospins are both T = 1/2, the isospin of the residual nucleus is also going to be T = 1/2, which is more likely to occur in general. If the reaction channels of (γ,p) and (γ,n) are open, they will be the major reaction processes. This isospin selection rule becomes important for light nuclei with N = Z in the excitation energy range below 20–30 MeV. In contrast, as the excitation energy of the target nucleus becomes higher, the level density of T = 1 becomes larger so that the emission of alpha (T = 0) is acceptable.

Most of evaporation models such as GEM do not take into account the isospin in the de-excitation process. Hence, we improved GEM in PHITS to consider the isospin selection rule for the giant resonance of light and Z = N nuclei, which are ⁶Li, ¹⁰B, ¹²C, ¹⁴N, and ¹⁶O, by introducing the excitation-energy-dependent branching ratio calculated by the shell model [11].

2.3. Quasideuteron disintegration process in the JAERI quantum molecular dynamics model

The quasideuteron disintegration process is the dominant process in photonuclear reactions between approximately 25 and 140 MeV. We therefore implemented the quasideuteron disintegration process in JQMD in order to extend its applicable photonuclear reaction energy up to 140 MeV.

2.3.1. Reaction cross section of quasideuteron disintegration process

JENDL/PD-2004 includes the total photonuclear reaction cross section below 140 MeV. In order to implement the quasideuteron process in PHITS, we needed to distinguish the contribution of the quasideuteron process from that of the giant resonance.

The quasideuteron disintegration cross section σ_qd(E) can be expressed using the free deuteron photodisintegration cross section ${\sigma }_{\mathit{d}}\left(E\right)$ [12,13]: (5) $${\sigma }_{\mathrm{qd}}\left(E\right)=\frac{L}{A}NZ{\sigma }_{\mathit{d}}\left(E\right)f\left(E\right),$$(5) where E stands for the incident photon energy, and L, A, N, Z, and f(E) are the Levinger parameter, mass number, neutron number, atomic number, and the Pauli-blocking function [14], respectively. In this study, ${\sigma }_{\mathit{d}}\left(E\right)$ was taken from [15] (6) $${\sigma }_{\mathit{d}}\left(E\right)=61.2\frac{{(E-2.226)}^{3/2}}{E^3}\left(\mathrm{mb}\right),$$(6) where E is expressed in MeV. For f(E), we adopted the equation developed by Chadwick et al. [14]: (7) $$\begin{array}{ccc}\hfill f\left(E\right)& =& 8.3714\times {10}^{-2}-9.8343\times {10}^{-3}E\hfill \\ +4.1222\times {10}^{-4}{E}^2-3.4762\times {10}^{-6}{E}^3\hfill \\ +9.3537\times {10}^{-9}{E}^4.\hfill \end{array}$$(7) The Levinger parameters for 68 nuclei whose cross section data are included in JENDL/PD-2004 were determined from the least-square fitting of their total cross sections assuming that there is no contribution from the giant resonance above 100 MeV. The results are shown in Figure 5.

Figure 5. The Levinger parameters obtained through least-square fitting the evaluated data. The fitted parameters are represented by circles. The dashed line shows the averaged value of the circles for Z > 20, that is, L = 7.88.

The values for the heavier nuclei of Z > 20 stay steady at around eight, though they vary widely in lighter nuclei. This is because the heavy nucleus can be considered as a particle assembly with less shell structure. We adopted a constant value of 7.88 as the Levinger parameter for the nuclei whose evaluated cross sections are not included in JENDL/PD-2004.

For example, the calculated σ_qd(E) for ¹²C is shown in Figure 6 in comparison with the total cross section in JENDL/PD-2004.

Figure 6. The photon-induced total cross section of ¹²C. The open rectangles are the evaluated data from JENDL/PD-2004, and the dashed line is the separated σ_qd(E).

Figure 6 shows that the σ_qd(E) is well separated from the total photonuclear reaction cross section.

2.3.2. JQMD

JQMD is used to describe the nuclear reactions from the perspective of the dynamics of the interacting nucleons. However, photons are not transported in the framework of JQMD. Therefore, we describe the quasideuteron disintegration reaction as the de-excitation process of the nucleus excited by the incident photon. The procedures to treat the quasideuteron disintegration process are (1) to randomly choose a proton–neutron pair in the target nucleus, (2) to transfer the incident photon's energy and momentum to the selected proton and neutron, (3) to judge whether the reaction is allowed or not by the Pauli exclusion principle, (4) to perform JQMD simulation until the switching time to the static stage, and (5) to simulate the static stage by GEM. If the generation of proton–neutron pairs is prohibited at step (3), the process will be brought back to step (1).

In the selection of a proton–neutron pair, we tried to pick a pair that is similar to a real deuteron in terms of their distance and relative moment. However, the results were almost the same as those obtained from the random selection. Therefore, we decided to choose a proton–neutron pair at random.

3. Results and discussion

We benchmarked the improved PHITS code by comparing the results with the experimental data, and the agreement was generally satisfactory. In this paper, typical examples of the benchmark are presented.

3.1. Verification of isospin selection rule in PHITS

The photonuclear reaction cross sections for the incident photon energy in the range 10–40 MeV were calculated to verify the isospin selection rule in PHITS. Figure 7 compares the calculated photonuclear reaction cross section of ¹²C with experimental data.

Figure 7. (γ,X) cross sections of ¹²C for an incident photon energy in the range 10–40 MeV. The experimental results of (γ,p) [16], (γ,n) [17], and (γ,3α) [18] are represented with the open triangles in each corresponding frame.

The calculated results were obtained from PHITS 2.64, considering and ignoring the isospin selection rule, respectively. As shown in Figure 7, the alpha-particle emission is suppressed owing to the isospin selection rule around the incident photon energy of 20 MeV. It is obvious that the results obtained without considering the isospin selection rule significantly underestimate the proton- and neutron-production cross section for the incident energy of approximately 20 MeV. On the other hand, PHITS 2.64 considering the isospin selection rule can reproduce the experimental proton- and neutron-production cross sections very well. As for the cross section of (γ,3α), PHITS 2.64 considering the isospin selection rule reproduces the peak and valley positions of the experimental data, though their absolute values are overestimated and underestimated, respectively. The cause of this discrepancy is not understood yet, and further studies are necessary for reproducing the experimental data.

3.2. Verification of quasideuteron disintegration process in PHITS

The double-differential cross section for the incident photon energy range 59.3–65.2 MeV was calculated to verify the quasideuteron disintegration process that is simulated in JQMD. Figure 8 shows the photon-induced proton-emission double-differential cross section of ⁴⁰Ca at 90° with respect to the photon beam.

Figure 8. (γ,p) double-differential cross section of ⁴⁰Ca for an incident photon energy in the range 59.3–65.2 MeV at an angle of 90° calculated by PHITS 2.64 in comparison with the measurement data obtained by Ryckbosch et al. [19]. The solid line represents the calculation output, while the triangles represent the experimental data.

PHITS 2.64 reasonably reproduces the experimental results in the overall proton-emission energy range. However, PHITS slightly overestimates in the higher emission energy region. This result may suggest that the proton and neutron generated by quasideuteron disintegration are less scattered and moderated in the QMD simulation.

3.3. Simulation of X-ray therapy

In order to estimate the impact of these improvements on practical radiological protection, we calculated the neutron average flux inside a water phantom irradiated by X-rays generated from W and Cu targets bombarded by 30-MV electrons. The spectrum of X-rays was also calculated by PHITS. The water phantom was assumed to be a 30-cm cube that represents the human body. The neutron average fluxes inside the water phantom calculated by PHITS 2.64 and 2.30 are depicted in Figure 9.

Figure 9. The neutron average flux inside the water phantom produced by a 30-MV electron-induced photon. The solid line represents the calculation results of PHITS 2.64, while the dashed line represents the calculation results of the older PHITS.

PHITS 2.64 gives about double the neutron flux mainly because the isospin selection rule is considered, which results in suppression of alpha emission and an increase in neutron emission. This tendency suggests that the neutron flux inside patients treated with high-energy X-ray therapy may be underestimated in the calculations that do not take into account the isospin selection rule. Though the neutron dose itself is small compared with the photon dose, neutrons have a higher relative biological effectiveness. Therefore, accurate estimation of the neutron dose is necessary in risk assessment of second primary cancerogenesis.

4. Conclusion

We improved the photonuclear reaction model of PHITS in terms of the following three points: (1) the total photonuclear reaction cross sections are replaced by the evaluated photonuclear reaction data contained in JENDL/PD-2004; (2) the evaporation model for the giant resonance of some light nuclei is modified considering the branching ratios calculated using the isospin selection rule; and (3) the quasideuteron disintegration process is implemented in the JQMD model of PHITS. The accuracy of the calculation was verified by comparison with the literature data. Owing to these improvements, PHITS can handle photonuclear reactions up to 140 MeV of incident photon energy.

All the improvements achieved in this study are included in the latest version of PHITS 2.64 (as of November 2013). The expansion of the photonuclear reaction to a higher incident photon energy will be the objective of a future study.