crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge pull request #1085 from amandalund/photon-documentation

Browse source 
Photon physics documentation
This commit is contained in:
Paul Romano 2018-10-06 10:09:34 -05:00 • committed by GitHub
commit 6dce5e3157
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
1 changed files with 663 additions and 30 deletions
Download patch file Download diff file Expand all files Collapse all files

693
docs/source/methods/photon_physics.rst
View file

@ -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