Enabling attractive-repulsive potentials in binary-collision-approximation monte-carlo codes for ion-surface interactions

Binary Collision Approximation (BCA) codes for ion-material interactions, such as SRIM, Tridyn, F-TRIDYN, and SDtrimSP, have historically been limited to screened Coulomb potentials even at low energies due to the difficulty in numerically solving the Distance of Closest Approach (DOCA) problem for attractive-repulsive potentials. Techniques such as direct n-body simulation or modifications to Newton’s method are either prohibitively costly or not guaranteed to work for all potentials. Advanced rootfinding techniques, such as companion matrix solvers, offer a solution. For many attractive-repulsive potentials, however, a companion matrix cannot be used directly, because there is no way to put the associated functions into a monomial basis form. A complementary technique is proxy rootfinding—by finding the best-fit polynomial approximant of a function, the zeros of the approximant can be guaranteed to be close to the zeros of the function. Using the Chebyshev basis and grid offers additional guarantees with regards to the quality of the approximation, the speed of convergence, and the avoidance of Runge’s phenomenon. By finding Chebyshev interpolants and using the Chebyshev-Frobenius companion matrix, the zeros of any real function on a bounded domain can be found. Here we show that using an Adaptive Chebyshev Proxy Rootfinder with Automatic Subdivision (ACPRAS) with appropriate scaling functions, numerical issues presented by attractive-repulsive potentials, including those of scale, can be handled. Using these techniques, we show that it is possible to include any physically reasonable interatomic potential in a BCA code, and to guarantee correctness of the resulting scattering angle calculations.


Introduction
The Binary Collision Approximation (BCA) is a popular, typically Monte Carlo, technique used for the numerical simulation of ion-matter interactions.Codes such as SRIM [1], TRIM [2], Tridyn [3,4], F-TRIDYN [5,6], MARLOWE [7], and SDtrimSP [8], are widely known implementations adopted for a variety of applications, from radiotherapy to nuclear fusion, and more generally for any problem involving the interaction of energetic ions with matter that consist of a series of atomic collisions.In most BCA codes, each binary collision is assumed to be perfectly elastic, with inelastic mechanisms calculated separately [9,10].The theory of elastic binary collisions between point particles has been well established since Rutherford [11].The key quantity of interest for scattering via a spherically symmetric potential is the scattering angle θ (see figure 1), which is responsible for determining nuclear stopping, the mechanics of reflection of energetic ions from surfaces, and the details of backward momentum transfer to surface atoms responsible for sputtering, processes that dominate plasma-material interactions in fusion.Scattering angles are found through the scattering integral, ( ) Equation (1) shows the scattering integral of two particles, where V(r) is a spherically symmetric interaction potential, p is the impact parameter, the perpendicular distance between the initial trajectory of the incident particle and the target particle, r, the variable of integration, is the distance between the two colliding particles in the center of mass frame, E r = E a M b /(M a + M b ) is the center of mass energy of the two particles a and b, and, importantly, r 0 , the lower bound of the scattering integral, is the distance of closest approach (DOCA).
Finding the value of the DOCA is a non-trivial pre-requisite to solving the classical scattering integral.If the value of the distance of closest approach is incorrect, the value of the scattering angle and thus nuclear stopping, reflection coefficients, implantation distributions, and sputtering yields will also be incorrect.If the value is sufficiently incorrect, such that in the computation of the scattering angle p 2 /r 2 + V(r)/E r > 1, the denominator of equation (1) can become complex, leading to unphysical scattering angles.For certain special cases, such as Coulombic scattering, it can be found analytically, but for most practical potentials it must be found numerically.Formally, the distance of closest approach is the distance r that is the outermost root of the distance of closest approach function g(r), For the special case of p = 0, that is, a head-on collision, it is apparent that the distance of closest approach is the distance at which the potential and kinetic energies are equal, as seen in equation (3).
The second term, p 2 /r 2 , in the DOCA function represents the rotational energy of the center of mass system.To solve for the distance of closest approach, Newton's method is most often used, as it is notably in the TRIM and Tridyn family of codes [2,12].Newton's method requires an initial guess r n=0 at iteration n = 0, and after that the numerical solution is found by repeated iteration of the standard Newton procedure, . For any screened Coulomb potential of the type V(r) ∼ e kr /r, Newton's method is sufficient.Screened Coulomb potentials decrease monotonically from + ∞ at r = 0 to zero as r approaches + ∞ .By examination of equation (2), this behavior guarantees that there is only one root to the DOCA function, and thus Newton's method is guaranteed to work with any non-singular initial guess.
For interaction potentials that are not as well-behaved as screened Coulomb potentials, however, Newton's method is not guaranteed to work [13].Specifically, attractive-repulsive potentials, which have a hard repulsive core as r approaches zero and an attractive regime starting a finite distance from the center of the potential, can have up to three zeros.For the classical scattering problem, the correct DOCA is specifically the outermost root of the DOCA function.This is because, for classical scattering, approaching atoms are assumed to start at infinity and approach the collision partner; the outermost root of the DOCA function will be the first distance at which the particles turn away from each other; for this reason, the DOCA is also known as the classical turning radius.Multiple zeros or multiplicity of zeros can lead Newton's method to converge to the wrong zero or fail to converge at all.For these potentials, which includes common attractive-repulsive potentials such as Lennard-Jones or Morse potentials, a different root-finding method that is robust to these problems must be used.Some examples of relevant interatomic potentials, as well as if they introduce multiple zeros into the distance of closest approach problem, are shown in table 1.
Interatomic potentials from table 1, including Coulomb, screened Coulomb, and attractive-repulsive potentials, are shown in figure 2. Some modern potentials are shown in figure 3, including potentials by Björkas et al [14], Wang et al [15], Marinica et al [16], and Derlet et al [17].Modern interaction potentials are usually derived from near-first-principle approaches, such as Density Functional Theory (DFT).
A previous attempt to include attractive-repulsive potentials in a BCA code was performed by Robinson in the MARLOWE code [18] for the special case of a Morse potential.However, the modified Newton's method outlined therein is not described in enough detail to recreate it, nor does it, making assumptions about the original implementation, appear to be guaranteed to work for any interaction potential.Other approaches have used direct n-body solvers to verify the results of the Frobenius companion matrix method for polynomial rootfinding [19] in the case of power-law potentials such as Lennard-Jones, but widespread application of attractive-repulsive potentials in BCA codes specifically is missing from the scientific literature.
In this work we present a general technique able to handle arbitrary attractive-repulsive potentials in BCA codes.The technique is based on an advanced strategy of rootfinding, named Adaptive Chebyshev Proxy Rootfinder with Automatic Subdivision (ACPRAS) [13].In the following sections, first we discuss scaling functions S(r) which can be applied to the DOCA function, giving practical examples of simple functions and Table 1.Some relevant interatomic potentials, including whether they introduce multiple zeros to the DOCA problem.For Lenz-Jensen, = y r 9.67 and the constants (a, b, c) = (0.3344, 0.485, 0.002 647).These potentials are included in RustBCA [20].
Lenz-Jensen V(r) = (1 techniques to remove singularities (section 2.1).Then we demonstrate the appearance of multiple roots in attractive-repulsive potentials, for the classical case of a Lennard-Jones potential (section 2.2).In section 3 we describe the ACPRAS technique applied to the DOCA problem in detail, and point the reader to a software implementation maintained via GitHub.In section 4 we report multiple numerical examples of implementation, showing how the ACPRAS method works for (1) simple Lennard-Jones type potentials, (2) modern interaction potentials, and (3) an extreme case of a transcendental function having a large number of zeros.Finally, we implemented ACPRAS in the open-source ion-material interactions code RustBCA [20] and tested it for several cases of practical interest, namely: (1) an ion beam of light isotopes (hydrogen and helium) impacting on a layered TiO2-Al-Si target, and (2) a hydrogen beam reflecting on a tungsten surface.We compared typical outputs of a BCA analysis, such as implantation profiles, damage production (Frenkel Pairs), and reflection coefficients, for universal screened Coulomb potentials versus modern attractive-repulsive potentials.
2. Singularities and root multiplicity in the DOCA problem 2.1.Singularities in the DOCA problem It can be seen by inspection of equation (2) that the DOCA function has multiple singularities, most notably as r approaches zero.In order to handle this function numerically, the singularities must be tamed.Strategies to avoid singularities exist.Transformations of the equation, such as squaring, will not modify its roots, since the domain is (0, ∞ ).Reversible changes of variables, such as r = 1/y 2 , can be applied while the roots of the original function remain recoverable.Applying these two strategies to the DOCA function can produce a singularity free equation, equation (4), ready for numerical analysis, Another strategy to improve the numerical resilience of a rootfinding technique is multiplying the function by a scaling function S(r).This can be done so long as S(r) has no roots on the domain of interest, in this case, (0, ∞ ).Some examples of valid scaling functions are ( ) ( ) = + S r r a n or S(r) = e − ar .Once a suitable set of transformations has been performed, the resulting form of the DOCA function can be used with rootfinding methods to find the distance of closest approach.

Multiple roots of attractive-repulsive potentials
To demonstrate the number of roots of an attractive-repulsive potential's DOCA function, here we consider a Lennard-Jones type potential.These results can, in principle, be extended to more general attractive-repulsive potentials.Define a Lennard-Jones type potential as an interatomic pair potential of the form, where r > 0, σ > 0, m > 1, n > 1, and where m > n and, typically, n divides m; that is, m = np.Potentials of this form are commonly used when a semi-realistic form of the interatomic potential is needed, but computational speed is more important than the exact details of the interaction.Limitations on the use and appropriateness of Lennard-Jones type potentials specifically is considered in [21].The DOCA function is the function whose largest root is the classical turning point of the binary collision problem.Other roots may exist, but these are in classically forbidden regions-areas where the particle must pass through an energy barrier larger than the incident energy to arrive there-and are not physically meaningful for classical scattering problems.This function, found here by squaring g(r), is This function is undefined at r = 0. To handle this function numerically, the singularity at zero must be tamed.For a Lennard-Jones type potential, this function becomes The singularity at zero can, assuming m 2, be handled by multiplying F(r) by r m , which only introduces zeros at r = 0.
If there is just one root to G(r), a simple method such as bisection could, in principle, be used.When there are multiple roots to G(r), this degree-m polynomial can be solved by using the Frobenius companion matrix method.In order to determine if this is necessary, it would be beneficial to know when G(r) has multiple real roots based on the input parameters, p and E r .Starting from G(r), we can solve the distance of closest approach function for the impact parameter p Examining p(r), it is apparent that r(p) has multiple zeros whenever p(r) has extrema.We can determine when those extrema exist by finding roots of the first derivative of p(r) This function has roots when the numerator is zero.Thus, we can reduce the problem of deciding whether there are multiple roots to a degree-m polynomial with nonzero coefficients of degree m, m − n, m − 2, and 0 to a problem in deciding whether the following polynomial with nonzero coefficients of degree m, m − n, and 0 has real roots, by making the transformation = z r 1 ( ) This is a 3-term degree-m polynomials with the coefficients A and B: In the form: For a typical Lennard-Jones potential, since n divides m, i.e. m = np, this can be further simplified with a substitution of u = z n : We have 3 cases where the number of real roots are immediately and easily decidable: p = 4, 3, and 2. Using the discriminants for degree-p polynomials where the highest-degree coefficient is 1, the coefficient of second highest degree is A, the coefficient of degree zero is B, and all other coefficients are zero: Taking the example of the Lennard-Jones 12-6 potential, where n = 6 and p = 2, we can find a simple condition on the center of mass energy for the existence of multiple roots of the distance of closest approach function: When the discriminant is less than zero, there is only one real root; this is our condition for the existence of multiple roots, and it can be found in terms of epsilon, When this simple condition is not met, there is exactly one real root to the distance of closest approach function, and a simple method such as bisection can be used to determine the value of the root.There are other cases that result in similarly simple expressions; for example, m = 9 n = 3, where n and m share the common factor of 3, yields the expression Even for energies above this value, however, Newton's method specifically may still fail.In figure 4, Newton's method and an advanced rootfinder are compared.Newton's method, even for energies above the energy threshold, here 0.1 and 1.0 eV, converges to the wrong root, especially notable at small impact parameters.Even when attractive-repulsive potentials can be guaranteed to have a single real root, the multiplicity of that root or the existence of complex roots can still lead to incorrect results for some rootfinders.Plotting the complex roots using a method such as the Delves and Lyness method [22][23][24], one can see some failure points more clearly.In particular, the Python 'pydelves' implementation [25] of the Delves and Lyness method was used for the complex rootfinding in this work.In figure 5, it is used to find all zeros of the DOCA function in the complex plane.Their real components are plotted as circles.Compared to the previous figure, This figure shows the failure of Newton's method by comparing solutions to the DOCA problem using Newton's method (lines) and the Chebyshev Proxy Rootfinder (filled circles).Here, R is the DOCA and p is the impact parameter.It is apparent that Newton's method often converges on the wrong real root or converges to a complex root with a small imaginary part.In contrast, the CPR rootfinder accurately determines all roots simultaneously, and, in the context of the DOCA problem, the outermost root can be directly computed to high accuracy.
figure 4, it is apparent that many places where Newton's method fails are those where multiple roots exist or where a complex root with a small imaginary component exists (small impact parameters and high energy).

Chebyshev interpolants as rootfinding proxies
One advanced rootfinding technique relies on the concept of proxy functions.In essence, a proxy function is built from some set of basis functions that, by some technique, approximates the true function for the purposes of rootfinding.Roots of the proxy function, which are simpler to find than roots of the true function by design, are thus proxies for the roots of the true function.Chebyshev polynomial proxies have been shown to be the optimal polynomial basis for proxy rootfinding, but other basis functions may be used [13].The theory and technique of Chebyshev Proxy Rootfinding for real roots is covered in detail by Boyd [13,[26][27][28][29]. Chebfun [30], an open-source MATLAB package that uses Chebyshev proxies for the manipulation of functions including rootfinding, is a popular example of Chebyshev proxies in action.
Chebyshev proxies for bounded functions can be created in a straightforward way, using the following ingredients.First, given the interval [a, b], and the function to be approximated, f (x), the Chebyshev grid points (equation (22)) and the values of ˜( ) f x , the Chebyshev interpolant of f (x) can be calculated.By using the Chebyshev grid, x k , Runge's pheonomenon can be avoided [31,32].The Chebyshev gridpoints can be calcualted in a straightforward way, Coefficients to the Nth order Chebyshev interpolant of f (x) on the Chebyshev grid can be found via the interpolation matrix, reported in [29], reproduced below in equation (23), and the original function values on the Chebyshev grid, as shown in equation (24).p i takes the value 2 if i = 0 or if i = N; it takes the value 1 otherwise.
Once the Chebyshev coefficients have been calculated, a convergence criterion can be constructed when performing adaptive interpolation to compare subsequent iterations.An error polynomial, δ 2N (x), can be calculated from the difference of two successive Chebyshev interpolants, ˜( ) where T i (x) is the ith .This plot shows the real parts of the complex roots of the DOCA function for the Lennard-Jones 12-6 function using Delves' method for complex rootfinding [22].Here, R is the distance of closest approach, a is the Bohr radius, p is the impact parameter, and epsilon is the reduced energy.
order Chebyshev basis polynomial, An associated error term can be found by summing the magnitudes of the coefficients of the error polynomial: Once the coefficients to a Chebyshev interpolant have been found, the Clenshaw-Curtis recurrence relation can be used to produce the value of the interpolant at any point x in [a, b].This recurrence relation is initialized with the following, Once the Clenshaw-Curits recurrence relation has been initialized, the value of the Chebyshev interpolant at a given x can be found with the iteration in table 2, repeated for all i in 1, 2,...,N + 1, where N is the order of the Chebyshev interpolant.Note that the following reproduction corrects a typo present in the original reference [13].Finally, the value of the interpolant at x is found from the final values of b 0 , b 3 , and the zeroth-order Chebyshev coefficient of the interpolant, To summarize, the coefficients of an Nth order Chebyshev interpolant ˜( ) f x N and its value at a given x can be calculated using the following procedure: 6.The value of the Chebyshev interpolant at the given x can then be found using b 0 , b 3 , and a 0 using equation (28) Following this procedure allows one to construct Chebyhsev proxies for the purpose of rootfinding.To find the roots of a given interpolant, a polynomial rootfinding technique must be used.Since polynomial rootfinding techniques are more robust than general transcendental rootfinding techniques, the use of proxies allows one to use a robust rootfinding method for solving transcendental equations.

Polynomial rootfinding with the frobenius companion matrix
One method of polynomial rootfinding, popularized by the computational linear algebra community, is known as the companion matrix method [33].This method solves for polynomial roots by finding the eigenvalues of a particular matrix constructed from the polynomial coefficients, termed the companion matrix.For every polynomial p(x), expressed in the monomial basis, Table 2. Algorithmic implementation of the Clenshaw-Curits relation.
the Frobenius companion matrix is defined as, This matrix is a sparse, nonsymmetric matrix whose eigenvalues are the roots of the polynomial p(x).This method of root-finding will find all roots, real and complex, simultaneously, at only the cost of finding the eigenvalues of an n × n matrix.For some potentials, such as Lennard-Jones, that can be reformulated into a monomial basis, this is the fastest root-finding method available guaranteed to find the correct root of the distance of closest approach function.

Chebyshev-frobenius companion matrix
A companion matrix method can also be used to directly find the roots of Chebyshev proxies.Since Chebyshev polynomials are not naturally expressed in the monomial basis, however, the form of the Chebyshev-Frobenius companion matrix is different from the monomial basis companion matrix, as seen in equation (31), To construct the Chebyshev-Frobenius matrix, one needs only the Chebyshev coefficients a i of the interpolant.The eigenvalues of this matrix will be the roots of the Chebyshev interpolant, and, in the case of proxies, the approximate roots of the original function.

Adaptive chebyshev proxy rootfinder with automatic subdivision
To practically use Chebyshev proxies to find the roots of real functions requires some extra machinery to handle issues related to the high dynamic range of attractive-repulsive potentials in high-energy classical scattering problems and their associated DOCA functions.These are adaptive interpolation, wherein the order of the Chebyhsev polynomial used for interpolation is doubled until it reaches a convergence criterion involving the magnitude of the error term described in equation (26), and automatic subdivision.Given a constraint on the maximum order of Chebyshev interpolant, an Adaptive Chebyshev Proxy Rootfinder with Automatic Subdivision (ACPRAS) can automatically determine a troublesome interval and subdivide it into two or more subintervals, each of which will more likely be easily fit by a single Chebyshev interpolant.It does so when the maximum order has been exceeded for a particular interval without meeting the convergence criterion.If an interval is split into two or more subintervals, the Chebyshev proxy procedure continues on each subinterval, and subsequently on any subsubintervals, and so on, until the convergence criterion is met.
A complete implementation of the ACPRAS method would proceed as follows, starting with the function f (x) whose roots are to be found, the interval on which to look for roots [a, b], an initial Chebyshev interpolant order N, a convergence criterion ò < error, and a maximum Chebyshev order N max : (i) For the interval [a, b], construct the Chebyshev gridpoints given a Chebyshev order N (equation ( 22)).
(ii) From the function evaluations on the grid, f (x k ), and the interpolation matrix I jk , calculate the Chebyshev coefficients a j (equations ( 23) and (24)).
(iii) Calculate the error term (equation ( 26)) between two successive interpolants (the initial interpolant being ˜( ) = f x 0 0 ) and check the convergence criterion.If it is met, proceed to step 4, and if not, return to step 1 with a doubled order 2N; for automatic subdivision, if 3.5.Implementation of ACPRAS for attractive-repulsive potentials An implementation of the ACPRAS method in the Rust programming language is available, titled rcpr [34].This crate, the term for Rust packages, exposes methods for Chebyshev Proxy rootfinding and monomial-form rootfinding.The rcpr crate has then been added to the ion-material interactions code RustBCA [20].For the first time, this enables a BCA code to include arbitrary interaction potentials, including cubic-spline piecewise potentials, Lennard-Jones potentials, and more, allowing a BCA code to remain accurate at lower energies or be compared directly to molecular dynamics codes; this latter possibility being mentioned in [35].Some features of attractive-repulsive potentials in the BCA are worthy of discussion.Particularly, attractiverepulsive potentials can exhibit an orbiting condition, where, for a given incident energy, a binary collision at a particular impact parameter can induce the particles to orbit each other.This condition, which involves a singularity (see figure 6), may lead to inaccurate scattering angles for energies and impact parameters near the orbiting condition [36].In order to overcome this difficulty, RustBCA includes multiple Gaussian quadrature schemes, including the arbitrary order Gauss-Mehler scheme for accurately resolving scattering angles near the orbiting condition.
Additionally, the distances of closest approach show the effect of choosing the outermost of the multiple turning radii (i.e.multiple roots of the DOCA function), as shown in figure 7. Essentially, this results in discontinuities in the distances of closest approach that other available BCA codes cannot reproduce.RustBCA, although slower when using ACRPAS instead of Newton's method and a high-order Gauss-Mehler quadrature instead of an approximation such as the MAGIC algorithm [1] implemented in the TRIM family of codes, nevertheless produces correct results.

Numerical results
In this section we report several examples of the ACPRAS technique.The methodology has been tested on attractive-repulsive potentials such as Lennard-Jones, and on modern Embedded Atom Model potentials obtained from DFT codes.We have also tested the ACPRAS method on extreme cases of transcendental functions with oscillatory behavior and a large number of zeros.The Python script used for these examples is open-source and available online [37].Finally, we implemented ACPRAS in the open-source Binary Collision Approximation code RustBCA [20], and analyzed the effect of an attractive-repulsive potential on two problems, (1) an ion beam of light isotopes (hydrogen and helium) impacting on a layered TiO2-Al-Si target, and (2) a hydrogen beam reflecting on a tungsten surface.

DOCA function with L-J potential
An example of ACPRAS applied to a DOCA function for a Lennard-Jones 12-6 potential is shown in figure 8.For this problem, the Lennard-Jones parameters were σ = 1 Angstrom and ò = 1 eV.Collision parameters for this problem were p=0.1 Angstrom and E r = 10 eV.A scaling function of the form S(r) = 1/(r + 1) 6 was used.Here, the distance of closest approach is correctly identified by automatically subdividing the interval around the zero point.For each specific problem, the maximum Chebyshev interpolant order and precision can be tuned to optimize the speed at which the algorithm reaches a solution.When the maximum order is low, the algorithm will prioritize making new subdivisions associated with small eigenvalue problems over solving fewer, larger eigenvalue problems.Automatic subdivision is a powerful tool, especially for problems where the dynamic range is high or difficult to predict in real time, such as in a BCA code.

Embedded Atom Model tungsten-tungsten potential
Applied to a much more difficult, piecewise-defined, cubic spline Embedded Atom Model tungsten-tungsten potential [17], ACPRAS performs very well, successfully finding all roots of the problem as seen in figure 9. Here, the collision parameters were p = 0.1 Angstrom and E r = 0.1 eV.This problem highlights the novel ability for ACPRAS to solve distance of closest approach problems for arbitrarily complex interaction potentials.By guaranteeing an accurate solution to the DOCA problem regardless of the complexity of the interaction, ACPRAS enables the use of advanced interaction potentials in the BCA, potentially significantly lowering the minimum energy of validity of the model and enabling direct comparisons to molecular dynamics codes [35].

Transcendental test function with arbitrary number of zeros
As an extreme case, we tested ACPRAS on a highly nonlinear transcendental function having many zeros, some very near each other, The function of equation (32), shown in figure 10, has both a periodic component and a 6th-degree polynomial component, both contributing to the large number of zero crossings.ACPRAS automatically split the domain into 20 intervals and correctly identifies all zeros (marked with stars in figure 10).Even though a physically realistic attractive-repulsive interatomic potential will never be as complex as equation (32), this test shows the robustness of the ACPRAS method for even an extreme rootfinding problem.

H and He beams implanting on layered TiO 2 -Al-Si target
This example is demonstrated in figure 11.Here, a 2 keV beam of hydrogen (red) and helium (fuchsia) is incident on a layered TiO2-Al-Si target.Hydrogen interacts with all other species through an arbitrarily defined  Lennard-Jones 6.5-6 potential, while all other interactions use the universal Kr-C screened Coulomb potential.The distinct nature of the attractive-repulsive trajectories is apparent in their wide lateral spread and more chaotic paths.Screened Coulomb potentials typically result in small angle scattering, but the same is not true of attractive-repulsive potentials in general.Many modern potentials are piecewise cubic spline potentials, fit using high-fidelity quantum simulations or to experiment.One such potential is reported in [14], modified from an earlier embedded atom potential to include correct high-energy interactions.This potential can be compared to the classical Kr-C potential for depth distributions of 2keV implantation of hydrogen on tungsten for the first time; such distributions are shown in figure 12.It is apparent that using an attractive-repulsive potential significantly changes the nature of  the elastic scattering for hydrogen on tungsten.Use of advanced, pair-specific potentials may yield BCA predictions significantly more accurate than those based entirely on universal screened Coulomb potentials.
Additionally, comparisons between traditionally used potentials such as Lenz-Jensen and Morse can be made.Lenz-Jensen is a potential particularly well suited for ion surface scattering.Morse is a potential that produces some realistic aspects of true interatomic potentials but is relatively simple to compute.In figure 13, a direct comparison between the potentials can be found in the implantation distributions and Frenkel-pair production distributions of hydrogen on a layered titanium dioxide, aluminum, and silicon surface.

Hydrogen reflection on tungsten
Finally, the Kr-C potential is compared to the realistic W-W/W-H potentials from Björkas and Wang respectively for reflected energy distributions.Hydrogen recycling is one of the key mechanisms at the plasmamaterial interface that can influence upstream plasma physics in fusion devices such as tokamaks.Prompt reflected velocity distributions of mostly entirely neutralized hydrogen can be calculated by the BCA, although, being difficult to measure directly, the accuracy of the distributions produced by the classical screened Coulomb potentials is not assured.In figure 14, a direct comparison of the reflected energy spectra for the Kr-C potential and W-W/W-H potentials is performed using RustBCA.It is apparent that the potential choice plays a large role in the shape of the reflected energy distribution-the peak energy in the realistic case is closer to the incident energy, and the amount of reflected hydrogen with less energy is drastically reduced.Reflection coefficients for the two simulations were 0.42 and 0.65 for Kr-C and W-W/W-H respectively-a 50% increase in the case of realistic potentials.Interestingly, the sputtering yields for the two simulations were 0.0020 and 0.020 respectively -a tenfold increase when the realistic potentials were used compared to Kr-C.

Conclusions
In this work we presented a numerical technique enabling realistic attractive-repulsive potentials in Binary Collision Approximation codes.The technique is based on an Adaptive Chebyshev Proxy Rootfinder with   Automatic Subdivision (ACPRAS), coupled to a proper re-scaling of the potential function.In practice, an arbitrary attractive-repulsive potential's associated distance of closest approach function is described in terms of its Chebyshev interpolants; then using the Chebyshev-Frobenius companion matrix, all of the roots, of the outermost is the distance of closest approach, can be found.Thanks to automatic subdivision and scaling functions, all numerical issues, including those of scale, can be handled.A guide to implementing ACPRAS in general and specifically for classical scattering problems has been reported in this work.
The technique was implemented in an open-source package written using the Rust programming language, hosted publicly on GitHub at [34].The ACPRAS technique was then implemented and tested in the RustBCA code, also publicly avaialble via GitHub [20].RustBCA, together with the ACPRAS algorithm, has allowed for the inclusion of realistic potentials in BCA codes for the first time.Using ACPRAS allows RustBCA to include any physically reasonable interatomic potential, guaranteeing correctness of the resulting scattering angle.
Examples showing the application of ACPRAS to problems involving ion-material interactions have been reported, providing evidence that attractive-repulsive potentials in BCA codes produce significantly different results than relying on historically used, universal screened Coulomb interaction potentials.Future studies using modern potentials, including potential direct comparisons to molecular dynamics, will show if the use of  attractive repulsive potentials in the BCA offers a compelling alternative molecular dynamics for low-energy implantation, reflection, and sputtering problems.

Figure 1 .
Figure1.Center of mass binary collision between two atoms.p is the imapct parameter, r 0 is the distance of closest approach, and θ is the center of mass scattering angle.

Figure 2 .
Figure 2. A selection of some historically important interatomic potentials.Apparent are the long-ranged nature of the Coulomb potential and the monotonicity of the screened Coulomb and Lenz-Jensen potentials (Kr-C, Born-Mayer).Note that this is shown in a shifted semilogarithmic form to exaggerate the short-range, high-energy behavior, with the zero position indicated by the dashed line.

Figure 3 .
Figure3.A selection of some advanced interatomic potentials.Bjorkas et al[14], a W-W potential, is a modified version of Derlet et al[17] that includes high-energy interactions.Marinica et al[16] is another low energy W-W potential, visibly unsuitable for high energy collisions (which should have a Coulomb like potential near the origin).Kr-C, the universal screened Coulomb potential is shown for comparison, as is the H-W potential from Wang et al[15].

Figure 4 .
Figure 4.This figure shows the failure of Newton's method by comparing solutions to the DOCA problem using Newton's method (lines) and the Chebyshev Proxy Rootfinder (filled circles).Here, R is the DOCA and p is the impact parameter.It is apparent that Newton's method often converges on the wrong real root or converges to a complex root with a small imaginary part.In contrast, the CPR rootfinder accurately determines all roots simultaneously, and, in the context of the DOCA problem, the outermost root can be directly computed to high accuracy.

Figure 5
Figure5.This plot shows the real parts of the complex roots of the DOCA function for the Lennard-Jones 12-6 function using Delves' method for complex rootfinding[22].Here, R is the distance of closest approach, a is the Bohr radius, p is the impact parameter, and epsilon is the reduced energy.

1 . 4 .
Calculate the Nth order Chebyshev gridpoints for the domain [a, b], x k , using equation (22) 2. Fill the values of the interpolation matrix I using equation (23) 3. Calculate the Chebyshev coefficients a j from the values of f (x) at the Chebyshev gridpoints x k and the interpolation matrix via equation (24) Initialize the Clenshaw-Curtis relation (equation (27)) 5. Calculate ˜( ) f x N by completing the algorithm in table 2 initalized in the previous step to find the values of b 0 and b 3 split [a, b] into two intervals [a, a + (b − a)/2] and [a + (b − a)/2, b] and return to step 1 with each new subinterval.(iv) From the Chebyshev coefficients a j , construct the Chebyshev-Frobenius companion matrix, ( ˜( )) A f x N (equation (31)).(v) Calculate the eigenvalues of the Chebyshev-Frobenius matrix, yielding the approximate roots of f (x) in [a, b]; finally, determine if the roots of the interpolant lie in the original domain [a, b]

Figure 6 .
Figure6.Scattering angles with a Lennard-Jones 6.5-6 potential found using ACPRAS.Epsilon is the reduced energy, p is the impact parameter, a the Bohr radius, and θ the scattering angle in the center of mass frame.Orbiting conditions are visible for ε = 5 × 0 −5 and 1 × 10 −4 at approximately 27 and 33 screening lengths respective.Here, the screening length used is the Lindhard screening length.

Figure 7 .
Figure 7. DOCA values for the Lennard-Jones 6.5-6 potential, found using ACPRAS.Epsilon is the reduced energy, p is the impact parameter, a the Bohr radius, and R the distance of closest approach.

Figure 8 .
Figure 8. ACPRAS applied to a distance of closest approach function for a 12-6 Lennard-Jones potential with σ = 1 Angstrom and ò = 1 eV, for p = 0.1 Angstrom and E r = 10 eV.The DOCA function is F(x) and the five domains are separated by plus sign markers.ACPRAS correctly identified the root to 6 significant digits,2 of which are shown.The order of each Chebyshev interpolant is shown in the legend.

Figure 9 .
Figure 9. ACPRAS applied to an advanced case of the distance of closest approach function.

Figure 10 .
Figure10.This plot shows a difficult function, equation(32), that has many zeros, some very close together, coming from both a periodic component and a polynomial component.ACPRAS was applied to find the roots of this function.Marked with stars are every zero ACPRAS identified.Subdivision borders are marked with plus signs, and each interpolant has a different color plotted on top of the original function in blue.ACPRAS identified all zeros to high precision, including the two very near zero, with 147 subdivisions.

Figure 12 .
Figure 12.Implanted ion depth distributions of hydrogen incident on a tungsten target.A comparison between an advanced cubic spline tungsten potential and the universal screened Coulomb Kr-C potential is shown, with the Kr-C potential leading to deeper implantation than the W-W/W-H potential.

Figure 11 .
Figure 11.Incident ion and displacet atom trajectories of 2 keV hydrogen (red) and helium (fuchsia) on a layered TiO2-Al-Si target.Titanium is marked in green, oxygen in purple, aluminum in black, and silicon in blue.The surface is at x = 0. Final positions of incident particles are marked with triangles.

Figure 13 .
Figure 13.Frenkel Pair (FP) production and implantation depth distributions for hydrogen on a layered titanium dioxide, aluminum, and silicon target.A comparison is shown between a Morse potential for the interaction between the incident ions and the target and a Lenz-Jensen potential.The titanium dioxide layer is 0.01 microns thick, the aluminum 0.04, and the silicon 1 micron.

Figure 14 .
Figure 14.Reflected hydrogen energy spectra from RustBCA simulations of 2 keV hydrogen noirmally incident on a tungsten target for the Kr-C potential and realistic W-W/W-H potentials.