1 Introduction
As a result of the acceleration or especially the deceleration of charged particles caused by the electrostatic field of atoms, electrons (or positrons) emit bremsstrahlung. Bremsstrahlung emission is of basic interest in many fields. Examples of applications where the bremsstrahlung interaction plays a vital role are, X-ray generators, radiation shielding, medical physics, and electron probe microanalysis [1]. In medical physics, the bremsstrahlung phenomenon is very important in both diagnostics radiology and radiotherapy because the fundamental characteristics of X-ray machine is described by the spectral distribution of photon beam originating from the target. Bremsstrahlung photons have much higher penetration than other charged particles and thus have large contribution in dose in deep regions of the irradiated samples, therefore an accurate Monte Carlo simulation of bremsstrahlung is required in a general purpose Monte Carlo code such as SuperMC [2].
SuperMC is a general purpose, easy-to-use Monte Carlo simulation program, developed by FDS Team [2-5], for modeling and simulating nuclear facilities. The released version of SuperMC can perform the simulation of neutrons, photons and coupled neutron and photon transport while Monte Carlo simulation of electron transport has been under development and will be released soon, this research work is part of the development of electron transport. The main technical features of SuperMC include hybrid MC-deterministic methods and the adoption of advanced information technologies, while the main usability features include automatic modeling of geometry and physics, visualization and virtual simulation and cloud computing services. SuperMC has been applied to fusion reactor FDS I[6], FDS II[7], FDS III[8], ITER [9-10] and other reactor studies[11] and is also designed for MC-based dose calculation engine of ARTS (Accurate/Advanced Radiotherapy System) [12].
Bremsstrahlung process is described by the differential cross section (DCS) in energy and direction of emitted photons and the direction of outgoing electrons. Atomic DCSs of bremsstrahlung emission by electrons are derived from conventional quantum electrodynamics methods and comprehensive review on the DCSs is described by Koch and Motz [13] and Tsai [14]. Efforts have been made to use accurate methods for sampling the photon energy and, especially, the angular distribution of bremsstrahlung photons, but still, there is need to use simple, efficient and accurate methods to sample energy and angular distribution. The 2BS formula of Koch and Motz was implemented by Alex F. Bielajew et al. [15] in EGS4 [16] for the angular distribution of bremsstrahlung photons. Acosta et al. [17] developed a parameterization with adjustable parameters determined from the fitting of shape functions provided by Kissel et al. [18] based on the calculations of partial wave methods by Tseng et al. [19]. Geant4 [20] uses a parameterization of Tsai [14] DDCS for sampling photon angular distribution, while Rodrigues et al. [21] implemented 2BN formula of Koch and Motz, in Geant4 as more accurate option for low energy region. EGSnrc [22] incorporated modified version of 2BS formula which converges to 2BN for the treatment of low energy region.
In this article, an algorithm is described for Monte Carlo simulation of electron bremsstrahlung with kinetic energy range of 1 keV–10 GeV and for Z = 1–100. We also describe the efficient and accurate methods by which the angular distribution and energy of bremsstrahlung photons are sampled. Photon energy is sampled according to most reliable Seltzer and Berger [23] scaled DCSs and for angular distribution of bremsstrahlung photons, a novel hybrid model is developed. In Section 2, DCSs used for sampling energy and angular distribution of emitted photons, and corresponding random sampling algorithms, are presented. Results of the newly developed MC simulation algorithm are presented in Section 3.
2 Monte Carlo simulation of electron bremsstrahlung
Bremsstrahlung process is described by the atomic differential cross section, in the direction θ and energy K of the emitted photons and the direction of the outgoing electrons. We consider the process in which electron of energy E (kinetic energy T) is accelerated or decelerated by the screened coulomb field of the target atom and emit bremsstrahlung photons of energy K in the interval [0, E] in the direction of polar angle θ. Angular deflections of the incident electron due to bremsstrahlung are already accounted for by electric scattering in SuperMC and assumed that the direction of electron is not altered by the emission of bremsstrahlung photon.
In conventional MC simulations, bremsstrahlung emission is modeled using analytical DCSs obtained from simple approximations and this was done to minimize the computer memory requirements. With the advancement in computer technology, however, it is now convenient to use the combination of numerical DCSs and analytical formulae, which are more accurate and reliable at present. In this study we used the models and numerical data which are considered as most reliable at present and their form suit for Monte Carlo simulation.
2.1. Sampling photon energy
In this section we describe the algorithm to sample photon energy, K, from scaled energy-loss total DCSs (electron-electron, electron-nucleus), χ(Z, T, w=K/T), tabulated by Seltzer and Berger [23] which are based on the combination of numerical partial wave results and high energy theory.
where, w is reduced photon energy,
For sampling photon energy, a simple and efficient composite sampling procedure is implemented. The distribution function
where Kc is threshold photon energy. As can be seen from Fig. 1, χ(Z, T, w) is relatively flat over the interval [Kc/T, 1] and can be used for rejection sampling while direct (inverse-transform) sampling method is used to sample w from the distribution
The algorithm is then
i. Calculate χmax (find maximum value of χ in K/T grid for given T and Z).
ii. Generate two random numbers r1 and r2 over the interval [0, 1].
iii. Sample w from the distribution w−1 as w = exp[ln(Kcr1/T)]
iv. If r2 ≤ χ(Z, T, w)/χmax, accept w; else, go to Step ii.
Efficiency and accuracy of this sampling algorithm is discussed in Section 3.1. A number of comparisons between measured and calculated scaled differential cross section in photon energy have been done [23]. These comparisons show that Seltzer and Berger’ scaled differential cross sections are the most reliable representation of bremsstrahlung at present.
2.2. Hybrid model for photon angular distribution
A new hybrid model for sampling the direction of bremsstrahlung photons has been developed which characterizes the explicit treatment of photon angular distribution for low and high energy incident electrons. High energy analytical formula of DDCS (differential in energy and angle of bremsstrahlung photon) described by Tsai (1974) [14] based on high energy theory of Davies, Bethe, Maximon and Olsen (DBMO) [24, 25] with emphasis on form factor screening corrections is used to sample the direction of emitted photons. However at low energies, 2BS formula of Koch and Motz and Tsai’s DDCS give rise to deviations which could reach 50 Degrees on the most probable angle [21]. This is of particular importance when simulating bremsstrahlung emission in very thin targets and addressing topics like medical imaging which works with X rays of about a few keV to a few hundred keV. Therefore an accurate treatment of the angular distribution of the bremsstrahlung photons, needs to be done for low energy incident electrons. Thus a new, accurate and fairly simple analytical distribution function has been developed for photon angular sampling with adjustable parameters determined from the fitting of the numerical values of shape functions tabulated by Kissel et al. [18].
2.2.1 High energy bremsstrahlung angular distribution
In this study, for high energy electrons, the angular distribution of emitted photons has been sampled accurately, according to the “exact form” of Tsai’s double differential cross section (DDCS) which includes atomic form factors. Thus the sampling is fully consistent with theory. It should be noted that the parameterization of the angular distribution instead of “exact form” may produce uncertainties and the accuracy of the parameterization needs to check against, many, possible electron and photon energies and materials. The double differential cross section, energy angle distribution of electron bremsstrahlung, described by Tsai [14] is
where y= K/E and u= θE /me (me = 0.511 MeV), α=1/137, Z is atomic number, and fc is coulomb correction. For the calculation of the function X, Tsai [14] provided simple expressions:
where a =184.15e−1/2Z−1/3/me, and a'=1194e−1/2Z−2/3/me are atomic parameters. The minimum momentum transfer used to calculate X is defined as tmin = {yme2(1+u2)/[2E(1−y)]} 2.
To sample the direction of bremsstrahlung photon from Tsai’s DDCS, the angular distribution function is defined as
where
As can be seen from Fig. 2, R(u) can be employed as rejection function and therefore a composite approach of direct sampling and rejection sampling is used to sample u from Eq. (5). Here u is directly sampled from Eq. (7) and then R(u) is used as rejection function. In order to employ rejection sampling, the location of the maximum,
Since true maxima is 3% (at the most) higher than the estimated maxima, therefore, this small error can be corrected by introducing an additional multiplicative factor to the estimated maxima, which is 1.03 (3%) for entire energy range.
The algorithm is then
i. Calculate Rmax from Eq. (8) for given Z, E and K.
ii. Generate two random numbers r1 and r2 over the interval [0, 1].
iii. Calculate normalization constant (Nu) of h(u) over [0,πE/me] as Nu= 1/arctan(πE/me)
iv. Sample u from h(u) as u = tan(r1/ Nu)
v. If r2 ≤ R(u)/Rmax, accept u; else, go to step ii.
The proposed sampling algorithm requires several logarithmic calculations in step (i) and in step (v) which makes it slightly slow, but holds high sampling efficiency.
2.2.2. Low energy bremsstrahlung angular distribution
Numerical values of the shape functions based on the partial wave methods [19], have been published by Kissel et al. [18] for 144 benchmark values of electron kinetic energy (1–500 keV), atomic number (Z=1–92) and reduced photon energy. For electron kinetic energies of <500 keV, the angular distribution is sampled from an analytical distribution function with adjustable parameters determined from the fitting of Kissel et al. [18] data. Kissel and coworkers have provided the analytical form of the shape functions in terms of Legendre polynomial, but their analytical expression does not suit for the sampling of polar angle θ. Therefore, in this study, a fairly simple parameterization of the partial wave double differential cross section with adjustable parameters is proposed as:
where θ = 0–π is the polar angle of bremsstrahlung photon; and A, B and C are adjustable parameters, which are determined by fitting to the numerical values of shape functions for 144 benchmark cases and stored in tables for rapid retrieves of their values. Some differences were found between the data and fits for lower photon energies, but still the fits are accurate enough for sampling angular distribution of bremsstrahlung photons. Comparisons between analytical function determined by Eq. (9) with the original shape functions tabulated by Kissel et al. [18] are shown in Fig. 4.
The best side of the proposed analytical form is that ρ has analytical integral and the corresponding cumulative probability distribution function has analytical inverse, therefore direct sampling procedure (inverse transform method) was applied to sample θ. The normalization constant NA is constant chosen to normalize the integral of ρ(θ) in the overall range of θ to unity.
Since the calculation of NA involves five exponential calculations therefore in order to make algorithm efficient, we calculated NA for 144 benchmarked values and stored in tables along with adjustable parameters for their rapid retrieves, for particular T, Z, and K/T. The polar angle is then sampled as
where r is random number uniformly distributed over [0, 1]. This algorithm requires one log(log) calculation and two exponential calculations in Eq. (11). As it uses direct sampling method, the algorithm is more efficient than the sampling method proposed by Acosta et al. [17] which uses four random number, on average, to sample each value of polar angle.
2.2.3. Comparison between Kissel shape functions, Tsai DDCS and fitting function
Figure 4 shows comparisons of the polar angle distribution of bremsstrahlung photons using Tsai DDCS, Kissel et al. data and the analytical distribution function, Eq. (9), with reduced photon energy, K/T = 0.6, for 50 keV, 100 keV and K/T = 0.8 for 500 keV electrons incident on Silver and Gold. For low energy about 50 keV and 100 keV there is a significant deviation in the most probable angle between Tsai’s DDCS and Kissel’s data/Analytical function (equation 9). With an increase in the electron energy all three approaches tend to overlap and present a good agreement. Since Kissel et al. data is very accurate for electron energies below 500 keV, therefore we recommended to use analytical distribution function given by Eq. (9) in this energy range. For energies of > 500 keV, Tsai DDCS is a good choice as it retains the high sampling efficiency and physics accuracy
3 Results and Discussion
3.1. Photon Energy, sampling efficiency and accuracy
3.1.1 Sampling Accuracy
The accuracy of the sampling algorithm is determined by comparing sampled photon energy distribution with interpolated values of scaled DCSs (Seltzer and Berger) [23]. Two examples are presented. In figure 5a, photon energy distribution is plotted for 100 keV electrons incident on Silver with Kc=10 keV. Similar comparison is plotted in Fig.5(b) except for 100 MeV electrons. The smooth curves are from equation (2) and histograms are sampled values of photon energies from MC simulation. Fig. 5 verified that the photon energy sampling algorithm works as expected. All plots in Fig. 5 are normalized to unity.
3.1.2. Sampling efficiency
The sampling efficiency, P, of the photon energy sampling algorithm is given in Table 1. The sampling efficiency is defined as the ratio of the accepted sampled values to the total number of trial. We have calculated the efficiency in terms of incident electron energy. Although the sampling efficiency severely depends on the threshold photon energy and also less on Z, the qualitative behavior remains the same for other values of Kc. We are therefore more interested in calculating P for different electron energies while keeping Kc = 1keV and Z = 43 (Ag).
| E (MeV) | 0.05 | 0.1 | 0.5 | 1 | 10 | 50 | 100 | 500 | 1000 | 10000 |
|---|---|---|---|---|---|---|---|---|---|---|
| P (%) | 90 | 84 | 79 | 80 | 87 | 91 | 92 | 94 | 95 | 95 |
From Table 1, the sampling efficiency in energy region of E < 100 keV and E > 1 MeV is higher than that in 100 keV–1 MeV, which is obvious from the shape of the energy spectrum, as shown in Fig. 1. Even for higher incident electron energies, the energy spectrum is flat except at end points of K/T.
3.2. Photon Angular Distribution, sampling accuracy and efficiency
This section concerns the verification that the angular sampling algorithms work as expected. The verification is performed by comparing the sampled angular distribution with the theoretical expressions. In Fig. 6, the sampled angular distribution from Tsai’ DDCS are plotted with the theoretical expression given in Eq. (5). The accuracy of angular sampling using analytical function in Eq. (9) is determined by fitting of Kissel’s data, as plotted in Fig. 7. The smooth curves are from theoretical expressions, and the histograms are sampled angular distribution from Monte Carlo simulation. All plots are normalized to unity. Error bars in Figs. 6 and 7 are uncertainty (±3σ) in each bin.
The sampling efficiency of our proposed angular sampling algorithm based on Tsai’s DDCS and angular sampling algorithm proposed by Bielajew et al. (1989) implemented in EGS5 [26] and Geant4 [20] are plotted in Fig. 8. The sampling efficiency is plotted for incident electron energies from 500 keV to 1000 MeV for silver and gold, and K/T= {0,0.5,1}. One can realize that the sampling efficiency based on Tsai’s DDCS is generally very good for both Z = 47 and Z = 79 and about 70% of average for K/T ~ 0,0.5 and reduces by few percent when K/T =1. The sampling efficiency for Z = 79 decreases, at most, to 55% at 1 MeV, only when K/T = 1 it retains 62% up to 1000 MeV. So, the sampling technique employed here is more efficient in the useful photon energy range (lower reduced photon energies) than the sampling technique used to sample the angular distribution of photons, reported by Alex F. Bielajew et al. (1989) [15] using 2BS formula of Koch and Motz [13] which is implemented in EGS5 [26] and Geant4 [20].
4 Conclusion
It can be concluded that the proposed algorithm provides accurate and efficient methods for sampling energy and direction of bremsstrahlung photons. The efficiency of sampling photon energy is more than 80% with this new algorithm. Hybrid model provides very accurate treatment of photon angular distribution by low as well as high energy incident electrons. Angular sampling which is based on Tsai DDCS is very efficient for useful photon energy range, compared to other available angular sampling techniques in high energy region (Alex F. Bielajew et al. 1989). Distribution function used for photon angular sampling using Tsai DDCS has contribution of atomic form factors and hence explains physics phenomena very well especially in high energy region. The parameterization of the shape functions given by Eq.(9) provides a very accurate fit to the Kissel’s benchmark data (Fig. 4) therefore simulated angular distribution is fully consistent with the most accurate and reliable theory at present and the sampling algorithm is more efficient than the previously developed method by Acosta et al. [24], which uses four random numbers, on average, to sample one photon energy.
It is worth mentioning that the approaches used for sampling energy and direction of bremsstrahlung photons can be used for a wide range of incident electron energies and materials and hence suitable for general purpose Monte Carlo simulations.
Monte Carlo simulation of bremsstrahlung emission by electrons
. Radiat Phys Chem, 2006, 75:1201-1219. DOI: 10.1016/j.radphyschem.2005.05.008.CAD-based Monte Carlo program for integrated simulation of nuclear system SuperMC
. Ann Nucl Energy, 2015, 82:161-168. DOI: 10.1016/j.anucene.2014.08.058.CAD-based interface programs for fusion neutron transport simulation
. Fusion Eng Des, 2009, 84:1987-1992. DOI: 10.1016/j.fusengdes.2008.12.041.Benchmarking of MCAM 4.0 with the ITER 3D model
. Fusion Eng Des, 2007, 82:2861-2866. DOI: 10.1016/j.fusengdes.2008.12.041.A Discrete ordinates nodal method for one-dimensional neutron transport calculation in curvilinear geometries
. Nucl Sci Eng, 1999, 133(3):350-357.Conceptual design of the fusion-driven subcritical system FDS-I
. Fusion Eng Des, 2006, 81, Part B:1305-1311. DOI: 10.1016/j.fusengdes.2005.10.015.Conceptual design of the China fusion power plant FDS-II
. Fusion Eng Des, 2008, 83(10-12):1683-1689. DOI: 10.1016/j.fusengdes.2008.06.048.Fusion-based hydrogen production reactor and its material selection
. J Nucl Mater, 2009, 386-388:122-126. DOI: 10.1016/j.jnucmat.2008.12.075.Conceptual design and testing strategy of a dual functional Lithium-Lead test blanket module in ITER and EAST
. Nucl Fusion, 2007, 47(11):1533-1539.Design analysis of the China dual-functional lithium lead (DFLL) test blanket module in ITER
. Fusion Eng Des, 2007, 82: 1893-1903. DOI: 10.1016/j.fusengdes.2007.08.012.Progress in fusion-driven hybrid system studies in China
. Fusion Eng Des, 2002, 63-64:73-80. DOI: 10.1016/S0920-3796(02)00239-9.Development of accurate/advanced radiotherapy treatment planning and quality assurance system (ARTS)
. Chinese Phys C, 2008, 32(Suppl. II):177-182.Bremsstrahlung cross section formulas and related data
. Rev Mod Phys, 1959, 31:920-955. DOI: 10.1103/RevModPhys.31.920.Pair production and bremsstrahlung of charged leptons
. Rev Mod Phys, 1974, 46:815-851. DOI: 10.1103/RevModPhys.46.815.The EGS4 code system
.Monte Carlo simulation of bresmsstrahlung emission by electrons
. Appl Phys Lett, 2002, 80:3228-3230. DOI: 10.1063/1.1473684.Shape functions for atomic field bremsstrahlung from electrons of kinetic energy 1-500 keV on selected neutral atoms 1≤Z≤92
. At Dta Nucl Data Tables, 1983, 28:381-460. DOI: 10.1016/0092-640X(83)90001-3.Electron bremsstrahlung angular distribution in the 1-500 keV energy range
. Phys Rev A, 1979, 19:187-195. DOI: 10.1103/PhysRevA.19.187.Geant4- a simulation toolkit
. Nucl Instrum Meth A, 2003, 506:250-303. DOI: 10.1016/S0168-9002(03)01368-8.Geant 4 applications and developments for medical physics experiments
. IEEE Trans Nucl Sci, 2004, 51:1412-1419. DOI: 10.1109/TNS.2004.832314.Bremsstrahlung energy spectra from electrons with kinetic energy 1 keV-10 GeV incident on screened nuclei and orbital electrons of neutral atoms with Z = 1-100
. At. Data Nucl. Data Tables, 1986, 35:345-418. DOI: 10.1016/0092-640X(86)90014-8.Theory of bremsstrahlung and pair production. II. Integral cross section for pair production
. Phys Rev, 1954, 93:788-795. DOI: 10.1103/PhysRev.93.788.Photon and electron polarization in high-energy bremsstrahlung and pair production with screening
. Phys Rev, 1989, 114:887-904. DOI: 10.1103/PhysRev.114.887.The EGS5 code system
.