Differential evolution optimization of Rutherford backscattering spectra Open Access

Rutherford backscattering spectrometry (RBS) is applied to analyze the elemental composition depth profile of thin films and near-surface regions. RBS is mainly recognized for the high sensitivity to heavy elements on a light substrate. The underlying physical phenomena that govern Rutherford backscattering spectrometry are well understood. The RBS spectrum can be predicted if the experimental conditions and the composition depth profile of the sample are known. However, the inverse problem, which is to determine the composition depth profile from an experimental RBS spectrum, is more challenging. The first challenge is that in principle, multiple composition depth profiles may lead to an identical RBS spectrum.¹ Thus, to fit an RBS spectrum, the analyst must restrict the solution space by defining the sample structure and the possible elements or by narrowing the solution space by imposing the consistent analysis of multiple spectra. The second challenge is that the sample composition parameters often appear as highly correlated in the optimization function.

There is a need for an algorithm that can find the sample parameters irrespective of the correlation. The potential of simulated annealing² in this respect was recognized.³ The simulated annealing algorithm has been successfully implemented^4,5 and is being used to analyze RBS spectra. In the early stage of optimization, simulated annealing performs well, while in the last stage, it appears to be less efficient. Therefore, the last stage of simulated annealing optimization is conventionally complemented with a local search routine. Artificial neural networks are attractive alternatives to aid with the analysis of Rutherford backscattering spectra.⁶ One of the main advantages of artificial neural networks is that the analysis is retrieved instantaneously, and thus, it is ideally suited to analyze large sets of data for which only a few parameters vary. Unfortunately, a new network must be constructed if the parameter search space changes or if the experimental conditions are varied. Besides, the use of artificial neural networks to analyze RBS data is found to add around 1% to the uncertainty budget.^6–8 After the analysis with an artificial neural network, one conventionally uses a traditional analysis code to verify the correctness and to improve the accuracy.^9,10

It is accepted that different optimization algorithms are better suited to solve different problems.¹¹ In particular, the differential evolution algorithm has marked a breakthrough in meta-heuristic optimization.¹² In the recent years, many more new and promising meta-heuristic algorithms have been proposed.^13–18 It is, therefore, challenging to identify the optimization algorithm that is best suited for a specific problem.^19–22 

In the case of Rutherford backscattering spectrometry, an exploratory study of different optimization strategies is reported.²³ However, a more detailed study to assess the performance of the various optimization algorithms is needed. In this work, we investigate the potential of differential evolution to fit Rutherford backscattering spectra. The code is named Ruthelde. We demonstrate that the implementation is robust, that it does require minimal human supervision, and that a very good precision and accuracy for the extracted parameters can be achieved.

The differential evolution (DE) algorithm was first proposed by Storn and Price in 1995 as a simple and efficient heuristic for global optimization over continuous spaces.¹² The authors show that the differential evolution algorithm is able (1) to handle non-differentiable, nonlinear, and multi-modal fitness functions, (2) is easily parallelizable, (3) is easy to control and set up, and (4) has good convergence properties.

The DE algorithm, as a population based search heuristics, at each generation $G$⁠, holds a population $SG$ of $N$ solutions. Each solution with index $j$⁠, also called an individual, refers to a $P$ dimensional parameter vector $XjG={x1,jG,x2,jG,…,xP,jG}$⁠. The elements of the parameter vector represent the fit parameters of a particular problem. The individuals of the initial generation (⁠$G=0$⁠) are chosen randomly from the entire parameter space. The optimization in differential evolution is realized as follows.

For each target vector $XjG∈SG$⁠, a so-called mutant vector $VjG$ is generated according to

where $F$ is the scaling factor (⁠$0≤F≤1$⁠) and $r1,r2,r3∈[1,2,…,N]$ point to randomly selected individuals of $SG$ with $r1≠r2≠r3$⁠.

Successively, in the crossover process, a new vector called trial vector $UjG={u1,jG,u2,jG,…,uP,jG}$ is derived from the target vector $XjG$ and the mutation vector $VjG$ by

with $rand(0,1)$ being a random number in $[0,1]$⁠. $CR$ is referred to as the crossover probability (⁠$0<CR<1$⁠) and represents the second input parameter of the DE algorithm.

In a final step, a comparison between the target vector and the trial vector is done according to the value of their fitness function $f(X)$⁠. The vector with the better (higher) fitness value is taken into the new generation according to

Thus, a new generation $SG+1={XjG+1:j=1,2,…,N}$ is produced. The average and the best fitness of population $G$ are defined as

It follows from Eq. (3) that $f¯G+1≥f¯G$ and that $fbestG+1≥fbestG$⁠.

The operations of mutation, crossover, and selection are continuously repeated until a predefined termination criterion is reached. The latter one could be a minimum value for average or best fitness, a fixed computation time, or a certain number of iterations. A characteristic that makes differential evolution work so successfully is that mutation is based on differences between individuals of the entire population. This allows each fit parameter to self-tune and gives an appropriate reduction in the magnitude as the optimization proceeds and convergence is approached. This functionality parallels the approach used successfully in evolution strategies where the mutation variances are self-tuning.²⁴ 

In its original form, the differential evolution algorithm has three fixed input parameters determining its performance: the population size $N$⁠, the scaling factor $F$⁠, and the crossover probability $CR$⁠. Over the years, several optimizations and derivations to differential evolution are proposed. Noteworthy are changes in the initialization,^25–27 modifications of the mutation,^28–30 variations in crossover^31–34 and selection,^35,36 and hybridization with other search algorithms.^37–39 Many efforts have also been spent to implement one or more DE parameters in an adaptive manner.^40–42 A comprehensive review on the development of differential evolution during the last two decades is given in Ref. 43.

The analysis of RBS spectra is commonly done by forward simulation and comparison of the simulated spectrum to the measured one. The fit parameters are defined as the input values that are used to simulate the spectrum. Common fit parameters in Rutherford backscattering analysis are the target model comprising the thicknesses of all layers as well as the elemental ratios in each layer and the experimental setup given by the detector calibration, detector resolution, and the detector solid angle in conjunction with the total applied charge.

The fitting of Rutherford backscattering spectra is typically achieved by minimizing the sum of squared residuals (⁠$χ2$⁠) between the simulated and the experimental spectra,

where $xiexp$ and $xisim$ represent the height of the experimental and simulated spectrum in channel $i$ and $Nc$ the total number of channels in the spectrum. Unfortunately, the sum of squared residuals is not normalized, and thus, it is less intuitive to the user since it depends on the shape and noise of the experimental spectra and on the number of channels. Therefore, we use a normalized fitness function, which is $χSG2/χ2$⁠, with $χSG2$ being the sum of squared residuals between the experimental spectrum and the spectrum smoothed by applying a Savitzky–Golay filter of width equal to the detector resolution. The later one can be regarded as the reference fit since due to its averaging nature, it generates a very smooth representation of the input data points. Thus, the fitness function $χSG2/χ2$ cannot overcome 1, and reaching any value close to 1 may be used as stopping criteria for the code.

The Ruthelde software is an implementation of differential evolution as described above by Eqs. (1)–(3) and without the introduction of adaptive parameters. Yet, as a surplus, we implemented a mechanism to retain the diversity in the population. This is done by introducing a threshold value $TR∈[0,1]$⁠. With every generation, $10%$ of the individuals in the population except the best individual is replaced by randomly initialized new ones if $f¯G/fbestG>TR$⁠.

The four parameters (⁠$N$⁠, $F$⁠, $CR$⁠, and $TR$⁠) that control the differential evolution algorithm are optimized by repeatedly fitting different RBS spectra and monitoring the performance. The performance for each set of four parameters is assessed through the steepness of the fitness improvement as a function of the number of generations and through the best final fitness after a fixed number of iterations. For each set of four parameters, the fitting is repeated 20 times with random seeds. The following parameters for the differential evolution algorithm are found to be versatile for fitting RBS spectra: $N=20$⁠, $F=0.6$⁠, $CR=0.9$⁠, and $TR=0.95$⁠. Based on our exploration, it appears that low values for $N$ (⁠$N<15$⁠) lead to a poor diversity within the population, which increases the chance of the algorithm being trapped in a local optimum. On the other hand, high values for $N$ (⁠$N>25$⁠) lead to a long computation time for each generation. Furthermore, we found that the differential evolution optimization performs well for a broad range of values for the scaling factor (⁠$F$⁠) and the crossover probability (⁠$CR$⁠), namely, $0.45<F<0.65$ and $0.65<CR<0.95$⁠. The threshold level (⁠$TR$⁠) is meant to retain a certain level of diversity within the population at the later stage of the fitting. The value of TR should be high enough (⁠$TR>0.9$⁠) not to introduce too much randomization in the early stage of optimization.

Irrespective of the efficiency of the search algorithm, the evaluated fit result can only be as good as the accuracy of the underlying spectrum simulation. If the physical processes are not correctly implemented, then the forward simulator will generate in-accurate spectra, and thus, the extracted fit parameters will be incorrect. Therefore, special care needs to be taken in the proper implementation of forward simulation.

The physics of Rutherford backscattering spectrometry is described with great accuracy by assuming elastic scattering with a central force field and by considering the gradual energy loss of the ions as they travel through the material.^44–46 We implemented the physics of RBS as described in the documentation of the Simnra code.^47,48 The electronic stopping and the nuclear stopping are implemented according to Ziegler–Biersack formalism.⁴⁹ The screening is calculated with the Andersen model.⁵⁰ The straggling is modeled as a Gaussian broadening as proposed by Bohr.⁴⁴ The propagation of straggling in thick layers is included. In the present implementation, neither the simulation of pileup nor the treatment of surface roughness is yet included. The integration of those and other features to advance the simulator is subject to ongoing code development.

To test the accuracy of the Rutherford backscattering implementation, we compare in Fig. 1 simulations based on the present code with simulations generated with Simnra. We adopt conditions that are close to experimental conditions that are used further on. The energy of the primary He$+$ ion beam is 1.52 MeV. The scattering angle is 170$°$⁠. Measured from the sample normal, the incidence angle is $−$11$°$ and the exit angle is 21$°$⁠. The detector solid angle is 0.42 msr. We assume a detector energy resolution of 15 keV and a bin-width of 1.5 keV/channel.

First, we consider the RBS spectrum for an Au on SiO$2$ on a Si substrate sample [Fig. 1(a)] to represent the very thin film case.⁵¹ The areal densities are $12×1015$ atoms/cm$2$ (approximately 2 nm) for the gold layer and $6.6×1017$ atoms/cm$2$ (⁠$≈$ 100 nm) for the silicon dioxide layer. The various signal edge positions in spectra generated with Ruthelde (black dots) and with Simnra (red line) agree to better than 100 eV. The analyses of the gold areal density in the two simulations using the surface approximation⁵² agree to better than 0.02%.

Second, in Fig. 1(b), we consider the Rutherford backscattering for a SrTiO$3$ layer of $625×1015$ atoms/cm$2$ (⁠$≈100$ nm) on a silicon substrate. The signals from backscattering on oxygen, on titanium, and on strontium are indicated in the figure. The intensities for the different elements obtained with Ruthelde and Simnra agree to better than 0.03%.

Third, Fig. 1(c) shows the Rutherford backscattering spectrum for an In$2$Ga$2$Zn$1$O$7$ layer of $1200×1015$ atoms/cm$2$ (⁠$≈140$ nm) on a silicon substrate. Again, one observes that both Ruthelde and Simnra agree very well. Note that the simulations possess the information about the contributions of the different elements even where they appear to overlap in the total spectrum. See, for example, the signals from backscattering on zinc, gallium, and indium indicated in the figure.

The differences between the output from Ruthelde and from Simnra are small and are similar to the differences between other state-of-the art implementations for RBS analysis.⁵¹ Therefore, we validate Ruthelde as a correct forward simulator for Rutherford backscattering. The three model cases are also selected as they represent different degrees of difficulty to fit corresponding experimental spectra as there is increasing overlap of the signals.

The simulations in Fig. 1 contain the element-specific probabilities to detect backscattered particles as a function of energy. We use a Monte Carlo code to generate a draw from the element-specific probability distributions to generate realistic test spectra. The accumulated charge is determined by the number of generated random numbers. For different synthetic spectra, the accumulated charge ranges between 1 and 1000 $μ$C. For each model system and for each chosen value of accumulated charge, we produce 30 independent synthetic spectra.

One such synthetic spectrum for each model system and assuming a detector solid angle of 0.42 msr is shown in Fig. 2. The synthetic spectra, shown as the black data points, correspond to an accumulated charge of 31.6 $μ$C for the Au [Fig. 2(a)] and the SrTiO$3$ [Fig. 2(c)] thin film systems. The synthetic spectrum for the In$2$Ga$2$Zn$1$O$7$ model system [Fig. 2(e)] is obtained for an accumulated charge of 56.2 $μ$C.

Additionally, we added another model system with a more complex, multi-layered target model, namely, 130 nm SrTiO$3$⁠, on top of 90 nm LaFeO$3$ on top of 145 nm SiO$2$ [Fig. 2(g)]. A grazing incidence geometry with an incident angle of 70$°$ and a scattering angle of 100$°$ was chosen. The solid angle and the applied charge are 1.0 msr and 25 $μ$C, respectively.

Note that the number of events for the different elements is monitored during the creation of the draw. For example, the number of backscattered particles related to the model system [Fig. 2(c)] is 34 666 events for strontium, 11 508 events for titanium, and 4064 events for oxygen. The knowledge about the intensity of the various contributions allows us to estimate their respective statistical uncertainties. Interestingly, the exact number of events for a specific contribution is different for the different independent synthetic spectra, even if the target accumulated charge is the same.

In this section, we investigate the performance of the differential evolution algorithm to find the elemental profile from a Rutherford backscattering spectrum.

The synthetic spectrum in Fig. 2(a) is analyzed with a two-layer model. In total, seven parameters need to be optimized: the areal densities of the two layers, the atomic ratio for the SiO$2$ layer, the detector gain (keV/channel), the detector offset, the detector resolution, and the charge solid-angle product (charge). A very good agreement between the synthetic spectrum and the fit after 500 generations (corresponding to 10 000 evaluations at the given population size $N=20$⁠) is observed in Fig. 2(a).

In this example, the search space for the Au thickness ranges from 0 to 90 nm. The thickness of the SiO$2$ layer is allowed to vary between 15 and 180 nm. In general, the parameter search space does not need to be narrowly confined. For instance, it is generally not needed to impose boundaries for the elemental ratios in a layer containing multiple elements. The boundaries for the layer thicknesses generally can range from 0 to 5 times the actual value. It is seen from Fig. 2 that the initial convergence is very efficient. In contrast, the boundaries for the detector energy calibration need to be defined more carefully. We have observed that a too wide search space leads to the algorithm being trapped in “non-physical” solutions. To our experience, the boundaries for the detector gain are best defined as plus or minus 10% of the expectation value.

Since the differential evolution algorithm is based on a stochastic process, the performance and the capabilities of the algorithm may be investigated by repeating the fitting multiple times and by inspecting the variance of the obtained parameters. The analysis of the spectra is repeated 32 times. The plot in Fig. 2(b) shows the relative standard deviation of a selection of fit parameters as a function of the number of evaluations. One observes that the relative standard deviation of the obtained parameters decreases as more evaluations are allowed and that some parameters converge faster and better as compared to others. For example, the detector energy calibration converges to a standard fitting error of less than 0.001% after 5000 evaluations. The charge solid-angle product has a slightly larger standard deviation. This can be understood since the charge solid-angle product depends both on the amplitude of the silicon signal in the spectrum as well as on the detector gain. Furthermore, a value for the areal density of oxygen and gold in the sample with a standard fitting error well below 0.01% is obtained after 7500 evaluations.

The synthetic spectrum for Sr$x$Ti$y$O$z$ in Fig. 2(c) is analyzed with a single-layer model. The seven parameters to be optimized are the areal densities of Sr, Ti, and O; the detector gain (keV/channel); the detector offset; the detector resolution; and the charge solid-angle product (charge). The areal density of the Sr$x$Ti$y$O$z$ layer is allowed to vary between 0 and $3×1018$ atoms/cm$2$⁠. The fit with a random initial guess and after 12 000 evaluations is shown as the red line. In Fig. 2(d), the convergence of the fit parameters as a function of the number of evaluations is shown. Remarkably, repeating the fit 30 times (with random starting conditions) results in values for the different fit parameters with a spread of less than 0.01%. Therefore, we demonstrate that differential evolution is a suitable algorithm to fit Rutherford backscattering spectra to determine the composition profile of compound layers with a high precision.

The synthetic spectrum for In$w$Ga$x$Zn$y$O$z$ in Fig. 2(e) is analyzed with a single-layer model. Eight parameters need to be optimized: the areal densities of In, Ga, Zn, and O; the detector gain (keV/channel); the detector offset; the detector resolution; and the charge solid-angle product (charge). The fit with a random initial guess and after 12 000 evaluations is shown as the red line. The convergence plot in Fig. 2(f) shows that repeated fits reproduce the same result with a spread of less than 0.01% for the detector gain, the charge, the indium and the oxygen areal densities, and with a spread of around 1% for zinc and gallium. The performance to fit most parameters is in line with the results for the other model systems. It can be understood that the refinement for gallium and zinc is less unique, as their signals are highly overlapping. In other words, the fit precision indicates how precisely a specific parameter can be extracted from a specific spectrum. The fit precision may be affected by signal overlap or by inaccuracies of the model.

Figure 2(g) shows the synthetic spectrum and the best fit for the multi-layered target model as described above (130 nm SrTiO$3$/90 nm LaFeO$3$/145 nm SiO$2$/Si). In total, 12 parameters need to be fitted by the algorithm. Figure 2(h) displays the convergence of the fit parameters. It is observed that the algorithm converges to a good fit despite the complex structure and many overlapping contributions in the spectrum structure. The example demonstrates that it is imperative for each new sample model to assess the convergence of the parameters as it is introduced here.

We demonstrated that a repeated refinement of a synthetic spectrum does produce parameter values with a very small dispersion. However, one must differentiate between the fit precision and the analysis uncertainty.

For example, the repeated analysis of the Sr$x$Ti$y$O$z$ spectrum in Fig. 2(c) yields areal densities for strontium, titanium, and oxygen that vary by 0.01% only. For the titanium areal density, the resulting value evaluates to $Ti=(126.324±0.007)×1015$ atoms/cm$2$ after 32 fits of one and the same spectrum. One must realize, however, that the spread on the values only reflects the fit precision, i.e., the capability of the algorithm to find this result.

In contrast, the artificial spectrum in Fig. 2(c) was derived from an input value of $Ti=125.0×1015$ atoms/cm$2$⁠. At first sight, it implies a significant deviation between the fit and fit precision on the one hand and the expectation value on the other hand. The reason for the apparent discrepancy is that the artificial spectrum was generated in a Monte Carlo process, which reflects the natural counting statistics. The two-sigma counting statistical uncertainty for a titanium signal of $Ti=125.0×1015$ atoms/cm$2$ (=11 508 events) is 215 events (⁠$=2×11508$⁠) or, equivalently, 1.9%. The latter value represents the minimum analysis uncertainty of the fit due to the counting statistics.

We plot in Fig. 3(a) as blue data points the titanium results for 32 independent synthetic spectra with nominally the same accumulated charge of 31.6 $μ$C and a detector solid angle of 0.42 msr. The horizontal axis indicates the number of titanium events in a specific synthetic spectrum. As expected from Poisson statistics, the different spectra contain a different amount of titanium events with a two-sigma of around 215 events. The vertical axis shows for each synthetic spectrum the refined areal density for titanium. The two-sigma fit error, as explained in Sec. VI, is plotted on the different data points. The fit error for one specific spectrum appears to be comparable to the size of the symbols. In contrast, the uncertainty to which the parameters of the sample can be refined is dominated by the statistical variations between the spectra. Thus, the vertical black error bar in Fig. 2(c) represents the true uncertainty on the analysis. Therefore, the correct presentation of the result of the analysis of the synthetic spectrum in Fig. 2(c) for the titanium areal density reads $Ti=(126.3±2.5)×1015$ atoms/cm$2$⁠.

Generally speaking, a proper way to estimate the uncertainty for a given RBS analysis is to generate a collection of independent synthetic spectra with a Monte Carlo method based on the best fitting model and to analyze the resulting spread of the analyzed parameters. All the uncertainties that are expected to contribute to the final result must be included during the creation of the collection of synthetic spectra.

While the above method allows one to estimate the uncertainty for a given experiment, the method can also be used to investigate the uncertainty as a function of the accumulated charge. We have created a collection of 32 independent synthetic spectra for different values of the accumulated charge. For each set of spectra, the uncertainty on the analyzed value is deduced as explained in Fig. 3(a). The dependency of the total uncertainty for the different fit parameters on the accumulated charge is plotted in Fig. 3(b). Observe that the average fit precision, indicated as the brown line in Fig. 3(b), is ten times smaller than the statistical uncertainties for the studied cases. The plot of the total analysis uncertainty as a function of the accumulated charge may be instructive for the experimenter to estimate the accumulated charge that is needed to reach a specific uncertainty. It is reassuring that the statistical variations in the case of non-overlapping signals derived with the Monte Carlo method closely match with the expected uncertainties from Poisson statistics.

Finally, we remark in Fig. 3(b) that the detector gain is determined with a good accuracy from the analysis of an experiment with even a modest accumulated charge. It is expected that one will take advantage of this fact to reliably calibrate and sum the various spectra of a multi-detector setup.⁵³ 

From Fig. 2, it is seen that it takes the algorithm between 2000 and 7000 evaluations to find a fit with $1%$ standard uncertainty. To give the reader a rough estimate, on a contemporary PC (i5 6267U), the computation time for a single spectrum simulation is about 55 ms. This leads to a refinement time between 2 and 7 min.

However, the time needed to generate a simulation can be substantially reduced to about 15 ms if the average isotope mass is used instead of evaluating the isotopes individually. Therefore, we implemented the fit routine in a way that the individual isotope contributions are neglected in the early stage of optimization and are automatically activated in the later stage. With this, the computation time is reduced, without loss of accuracy, to a period of 1–3 min.

The performance of Ruthelde to fit Rutherford backscattering spectra is comparable to the one reported by Silva et al. for Multi-Simnra using adaptive differential evolution with crossover rate repair.²³ A quantitative comparison cannot be made though since different target models are used for benchmarking and since the actual execution time depends on the implementation of the forward simulator.

The Ruthelde software can be used either via its graphical user interface or via the command prompt. The latter enables one to initiate multiple instances of the program in separate threads. Therefore, it provides a powerful tool for the unattended fitting of Rutherford backscattering spectra, as illustrated in Secs. IX.

We investigate the performance of Ruthelde to analyze a large set of experimental Rutherford backscattering spectra obtained on strontium-titanate thin films. Strontium-titanate is an important high-$κ$ material for the semiconductor technology,⁵⁴ but it is known that the electrical and optical properties depend critically on the stoichiometry and the thickness.^55,56

The studied system consists of an Sr$x$Ti$y$O$z$ layer, approximately 75 nm thick, grown with molecular beam epitaxy (MBE) in a Riber MBE49OX on a 200 mm diameter (001)-oriented silicon substrate. The strontium is evaporated from a double filament effusion cell, and the titanium is evaporated from an electron beam evaporator. An excellent stability of the titanium beam during the deposition is obtained using a mass spectrometer signal as a feedback loop for the electron beam evaporator power. The silicon substrate is kept fixed during the deposition (not rotating), which is expected to intentionally lead to a composition gradient across the wafer both for strontium and for titanium. Molecular oxygen is fed into the deposition system to a pressure of 0.13 mPa.⁵⁷ The introduction of oxygen into the vacuum chamber, besides allowing the Sr$x$Ti$y$O$z$ compound to form, is also known to change the growth rate for titanium and strontium. Since the Ti flux is kept constant by the feedback loop, the oxygen will affect the strontium growth rate, and as a result, it is expected that the Sr/Ti ratio is not constant as a function of depth. The Rutherford backscattering spectra are acquired with the experimental conditions as described in Sec. IV.

As opposed to the analysis of the synthetically generated spectrum, here, the Sr$x$Ti$y$O$z$ layer is modeled as if it is composed of six Sr$x$Ti$y$O$z$ layers, each layer with an independently variable concentration. Therefore, it allows one to test the ability to fit a continuous profile.⁵⁸ In total, there are 17 degrees of freedom. The quality of the fit for the measured Rutherford backscattering spectra is evident from Figs. 4(a)–4(c). Note also the sloped Sr signal as compared to the one in Fig. 2(c). This is properly accounted for by the fit (red lines) by modeling the concentration gradient as a function of depth.

A similar analysis is done for the Rutherford backscattering spectra obtained at other locations along the direction of the gradient on the wafer. Notably, the same input file is used for the differential evolution algorithm to fit the various spectra. Also, no initial guess is provided to the algorithm, and the spectra are analyzed without human intervention.

The data points in Fig. 4(d) represent the areal density for the elemental constituents as a function of the position on the wafer. One may observe that the titanium and the strontium areal densities have an opposite trend as a function of the position on the wafer. This is as expected from the deposition conditions. One may also observe that the oxygen areal density decreases with a decreasing relative concentration of titanium. The results demonstrate that it becomes possible to study the relation between the composition of an Sr$x$Ti$y$O$z$ layer and its opto-electronic properties. More generally, the results convincingly demonstrate that the differential evolution algorithm is capable of refining the areal densities and elemental depth profiles as a function of the position on the wafer from the respective Rutherford backscattering spectra.

The two-sigma fit errors plotted in Fig. 4(d) are obtained by repeating the analyses 30 times. While the fit errors for titanium and strontium are smaller than the data points, the fit errors for oxygen are clearly visible. As an example, for the spectrum near the center of the wafer, the two-sigma fit error on the areal density is 0.8% for titanium and strontium and 3.8% for oxygen. Interestingly, the fit error obtained on the experimental spectra exceeds the fit error obtained on the synthetic spectra. This may be attributed to (a) an imperfect assumption of the sample model or experimental conditions or (b) additional noise in the experimental spectrum due to limited differential non-linearity of the detector electronics. Figuring out which effect contributes to what extent to the increased fit error is subject to ongoing investigations. Furthermore, the fit error can be compared with the counting statistical two-sigma uncertainty on the areal density, which is, for the given example, 0.9% for titanium, 0.5% for strontium, and 1.6% for oxygen. Thus, in the present case, the fit error is comparable to or greater than the statistical uncertainty. In general, one must determine the fit error and the statistical uncertainty independently and report a total analysis uncertainty that covers both contributions by adding the variances.

Rutherford backscattering spectrometry is useful to analyze the presence of heavy elements near the surface of a lighter substrate. The difficulty lies in modeling the measured data since there is a high degree of correlation between the various parameters.

We demonstrate that differential evolution is a suitable algorithm to aid with the analysis of Rutherford backscattering spectra. Our work shows that the method is very robust and that very accurate and very precise results can be obtained. Besides, we present a method to estimate the fit precision on the one hand and the uncertainty that derives from the counting statistics on the other hand. The reported total uncertainty should account for both sources of uncertainty.

The capabilities of differential evolution to aid with the analysis of Rutherford backscattering spectra are demonstrated on different model systems. Besides, it is also demonstrated on a practical use case: a large collection of experimental spectra. With this use case, it is also shown that the method can refine a depth-dependent concentration gradient.

This project has received funding from the European Union’s Horizon 2020 research and innovation program under Grant Agreement No. 824096. The support from Johan Desmet, Praveen Dara, and Michiel Jordens at IMEC is greatly acknowledged. The work by C.M. has been supported by the ERC NOTICE GA 864483. The work by N.C. has been supported by the Research Foundation—Flanders (FWO) (Grant No. 1S45421N).

The authors have no conflicts to disclose.

René Heller: Conceptualization (lead); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (lead); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Nico Klingner: Data curation (equal); Formal analysis (equal); Investigation (equal); Resources (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Niels Claessens: Formal analysis (equal); Investigation (equal); Resources (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Clement Merckling: Data curation (equal); Formal analysis (equal); Investigation (equal); Resources (equal); Writing – review & editing (equal). Johan Meersschaut: Data curation (equal); Formal analysis (equal); Investigation (equal); Resources (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The Ruthelde code is implemented in the Java programming language, and thus, its usage is operating system independent. Besides the optimization algorithm presented here, the software is also equipped with a responsive user interface that includes features, such as stopping calculation, penetration depth plotting, and scattering kinematics calculation.

The source code is available as an open source under the GNU General Public License Version 2 on GitHub.⁵⁹ A pre-compiled executable is included in the repository as well. The synthetic RBS spectra used for this study are available as an open source under the Creative Commons Attribution 4.0 International license on RODARE.⁶⁰