source: trunk/documents/UserDoc/DocBookUsersGuides/PhysicsReferenceManual/latex/electromagnetic/standard/nuclearrec.tex @ 1345

\section[Ion Scattering]{Ion Scattering}

The necessity of accurately computing the characteristics of interatomic
scattering arises in many disciplines in which energetic ions pass
through materials. Traditionally, solutions to this problem not involving
hadronic interactions have been dominated by the multiple scattering,
which
is reasonably successful, but not very flexible. In particular, it
is relatively difficult to introduce into such a system a particular
screening function which has been measured for a specific atomic pair,
rather than the universal functions which are applied. 
 In many problems
of current interest, such as the behavior of semiconductor device
physics in a space environment, nuclear reactions, particle showers,
and other effects are critically important in modeling the full details 
of ion transport. The process $G4ScreenedNuclearRecoil$ provides 
simulation of ion elastic scattering \cite{ion_scat_model}.

\subsection{Method}

The method used in this computation is a variant of a subset of the
method described in Ref.\cite{MendenhallWellerXSection}.
A very short recap of the basic material is included here. The scattering
of two atoms from each other is assumed to be a completely classical
process, subject to an interatomic potential described by a potential
function\begin{equation}
V(r)=\frac{Z_{1}Z_{2}e^{2}}{r}\phi\left(\frac{r}{a}\right)\label{VR_eqn}\end{equation}
where $Z_{1}$ and $Z_{2}$ are the nuclear proton numbers, $e^{2}$
is the electromagnetic coupling constant ($q_{e}^{2}/4\pi\epsilon_{0}$
in SI units), $r$ is the inter-nuclear separation, $\phi$ is the
screening function describing the effect of electronic screening of
the bare nuclear charges, and $a$ is a characteristic length scale
for this screening. In most cases, $\phi$ is a universal function
used for all ion pairs, and the value of $a$ is an appropriately
adjusted length to give reasonably accurate scattering behavior. In
the method described here, there is no particular need for a universal
function $\phi$, since the method is capable of directly solving
the problem for most physically plausible screening functions. It
is still useful to define a typical screening length $a$ in the calculation
described below, to keep the equations in a form directly comparable
with our previous work even though, in the end, the actual value is
irrelevant as long as the final function $\phi(r)$ is correct. From
this potential $V(r)$ one can then compute the classical scattering
angle from the reduced center-of-mass energy $\varepsilon\equiv E_{c}a/Z_{1}Z_{2}e^{2}$
(where $E_{c}$ is the kinetic energy in the center-of-mass frame)
and reduced impact parameter $\beta\equiv b/a$\begin{equation}
\theta_{c}=\pi-2\beta\int_{x_{{\scriptscriptstyle 0}}}^{\infty}f(z)\, dz/z^{2}\label{theta_eqn}\end{equation}
where\begin{equation}
f(z)=\left(1-\frac{\phi(z)}{z\,\varepsilon}-\frac{\beta^{2}}{z^{2}}\right)^{-1/2}\label{fz_eqn}\end{equation}
and $x_{{\scriptscriptstyle 0}}$ is the reduced classical turning
radius for the given $\varepsilon$ and $\beta$. 

The problem, then, is reduced to the efficient computation of this
scattering integral. In our previous work, a great deal of analytical
effort was included to proceed from the scattering integral to a full
differential cross section calculation, but for application in a Monte-Carlo
code, the scattering integral $\theta_{c}(Z_{1},\, Z_{2},\, E_{c},\, b)$
and an estimated total cross section $\sigma_{{\scriptscriptstyle 0}}(Z_{1},\, Z_{2},\, E_{c})$
are all that is needed. Thus, we can skip algorithmically forward
in the original paper to equations 15-18 and the surrounding discussion
to compute the reduced distance of closest approach $x_{{\scriptscriptstyle 0}}$.
This computation follows that in the previous work exactly, and will
not be reintroduced here.

For the sake of ultimate accuracy in this algorithm, and due to the
relatively low computational cost of so doing, we compute the actual
scattering integral (as described in equations 19-21 of \cite{MendenhallWellerXSection})
using a Lobatto quadrature of order 6, instead of the 4th order method
previously described. This results in the integration accuracy exceeding
that of any available interatomic potentials in the range of energies
above those at which molecular structure effects dominate, and should
allow for future improvements in that area. The integral $\alpha$
then becomes (following the notation of the previous paper)
\begin{equation}
\alpha\approx\frac{1+\lambda_{{\scriptscriptstyle 0}}}{30}+\sum_{i=1}^{4}w'_{i}\, f\left(\frac{x_{{\scriptscriptstyle 0}}}{q_{i}}\right)\label{alpha_eq}\end{equation}
where 
\begin{equation}
\lambda_{{\scriptscriptstyle 0}}=\left(\frac{1}{2}+\frac{\beta^{2}}{2\, x_{{\scriptscriptstyle 0}}^{2}}-\frac{\phi'(x_{{\scriptscriptstyle 0}})}{2\,\varepsilon}\right)^{-1/2}\label{lambda_eqn}\end{equation}
\\
$w'_{i}\in${[}0.03472124, 0.1476903, 0.23485003, 0.1860249] \\
$q_{i}\in${[}0.9830235, 0.8465224, 0.5323531, 0.18347974]\\
Then \\
\begin{equation}
\theta_{c}=\pi-\frac{\pi\beta\alpha}{x_{{\scriptscriptstyle 0}}}\label{thetac_alpha_eq}\end{equation}


The other quantity required to implement a scattering process 
is the total scattering cross section $\sigma_{{\scriptscriptstyle 0}}$
for a given incident ion and a material through which the ion is propagating.
This value requires special consideration for a process such as screened
scattering. In the limiting case that the screening function is unity,
which corresponds to Rutherford scattering, the total cross section
is infinite. For various screening functions, the total cross section
may or may not be finite. However, one must ask what the intent of
defining a total cross section is, and determine from that how to
define it. 

In Geant4, the total cross section is used to determine a mean-free-path
$l_{\mu}$ which is used in turn to generate random transport distances
between discrete scattering events for a particle. In reality, where
an ion is propagating through, for example, a solid material, scattering
is not a discrete process but is continuous. However, it is a useful,
and highly accurate, simplification to reduce such scattering to a
series of discrete events, by defining some minimum energy transfer
of interest, and setting the mean free path to be the path over which
statistically one such minimal transfer has occurred. This approach
is identical to the approach developed for the original TRIM code
\cite{TRIM1980}. As long as the minimal interesting energy transfer
is set small enough that the cumulative effect of all transfers smaller
than that is negligible, the approximation is valid. As long as the
impact parameter selection is adjusted to be consistent with the selected
value of $l_{\mu}$, the physical result isn't particularly sensitive
to the value chosen.

Noting, then, that the actual physical result isn't very sensitive
to the selection of $l_{\mu},$ one can be relatively free about defining
the cross section $\sigma_{{\scriptscriptstyle 0}}$ from which $l_{\mu}$
is computed. The choice used for this implementation is fairly simple.
Define a physical cutoff energy $E_{min}$ which is the smallest energy
transfer to be included in the calculation. Then, for a given incident
particle with atomic number $Z_{1}$, mass $m_{1}$, and lab energy
$E_{inc}$, and a target atom with atomic number $Z_{2}$ and mass
$m_{2}$, compute the scattering angle $\theta_{c}$ which will transfer
this much energy to the target from the solution of\begin{equation}
E_{min}=E_{inc}\,\frac{4\, m_{1}\, m_{2}}{(m_{1}+m_{2})^{2}}\,\sin^{2}\frac{\theta_{c}}{2}\label{recoil_kinematics}\end{equation}.
Then, noting that $\alpha$ from eq. \ref{alpha_eq} is a number very close to unity, one can solve for an approximate 
impact parameter $b$ with a single root-finding operation to find the classical turning point.
Then, define the total cross section to be $\sigma_{{\scriptscriptstyle 0}}=\pi b^{2}$,
the area of the disk inside of which the passage of an ion will cause
at least the minimum interesting energy transfer. Because this process
is relatively expensive, and the result is needed extremely frequently,
the values of $\sigma_{{\scriptscriptstyle 0}}(E_{inc})$ are precomputed
for each pairing of incident ion and target atom, and the results
cached in a cubic-spline interpolation table. However, since the actual result isn't very critical, the
cached results can be stored in a very coarsely sampled table without
degrading the calculation at all, as long as the values of the $l_{\mu}$
used in the impact parameter selection are rigorously consistent with
this table.

The final necessary piece of the scattering integral calculation is
the statistical selection of the impact parameter $b$ to be used
in each scattering event. This selection is done following the original
algorithm from TRIM, where the cumulative probability distribution
for impact parameters is \begin{equation}
P(b)=1-\exp\left(\frac{-\pi\, b^{2}}{\sigma_{{\scriptscriptstyle 0}}}\right)\label{cum_prob}\end{equation}
where $N\,\sigma_{{\scriptscriptstyle 0}}\equiv1/l_{\mu}$ where $N$
is the total number density of scattering centers in the target material
and $l_{\mu}$ is the mean free path computed in the conventional
way. To produce this distribution from a uniform random variate $r$
on (0,1], the necessary function is\begin{equation}
b=\sqrt{\frac{-\log r}{\pi\, N\, l_{\mu}}}\label{sampling_func}\end{equation}
This choice of sampling function does have the one peculiarity that
it can produce values of the impact parameter which are larger than
the impact parameter which results in the cutoff energy transfer,
as discussed above in the section on the total cross section, with
probability $1/e$. When this occurs, the scattering event is not
processed further, since the energy transfer is below threshold. For
this reason, impact parameter selection is carried out very early
in the algorithm, so the effort spent on uninteresting events is minimized.

The above choice of impact sampling is modified when the mean-free-path
is very short. If $\sigma_{{\scriptscriptstyle 0}}>\pi\left(\frac{l}{2}\right)^{2}$
where $l$ is the approximate lattice constant of the material, as
defined by $l=N^{-1/3}$, the sampling is replaced by uniform sampling
on a disk of radius $l/2$, so that\begin{equation}
b=\frac{l}{2}\sqrt{r}\label{flat_sampling}\end{equation}
This takes into account that impact parameters larger than half the
lattice spacing do not occur, since then one is closer to the adjacent
atom. This also derives from TRIM. 

One extra feature is included in our model, to accelerate the production
of relatively rare events such as high-angle scattering. This feature
is a cross-section scaling algorithm, which allows the user access
to an unphysical control of the algorithm which arbitrarily scales
the cross-sections for a selected fraction of interactions. This is
implemented as a two-parameter adjustment to the central algorithm.
The first parameter is a selection frequency $f_{h}$ which sets what
fraction of the interactions will be modified. The second parameter
is the scaling factor for the cross-section. This is implemented by,
for a fraction $f_{h}$ of interactions, scaling the impact parameter
by $b'=b/\sqrt{scale}$. This feature, if used with care so that it
does not provide excess multiple-scattering, can provide between 10
and 100-fold improvements to event rates. If used without checking
the validity by comparing to un-adjusted scattering computations,
it can also provide utter nonsense. 


\subsection{Implementation Details}

\label{Lobatto}The coefficients for the summation to approximate
the integral for $\alpha$ in eq.(\ref{alpha_eq}) are derived from
the values in Abramowitz \& Stegun \cite{AbramowitzStegunLobatto},
altered to make the change-of-variable used for this integral. There
are two basic steps to the transformation. First, since the provided
abscissas $x_{i}$ and weights $w_{i}$ are for integration on {[}-1,1],
with only one half of the values provided, and in this work the integration
is being carried out on {[}0,1], the abscissas are transformed as:\begin{equation}
y_{i}\in\left\{ \frac{1\mp x_{i}}{2}\right\} \label{as_yxform}\end{equation}
 Then, the primary change-of-variable is applied resulting in:\begin{eqnarray}
q_{i} & = & \cos\frac{\pi\, y_{i}}{2}\label{q_xform}\\
w'_{i} & = & \frac{w_{i}}{2}\sin\frac{\pi\, y_{i}}{2}\label{wprime_xform}\end{eqnarray}
 except for the first coefficient $w'_{1}$where the $\sin()$ part
of the weight is taken into the limit of $\lambda_{{\scriptscriptstyle 0}}$
as described in eq.(\ref{lambda_eqn}). This value is just $w'_{1}=w_{1}/2$.

\subsection{Status of this document}
 06.12.07  created by V. Ivanchenko from paper of M.H.~Mendenhall and R.A.~Weller\\
06.12.07  further edited by M. Mendenhall to bring contents of paper up-to-date with current implementation.

\begin{latexonly}

\begin{thebibliography}{99}

\bibitem{ion_scat_model}
M.H.~Mendenhall, R.A.~Weller, An algorithm for computing screened Coulomb
scattering in Geant4, {\em Nucl. Instr. Meth. B 227 (2005) 420.}

\bibitem{MendenhallWellerXSection}
M.H.~Mendenhall, R.A.~Weller, Algorithms for the rapid computation of
  classical cross sections for screened coulomb collisions, 
{\em Nucl. Instr. Meth. in Physics Res. B58 (1991) 11.}

\bibitem{TRIM1980}
J.P.~Biersack, L.G.~Haggmark, A Monte Carlo computer program for the
  transport of energetic ions in amorphous targets, 
{\em Nucl. Instr. Meth. in Physics Res. 174 (1980) 257.}

\bibitem{AbramowitzStegunLobatto}
M.~Abramowitz, I.~Stegun (Eds.), Handbook of Mathematical Functions, 
{\em Dover, New York, 1965, pp. 888, 920.}

\end{thebibliography}

\end{latexonly}

\begin{htmlonly}

\subsection{Bibliography}

\begin{enumerate}

\item
M.H.~Mendenhall, R.A.~Weller, An algorithm for computing screened Coulomb
scattering in Geant4, {\em Nucl. Instr. Meth. B 227 (2005) 420.}

\item
M.H.~Mendenhall, R.A.~Weller, Algorithms for the rapid computation of
  classical cross sections for screened coulomb collisions, 
{\em Nucl. Instr. Meth. in Physics Res. B58 (1991) 11.}

\item
J.P.~Biersack, L.G.~Haggmark, A Monte Carlo computer program for the
  transport of energetic ions in amorphous targets, 
{\em Nucl. Instr. Meth. in Physics Res. 174 (1980) 257.}

\item
M.~Abramowitz, I.~Stegun (Eds.), Handbook of Mathematical Functions, 
{\em Dover, New York, 1965, pp. 888, 920.}

\end{enumerate}

\end{htmlonly}