Self-consistent ion beam analysis:
an approach by multi-objective optimization

MultiSIMNRA uses an objective function based on the calculation of quadratic differences ($\chi^{2}$ function). However, the simple combination of the standard $\chi^{2}$ function for each spectrum can result in a biased analysis if there are unaccounted model deviations since spectra with higher counting statistics receive a higher weight in the minimization process [29]. In addition, in numerical experiments it has also been observed that an adjusted weighting of the individual spectra results in faster convergence. Differences in counting statistics can occur either because of different experimental conditions (for example different solid angles for different detectors or different accumulated charges in different measurements) or because of different cross-sections involved in the physical processes [30, 31]. We compensate this unbalance by normalizing the $\chi^{2}$ function by its expected value for each spectrum.

where $max(M_{i},2)$ is used to avoid the null uncertainties.

The $\chi^{2}$ definition for minimization using MultiSIMNRA is:

where $M_{i}$ and $T_{i}$ are respectively the calculated and experimental value of counts in channel $i$ of a spectrum $s$, $S$ is the total number of spectra, $N_{s}$ is the number of channels in the fitting region of the spectrum $s$, $P_{s}$ is the number of fitting parameters in the setup configuration for the measurement of the spectrum $s$, and $P_{p}$ is the number of fitting parameter in the sample depth profile (note that $DoF=N_{s}-P_{s}-P_{p}$ for each spectrum). MultiSIMNRA also enables the user to change the parameter $a_{s}$ (which is equal to one by default) which is a parameter that serves to provide a intentional unbalancing to the contribution of a certain spectrum to the total $\chi^{2}$.

Our choice of an objective function based on reduced $\chi^{2}$ instead of the $\chi^{2}$ itself can be justified since each spectrum has a different expected value for the latter. Thus, each spectrum contributes differently to the global value of an objective function based on the $\chi^{2}$ calculation. This difference affects the balance of the objective function attributing different statistical weights to each spectra [30]. However, an objective function based on the reduced form of the $\chi^{2}$, e.g. normalized by the number of degrees-of-freedom, ensures that all the spectra have the same expected value (equals to unity) and consequently the same contribution to the global objective function. Thus, the use of the mean reduced $\chi^{2}$ as the objective function ensures that all spectra in the analysis have the same statistical weights in the minimization process, and that the expected value of Eq. 2 is also equal to one, independently on the number of spectra in the joint analysis.

Formally, this is a multi-objective optimization problem solved by the weighted-sum method [29]. By adopting the weights to normalize the individual objective functions we equalize the scale of comparison provided by the global objective function. On this terms, the $a_{s}$ of eq. 2 can be interpreted as a way to express preferences in the analysis, whilst leaving an arbitrary weight ($a_{s}$) to the analyst to be used as an ’ad-hoc’ quality parameter. A higher $a_{s}$ represents a higher preference [30].

This problem was addressed differently by Barradas et. al[32], who introduced a weight for each partial $\chi^{2}$ in the total $\chi^{2}$ calculation, in an attempt to correct for this bias. The basic idea is that, by assuming an objective function proportional to the squared difference of the experimental data and the simulated model, it enables the calculation of its expected value for each spectrum. This can be done as:

The authors argued that in a less-than-perfect fit, the quadratic difference grows faster than proportional to $I_{exp}$, so they introduced a weight for each spectrum of $I_{exp}^{1.5}$ in an attempt to balance the objective function. It is worth to point out that, by doing so, the authors introduced a biased term that privileges spectra with a lower counting statistics, since the punishment in the objective function is stronger for spectra with higher integrated counts than for spectra with lower integrated counts (when compared to the unbiased weight of $I_{exp}$). The authors also included an extra term as a penalty for the increment of the number of layers ($n_{lay}$) in the sample description [33]. In this case, the combination of the partial $\chi^{2}$ for a balanced objective function is as follows:

The NDF’s objective function is also implemented in MultiSIMNRA for comparison reasons.

Other objective functions may be available in NDF, mainly for the Bayesian inference method of uncertainty estimation [34]. We refer to eq. 4 as an alternative example, and as the only version published until now for the NDF’s objective function.

One should keep in mind that due to noise, imperfections in physical modeling, measuring artifacts, etc. it may not exist a unique set of parameters that minimizes the global and all the individual objective functions simultaneously. The Pareto optimal condition defines there is a set of parameters that minimize the global objective function, even if not presenting an optimized solution for at least one of the single objective functions. Therefore, the Pareto optimal is the best solution possible within the average agreement of all measurements, and the choice of the weights may affect the Pareto optimal [29].

Perhaps, this is the most critical argument in favor of adopting the weighted-sum method for this multi-objective optimization. Suppose all different spectra have about the same weight when combining the different measurements. In that case, the unknown or unaccounted systematic errors are averaged in the analysis, and no particular deviation of the forward model is preferred over the other. Otherwise, using an unbalanced objective function, spectra with higher statistics may dominate the analysis, even though they may have more flawed theoretical understanding, while neglecting those with fewer counts but possibly better modeled.

At the end of the minimization process, our code can also perform a statistically sound evaluation of the uncertainty of the fitted parameters. This evaluation is done using a Monte-Carlo calculation to explore the values of $\chi^{2}$ in the neighborhood of all fitted parameters. For the determination of a confidence region with $1-\gamma$ confidence level, random variations of the parameters are accepted only if the changes in the $\chi^{2}$ are lower than a critical value ($\chi^{2}_{\gamma}$) determined by the $\chi^{2}$ probability density function and the confidence level ($1-\gamma$) [35]. The accepted values are stored in a buffer file to account for the statistics in the near region of the Pareto optimal. The confidence level is determined by the user, and considered as 68.3% (one standard deviation) by default.

The critical value $\chi^{2}_{\gamma}$ is defined by the equation:

where $f(z;n)$ is the $\chi^{2}$ probability density function, $z$ is the set of parameters and $n$ is the number of parameters. That is:

When the buffer file has enough statistics, the data are used to calculate the variances and the covariances, as well as, the confidence limits for the depth profile considering the uncertainties of all parameters. The size of the buffer is user defined (3000 is set by default, but other values are also available, like the value of 40000 as recommended by the Guide of Uncertainties and Measurements in its appendix on Monte-Carlo calculation of uncertainties).

Interpretation of a set of experimental measurements requires to choose a model (ie. the depth discretization of the sample) which will fit the data to some extent out of a family of possible alternative models (like different concentration profiles or layer settings). The tempting choice of selecting the model providing the best fit (lowest $\chi^{2}$) is likely to over-fit the experimental data [36].

In MultiSIMNRA the number of layers, and consequently the number of fitting parameters, is defined by the user. Layers are not automatically created, albeit there are tools to slice, scale, merge and insert layers, making the layer management easier. However, it is known that an exaggerated number of fitting parameters leads to a reduction of the objective function, but it also introduces noise (artifacts) in the final solution by an over-fit of the experimental data.

A principled guideline - automatically incorporating Ockham’s razor - is given by a Bayesian model selection approach [37], which balances the goodness-of-fits of the models (or more precisely: the posterior probability distribution of the model parameters) with their complexities. The relevant quantity for this is the so called evidence [38] or marginal likelihood.

Unfortunately, despite recent algorithmic progress in computing this quantity (e.g. using Nested sampling [38, 39]) it is still time consuming and cumbersome to compute since the plain Markov Chain Monte Carlo (MCMC) is too inefficient. For this reason, several easier accessible approximations have been suggested. In the limit of many data and the (in that case reasonable) assumption of a multivariate Gaussian posterior distribution, the Akaike information criterion (AIC) can be derived. It depends - besides $\chi^{2}$ - only on the number of parameters $n$ and the number of data-points $N$ and is given by [40]:

The $AIC$ estimator is a asymptotically unbiased quantitative measure of Occam’s razor. The minimum value of the $AIC$ coefficient represents the simplest description of the sample with no loss of information. Note that eq. 6 is proportional to the number of fitting parameters thus a reduction of $\chi^{2}$ is not necessarily followed by a reduction of the $AIC$ coefficient. Hence, a reduction of the $\chi^{2}$ function not followed by a reduction of the $AIC$ estimator, represents an over-fitting of the experimental data.

Comparing the $AIC$ estimator calculated for different description of the sample enable the quantitative comparison of different fitted models (differences regarding the number of fitting parameters) in order to judge the level of empirical support of each model (see table 1).

At the time of this publication four algorithms for minimizing the objective function were implemented and are available as an option for the user. The first option is the Nelder-Mead’s Simplex (NMS) [41, 42], a direct search method, computationally simple which does not require calculation of derivatives. This is the same algorithm implemented in the fitting tools of SIMNRA. It is known that NMS may become less efficient for large number of fitting parameters and that it eventually stagnates in local minima. To correct for these issues, in our implementation, adaptive Simplex-parameters was implemented [43] and also oriented restarts [44], which are two features that aim to improve the effectiveness of the NMS algorithm. The second algorithm is a modified version of the Levenberg-Marquardt algorithm (LMA), which is a robust and fast deterministic optimization algorithm [45], but very sensitive to the initial parameters and does not guarantee global minimization. The third method is the evolutionary annealing-simplex algorithm [46], which is a probabilistic heuristic global optimization technique that joins ideas from different methodological approaches, where a generalized downhill simplex methodology is coupled with a simulated annealing procedure. The fourth method is a modified version of the evolutionary algorithm, Differential Evolution (DE), called Adaptive Differential Evolution with crossover rate repair (RCR-JADE) [47]. An interesting characteristic of DE is its robustness to the initial guess, but it is sensitive to the search-space boundaries.

An initial guess for the values of concentrations of elements in the sample are informed by the user to the MultiSIMNRA software in the form of a table, together with its range limits. The program also includes some tools to handle the sample description (like scale layer, slice layer, remove layer, etc.) that make the user interaction easier.

A performance test was done to compare the convergence time of all four optimization algorithms. The test consisted in using the MultiSIMNRA code to fit a set of four simulated spectra of an idealized sample (with Poisson noise added) in different configurations of energy and scattering angle. The sample consisted of a three layered model with a total of seven fitting parameters. The initial guess to start the fitting procedure was the same for all algorithms and chosen to be very far from the known solution and no human intervention was made during the optimization process. All the three algorithms converged to the known solution (with the exception of the Levenberg–Marquardt algorithm that slowly converged to a configuration just close to the known solution, probably due to the derivative instabilities generated by the inserted noise). The time evolution of the objective function for each optimization algorithm is presented in fig. 1(a).

It is worth to mention that, considering the artificial noise added to the simulated spectra, a good solution was obtained for three of the four implemented algorithms in a convergence time of less than 2 minutes. Starting with a good initial guess this time can be drastically reduced.

Another performance test was executed to compare the convergence time of the four optimization algorithms in more complex analysis. This time, the model included 32 free parameters (26 for the target description and 6 for the experimental setup). Fig. 1(b) shows that all the four algorithms converged to the minimum, but Leveberg-Marquart persists as the less effective algorithm. Nealder-Mead’s simplex initiated the convergence with a slight delay compared to differential evolution and the evolutionary annealing simplex. Quickly converged so that the three algorithms converged in approximately the same total convergence time.

All performance tests were done on a desktop with the last generation of x86-64 processor.

The basic motivation behind the self-consistent approach is to integrate the maximum of information obtained from different ion beam techniques. In this sense, the use of prior information in the analysis is very important and can considerably increase the quality of the solution. Butler [11] proposed a three steps approach in the analysis to reduce ambiguities: i - the use of multiple detectors; ii - the information of chemical formulas; iii - the information on diffusion profiles (what Butler called thermodynamics information). We can summarize these three steps into two more generic steps: i - the combination of multiple measurements (or the Alkemade statement [12]); ii - the use of prior information. The first is a generalization of the first Butler’s rule, and the latter is a generalization of the second and the third Butler’s rules combined: any information on the chemical state of the elements in the sample, and on the diffusion profile can be considered prior information, as is the information on the detector electronics and experimental configuration.

Prior information is used in MultiSIMNRA by the definition of constraints through the definition of Links. A Link is a feature of MultiSIMNRA that defines the relationship between the fitting parameters. The definition of Links is important not just because of the reduction of the convergence time by the direct reduction of fitting parameters, but also to establish boundary conditions and the reduction of degrees of freedom for the fitting process, taking advantage of a prior knowledge on certain parameters. The definition of Links tying the fitting parameters increases the accuracy of the analysis and is a way to deal with problems releated to the curse of dimensionality [48].

The accuracy of the simulation codes and of fundamental physics models are crucial not just for the correctness of the obtained results but also for the algorithm convergence in the objective function optimization. Deviations in the modeling of physical phenomena become more evident when combining the information contained in different spectra due to inconsistencies in the simulation of the different conditions [49]. All this may lead to biased result or misinterpretation of data. Thus, the accuracy of the codes must be continuously confronted with experimental data. In a software inter-comparison promoted by IAEA [50, 51], it was demonstrated that the main differences on the simulation results of the most used codes lies on differences in the databases of fundamental physics data, such as stopping powers, cross-sections, straggling and multiple or plural scattering. Consequently, the quality of fundamental physical data are crucial for a self-consistent analysis of multiple measurements.

According to reference [52], the semi-empirical approach adopted in SRIM is still the best approximation of the available experimental data for elements with $1\leq Z\leq 92$, even though, according to the SRIM changelog, no major updates on its database were made since the 2003 version. Merit should be given to some theoretical developments [53] but regarding the level of agreement to experimental data aiming at accurate IBA calculations and computational performance, the semi-empirical approach of SRIM is still unbeaten. Thus, new accurate and traceable stopping power measurements and continuous updates of stopping power databases are highly necessary.

The situation for the cross-sections database as provided by SigmaCalc [7] is good concerning its continuous update, and the efforts to increase the coverage of elements and the intensive tests by means of benchmarking experiments [54]. However, concerning the uncertainties and traceability, the database is still scarse [55], and restricted to a single element (Carbon). Besides, no uncertainty analysis is made on the energy reference for such database.

MultiSIMNRA has the option to use SRIM [56] stopping powers and SigmaCalc [7] evaluated cross-sections for all spectra simulations, ensuring consistency in all physical models. Despite the good internal database for stopping powers of SIMNRA, the use of SRIM data is favored for traceability purposes, since its results are intensively confronted to experimental data and vice-versa [57]. For the cross-sections, the evaluated data provided by SigmaCalc can be obtained in [58, 59] and incorporated in the SIMNRA reaction list. Starting with the new SIMNRA version [8], a built-in version of SigmaCalc is implemented and the download from an external source is no longer necessary.

The use of experimental cross sections in situations where there are no evaluated data may rise the problem of how to combine different accuracies in the fundamental databases to calculate the statistical weights in the objective function (cross-sections of different elements or scattering angles, for instance). This is a topic not yet adequately addressed in the literature and therefore not considered in MultiSIMNRA.

To illustrate how the ambiguities inherent to an RBS analysis can be solved, Butler [11] proposed an exercise based on the fit of simulated data for an hypothetical sample of metallurgical interest. The sample contains a depth profile of the elements O, Al, Ni, and Cr (see fig. 3). What is important about this exercise is that the RBS simulation provides a spectrum with overlapping signals from different elements, that opens space for ambiguities in the interpretation of the data. In this sense, Butler also proposed three steps in the analysis to avoid the ambiguities and to extend the applicability of ion beam techniques (see section on the use of prior information). This test was also successfully accomplished by Barradas using the NDF software [33] following the three Butler’s steps. Here, we adopted a different approach: For a more realistic simulation, a Poisson noise was added to the data and, instead of using the chemical state information, we followed the Alkemade [12] statement (at least $N-1$ spectra may be necessary to solve $N$ variables). Thus, without any previous knowledge of the sample, all information was extracted only by the self-consistent analysis.

Fig.  2 shows the simulated data (continuous line) of Butler’s hypothetical sample at four different experimental conditions with the addition of Poisson noise, and the superposition of the fitted curves obtained by MultiSIMNRA (discontinuous lines). The spectrum in fig.  2-(a) was generated at the same conditions as in the Butler’s original paper: incident helium ion energy, $2000$ keV; detector angle, $170^{\circ}$; calibration, $5$ keV/ch; intercept, $25$ keV, and resolution, $15$ keV. It was calculated for a detector with $1$ msr solid angle and $10$ $\mu$C integrated charge. The spectrum in fig.  2-(b) was generated at the same conditions as fig.  2-(a) but with the detector placed at $120^{\circ}$ scattering angle. Data in fig.  2-(c) and (d) were generated in the same conditions as fig.  2-(a) and (b), respectively, but with $3050$ keV incident helium ion energy to take advantage of the resonant non-Rutherford scattering cross-section for the ¹⁶O($\alpha$,$\alpha$)¹⁶O reaction. The cross-section used for the ¹⁶O($\alpha$,$\alpha$)¹⁶O reaction was calculated using SigmaCalc  [7].

The depth profiles obtained by MultiSIMNRA are presented in fig. 3. The dotted lines are the data used to generate all spectra presented in fig. 2. The continuous lines are the profiles obtained using MultiSIMNRA and the dashed painted region gives the confidence limits for the result. It is possible to observe the good agreement between the original profile and the fitted model.

One aspect worth mentioning is that uncertainty bands represent all possible profiles that passed the inference test for uncertainty estimation (see the section on uncertainty calculation). Thus the profiles and uncertainty bands are highly dependent on the number of layers of the model (see the section on model selection based on Ockham’s razor). Therefore, it is not easy to represent the variances and covariances between all fitted parameters graphically. However, MultiSIMNRA provides the correlation matrix and total amounts of the elements in the final report. The latter presents minor uncertainties compared to the uncertainty bands because it considers the correlations properly.

To demonstrate the advantages of joint processing of spectra obtained by multiple ion beam techniques, we present the analysis of a thin-film iron samples using heavy-ion time-of-flight RBS data and time-of-flight ERD data. The sample is a thin iron film with 0.7at.% of tungsten, used as a model system for the EUROFER low activation alloy intended to be used in future fusion power plants. The film was sputter deposited on silicon to investigate the preferential sputtering of iron and the tungsten accumulation at the surface during low energy deuterium irradiation [24]. Thus, top layers’ depth profile characterization in a nanometer range is of high importance to this study, either for heavier (Fe and W) or for lighter (C and O) elements.

The tungsten accumulation at the surface of the iron films is of high interest in the design and construction of fusion reactors. For our purpose, it presents a stringent test of multiple spectra analysis by MultiSIMNRA, using two different scattering techniques in two different laboratories: the ToF-RBS data was obtained at the Max-Planck Institut für Plasmaphysik in Garching, Germany [24], and the ToF-ERD data was obtained at the Ruder Boscovic Institute in Zagreb, Croatia [60]. Fig. 4 shows experimental data and the simulated spectra for the sample’s adjusted model.

The top plots in fig. 4 are the ToF-RBS data: energy spectra of scattered 11.5 MeV silicon ions with incidence perpendicular to the surface on the left, and scattered 5.0 MeV silicon ions with an incidence angle of 60^∘ on the right. In both cases, the detector was placed at 150^∘ scattering angle [24]. The bottom plots show the ToF-ERD data for each recoiled ion species: C, O, and Fe, respectively. The incident beam was 20.0 MeV iodine ions at an incidence angle of 70^∘. The detector was placed at 37.5^∘ scattering angle. The scattering angle is always in respect to the incident beam [60].

The only assumptions used are the full oxidation of the iron and of the tungsten at the surface. These assumptions were included in the analysis as prior information, and imposed as a constraint by setting a Link (see section on the use of prior information) between the elements, considering Fe₂O₃ and WO chemical states at the surface. The final result of the depth profile model for elemental concentrations is finally presented in fig. 5 and compared with sputter-XPS data (the two discontinuous lines, differing only on the sputtering rate to correct the depth scale).

The excellent agreement between the self-consistent analysis of the ion beam methods and the sputter-XPS measurement is evident. The relative higher tungsten concentration in the XPS analysis is probably a sputtering enhancement caused by preferential sputtering of iron due to the argon primary beam (10 keV energy and incidence angle of 20^∘), hence an XPS measurement artifact. This points to another positive feature of the ion beam methods: direct assessment of depth profile information with minimum sample modification. For the sputter-XPS data, the depth scale was obtained using sputtering rates as calculated by SDTrim.SP [61] (version 6.05) SRIM [56] (version 2013).