© S. Chen et al., Published by EDP Sciences, 2025

1. Introduction

Monte Carlo methods, are known for their capacity to realistically depict stochastic phenomena, facilitate the simulation of intricate physical experiments, exhibit minimal geometric constraints, and harness parallel computing capabilities, thereby experiencing a utilization surge within the fields of nuclear technology [1, 2]. The photon and electron coupling transport simulation using Monte Carlo method has been successfully applied in nuclear technology applications such as industrial irradiation, non-destructive testing, nuclear logging, nuclear medicine and the development of nuclear instruments. In industrial irradiation, Monte Carlo particle transport simulations are employed for the design and verification of dose detection system, as well as the optimization of shielding designs for irradiation equipment [3, 4]. In the non-destructive testing, Monte Carlo particle transport simulations can be used to optimize collimators, radiation shielding, imaging quality [5, 6], etc. For example, in X-ray imaging, the thickness of the intensifying screen is optimized to improve imaging quality. In the design and parameter calibration of detectors, the use of Monte Carlo particle transport simulation can significantly improve the economic efficiency of related work [7, 8]. For example, simulation of detectors with different structures and sizes can reduce the number of experiments. In nuclear logging, the development of equipment and the production of interpretation chart require particle transport to improve the performance [9, 10]. In nuclear medicine, Monte Carlo particle transport methods are applied in the production of medical isotopes, dose evaluation, and the development of nuclear medicine equipment [11, 12].

There are several key issues about the Monte Carlo particle transport method in nuclear technology applications. In radiation therapy, it is necessary to use medical images with DICOM (Digital Imaging and Communications in Medicine) format to generate the computational model for particle transport simulation, enabling the evaluation and verification of dose distribution in different tissues and organs. Simulations of multiple particles coupling transport are also necessary for electron accelerators and detectors, especially photon and electron coupling transport, allowing for the accurate simulation of the impact of fine structures on key parameters. In detector applications, deep penetration issues are often involved, and the variance reducing of pulse height tally can significantly improve computational efficiency. R&D have been conducted on the above issues based on the TopMC (Multi-functional Program for Neutronics Calculation, Nuclear Design and Safety Evaluation) to achieve widespread applications, including the design and application of nuclear detection and acceleration, radiation processing and nuclear medicine.

TopMC is the totally updated and extended version of SuperMC, which has been developing by the FDS Consortium for more than 20 years [13–15]. SuperMC focuses on the transport of uncharged particles, including neutron and photon, as well as the depletion, activation and shutdown dose calculations. It is characterized by high efficiency and high fidelity neutronics calculation, accurate modelling, visualized intelligent analysis and virtual simulation. These features have made SuperMC widely applicable in major nuclear engineering projects, especially in nuclear energy systems, including commercial nuclear power plants and advanced nuclear energy systems. Building on this foundation, TopMC 1.0 introduces several methods and functions to address the above key issues in nuclear technology applications by using Monte Carlo transport simulation, including the electron transport, DICOM images based voxel model, photon-electron coupling transport, pulse height tally with variance reduction. Compared to SuperMC, TopMC expands its functionality and user-friendly features, enhancing its performance and broadening application scope.

2. Methods and techniques

To expand the application of TopMC in nuclear technologies, the following key technologies have been developed in the modeling and analysis modules of TopMC: DICOM images based voxel modeling and particle tracking, electron transport and its coupling with photon and neutron, and pulse height tally with variance reduction. The details are as follows.

2.1. DICOM images based voxel model

The DCMTK (DICOM Toolkit) toolkit is used to construct and convert the DICOM image files, which is a collection of libraries and applications implementing large parts of the DICOM standard. The pixel, pitch and thickness of CT (Computed Tomography) images can be obtained, and the geometry model can be generated by the pre-process module of TopMC. To increase the efficiency of calculations and analysis, the voxel model can be coarsened by a specified proportion, and the CT values would be recalculated based on the new voxel model and a volume averaging method (see Fig. 1).

Fig. 1. CT image with different pixel visualized by the pre-process module of TopMC. (a) pixel 512 × 512, (b) pixel 128 × 128.

The range of CT values corresponds to different tissues and organs. According to the ICRU (International Commission on Radiation Units and Measurements) report, the material density and nuclide composition can be derived from the CT values [16], allowing for the automatic generation of the material of each voxel. In addition, RT (Radiation Therapy) structure files can be imported to accurately define tumours and organs instead of CT values.

Particles frequently traversing each voxel in the model can lead to low computational efficiency. To address this, the Delta-Tracking particle tracing method is employed. The transport length is sampled using a fictitious cross-section according to the following formula. Rejection sampling is then used to determine whether the particle undergoes a physical reaction or a fictitious reaction.

$$\begin{array}{c}\hfill \begin{array}{c}{\mathrm{\Sigma }}_{t,max}=\left\{{\mathrm{\Sigma }}_{t,1},{\mathrm{\Sigma }}_{t,2},\dots ,{\mathrm{\Sigma }}_{t,M}\right\}\hfill \\ {\displaystyle L=-\frac{ln\xi }{{\mathrm{\Sigma }}_{t,max}}}\hfill \end{array}\end{array}$$

where Σ_t, max is the fictitious cross-section, Σ_t, i is the macro total cross-section of material i, and ξ is a random number.

Moreover, in the presence of strong absorber, the Delta-Tracking method can experience reduced efficiency due to a large number of fictitious reactions. To address this, a coupled tracking method that combines Ray-Tracking and Delta-Tracking has been developed. This method uses Ray-Tracking when the following condition is satisfied; otherwise, it uses Delta-Tracking method. Here, ε is a constant:

$$\begin{array}{c}\hfill {\displaystyle \frac{{\mathrm{\Sigma }}_{t,i}}{{\mathrm{\Sigma }}_{t,max}}<1-\epsilon }\end{array}$$

Finally, to further accelerate the search for the voxel in which a particle is located, a quicksort algorithm is used to create an index between coordinates and voxel positions. By combining absolute and relative search methods, the voxel containing the particle can be quickly identified.

2.2. Photon-electron coupling transport

To address the characteristics of high collision frequency and low energy loss per collision in electronic transport, the condensed history method is adopted, which divides a random walk of the electron into several electron steps, and multiple collisions in each electron step are treated as one collision. Here, the energy of electrons can be used to divide the electron walk history into steps.

$$\{\begin{array}{c}\mathrm{\Delta }{s}^{\prime }=-{\int }_{E_{n-1}}^{E_n}\frac{1}{-\mathrm{d}E/\mathrm{d}s}\mathrm{d}E\\ \frac{E_n}{E_{n-1}}=k\end{array}$$

where Δs^′ is the electron step length, E_n is the energy of step n, −dE/ds is stopping power, and k is a constant.

The angular deflection in each electron step is calculated according to the Goudsmit-Saunderson theory [17],

$$\begin{array}{c}\hfill \begin{array}{c}F(s,\mu )=\sum _{l=0}^{\infty }\left({\displaystyle l+\frac{l}{2}}\right)exp(-s{G}_l)P_l\left(\mu \right)\hfill \\ G_l=2\pi N{\int }_{-1}^1{\displaystyle \frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }}}\left[1-P_l\left(\mu \right)\right]\mathrm{d}\mu \hfill \end{array}\end{array}$$

where s is the electron step length, μ is the cosine of the scattering angle, $P_l(\mu )$ is the Legendre polynomials, N is the nucleon density, and dσ/dΩ is the microscopic cross-section.

The energy loss of each step will fluctuate and is calculated by following formula,

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \frac{\mathrm{d}E}{\mathrm{d}s}=NZ\frac{2\pi e^4}{m_e{v}^2}\left[{\displaystyle ln\frac{E^2\left(\tau +2\right)}{2I^2}+f^{-}\left(\tau ,{\epsilon }_m\right)-\delta }\right]}\hfill \\ f^{-}\left(\tau ,{\epsilon }_m\right)=-{\beta }^2+\left(1-ln2\right)+\left({\displaystyle \frac{1}{8}+ln2}\right){\left({\displaystyle \frac{\tau }{\tau +1}}\right)}^2\hfill \end{array}\end{array}$$

where Z is the average atomic number, I is average ionization potential, τ is electron kinetic energy per unit mas, δ is density correction, β = v/c, v is electron velocity, and c is light speed. Furthermore, the convolution of the Landau distribution and the Gaussian distribution is used to perform energy straggling to correct the energy after multiple scattering,

$$\begin{array}{c}\hfill {\displaystyle f\left(s,\mathrm{\Delta }E\right)=\frac{1}{\sqrt{2\pi }\sigma }{\int }_{-\infty }^{+\infty }{\displaystyle f\left(s,\mathrm{\Delta }{E}^{{}^{\prime }}\right)exp(\frac{\mathrm{\Delta }E-\mathrm{\Delta }{E}^{{}^{\prime }}}{}\mathrm{2\sigma }}d\mathrm{\Delta }{E}^{{}^{\prime }}}\end{array}$$

where σ is the variance. Here, the electron energy loss due to the bremsstrahlung photon production in each step is not included in the above formula, and it will be subtracted in the bremsstrahlung process.

Electrons can emit secondary electrons through ionizations and excitation reactions, and emit photons through bremsstrahlung, transitions and positron annihilation. The photoelectric effect, Compton scattering and electron pair effect of photons can produce electrons. In the simulation of photon and electron coupling transport, if electron transport is turned off during photon transport, thick target bremsstrahlung is used to deposit electron energy in situ to produce photons; if photon transport is turned off during electron transport, photon energy is also deposited in situ. Therefore, performing the photon and electron coupling transport simulation is essential to get precise results.

2.3. Pulse height tally with variance reduction

The research and application of gamma detectors require pulse height tally to calculate parameters such as detection efficiency and response function, and application scenarios involve deep penetration problems. To improve calculation efficiency, pulse height tally with variance reduction had been developed [18]. This method involves constructing a track history tree, where the nodes of the tree are divided into reduced variance nodes and physical nodes. Each node represents a nuclear reaction event, a variance reduction event, or the termination event of a particle. The child nodes under each node are secondary particles generated by an event, and the parent and child nodes are connected through branches. The branches connecting nodes represent the tracks generated after a reaction or a variance reduction event, and their number corresponds to the number of particles.

For each event that occurs, a child node is generated under the current node. If the event is a real physical reaction, the child node is a physical node; otherwise, it is a variance reduction node. The weight changes are recorded before and after the event process in the current branch. Additionally, the energy changes of particles in each counting cell for the current branch are also recorded.

After each source particle has been transported, the track history tree can be backtracked from the last node to the first node, and be merged to obtain the pulse height contribution of the source particle. Figure 2 shows a typical node, and illustrates the method of energy deposition merging.

Fig. 2. Node of pulse height tally with variance reduction technology.

The two branches in this node represent the generation of two secondary particles through an event. The weight of the event process is changed to W₁ and W₂, and the energies deposited by secondary particles before the event is E₁ and E₂. The energy deposition of secondary particles generated by the first particle is e_p1, and the corresponding weight is W_p1; the energy depositions of secondary particles generated by the second particle are E_p21 and E_p22, with corresponding weights of W_p21 and W_p22. If the current node is a physical node, two sets of energy deposits are generated after merged, and the energy depositions and weights are shown in the following formula.

$$\begin{array}{c}\hfill \{\begin{array}{c}{E}_1+E_2+E_{p1}+E_{p21},W_1\times W_2\times W_{p1}\times W_{p21}\hfill \\ E_1+E_2+E_{p1}+E_{p22},W_1\times W_2\times W_{p1}\times W_{p22}.\hfill \end{array}\end{array}$$

If the node is a variance reduction node, three sets of energy deposition and weight will be generated.

$$\begin{array}{c}\hfill \{\begin{array}{c}{E}_1+E_{p1},W_1\ast W_{p1}\hfill \\ E_2+E_{p21},W_2\ast W_{p21}\hfill \\ E_2+E_{p22},W_2\ast W_{p22}.\hfill \end{array}\end{array}$$

For a termination node in the tree, its two child nodes will contain no data, then the calculations are simplified as following:

$$\begin{array}{cc}& \text{Physical node:}{E}_1+E_2,W_1+W_2\hfill \\ \hfill & \text{Variance reduction node:}\{\begin{array}{c}{E}_1\hfill \\ E_2\hfill \end{array},\{\begin{array}{c}{W}_1\hfill \\ W_2.\hfill \end{array}\hfill \end{array}$$

By sequentially merging the data of each node, the energy deposition and weight generated by the current particle can be obtained. This allows for the efficient calculation of the final pulse height tally using variance reduction techniques.

3. Verification and application

TopMC/SuperMC has been verified and validated by more than 2000 benchmark models and experiments, such as International Criticality Safety Benchmark Evaluation Project (ICSBEP), International Reactor Physics Benchmark Experiment (IRPhEP) and Shielding Integral Benchmark Archive and Database (SINBAD); pressurized water reactors BEAVRS [19], HM, and TCA; fast reactors IAEA-BN600 and IAEA-ADS; and fusion reactors ITER benchmark model [20–23], HCPB DEMO [24, 25] and FDS-II [26, 27]. To ascertain its validity, the performance of TopMC/SuperMC was assessed through series of tests, which include different physics processes (bremsstrahlung, scattering, ionization, and positron annihilation), sources (point, surface and volume), tallies (surface current and flux, cell flux and pulse height), energy spectra, variance reduction methods (energy splitting and roulette, weight window, dxtran sphere and exponent), and coupling transport (photon-electron and neutron-photon-electron). Moreover, a diverse array of application studies encompassing BNCT (Boron Neutron Capture Therapy) therapy, nuclear well-logging, gamma radiation detection systems, and electron accelerator simulations are performed. These validations underscored both the precision and efficacy of TopMC across various facets of nuclear technology applications.

3.1. Dose evaluation of BNCT

First, neutron beam irradiation tests are performed using both a homogeneous water phantom or water-equivalent phantom and a heterogeneous brain phantom. The thermal/fast neutron flux densities and photon absorbed doses at different positions are calculated and compared to verify the dose evaluation, and results deviations are less than 7%.

Afterwards, a clinical case is used to perform dose evaluation of BNCT. The CT images consist of 96 layers, each with a size of 30 cm × 30 cm × 28.8 cm, and a resolution of 512 × 512 pixels. A circular neutron source with a diameter of 6 cm is set above the patient’s head. The boron concentration in the blood is assumed to be 1500 ppm, and the ratio of boron concentration in the blood, tissue and tumour is 1:1:1.5. Using the conversion function of TopMC, the calculation model of TopMC and the reference program are generated automatically from the CT images and RT-Structure file. The biological doses in the human body are calculated in (Fig. 3). Results show a good agreement between TopMC and the reference program.

Fig. 3. BNCT simulation results. (a) TopMC, (b) reference program.

3.2. Nuclear logging

Rays can interact with geological materials in various ways, producing secondary gamma rays. The intensity and energy of gamma rays are related to the deceleration and absorption properties of the formation, thus the formation parameters, such as density and saturation, can be inverted from the detected gamma ray information.

The logging model calculated in this contribution is shown in Figure 4, which includes the well, water, rock, shield, NaI detectors, and a cesium source. The well is a cylinder with a diameter of 20 cm and a height of 60 cm; the diameter of the rock layer is 100 cm; the radius of the shield is 3.7 cm and the height is 36 cm; the radiuses of NaI detectors are 2 cm and 1 cm respectively, with heights of 5 cm and 2 cm. The energy of photon source is 0.662 MeV.

Fig. 4. Model of nuclear logging. (a) Sectional view, (b) logging instrument.

The pulse height spectra of two detectors are calculated, and the normalized pulse height spectra are shown in Figure 5. Detector A is closer to the source. As seen in the figure, the pulse height spectrum calculated by TopMC is in good agreement with that of the reference program, with relative errors less than 2%. In the range with large energy deposition values, the relative error between TopMC and the reference program increases due to statistical fluctuations. However, due to the small absolute values of these calculation results, it does not significantly affect the simulation of parameters, such as detector detection efficiency.

Fig. 5. Energy distribution of pulses created in the detector. (a) Detector A, (b) detector B.

3.3. Gamma detector response

A high purity germanium detector is modelled as shown in Figure 6. The manufacturer specified the crystal diameter as 58.9 mm, the length as 54.3 mm, the thickness of the aluminum layer and the inactive germanium layer as 1.27 mm and 0.7 mm, respectively.

Fig. 6. Model of HPGe detector.

The pulse height spectrum and response functions in the detector are shown in Figures 7 and 8 respectively, for the 0.662 MeV and 10 MeV γ incidents. It can be seen that the results of TopMC are in good agreement with the reference values.

Fig. 7. Deposition spectrum of gamma-ray for HPGe detector. (a) 0.662 MeV (b) 10 MeV.

Fig. 8. Response function of gamma-ray for HPGe detector. (a) 0.662MeV, (b) 10 MeV.

3.4. Electron accelerator head design

The design of electron accelerators is crucial for the development of radiotherapy equipment. The model of an accelerator head is shown in Figure 9, which consists of a primary collimator, a flattening filter, an ionization chamber, upper and lower jaws. The primary collimator and jaws are mainly composed of tungsten, while the flattening filter is mainly composed of lead. A water phantom is placed 1 meter away from the target, and 2 sets of detectors are positioned in the water phantom. The angle between the 2 sets of detectors are 45°, and the detectors are balls with a radius of 1 cm. Figure 10 shows the energy spectrum of the photon source.

Fig. 9. Model of accelerator.

Fig. 10. Energy spectrum of X-ray.

The energy depositions in detectors are calculated using 50 million particle history. Figure 11 shows the calculation results of the two sets of detectors. The first set of detectors is positioned facing the opening position of the tungsten door, and the second set of detectors is at a 45° angle relative to the first set of detectors. The results exhibit symmetrical characteristics about the central axis, with good uniformity within a 2 cm range at the central position. There is a slight decrease in energy deposition at the central position, which is due to the thicker central position of the flattening filter.

Fig. 11. Deposition energy of detectors.

3.5. Radiation therapy

Microspheres are a type of granule used in brachytherapy, which have a wide range of diagnostic, therapeutic, and research applications. In this case, 120 microspheres are used (the red part in Fig. 12a), each with a radius of 0.03 cm, and are located within a cylinder. The cylinder has a height of 3.018 cm and a radius of 0.3 cm. Detectors are placed at distances of 0 cm, 0.075 cm, 0.15 cm, and 0.3 cm from the surface of the cylinder (the blue part in Fig. 12b), with 56 detectors at each distance. These 56 detectors are uniformly distributed in 8 radial and 7 axial directions, with a spacing of 0.15 cm between detectors in the axial direction. The radius of each detector is 0.01 cm. The mass density of the microspheres is 3 g/cm³, and they are composed of 10.5% praseodymium, 63.2% oxygen, 10.5% aluminum, and 15.8% silicon, with the remaining being water. Electrons are released from the 120 microspheres with equal probability.

Fig. 12. Model of radioactive microspheres. (a) Side view, (b) top view.

The energy deposition in each detector was calculated. For comparison, multiple detectors aligned in a straight line in the radial direction, as well as detectors at axial distances of 0 cm and 0.3 cm from the surface of the cylinder, were selected. The relative errors between TopMC and the reference program are consistently within 1%, and the overall trends are consistent (Tabs. 1–3).

4. Conclusion

A series of methods and techniques have been developed in TopMC, including DICOM images based voxel model, photon-electron coupling transport, pulse height tally with variance reduction, to facilitate the broad application in nuclear technology. The correctness and effectiveness have been verified through cases studies in different aspects of nuclear technology applications.

Looking ahead, the enhancement of TopMC will focus on expanding its functionality and improving user-friendly features to broaden its application scope and enhance user satisfaction. The extension of the energy spectrum coverage through the implementation of a Class-II condensed history algorithm for electron transport is envisaged, aiming at meeting the precise demands of biomedical dosimetry, detection simulations, and other applications. New capabilities such as large perturbation, continue-run, sensitivity and uncertainty analysis will be added to improve the overall efficiency.

Besides, advanced Artificial Intelligence (AI) algorithms are being studied and integrated in TopMC. By coupling with advanced AI algorithms, TopMC will focus on the intelligent design and prediction, semantics based modeling and visual analysis of nuclear systems, aiming to revolutionize nuclear design and safety evaluation. A new AI-powered version can be anticipated to be released in the upcoming year.

Acknowledgments

We thank other FDS consortium members, and authors and contributors of references cited in this work.

Funding

This work was supported by the SuperSafety Science & Technology Co., Ltd., FDS consortium (YF-AC1.3-TopMC1.0).

Conflicts of interest

The authors declare that they have no competing interests to report.

Data availability statement

No data are associated with this article.

Author contribution statement

The study conception and design is contributed by Prof. Yican Wu. The first draft of the manuscript was written by Shanqi Chen and Jieqiong Jiang, all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

References

Shanqi Chen is currently senior engineer of SuperSafety Science & Technology Co., Ltd. His research interests include nuclear safety and nuclear technology application.

BinWu is currently chief engineer of SuperSafety Science & Technology Co., Ltd. His research interests include advanced nuclear system and nuclear technology cross-application.

Zhengyun Dong is currently senior engineer of SuperSafety Science & Technology Co., Ltd. His research interests include neutron transport theory and methods, advanced nuclear system.

Shaoheng Zhou is currently deputy chief engineer of International Academy of Neutron Science. His research interests include neutronics simulation and analysis.

Guangyao Sun is currently senior researcher of International Academy of Neutron Science. His research interests include neutron transport theory and methods, advanced nuclear system and nuclear technology application.

Li He is the doctoral student of University of Science and Technology of China. Her research interests include the application of AI technology in the field of nuclear energy and nuclear technology.

Dong Yao is doctoral student of University of Science and Technology of China. His research interests include the methods of fast Monte Carlo and the applications of dose evaluation in the radiotherapy.

Shengpeng Yu is currently senior researcher of International Academy of Neutron Science. His research interests include neutron transport theory and methods, advanced nuclear system and nuclear technology application.

Quan Gan is currently senior researcher of International Academy of Neutron Science. His research interests include neutron transport theory and methods, advanced nuclear system and nuclear technology application.

Lijuan Hao is currently senior researcher of International Academy of Neutron Science. Her research interests include neutron transport theory and methods, advanced nuclear system and nuclear technology application.

Jing Song is currently senior researcher of International Academy of Neutron Science. Her research interests include neutron transport theory and methods, advanced nuclear system and nuclear technology application.

Pengcheng Long is currently senior researcher of FDS Consortium. His research interests include neutron transport theory and methods, neutronics simulation and analysis.

Yazhou Li is currently senior researcher of FDS Consortium. His research interests include neutron transport theory and methods, reactor safety assessment and energy policy.

Jieqiong Jiang is currently senior researcher of FDS Consortium. Her research interests include neutron transport theory and methods, advanced nuclear system and nuclear technology application.

Fang Wang is currently senior researcher of FDS Consortium. Her research interests include neutron transport theory and methods, advanced nuclear system and nuclear technology application.

Liqin Hu is currently senior researcher of FDS Consortium. Her research interests include neutron transport theory and methods, advanced nuclear system and nuclear technology application.

Yican Wu has been working on theory, modeling, experiments, design, and analysis of advanced nuclear system, and is involved in both basic and applied research mainly concentrating on the areas of nuclear science and engineering, radiation medical physics and technology, computer simulation and software engineering, and other interdisciplinary research.

All Tables

All Figures

	Fig. 1. CT image with different pixel visualized by the pre-process module of TopMC. (a) pixel 512 × 512, (b) pixel 128 × 128.
In the text

	Fig. 2. Node of pulse height tally with variance reduction technology.
In the text

	Fig. 3. BNCT simulation results. (a) TopMC, (b) reference program.
In the text

	Fig. 4. Model of nuclear logging. (a) Sectional view, (b) logging instrument.
In the text

	Fig. 5. Energy distribution of pulses created in the detector. (a) Detector A, (b) detector B.
In the text

	Fig. 6. Model of HPGe detector.
In the text

	Fig. 7. Deposition spectrum of gamma-ray for HPGe detector. (a) 0.662 MeV (b) 10 MeV.
In the text

	Fig. 8. Response function of gamma-ray for HPGe detector. (a) 0.662MeV, (b) 10 MeV.
In the text

	Fig. 9. Model of accelerator.
In the text

	Fig. 10. Energy spectrum of X-ray.
In the text

	Fig. 11. Deposition energy of detectors.
In the text

	Fig. 12. Model of radioactive microspheres. (a) Side view, (b) top view.
In the text