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 [MW05]. This process is available with extended electromagnetic example TestEm7.
Method
The method used in this computation is a variant of a subset of the method described in Ref.[MW91]. 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
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\)
where
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 [MW91]) 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)
where
Then
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 [BH80]. 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
Then, noting that \(\alpha\) from Eq.(61) 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
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
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
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.
Implementation Details
The coefficients for the summation to approximate the integral for \(\alpha\) in Eq.(61) are derived from the values in Abramowitz & Stegun [MA65], 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:
Then, the primary change-of-variable is applied resulting in:
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.(62). This value is just \(w'_{1}=w_{1}/2\).
Bibliography
- BH80
J.P. Biersack and L.G. Haggmark. A monte carlo computer program for the transport of energetic ions in amorphous targets. Nucl. Instr. and Meth. in Phys. Research, 174():257, 1980.
- MA65
I. Stegun (Eds.) M. Abramowitz. Handbook of Mathematical Functions. Dover, New York, edition, 1965.
- MW91(1,2)
M.H. Mendenhall and R.A. Weller. Algorithms for the rapid computation of classical cross sections for screened coulomb collisions. Nucl. Instr. and Meth. in Phys. Research B, 58():11, 1991.
- MW05
M.H. Mendenhall and R.A. Weller. An algorithm for computing screened coulomb scattering in geant4. Nuclear Instruments and Methods in Physics Research B, 227():420, 2005.