diff --git a/docs/source/methods/photon_physics.rst b/docs/source/methods/photon_physics.rst index 4d58ee0db2..028a71de72 100644 --- a/docs/source/methods/photon_physics.rst +++ b/docs/source/methods/photon_physics.rst @@ -85,26 +85,26 @@ momentum transfer is traditionally expressed as .. math:: :label: momentum-transfer - x = \kappa \alpha \sqrt{1 - \mu} + x = a k \sqrt{1 - \mu} -where :math:`\alpha` is the ratio of the photon energy to the electron rest -mass, and the coefficient :math:`\kappa` can be shown to be +where :math:`k` is the ratio of the photon energy to the electron rest +mass, and the coefficient :math:`a` can be shown to be .. math:: - :label: kappa + :label: omega - \kappa = \frac{m_e c^2}{\sqrt{2}hc} \approx 29.14329, + a = \frac{m_e c^2}{\sqrt{2}hc} \approx 29.14329~\unicode{x212B}, where :math:`m_e` is the mass of the electron, :math:`c` is the speed of light in a vacuum, and :math:`h` is Planck's constant. Using :eq:`momentum-transfer`, -we have :math:`\mu = 1 - [x/(\kappa\alpha)]^2` and :math:`d\mu/dx^2 = --1/(\kappa\alpha)^2`. The probability density in :math:`x^2` is +we have :math:`\mu = 1 - [x/(ak)]^2` and :math:`d\mu/dx^2 = +-1/(ak)^2`. The probability density in :math:`x^2` is .. math:: :label: coherent-pdf-x2 p(x^2) dx^2 = p(\mu) \left | \frac{d\mu}{dx^2} \right | dx^2 = \frac{2\pi - r_e^2 A(\bar{x}^2,Z)}{(\kappa\alpha)^2 \sigma(E)} \left ( + r_e^2 A(\bar{x}^2,Z)}{(ak)^2 \sigma(E)} \left ( \frac{1 + \mu^2}{2} \right ) \left ( \frac{F(x, Z)^2}{A(\bar{x}^2, Z)} \right ) dx^2 where :math:`\bar{x}` is the maximum value of :math:`x` that occurs for @@ -113,7 +113,7 @@ where :math:`\bar{x}` is the maximum value of :math:`x` that occurs for .. math:: :label: xmax - \bar{x} = \kappa \alpha \sqrt{2} = \frac{m_e c^2}{hc} \alpha, + \bar{x} = a k \sqrt{2} = \frac{m_e c^2}{hc} k, and :math:`A(x^2, Z)` is the integral of the square of the form factor: @@ -154,6 +154,8 @@ section. The complete algorithm is as follows: 6. If :math:`\xi_2 < (1 + \mu^2)/2`, accept :math:`\mu`. Otherwise, repeat the sampling at step 3. +.. _incoherent-sampling: + Incoherent (Compton) Scattering ------------------------------- @@ -166,12 +168,11 @@ the two authors who discovered it: .. math:: :label: klein-nishina - \frac{d\sigma_{KN}}{d\mu} = \pi r_e^2 \left ( \frac{\alpha'}{\alpha} \right - )^2 \left [ \frac{\alpha'}{\alpha} + \frac{\alpha}{\alpha'} + \mu^2 - 1 - \right ] + \frac{d\sigma_{KN}}{d\mu} = \pi r_e^2 \left ( \frac{k'}{k} \right)^2 \left + [ \frac{k'}{k} + \frac{k}{k'} + \mu^2 - 1 \right ] -where :math:`\alpha` and :math:`\alpha'` are the ratios of the incoming and -exiting photon energies to the electron rest mass energy equivalent (0.511 MeV), +where :math:`k` and :math:`k'` are the ratios of the incoming and exiting +photon energies to the electron rest mass energy equivalent (0.511 MeV), respectively. Although it appears that the outgoing energy and angle are separate, there is actually a one-to-one relationship between them such that only one needs to be sampled: @@ -179,32 +180,31 @@ only one needs to be sampled: .. math:: :label: compton-energy-angle - \alpha' = \frac{\alpha}{1 + \alpha(1 - \mu)}. + k' = \frac{k}{1 + k(1 - \mu)}. -Note that when :math:`\alpha'/\alpha` goes to one, i.e., scattering is elastic, -the Klein-Nishina cross section becomes identical to the Thomson cross -section. In general though, the scattering is inelastic and is known as Compton -scattering. When a photon interacts with a bound electron in an atom, the -Klein-Nishina formula must be modified to account for the binding effects. As in -the case of coherent scattering, this is done by means of a form factor. The -differential cross section for incoherent scattering is given by +Note that when :math:`k'/k` goes to one, i.e., scattering is elastic, the +Klein-Nishina cross section becomes identical to the Thomson cross section. In +general though, the scattering is inelastic and is known as Compton scattering. +When a photon interacts with a bound electron in an atom, the Klein-Nishina +formula must be modified to account for the binding effects. As in the case of +coherent scattering, this is done by means of a form factor. The differential +cross section for incoherent scattering is given by .. math:: :label: incoherent-xs \frac{d\sigma}{d\mu} = \frac{d\sigma_{KN}}{d\mu} S(x,Z) = \pi r_e^2 \left ( - \frac{\alpha'}{\alpha} \right )^2 \left [ \frac{\alpha'}{\alpha} + - \frac{\alpha}{\alpha'} + \mu^2 - 1 \right ] S(x,Z) + \frac{k'}{k} \right )^2 \left [ \frac{k'}{k} + \frac{k}{k'} + \mu^2 - 1 + \right ] S(x,Z) where :math:`S(x,Z)` is the form factor. The approach in OpenMC is to first sample the Klein-Nishina cross section and then perform rejection sampling on the form factor. As in other codes, `Kahn's rejection method`_ is used for -:math:`\alpha < 3` and a direct method by Koblinger_ is used for :math:`\alpha -\ge 3`. The complete algorithm is as follows: +:math:`k < 3` and a direct method by Koblinger_ is used for :math:`k \ge 3`. +The complete algorithm is as follows: -1. If :math:`\alpha < 3`, sample :math:`\mu` from the Klein-Nishina cross - section using Kahn's rejection method. Otherwise, use Koblinger's direct - method. +1. If :math:`k < 3`, sample :math:`\mu` from the Klein-Nishina cross section + using Kahn's rejection method. Otherwise, use Koblinger's direct method. 2. Calculate :math:`x` and :math:`\bar{x}` using :eq:`momentum-transfer` and :eq:`xmax`, respectively. @@ -215,19 +215,370 @@ the form factor. As in other codes, `Kahn's rejection method`_ is used for Doppler Energy Broadening +++++++++++++++++++++++++ -LA-UR-04-0487_ and LA-UR-04-0488_ +Bound electrons are not at rest but have a momentum distribution that will +cause the energy of the scattered photon to be Doppler broadened. More tightly +bound electrons have a wider momentum distribution, so the energy spectrum of +photons scattering off inner shell electrons will be broadened the most. +In addition, scattering from bound electrons places a limit on the maximum +scattered photon energy: + +.. math:: + :label: max-energy-out + + E'_{\text{max}} = E - E_{b,i}, + +where :math:`E_{b,i}` is the binding energy of the :math:`i`-th subshell. + +Compton profiles :math:`J_i(p_z)` are used to account for the binding effects. +The quantity :math:`p_z = {\bf p} \cdot {\bf q}/q` is the projection of the +initial electron momentum on :math:`{\bf q}`, where the scattering vector +:math:`{\bf q} = {\bf p} - {\bf p'}` is the momentum gained by the photon, +:math:`{\bf p}` is the initial momentum of the electron, and :math:`{\bf p'}` +is the momentum of the scattered electron. Applying the conservation of energy +and momentum, :math:`p_z` can be written in terms of the photon energy and +scattering angle: + +.. math:: + :label: pz + + p_z = \frac{E - E' - EE'(1 - \mu)/(m_e c^2)}{-\alpha \sqrt{E^2 + E'^2 - + 2EE'\mu}}, + +where :math:`\alpha` is the fine structure constant. The maximum momentum +transferred, :math:`p_{z,\text{max}}`, can be calculated from :eq:`pz` using +:math:`E' = E'_{\text{max}}`. The Compton profile of the :math:`i`-th electron +subshell is defined as + +.. math:: + :label: compton-profile + + J_i(p_z) = \int \int \rho_i({\bf p}) dp_x dp_y, + +where :math:`\rho_i({\bf p})` is the initial electron momentum distribution. +:math:`J_i(p_z)` can be interpreted as the probability density function of +:math:`p_z`. + +The Doppler broadened energy of the Compton-scattered photon can be sampled by +selecting an electron shell, sampling a value of :math:`p_z` using the Compton +profile, and calculating the scattered photon energy. The theory and methods +used to do this are described in detail in LA-UR-04-0487_ and LA-UR-04-0488_. +The sampling algorithm is summarized below: + +1. Sample :math:`\mu` from :eq:`incoherent-xs` using the algorithm described in + :ref:`incoherent-sampling`. + +2. Sample the electron subshell :math:`i` using the number of electrons per + shell as the probability mass function. + +3. Sample :math:`p_z` using :math:`J_i(p_z)` as the PDF. + +4. Calculate :math:`E'` by solving :eq:`pz` for :math:`E'` using the sampled + value of :math:`p_z`. + +5. If :math:`p_z < p_{z,\text{max}}` for shell :math:`i`, accept :math:`E'`. + Otherwise repeat from step 2. Compton Electrons +++++++++++++++++ +Because the Compton-scattered photons can transfer a large fraction of their +energy to the kinetic energy of the recoil electron, which may in turn go on to +lose its energy as bremsstrahlung radiation, it is necessary to accurately +model the angular and energy distributions of Compton electrons. The energy of +the recoil electron ejected from the :math:`i`-th subshell is given by + +.. math:: + :label: compton-electron-energy + + E_{-} = E - E' - E_{b,i}. + +The direction of the electron is assumed to be in the direction of the momentum +transfer, with the cosine of the polar angle given by + +.. math:: + :label: compton-electron-mu + + \mu_{-} = \frac{E - E'\mu}{\sqrt{E^2 +E'^2 - 2EE'\mu}} + +and the azimuthal angle :math:`\phi_{-} = \phi + \pi`, where :math:`\phi` is +the azimuthal angle of the photon. The vacancy left by the ejected electron is +filled through atomic relaxation. Photoelectric Effect -------------------- +In the photoelectric effect, the incident photon is absorbed by an atomic +electron, which is then emitted from the :math:`i`-th shell with kinetic energy + +.. math:: + :label: photoelectron-kinetic-energy + + E_{-} = E - E_{b,i}. + +Photoelectric emission is only possible when the photon energy exceeds the +binding energy of the shell. These binding energies are often referred to as +edge energies because the otherwise continuously decreasing cross section has +discontinuities at these points, creating the characteristic sawtooth shape. +The photoelectric effect dominates at low energies and is more important for +heavier elements. + +When simulating the photoelectric effect, the first step is to sample the +electron shell. The shell :math:`i` where the ionization occurs can be +considered a discrete random variable with probability mass function + +.. math:: + :label: photoelectron-shell-pdf + + p_i = \frac{\sigma_{\text{pe},i}}{\sigma_{\text{pe}}}, + +where :math:`\sigma_{\text{pe},i}` is the cross section of the :math:`i`-th +shell, and the total photoelectric cross section of the atom, +:math:`\sigma_{\text{pe}}`, is the sum over the shell cross sections. Once the +shell has been sampled, the energy of the photoelectron is calculated using +:eq:`photoelectron-kinetic-energy`. + +To determine the direction of the photoelectron, we implement the method +described in Kaltiaisenaho_, which models the angular distribution of the +photoelectrons using the K-shell cross section derived by Sauter (K-shell +electrons are the most tightly bound, and they contribute the most to +:math:`\sigma_{\text{pe}}`). The non-relativistic Sauter distribution for +unpolarized photons can be approximated as + +.. math:: + :label: sauter + + \frac{d\sigma_{\text{pe}}}{d\mu_{-}} \propto + \frac{1 - \mu_{-}^2}{(1 - \beta_{-} \mu_{-})^4}, + +where :math:`\beta_{-}` is the ratio of the velocity of the electron to the +speed of light, + +.. math:: + :label: beta-2 + + \beta_{-} = \frac{\sqrt{(E_{-}(E_{-} + 2m_e c^2)}}{E_{-} + m_e c^2}. + +To sample :math:`\mu_{-}` from the Sauter distribution, we first express +:eq:`sauter` in the form: + +.. math:: + :label: photoelectron-mu-pdf + + f(\mu_{-}) = \frac{3}{2} \psi(\mu_{-}) g(\mu_{-}), + +where + +.. math:: + :label: mu-pdf-factors + + \psi(\mu_{-}) &= \frac{(1 - \beta_{-}^2)(1 - \mu_{-}^2)}{(1 - + \beta_{-}\mu_{-})^2}, \\ + g(\mu_{-}) &= \frac{1 - \beta_{-}^2}{2 (1 - \beta_{-}\mu_{-})^2}. + +In the interval :math:`[-1, 1]`, :math:`g(\mu_{-})` is a normalized PDF and +:math:`\psi(\mu_{-})` satisfies the condition :math:`0 < \psi(\mu_{-}) < 1`. +The following algorithm can now be used to sample :math:`\mu_{-}`: + +1. Using the inverse transform method, sample :math:`\mu_{-}` from + :math:`g(\mu_{-})` using the sampling formula + + .. math:: + + \mu_{-} = \frac{2\xi_1 + \beta_{-} - 1}{2\beta_{-}\xi_1 - \beta_{-} + 1}. + +2. If :math:`\xi_2 \le \psi(\mu_{-})`, accept :math:`\mu_{-}`. Otherwise, + repeat the sampling from step 1. + +The azimuthal angle is sampled uniformly on :math:`[0, 2\pi)`. + +The atom is left in an excited state with a vacancy in the :math:`i`-th shell +and decays to its ground state through a cascade of transitions that produce +fluorescent photons and Auger electrons. Pair Production --------------- +In electron-positron pair production, a photon is absorbed in the vicinity of +an atomic nucleus or an electron and an electron and positron are created. Pair +production is the dominant interaction with matter at high photon energies and +is more important for high-Z elements. When it takes place in the field of a +nucleus, energy is essentially conserved among the incident photon and the +resulting charged particles. Therefore, in order for pair production to occur, +the photon energy must be greater than the sum of the rest mass energies of the +electron and positron, i.e., :math:`E_{\text{threshold,pp}} = 2 m_e c^2 = +1.022` MeV. + +The photon can also interact in the field of an atomic electron. This process +is referred to as "triplet production" because the target electron is ejected +from the atom and three charged particles emerge from the interaction. In this +case, the recoiling electron also absorbs some energy, so the energy threshold +for triplet production is greater than that of pair production from atomic +nuclei, with :math:`E_{\text{threshold,tp}} = 4 m_e c^2 = 2.044` MeV. The ratio +of the triplet production cross section to the pair production cross section is +approximately 1/Z, so triplet production becomes increasingly unimportant for +high-Z elements. Though it can be significant in lighter elements, the momentum +of the recoil electron becomes negligible in the energy regime where pair +production dominates. For our purposes, it is a good approximation to treat +triplet production as pair production and only simulate the electron-positron +pair. + +Accurately modeling the creation of electron-positron pair is important because +the charged particles can go on to lose much of their energy as bremsstrahlung +radiation, and the subsequent annihilation of the positron with an electron +produces two additional photons. We sample the energy and direction of the +charged particles using a semiempirical model described in Salvat_. The +Bethe-Heitler differential cross section, given by + +.. math:: + :label: bethe-heitler + + \frac{d\sigma_{\text{pp}}}{d\epsilon} = \alpha r_e^2 Z^2 + \left[ (\epsilon^2 + (1-\epsilon)^2) (\Phi_1 - 4f_C) + + \frac{2}{3}\epsilon(1-\epsilon)(\Phi_2 - 4f_C) \right], + +is used as a starting point, where :math:`\alpha` is the fine structure +constant, :math:`f_C` is the Coulomb correction function, :math:`\Phi_1` and +:math:`\Phi_2` are screening functions, and :math:`\epsilon = (E_{-} + m_e +c^2)/E` is the electron reduced energy (i.e., the fraction of the photon energy +given to the electron). :math:`\epsilon` can take values between +:math:`\epsilon_{\text{min}} = k^{-1}` (when the kinetic energy of the electron +is zero) and :math:`\epsilon_{\text{max}} = 1 - k^{-1}` (when the kinetic +energy of the positron is zero). + +The Coulomb correction, given by + +.. math:: + :label: coulomb-correction + + f_C = \alpha^{2}Z^{2} \big[&(1 + \alpha^{2}Z^{2})^{-1} + 0.202059 + - 0.03693\alpha^{2}Z^{2} + 0.00835\alpha^{4}Z^{4} \\ + &- 0.00201\alpha^{6}Z^{6} + 0.00049\alpha^{8}Z^{8} + - 0.00012\alpha^{10}Z^{10} + 0.00003\alpha^{12}Z^{12}\big] + +is introduced to correct for the fact that the Bethe-Heitler differential cross +section was derived using the Born approximation, which treats the Coulomb +interaction as a small perturbation. + +The screening functions :math:`\Phi_1` and :math:`\Phi_2` account for the +screening of the Coulomb field of the atomic nucleus by outer electrons. Since +they are given by integrals which include the atomic form factor, they must be +computed numerically for a realistic form factor. However, by assuming +exponential screening and using a simplified form factor, analytical +approximations of the screening functions can be derived: + +.. math:: + :label: screening-functions + + \Phi_1 &= 2 - 2\ln(1 + b^2) - 4b\arctan(b^{-1}) + 4\ln(Rm_{e}c/\hbar) \\ + \Phi_2 &= \frac{4}{3} - 2\ln(1 + b^2) + 2b^2 \left[ 4 - 4b\arctan(b^{-1}) + - 3\ln(1 + b^{-2}) \right] + 4\ln(Rm_{e}c/\hbar) + +where + +.. math:: + :label: b + + b = \frac{Rm_{e}c}{2k\epsilon(1 - \epsilon)\hbar}. + +and :math:`R` is the screening radius. + +The differential cross section in :eq:`bethe-heitler` with the approximations +described above will not be accurate at low energies: the lower boundary of +:math:`\epsilon` will be shifted above :math:`\epsilon_{\text{min}}` and the +upper boundary of :math:`\epsilon` will be shifted below +:math:`\epsilon_{\text{max}}`. To offset this behavior, a correcting factor +:math:`F_0(k, Z)` is used: + +.. math:: + :label: correcting-factor + + F_0(k, Z) =~& (0.1774 + 12.10\alpha Z - 11.18\alpha^{2}Z^{2})(2/k)^{1/2} \\ + &+ (8.523 + 73.26\alpha Z - 44.41\alpha^{2}Z^{2})(2/k) \\ + &- (13.52 + 121.1\alpha Z - 96.41\alpha^{2}Z^{2})(2/k)^{3/2} \\ + &+ (8.946 + 62.05\alpha Z - 63.41\alpha^{2}Z^{2})(2/k)^{2}. + +To aid sampling, the differential cross section used to sample :math:`\epsilon` +(minus the normalization constant) can now be expressed in the form + +.. math:: + :label: pp-pdf + + \frac{d\sigma_{\text{pp}}}{d\epsilon} = + u_1 \frac{\phi_1(\epsilon)}{\phi_1(1/2)} \pi_1(\epsilon) + + u_2 \frac{\phi_2(\epsilon)}{\phi_2(1/2)} \pi_2(\epsilon) + +where + +.. math:: + :label: u + + u_1 &= \frac{2}{3} \left(\frac{1}{2} - \frac{1}{k}\right)^2 \phi_1(1/2), \\ + u_2 &= \phi_2(1/2), + +.. math:: + :label: phi + + \phi_1(\epsilon) &= \frac{1}{2}(3\Phi_1 - \Phi_2) - 4f_{C}(Z) + F_0(k, Z), \\ + \phi_2(\epsilon) &= \frac{1}{4}(3\Phi_1 + \Phi_2) - 4f_{C}(Z) + F_0(k, Z), + +and + +.. math:: + :label: pi + + \pi_1(\epsilon) &= \frac{3}{2} \left(\frac{1}{2} - \frac{1}{k}\right)^{-3} + \left(\frac{1}{2} - \epsilon\right)^2, \\ + \pi_2(\epsilon) &= \frac{1}{2} \left(\frac{1}{2} - \frac{1}{k}\right)^{-1}. + +The functions in :eq:`phi` are non-negative and maximum at :math:`\epsilon = +1/2`. In the interval :math:`(\epsilon_{\text{min}}, \epsilon_{\text{max}})`, +the functions in :eq:`pi` are normalized PDFs and +:math:`\phi_i(\epsilon)/\phi_i(1/2)` satisfies the condition :math:`0 < +\phi_i(\epsilon)/\phi_i(1/2) < 1`. The following algorithm can now be used to +sample the reduced electron energy :math:`\epsilon`: + +1. Sample :math:`i` according to the point probabilities + :math:`p(i=1) = u_1/(u_1 + u_2)` and :math:`p(i=2) = u_2/(u_1 + u_2)`. + +2. Using the inverse transform method, sample :math:`\epsilon` from + :math:`\pi_i(\epsilon)` using the sampling formula + + .. math:: + + \epsilon &= \frac{1}{2} + \left(\frac{1}{2} - \frac{1}{k}\right) + (2\xi_1 - 1)^{1/3} ~~~~&\text{if}~~ i = 1 \\ + \epsilon &= \frac{1}{k} + \left(\frac{1}{2} - + \frac{1}{k}\right) 2\xi_1 ~~~~&\text{if}~~ i = 2. + +3. If :math:`\xi_2 \le \phi_i(\epsilon)/\phi_i(1/2)`, accept + :math:`\epsilon`. Otherwise, repeat the sampling from step 1. + +Because charged particles have a much smaller range than the mean free path of +photons and because they immediately undergo multiple scattering events which +randomize their direction, it is sufficient to use a simplified model to sample +the direction of the electron and positron. The cosines of the polar angles are +sampled using the leading order term of the Sauter–Gluckstern–Hull +distribution, + +.. math:: + :label: sauter–gluckstern–hull + + p(\mu_{\pm}) = C(1 - \beta_{\pm}\mu_{\pm})^{-2}, + +where :math:`C` is a normalization constant and :math:`\beta_{\pm}` is the +ratio of the velocity of the charged particle to the speed of light given in +:eq:`beta-2`. + +The inverse transform method is used to sample :math:`\mu_{-}` and +:math:`\mu_{+}` from :eq:`sauter–gluckstern–hull`, using the sampling formula + +.. math:: + :label: sample-mu + + \mu_{\pm} = \frac{2\xi - 1 + \beta_{\pm}}{(2\xi - 1)\beta_{\pm} + 1}. + +The azimuthal angles for the electron and positron are sampled independently +and uniformly on :math:`[0, 2\pi)`. ------------------- Secondary Processes @@ -248,19 +599,297 @@ additional photons. Atomic Relaxation ----------------- +When an electron is ejected from an atom and a vacancy is left in an inner +shell, an electron from a higher energy level will fill the vacancy. This +results in either a radiative transition, in which a photon with a +characteristic energy (fluorescence photon) is emitted, or non-radiative +transition, in which an electron from a shell that is farther out (Auger +electron) is emitted. If a non-radiative transition occurs, the new vacancy is +filled in the same manner, and as the process repeats a shower of photons and +electrons can be produced. + +The energy of a fluorescence photon is the equal to the energy difference +between the transition states, i.e., + +.. math:: + :label: fluorescence-photon-energy + + E = E_{b,v} - E_{b,i}, + +where :math:`E_{b,v}` is the binding energy of the vacancy shell and +:math:`E_{b,i}` is the binding energy of the shell from which the electron +transitioned. The energy of an Auger electron is given by + +.. math:: + :label: auger-electron-energy + + E_{-} = E_{b,v} - E_{b,i} - E_{b,a}, + +where :math:`E_{b,a}` is the binding energy of the shell from which the Auger +electron is emitted. While Auger electrons are low-energy so their range and +bremsstrahlung yield is small, fluorescence photons can travel far before +depositing their energy, so the relaxation process should be modeled in detail. + +Transition energies and probabilities are needed for each subshell to simulate +atomic relaxation. Starting with the initial shell vacancy, the following +recursive algorithm is used to fill vacancies and create fluorescence photons +and Auger electrons: + +1. If there are no transitions for the vacancy shell, create a fluorescence + photon assuming it is from a captured free electron and terminate. + +2. Sample a transition using the transition probabilities for the vacancy + shell as the probability mass function. + +3. Create either a fluorescence photon or Auger electron, sampling the + direction of the particle isotropically. + +4. If a non-radiative transition occurred, repeat from step 1 for the vacancy + left by the emitted Auger electron. + +5. Repeat from step 1 for vacancy left by the transition electron. Electron-Positron Annihilation ------------------------------ +When a positron collides with an electron, both particles are annihilated and +generally two photons with equal energy are created. If the kinetic energy of +the positron is high enough, the two photons can have different energies, and +the higher-energy photon is emitted preferentially in the direction of flight +of the positron. It is also possible to produce a single photon if the +interaction occurs with a bound electron, and in some cases three (or, rarely, +even more) photons can be emitted. However, the annihilation cross section is +largest for low-energy positrons, and as the positron energy decreases, the +angular distribution of the emitted photons becomes isotropic. + +In OpenMC, we assume the most likely case in which a low-energy positron (which +has already lost most of its energy to bremsstrahlung radiation) interacts with +an electron which is free and at rest. Two photons with energy equal to the +electron rest mass energy :math:`m_e c^2 = 0.511` MeV are emitted isotropically +in opposite directions. Bremsstrahlung -------------- +When a charged particle is decelerated in the field of an atom, some of its +kinetic energy is converted into electromagnetic radiation known as +bremsstrahlung, or 'braking radiation'. In each event, an electron or positron +with kinetic energy :math:`T` generates a photon with an energy :math:`E` +between :math:`0` and :math:`T`. Bremsstrahlung is described by a cross section +that is differential in photon energy, in the direction of the emitted photon, +and in the final direction of the charged particle. However, in Monte Carlo +simulations it is typical to integrate over the angular variables to obtain a +single differential cross section with respect to photon energy, which is often +expressed in the form + +.. math:: + :label: bremsstrahlung-dcs + + \frac{d\sigma_{\text{br}}}{dE} = \frac{Z^2}{\beta^2} \frac{1}{E} + \chi(Z, T, \kappa), + +where :math:`\kappa = E/T` is the reduced photon energy and :math:`\chi(Z, T, +\kappa)` is the scaled bremsstrahlung cross section, which is experimentally +measured. + +Because electrons are attracted to atomic nuclei whereas positrons are +repulsed, the cross section for positrons is smaller, though it approaches that +of electrons in the high energy limit. To obtain the positron cross section, we +multiply :eq:`bremsstrahlung-dcs` by the :math:`\kappa`-independent factor used +in Salvat_, + +.. math:: + :label: positron-factor + + F_{\text{p}}(Z,T) = + & 1 - \text{exp}(-1.2359\times 10^{-1}t + 6.1274\times 10^{-2}t^2 - 3.1516\times 10^{-2}t^3 \\ + & + 7.7446\times 10^{-3}t^4 - 1.0595\times 10^{-3}t^5 + 7.0568\times 10^{-5}t^6 \\ + & - 1.8080\times 10^{-6}t^7), + +where + +.. math:: + :label: positron-factor-t + + t = \ln\left(1 + \frac{10^6}{Z^2}\frac{T}{\text{m}_\text{e}c^2} \right). + +:math:`F_{\text{p}}(Z,T)` is the ratio of the radiative stopping powers for +positrons and electrons. Stopping power describes the average energy loss per +unit path length of a charged particle as it passes through matter: + +.. math:: + :label: stopping-power + + -\frac{dT}{ds} = n \int E \frac{d\sigma}{dE} dE \equiv S(T), + +where :math:`n` is the number density of the material and :math:`d\sigma/dE` is +the cross section differential in energy loss. The total stopping power +:math:`S(T)` can be separated into two components: the radiative stopping +power :math:`S_{\text{rad}}(T)`, which refers to energy loss due to +bremsstrahlung, and the collision stopping power :math:`S_{\text{col}}(T)`, +which refers to the energy loss due to inelastic collisions with bound +electrons in the material that result in ionization and excitation. To obtain +the radiative stopping power for positrons, the radiative stopping power for +electrons is multiplied by :eq:`positron-factor`. Currently, the collision +stopping power for electrons is also used for positrons. + +While the models for photon interactions with matter described above can safely +assume interactions occur with free atoms, sampling the target atom based on +the macroscopic cross sections, molecular effects cannot necessarily be +disregarded for charged particle treatment. For compounds and mixtures, the +bremsstrahlung cross section is calculated using Bragg's additivity rule as + +.. math:: + :label: material-bremsstrahlung-dcs + + \frac{d\sigma_{\text{br}}}{dE} = \frac{1}{\beta^2 E} \sum_i \gamma_i Z^2_i + \chi(Z_i, T, \kappa), + +where the sum is over the constituent elements and :math:`\gamma_i` is the +atomic fraction of the :math:`i`-th element. Similarly, the radiative stopping +power is calculated using Bragg's additivity rule as + +.. math:: + :label: material-radiative-stopping-power + + S_{\text{rad}}(T) = \sum_i w_i S_{\text{rad},i}(T), + +where :math:`w_i` is the mass fraction of the :math:`i`-th element. The +collision stopping power, however, is a function of certain quantities such as +the mean excitation energy :math:`I` and the density effect correction +:math:`\delta_F` that depend on molecular properties. These quantities cannot +simply be summed over constituent elements in a compound, but should instead be +calculated for the material. Currently, we use Bragg's additivity rule to +calculate the collision stopping power as well, but this is not a good +approximation and should be fixed in the future. + .. _ttb: Thick-Target Bremsstrahlung Approximation +++++++++++++++++++++++++++++++++++++++++ +Since charged particles lose their energy on a much shorter distance scale than +neutral particles, not much error should be introduced by neglecting to +transport electrons. However, the bremsstrahlung emitted from high energy +electrons and positrons can travel far from the interaction site. Thus, even +without a full electron transport mode it is necessary to model bremsstrahlung. +We use a thick-target bremsstrahlung (TTB) approximation based on the models in +Salvat_ and Kaltiaisenaho_ for generating bremsstrahlung photons, which assumes +the charged particle loses all its energy in a single homogeneous material +region. + +To model bremsstrahlung using the TTB approximation, we need to know the number +of photons emitted by the charged particle and the energy distribution of the +photons. These quantities can be calculated using the continuous slowing down +approximation (CSDA). The CSDA assumes charged particles lose energy +continuously along their trajectory with a rate of energy loss equal to the +total stopping power, ignoring fluctuations in the energy loss. The +approximation is useful for expressing average quantities that describe how +charged particles slow down in matter. For example, the CSDA range approximates +the average path length a charged particle travels as it slows to rest: + +.. math:: + :label: csda-range + + R(T) = \int^T_0 \frac{dT'}{S(T')}. + +Actual path lengths will fluctuate around :math:`R(T)`. The average number of +photons emitted per unit path length is given by the inverse bremsstrahlung +mean free path: + +.. math:: + :label: inverse-bremsstrahlung-mfp + + \lambda_{\text{br}}^{-1}(T,E_{\text{cut}}) + = n\int_{E_{\text{cut}}}^T\frac{d\sigma_{\text{br}}}{dE}dE + = n\frac{Z^2}{\beta^2}\int_{\kappa_{\text{cut}}}^1\frac{1}{\kappa} + \chi(Z,T,\kappa)d\kappa. + +The lower limit of the integral in :eq:`inverse-bremsstrahlung-mfp` is non-zero +because the bremsstrahlung differential cross section diverges for small photon +energies but is finite for photon energies above some cutoff energy +:math:`E_{\text{cut}}`. The mean free path +:math:`\lambda_{\text{br}}^{-1}(T,E_{\text{cut}})` is used to calculate the +photon number yield, defined as the average number of photons emitted with +energy greater than :math:`E_{\text{cut}}` as the charged particle slows down +from energy :math:`T` to :math:`E_{\text{cut}}`. The photon number yield is +given by + +.. math:: + :label: photon-number-yield + + Y(T,E_{\text{cut}}) = \int^{R(T)}_{R(E_{\text{cut}})} + \lambda_{\text{br}}^{-1}(T',E_{\text{cut}})ds = \int_{E_{\text{cut}}}^T + \frac{\lambda_{\text{br}}^{-1}(T',E_{\text{cut}})}{S(T')}dT'. + +:math:`Y(T,E_{\text{cut}})` can be used to construct the energy spectrum of +bremsstrahlung photons: the number of photons created with energy between +:math:`E_1` and :math:`E_2` by a charged particle with initial kinetic energy +:math:`T` as it comes to rest is given by :math:`Y(T,E_1) - Y(T,E_2)`. + +To simulate the emission of bremsstrahlung photons, the total stopping power +and bremsstrahlung differential cross section for positrons and electrons must +be calculated for a given material using :eq:`material-bremsstrahlung-dcs` and +:eq:`material-radiative-stopping-power`. These quantities are used to build the +tabulated bremsstrahlung energy PDF and CDF for that material for each incident +energy :math:`T_k` on the energy grid. The following algorithm is then applied +to sample the photon energies: + +1. For an incident charged particle with energy :math:`T`, sample the number of + emitted photons as + + .. math:: + + N = \lfloor Y(T,E_{\text{cut}}) + \xi_1 \rfloor. + +2. Rather than interpolate the PDF between indices :math:`k` and :math:`k+1` + for which :math:`T_k < T < T_{k+1}`, which is computationally expensive, use + the composition method and sample from the PDF at either :math:`k` or + :math:`k+1`. Using linear interpolation on a logarithmic scale, the PDF can + be expressed as + + .. math:: + + p_{\text{br}}(T,E) = \pi_k p_{\text{br}}(T_k,E) + \pi_{k+1} + p_{\text{br}}(T_{k+1},E), + + where the interpolation weights are + + .. math:: + + \pi_k = \frac{\ln T_{k+1} - \ln T}{\ln T_{k+1} - \ln T_k},~~~ + \pi_{k+1} = \frac{\ln T - \ln T_k}{\ln T_{k+1} - \ln T_k}. + + Sample either the index :math:`i = k` or :math:`i = k+1` according to the + point probabilities :math:`\pi_{k}` and :math:`\pi_{k+1}`. + +3. Determine the maximum value of the CDF :math:`P_{\text{br,max}}`. + +3. Sample the photon energies using the inverse transform method with the + tabulated CDF :math:`P_{\text{br}}(T_i, E)` i.e., + + .. math:: + + E = E_j \left[ (1 + a_j) \frac{\xi_2 P_{\text{br,max}} - + P_{\text{br}}(T_i, E_j)} {E_j p_{\text{br}}(T_i, E_j)} + 1 + \right]^{\frac{1}{1 + a_j}} + + where the interpolation factor :math:`a_j` is given by + + .. math:: + + a_j = \frac{\ln p_{\text{br}}(T_i,E_{j+1}) - \ln p_{\text{br}}(T_i,E_j)} + {\ln E_{j+1} - \ln E_j} + + and :math:`P_{\text{br}}(T_i, E_j) \le \xi_2 P_{\text{br,max}} \le + P_{\text{br}}(T_i, E_{j+1})`. + +We ignore the range of the electron or positron, i.e., the bremsstrahlung +photons are produced in the same location that the charged particle was +created. The direction of the photons is assumed to be the same as the +direction of the incident charged particle, which is a reasonable approximation +at higher energies when the bremsstrahlung radiation is emitted at small +angles. .. _Koblinger: https://doi.org/10.13182/NSE75-A26663 @@ -273,3 +902,7 @@ Thick-Target Bremsstrahlung Approximation .. _LA-UR-04-0487: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-04-0487.pdf .. _LA-UR-04-0488: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-04-0488.pdf + +.. _Kaltiaisenaho: https://aaltodoc.aalto.fi/bitstream/handle/123456789/21004/master_Kaltiaisenaho_Toni_2016.pdf + +.. _Salvat: http://www.oecd-nea.org/globalsearch/download.php?doc=77434