Recent improvements of the particle and heavy ion transport code system – PHITS version 3.33

ABSTRACT

The Particle and Heavy Ion Transport code System (PHITS) is a general-purpose Monte Carlo radiation transport code that can simulate the behavior of most particle species with energies up to 1 TeV (per nucleon for ions). Its new version, PHITS3.33, was recently developed and released to the public. In the new version, the compatibility with nuclear data libraries and the algorithm of the track-structure modes have been improved, and they are recommended to be used for certain simulation conditions. Some utility functions and software have been developed and integrated into the new PHITS package, such as PHITS Interactive Geometry viewer in 3D (PHIG-3D) and RadioTherapy packaged based on PHITS (RT-PHITS). With these upgraded features, PHITS can be applied in a wide diversity of fields – beyond traditional nuclear engineering domains – including cosmic-ray, environmental, medical, life, and material sciences. In this paper, we summarize the upgraded features of PHITS3.33 with respect to the physics models, utility functions, and application software introduced since the release of PHITS3.02 in 2017.

KEYWORDS:

• PHITS
• Monte Carlo
• particle transport simulation
• nuclear reaction
• atomic interaction
• nuclear data
• track-structure simulation

1. Introduction

The Particle and Heavy Ion Transport code System, PHITS [1], is a general-purpose Monte Carlo radiation transport code that can simulate the behavior of most particle species with energies up to 1 TeV (per nucleon for ions). Before 2019, the code was distributed through the Research organization for Information Science and Technology (RIST), the Data Bank of the Organization for Economic Co-operation and Development’s Nuclear Energy Agency (OECD/NEA DB), and the Radiation Safety Information Computational Center (RSICC). Starting from 2019, Japan Atomic Energy Agency (JAEA) directly provides the code free of charge upon requests, following the Japanese export control rules. Since then, more than 4,000 users from 68 countries have been registered to the code, and the total number of users is expected to be over 10,000, including those who obtained PHITS through RIST, OECD/NEA DB, and RSICC before 2019.

In 2023, we released a new version of the code, PHITS3.33. The basic functionalities of PHITS implemented before the release of version 3.02 in 2017 were described in our previous papers [1–3]. In addition, the comprehensive benchmark studies of the recent versions of PHITS were performed and reported in the references [4,5], including another paper in this special issue [6]. Thus, this paper focuses on summarizing the upgrade features of PHITS introduced after version 3.02 in terms of physics models, utility functions, and embedded application software.

2. Upgrades in physics models

Figure 1 shows the physics models and data libraries recommended for use in PHITS3.33 to simulate nuclear and atomic collisions. They are the same as those recommended in PHITS3.02 except for the introduction of JENDL-5 [7] for neutrons and protons up to 200 MeV for certain applications and KURBUC [8,9]/ITSART [10] for proton and ion track-structure simulations. The details of the upgrade features made on the physics models in PHITS after version 3.02 are listed below.

Figure 1. Map of the physics models and the data libraries recommended for use in PHITS3.33 to simulate nuclear and atomic collisions.

2.1. Improvements on handling nuclear data libraries

Functions to read nuclear data libraries for deuterons, α particles, and photons have been implemented. In addition, the JENDL-5 deuteron sublibrary, which is mostly based on JENDL/DEU-2020 [11], has been included in the PHITS package. The database of the neutron production yields from d-D reaction written in the PHITS-readable format (frag data) has also been prepared [12]. Furthermore, treatment of thermal neutron scattering libraries, including phonon resonances, has been improved, with the support for continuous energy and mixed elastic libraries. These new features enable reliable design of neutron sources based on deuteron beams using PHITS.

In addition, JENDL-5 for protons and neutrons up to 200 MeV for 30 selected nuclei that are widely used in the particle transport simulations has been included in the PHITS package. Their high-energy parts are mostly based on JENDL-4.0/HE [13,14], a high-energy extended version of JENDL-4.0 [15], whose accuracy was well validated for neutron production and penetration calculations [5]. Owing to these features, JENDL-5 is now recommended for shielding calculations and neutron-source designs using PHITS. Note that nuclear data libraries with energies over 20 MeV cannot be used as an event generator, so the conventional nuclear reaction models, INCL4.6 [16] coupled with GEM [17], are still recommended for uses of other applications such as medical physics and detector designs. The complete dataset of JENDL-5 written in the PHITS-readable ACE format is available from the JAEA website [18].

For investigating the impact of this change on shielding calculations, we examined the shielding profiles of concrete against high-energy neutrons generated by p-Li and p-Cu reactions using PHITS ver. 3.02 and 3.33. The default nuclear data library and reaction models such as JENDL-4.0 and INCL4.6 were employed in PHITS3.02, while JENDL-5 was selected in PHITS3.33 for simulating neutron- and proton-induced nuclear reactions below 200 MeV by following our new recommendation for shielding applications. In the PHITS simulations, sequential cylinders of a thick Li or Cu target and a concrete shield with a radius of 100 cm were irradiated by 100 MeV protons, and the energy spectrum of neutrons incident to the concrete shield and the effective doses inside it were calculated using [t-cross] and [t-track] tallies, respectively. The thicknesses of the target and the concrete shield were set to 15 and 400 g/cm², respectively. The incident protons were fully stopped in the target. The history numbers for all simulations were fixed to 10⁸, and the variance reduction techniques and the parallelization protocols were not used in the simulations.

The calculated results are shown in Figures 2 and 3. The statistical uncertainties are negligibly small – less than a few percent in most cases. For the Li target, the neutron spectrum obtained from PHITS3.33 is much harder than that from PHITS3.02 because the quasi-mono-energetic peak of neutrons generated by p-⁷Li reaction can be reproduced by JENDL-5. The tendency is opposite in the case of Cu target; the neutron spectrum becomes slightly softer when PHITS3.33 was adopted in the simulation. Consequently, PHITS3.33 gave lower and higher gradients of the depth-dose curves for Li and Cu targets, respectively, compared with PHITS3.02, though the difference between the versions is trivial in the Li target case. These tendencies are consistent with the benchmark results published in our previous paper using JENDL-4.0/HE [5], suggesting that the appropriate shielding thickness for high-energy accelerator facilities tends to be thinner when the version of PHITS employed in the shielding calculation is updated from 3.02 to 3.33.

Figure 2. Neutron spectra incident to the concrete shield placed behind a thick (a) li or (b) Cu target irradiated by 100 MeV protons calculated by PHITS ver. 3.02 and 3.33.

Figure 3. Effective doses inside the concrete shield placed behind a thick (a) li or (b) Cu target irradiated by 100 MeV protons calculated by PHITS ver. 3.02 and 3.33.

It should be noted that the computational time is reduced by approximately one-third when the new version was employed in these shielding simulations; they are approximately 18,400 and 6,200 seconds for versions 3.02 and 3.33, respectively, in the case of the Cu target using a conventional Windows PC with a 3.6 GHz CPU. The shorter computational time of the new version is owing to the adoption of the nuclear data library (JENDL-5) for simulating proton and neutron-induced reactions below 200 MeV, which is much faster than the conventional nuclear reaction models (INCL4.6 + GEM).

2.2. Improvements on nuclear reaction models

A model for simulating charge-current and neutral-current neutrino interactions with ¹H, ²H, and electrons up to 150 MeV has been developed and implemented in PHITS, based on the cross sections evaluated by Strumia et al. [19] and Nakamura et al. [20]. This model allows the users to simulate the response of D₂O or DHO water-tank neutrino detectors. Note that 150 MeV as the energy upper limit is sufficient for solar and reactor neutrinos.

A new model to describe the fission probability of sub-actinide interacting with MeV-range radiation has been implemented in PHITS [21]. This model is used in tandem with INCL4.6 [16] and GEM [17] to consider the competition between evaporation and fission.

The gamma de-excitation model, EBITEM [22], has been upgraded to version 2. The details of the upgrade features are described in another paper in this special issue [23].

2.3. Improvements on atomic interaction models

KURBUC [8,9], a well-verified track-structure code applicable to protons up to 300 MeV and carbon ions up to 10 MeV/n in water, has been incorporated into PHITS [24]. Examples of the carbon-ion trajectories in water simulated using the continuous slowing-down approximation based on ATIMA [25] and track-structure mode based on KURBUC in PHITS are shown in Figure 4. In addition, a new track-structure code applicable to arbitrary ions and materials, ITSART [10], was developed and implemented. With these new features together with the electron track structure (ETS) mode [26] implemented before version 3.02, PHITS has become a unique tool that can handle both microscopic and macroscopic particle transport simulations in any material, though the accuracy of the microscopic simulation has not been fully verified except for water [27]. Note that the track-structure modes are not the default in PHITS because they are too time-consuming to be utilized in macroscopic particle transport simulations.

Figure 4. Trajectories of carbon ions with 5 MeV/n in water simulated using the continuous slowing-down approximation, based on ATIMA (upper), and track-structure mode based on KURBUC (lower) in PHITS.

A model for simulating the photon-induced muon pair production has been developed and implemented in PHITS [28]. In this model, the differential cross sections of the reaction are calculated by the Born approximation based on Tsai [29]. The reproduction of this process is of great importance for shielding calculations of high-energy electron accelerators.

A function to read the stopping power of each material from a user-supplied table has been implemented. As for the sample data, the stopping power of protons and α particles in various materials provided by the National Institute of Standards and Technology (NIST) database [30] as well as those of protons and carbon ions in liquid water, graphite, and dry air given in the International Commission on Radiation Units and Measurements (ICRU) Report 90 [31] were included in the PHITS package. In addition, the database of density correction factors used in the evaluation of ICRU Report 90 was prepared for the EGS mode [32].

3. Upgrades in utility functions

This section summarizes the major upgrade features included after version 3.02 in terms of the utility functions for extending the capabilities of PHITS.

3.1. Improvements on tallies

A function to analyze the multiple tally results obtained from the same or different PHITS simulations, so-called anatally, has been developed and implemented. Using this function, the systematic uncertainties due to the fluctuation of a certain parameter in geometry, source, or reaction model can be deduced [33]. Furthermore, the users can provide their own program with this function for analyzing the tally results in their specific way. With this feature, the photon-isoeffective doses can be directly calculated with PHITS by providing a user-defined program for analyzing the results of the [t-deposit], [t-let], and [t-sed] tallies based on the stochastic microdosimetric kinetic model [34]. This program is included in the PHITS package as an example of the user-defined program for the anatally function [35].

A new option in the [t-dpa] tally for calculating displacements-per-atom, so-called DPA, based on the athermal-recombination-corrected (arc) model [36] has been introduced. This feature greatly improves the reliability in estimating the operating lifetimes of materials in radiation environments because the arc model can consider the defect production efficiency for each material, which was ignored in the standard model adopted in the original [t-dpa] tally [37].

A new option for the microdosimetric function, i.e., the [t-sed] tally, has been developed using the track-structure algorithms in PHITS [38]. With this option, PHITS can calculate the number of ionizations and electronic excitations occurring in microscopic water targets distributed in macroscopic matter. This feature allows for direct estimation of the DNA damage yields created in the whole human body. This capability is unique to PHITS because other codes explicitly distinguish between the track-structure and condensed-history simulations so that it is impractical to perform the DNA-scale microscopic simulation in the macroscopic matter.

The decay chain analysis code bundled with PHITS, DCHAIN-SP [39], together with its associated tally in PHITS, [t-dchain], have been modernized with updated data libraries and new features, such as the support for tetrahedral and three-dimensional grid mesh geometries. The modernized version, named DCHAN-PHITS [40], allows easy transfer of the resulting decay radiation spectra directly into a secondary PHITS simulation to obtain the corresponding decay doses from emitted α and β particles, and γ-rays.

3.2. Improvements on source-generation functions

New options in the [source] section have been introduced for reproducing the cosmic-ray environments in free space, on typical low-earth orbits (LEO), and anywhere in the Earth’s atmosphere. They are called the free-space, the LEO, and the terrestrial modes, respectively. In the free-space mode, the fluxes of the galactic cosmic-rays (GCR) and the solar energetic particles (SEP) are determined from the models proposed by Matthiä et al. [41] and Tylka et al [42], respectively. In the LEO mode, the influence of the Earth’s magnetosphere on the GCR and SEP fluxes can be considered using the pre-calculated geomagnetic transmission functions [43], while the fluxes of trapped protons are determined from the AP8 model [44] implemented in SPENVIS [45]. In the terrestrial mode, the PARMA model [46,47] is adopted for calculating the GCR fluxes, including their secondary particles generated through the extensive air showers in the atmosphere. The complicated zenith angle distributions of the terrestrial cosmic-ray fluxes are also considered in the mode.

3.3. Other improvements

Coupled calculations with thermal analysis codes, such as ANSYS® Fluent, have become feasible. For this purpose, we improved PHITS to be capable of reading the NASTRAN bulk data and outputting the tally results in the OpenFOAM field data, which are general formats for representing the tetrahedral mesh geometry used in thermal analysis codes.

A function to read electro-magnetic field maps written in xyz or r-z grids has been implemented. With this implementation, PHITS can trace the particle trajectories not only in idealized electro-magnetic fields such as dipole and quadrupole magnets but also in complex fields with arbitrary direction and strength.

A function for calculating the response functions of organic scintillators based on SCINFUL-QMD [48] has been implemented [49]. For this implementation, we updated the total, elastic, and inelastic cross-section data of the hydrogen and carbon nuclei for neutrons and developed a new scorer to analyze the light outputs from a scintillator. In addition to detector response modelling, the updated cross sections are also useful for estimating the relative biological effectiveness (RBE) for neutron exposure, because the cross sections of some updated channels, such as the n(12C, 3α)n’ reaction, are hardly reproduced by conventional nuclear reaction models but are critical for the biological calculations.

A function to add user-defined particles with specific cross sections and decay schemes into the list of transport particles in PHITS has been implemented. The energy and angular distributions of particles generated in decay can be arbitrarily set by the user. Input files and source codes dedicated to this function are encouraged to share in GitHub for promoting the use of PHITS in the research field of elementary particle physics.

A continuous energy adjoint transport algorithm for photons has been developed and implemented in PHITS [50]. This algorithm can improve the efficiency of photon transport problems where the spatial size of the tally is relatively small compared to the source. Such problems are commonly seen in environmental radioactivity calculations.

4. Upgrades in embedded application software

The PHITS package contains various application software that help users to create input files and analyze/visualize the tally results. The following are the major application software included after version 3.02.

PHITS Interactive Geometry viewer in 3D (PHIG-3D) has been developed based on Gxsview [51]. It can read a PHITS input file and visualize the 3-dimensional (3D) view of its geometry by interactively changing the camera location and angle. The cross-sectional view can also be drawn by specifying a certain cutting surface, as shown in Figure 5. The geometry errors, i.e., double defined and undefined regions, can be detected and depicted in the geometry view. Note that PHITS itself has a native function for drawing the geometry in 3D with the [t-3dshow] tally, but its camera location and angle cannot be controlled interactively and users must modify the tally and re-execute PHITS to change the view.

Figure 5. Sample screen of PHIG-3D.

RadioTherapy package based on PHITS (RT-PHITS) has been developed for applying PHITS to medical physics using the Digital Imaging and Communications in Medicine (DICOM) dataset [52,53]. It consists of four major modules named CT2PHITS, PET2PHITS, Plan2PHITS, and PHITS2DICOM. The first two modules convert the CT and PET/SPECT images to the PHITS input format for reproducing the patient’s body geometry and the radionuclide distribution inside it, respectively. Plan2PHITS enables reproduction of the beam configuration for each patient treated by carbon-ion therapy based on the RT-Plan data. PHITS2DICOM converts the [t-deposit] tally outputs to the RT-Dose data for visualizing and analyzing the 3D dose distribution calculated by PHITS using commercial DICOM software. RT-PHITS is equipped with a graphical user interface (GUI) as shown in Figure 6, which allows users to interactively control its parameters, such as the area of data conversion.

Figure 6. Sample screen of RT-PHITS.

A program for calculating the DNA damage yields using the patterns of ionizations and excitations obtained from the track-structure simulation of PHITS has been developed. The calculated results can be outputted in the standard DNA damage (SDD) data format [54]. The accuracy of this program for estimating the yields of single-, double-strand breaks (SSB, DSB) and complex DSB for electron and proton irradiations has been well validated [55,56].

A program for converting data files in the Experimental Nuclear Reaction Data (EXFOR) database operated by IAEA [57] to the frag data format in PHITS has been developed. It is named Exfor2frag and is included in the PHITS package. The types of cross sections that can be converted to the frag data format are non-elastic cross sections, nuclear production cross sections, and differential cross sections of nuclear production cross sections for particle emission energies and emission angles.

5. Conclusions and future work

We upgraded various aspects of the PHITS code and released the new version as PHITS3.33. The new version is available through not only JAEA but also RIST, OECD/NEA DB, and RSICC, depending on the country where the users’ institutes are located. With these upgraded features, PHITS can be applied in a wide diversity of fields – beyond traditional nuclear engineering domains – including cosmic-ray, environmental, medical, life, and material sciences. More improvements and developments are work-in-progress. The improvement and verification of the track-structure modality are in progress, particularly for materials other than water [58,59]. A text editor specialized for making PHITS input files with user-assistance functions is under development and to be released in the near future.

Acknowledgments

We are grateful for Mr. K. Sakamoto, Mr. T. Imaki, Mr. T. Kitamura, and Mr. T. Miura of Computational Science & e-Systems (CCSE) of JAEA, Dr. N. Furutachi and Mr. A. Murofushi of RIST, and Mr. T. Sugita and Mr. T. Okamoto of NAIS Co. Ltd. for their support in programming PHITS, Dr. S. Ohnishi of National Maritime Research Institute (NMRI) for developing PHIG-3D, Dr. José Ignacio Marquez of European Spallation Source (ESS) for his support on improving the function to read thermal scattering libraries, Dr. T. Nishitani of Nagoya University for providing the database of the neutron production yields from d-D reaction, Dr. C. Konno and Dr. K. Tada of Research Group for Reactor Physics and Thermal Hydraulics Technology of JAEA and Dr. S. Nakayama of Nuclear Data Center of JAEA for providing JENDL-5 written in the PHITS-readable ACE format, Dr. D. Satoh of Research Group for Radiation Transport Analysis of JAEA for implementing the algorithm and database of SCINFUL-QMD, Dr. Y. Namito and Dr. H. Hirayama of KEK and Dr. M. Shimizu of AIST for providing the database of density correction factors used in the EGS mode, Dr. A. Malins of CCSE of JAEA for implementing the continuous energy adjoint transport algorithm, Dr. H. Iwamoto of Research Group for Nuclear Transmutation System of JAEA for his support on improving the high-energy fission model, Dr. T. Liamsuwan of Chulabhorn Royal Academy and Dr. H. Nikjoo for their support on implementing the algorithm of KURBUC, and Dr. Y. Koba, Dr. S. Yonai, Dr. S. Matsumoto, and Dr. W. Chang of National Institutes for Quantum Science and Technology, and Dr. T. Watabe and Mr. H. Sasaki of Osaka University for their support in developing RT-PHITS.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was partially supported by Japan Society for the Promotion of Science (JSPS), KAKENHI, Grant Numbers 16H02959, 18K14159, 18K19899, 21H02872, 22H03744, 26790072.