Nuclear Inst. and Methods in Physics Research, B

Semi-empirical formulae like the Yamamura model provide a quick reference of sputter yields for applications such as sputter depth profiling and secondary ion mass spectrometry. Fit parameters in such models are prone to errors which can propagate into the prediction of sputter yields. We compare experimental sputter yields to predictions of the Yamamura model using a Bayesian Markov Chain Monte Carlo (MCMC) algorithm. The model parameters Q (linear scaling) and s (power-law scaling) are explored. The results from MCMC are then compared to propositions of Seah (Surface and Interface Analysis, 2005) and extended to a collection of target materials by fitting simulated yields for argon and neon using TRIDYN. Q was found to be proportional to the threshold energy and a simple relation is proposed. The simplicity notwithstanding, Q is speculated to be dependent on a multitude of parameters such as density, energy transfer and orbital filling.


Introduction
Ballistic processes involved in energetic ion-target surface interactions find applications in thin film physical-vapour deposition [1], space [2] and fusion [3] research. They are also widely exploited in metrology applications such as for sputter depth profiling [4], lowenergy ion scattering [5,6], and secondary ion mass spectroscopy. The interactions in each case lead to a removal of target material (sputtering). While this may be beneficial to some applications and detrimental to others, in either case, the quantification of sputtering efficiency is strongly desired. Errors in predictions could lead to a range of issues, starting with minor ones, such as deposition of thinner (or thicker) films which are relatively simple to correct. However, improper sputter efficiencies can also lead to major issues such as errors in quantification of material and sample depth in metrology applications, or to device failure during operation of materials exposed to plasmas [7].
The quantification of sputter yields has improved with growth of understanding of the nature of processes involved in ion-target interactions. Models accounting for energy losses from an ion into the metal via nuclear [8,9] and electronic [10,11] channels were developed along with inter-atomic potentials [9]. Sigmund first formulated a semiempirical predictive equation for sputter yields [12]. Improvements to the model led to greater precision in predictability and extension of usage to low ion energies near the sputter threshold. Modern technology facilitated the use of computational methods like Monte to be effectively constant which is unphysical. The data for neon in Seah's work also showed random scatter and could not be described in a consistent manner as for argon.
In this work, we look at the Yamamura model in an effort to better understand the physicality of the parameter Q and its relation to other ion-target parameters and find a generality for incident ion species. We first study correlations between fit parameters by comparing the model to experimentally determined sputter yields using a Bayesian parameter space searching algorithm. Subsequently, we compare the values for the parameter Q (linear scaling factor) obtained from this method, with those proposed by Seah [19] for argon. We show that upon observing correlations between Q and various ion-target relevant parameters, we find consistent structure in Q as a function of the energy transfer parameter for datasets of both argon and neon. Furthermore, using consistent datasets obtained from Monte Carlo simulations, we show a linear dependence between Q and E th for the two ion species.

Yamamura Model
Following the initial proposition by Sigmund [12] for predicting the sputter yield of a target irradiated by ions with energy >1 keV, many developments to his propositions have evolved [15,16,21,22] for ion energies near the sputter threshold. One of the most extensively used semi-empirical models which we shall discuss in this section is the one developed by Yamamura [15]. It builds upon Sigmund's model and additionally accounts for electronic energy loss as well as a sputter threshold, and is defined as: where Y is the sputter yield, which varies as a function of the incident ion energy, E (in eV). M 1 and M 2 are the masses of the ion and target atoms (in amu), respectively. is the reduced energy. * is an energy independent function of the ratio of ion and target masses. S is the nuclear stopping power of the target (in eV Å 2 atom −1 ), E ℎ is the sputter threshold. U is the surface binding energy, approximated by the sublimation energy of the target material. The factor 0.042 has units of Å −2 . is a parameter that factors in the contribution of reflected ions to the recoil cascade and takes the form: where W is taken from tabulated values (given for 32 elements) or is otherwise approximated as W = 0.35× S . k is the correction factor for the electronic stopping power, set to 1. A detailed treatise on the definitions and functional forms of parameters * , S (E) and can be found in [15]. Q(Z 2 ), s and E ℎ are parameters that take values given by Yamamura in the form of lookup tables, with E ℎ being calculated as: and, is the energy transfer factor for elastic collisions. Unlike ℎ , the parameters Q and s do not have an empirical expression and their values need to be obtained for each target element. The parameter s is tabulated in the paper by Yamamura [15] and varies for each target element with values of either 2.5 or 2.8. Q is considered to be dependent on the target only and is a function of its atomic number. However, no functional form was described in the work. Instead, Q values were tabulated for 34 target materials by fitting experimental (and simulated) sputter yield data and using Q as a fit parameter. Where data was unavailable (42 target elements), Q was set to 1. It was shown by Seah [19], that this approximation of Q leads to errors in prediction of sputter yields by at least 20% for targets where Q was not tabulated (set to unity). These errors can be significant for applications such as sputter depth profiling, low energy ion scattering or secondary ion mass spectroscopy. In the following section, we shall look at the corrections proposed in [19] for the parameter Q for argon ions.

Corrections to Q: Q eff
In order to reduce the errors arising from unknown values of Q, the dependence of Q on the target material was explored by Seah [19,20]. There, Q was considered to have a dependence (and dimensions) of atomic density. To account for this dependence on the target density, an average inter-atomic spacing, r, was defined as: where is Avogadro's number and is the density of the target (kg/m 3 ). Q, now termed Q eff , is then allowed to scale with the mass and density of the target as: where a, b, c, d, e, and f are fitting coefficients -with values: 0.0202, 19.0, 14.6, 0.0166, 5, and 50, respectively for argon ion bombardment -that include any dependence of Q on the ion specie. Albeit artificially (due to the coefficients), Q eff now scales only with the target parameters (mass and density). Eq. (6) reduces the number of tabulated values from 77 in total (one for each target) to 6 per ion specie, while decreasing the relative error between predictions and experimental data to 9% on average. The Gaussian terms are used to describe the fluctuations in Q as a function of target mass and provide a good approximation for elements with M 2 ≤ 100 for argon ion bombardment. In [19,20], the coefficients in Eq. (6) were tabulated for argon for the Yamamura model. However, similar treatment was not performed for neon due to possible correlations of Q eff with s and large Q eff r 3 oscillations at low target masses. We shall, in the following section, describe a probabilistic method for searching the parameter space while we simultaneously observe correlations between model parameters and quantify uncertainty in the parameters of the Yamamura model.

TRIDYN simulations
A comparison of the models to experiments is not always possible to due the lack of experimental sputter yields. To circumvent this deficiency, we rely on TRIDYN simulations of sputter yields for a wide range of elements that are solid at standard temperature and pressure, as a means to understand the variation in an ideal case.
Simulation codes such as TRIDYN are based on the binary collision approximation and are subject to variations in output depending on the material inputs such as the surface binding energy as well as incident ion energies. This can lead to discrepancies between simulation output and experiments, especially for reactive ion species which form compounds with the target elements or exhibit strong diffusion effects. Despite these drawbacks, we have shown TRIDYN to be effective for reactive ion species interacting at energies near the sputter threshold [23,24]. Nevertheless, to minimize errors from dynamic compositional changes and its effects on surface binding of atoms, especially where experimental data is lacking for comparison, we limit ourselves in this work to inert gas ions incident normal to the target surface.
Simulations were performed in the dynamic mode for argon and neon ions bombarding elements with atomic number (Z) between 3 ≤ ≤ 98. S for the target was set to the sublimation energy in each case and the simulations were performed until a fluence of 1×10 17 ions/cm 2 . A maximum of 5% of gas species were allowed to be implanted and retained within the target.

Bayesian MCMC
In order to assess the best-fit parameters, their correlations with each other and their uncertainties, Bayesian analysis is a powerful tool. The power of Bayesian analysis to map the correlations for a large set of fit parameters has been demonstrated for a wide variety of fields in physics [25] and astronomy [26]. Bayesian analysis relies on probabilistic definitions based on Bayes' rule [27] which is given as: where ( | , ) is the posterior probability density function (PDF) which describes the probability for a parameter set = (where each represents a sputter yield datum for a given ion-target combination) for a model ( ) [25,28]. Here, are the number of data points and the number of model parameters. ( | ) (read as given M) encodes the prior knowledge of the model parameters the probability for the parameters or our confidence in their values, and is called the prior PDF. ( | , ) is the likelihood function () which is defined later. The denominator, ( | ) is often termed as the Bayesian evidence, which refers to the validity of the model in predicting the data, and for practical purposes serves as a normalizing factor in Eq. (7).
Bayesian Markov Chain Monte Carlo (MCMC) allows for a means of estimating the posterior probability defined by Bayes' theorem by randomly sampling the PDFs for the posterior. We assume no a-priori knowledge of the mean and distribution of the model parameters and for simplicity, use a uniform PDF as priors for each of the model parameters. When uniform priors are used, the posterior PDF is proportional to the likelihood function as: The random sampling is performed by using an affine invariant ensemble sampler, emcee [29]. The affine invariant sampler is advantageous over other common sampling methods [30,31] as it is unperturbed by affine transformations (rotations, stretching) of the parameter space. It implies that the sampler is not sensitive to the scales of the parameters involved. The sampling of the posterior PDF occurs through assignment of a set (ensemble) of walkers. The walkers are kdimensional vectors which randomly move across the parameter space. Their movement is governed by the likelihood function, , which for emcee is defined as a log-likelihood: where Y is the experimental yield with uncertainty e , and the subscript i denotes a datum for n data points in a dataset D. Y yama is the yield predicted by the model (Eq. (1)). The walkers move in a two-step process, first by selecting a point to jump to, and then deciding whether the jump occurs. The second step helps to map out the contours in the parameter space and the walkers cluster around locations of highest probability. The walker movement process is repeated thousands of times for each walker until the entire posterior PDF is mapped and the likelihood is maximized. As a result, the posterior PDF of each model parameter is populated by the walkers and correlations between parameters can be evaluated from the PDF maps. As a consequence of this analysis, the expectations values, maximum likelihoods for parameters and the confidence in their value can be readily obtained.
In order for defining the problem for such Bayesian analyses, a prior PDF is important as it encodes the knowledge of the system known a-priori. This is a subject of importance and ongoing research for a generalized choice of priors that is ''objective'' [32]. Priors such as Uniform priors are commonly used, and Jeffery's prior [33] is a general purpose prior but more complex in definition. Here, we use uniform priors for simplicity. More complex priors may be considered for future works.
For purposes of our analysis, Eq. (1) was used as the model (M). The parameter space was the set of free-parameters: Q, s, E th . All MCMC calculations were performed using an ensemble of 50 walkers for a step count of 25000 steps. The first 5000 steps were discarded (burn-in) to negate any bias introduced by the choice of starting positions of the ensemble. The priors were chosen to be uniformly distributed with the starting value equal to the tabulated values (for Q and s) [15] or as calculated using Eq. (3) for the sputter threshold (see Table 1). The range for Q was allowed to vary over an order of magnitude (on each side) from the tabulated value while the ℎ was varied from 1 eV to 60 eV (the range of variation of thresholds for most metals). For s, the range was varied from 0 to 4. The choice of the lower bound is the extreme case where there is no threshold effect and the model is reduced to the initial formulation by Sigmund [12]. For the upper bound, literature evidence [15,16] suggests optimal values of 2.5 or 2.8 depending on target material. It is possible that s can vary beyond the tabulations and a higher limit was chosen for exploring probable values.
Following every MCMC calculation, we present the expectation values of the marginalized PDFs for the parameters which we shall refer to as 'best-fit' estimates. These values are then used in the Yamamura model to depict the theoretical sputter yield curve or for further analysis.

Experimental method
Sputter yields were obtained experimentally for transition metal films of molybdenum (Mo), ruthenium (Ru), palladium (Pd) and tungsten (W) by measuring mass loss from a quartz crystal micro-balance (QCM). Elemental films of 400 nm thickness were deposited on QCM crystal substrates using a DC magnetron sputter deposition technique. The deposition facility was baked out to a base pressure of 1 × 10 −8 mbar with the dominant background specie being water vapour. The metal targets were pre-sputtered for 20 min to remove oxides and contaminants from target surface. Sputter deposition was carried out under an argon feed gas at a gas flow of 30 sccm at a constant power of 400-500 W. The thickness during sputter deposition was monitored using the frequency response of the QCMs.
For the sputter yield experiments, argon and neon ions were generated using a 15 cm DC Kaufman ion source (Veeco Instruments). The ion flux was monitored using a Faraday cup, simultaneously along with the mass loss. The ion energies were pre-characterized by the Faraday cup (FC) which also served as a retarding field energy analyser. Ion energies for obtaining sputter yield data were varied in random order between 50 eV and 300 eV to avoid systematic errors. QCM frequency responses were converted to thickness estimates using the Zmatch method [34]. Sputter yields were calculated using the thickness difference before and after each exposure, and the irradiated fluence measured simultaneously using the FC. Details of the experimental setup and data processing can be found in Refs. [34,35].

Experimental sputter yields and model parameters: Ar +
Using Ar + as the bombarding ion specie, sputter yields were experimentally determined by measuring changes to the QCM frequency response. The results are shown in Fig. 1. The data was compared to the P. Phadke et al. predictions of the Yamamura model, replacing Q with Q eff values from Eqs. (5) and (6) using fit parameters (a-f) from [19]. In addition, the experimental data was used to determine posterior probabilities for the parameters: Q, s and E ℎ from the Yamamura model which then can be compared with values from Eqs. (5) and (6). Availability and ion energy range of literature data on yields differ between elements Ru, Mo, Pd and W. In order to prevent varying credibility in the best-fit values from limited reference datasets, reference measurements were not used as inputs into the MCMC algorithm. We shall see in a later section, the comparison of the obtained posterior probabilities in comparison to literature datasets.
Experimentally, the sputter yields for molybdenum, ruthenium and tungsten are lower than reports from literature values obtained from experiments from bulk samples [36,37] or e-beam deposited films [34]. This is possibly due to our sputter deposited films being denser (due to energetics of the deposition process) and more flat on a nanometre scale roughness, compared to materials previously studied in literature. However, this hypothesis is inconsistent as can be observed in the case of palladium films ( Fig. 1(d) and as we shall see, in the case of Ne + ion bombardment. The dependence of sputter yields on structure, density and roughness and the interplay between them, therefore, appears to be a complex function of these parameters and requires dedicated experiments. Differences arising from these parameters would exist for every dataset (measured in this work and in literature) and precise understanding of this interplay is beyond the scope of this paper. We shall instead focus on a global behaviour of sputter yields. For this, we use the ensemble sampler described in Section 3 for our dataset as well as a collection of available literature data (see Section 5.3).
An example of the posterior probabilities generated from the MCMC algorithm for palladium is shown in Fig. 2. Similar plots which were generated for the other metals studied are shown in Supplementary Information. The MCMC maps depict a negative correlation for E th with both Q and s. However, Q and s are positively correlated with Q straying to a 12% larger value from the Q eff calculated in [19], and 50% larger than the original tabulation of Q in [15]. The threshold for palladium sputtering by argon is found to be 80% higher than that estimated by Eq. (3). This relatively high value for E th arises partly due to the lack of large sets of data below ∼50 eV.
The power of the Bayesian method lies in the updating of posterior probabilities via better likelihood and prior estimates which are facilitated by availability of datasets. We shall compare the 'best fits' from our dataset with a compendium of literature values in Section 5.3. Updated Yamamura curves using the expectation values ('best fits') from the MCMC sampling are shown in Fig. 1 along with a shaded region highlighting the 1 variations in the yields from the parameter estimates. As an aside, an area within the 50th quantile region (contours from Fig. 2) could be used but in most cases the 1 regions are larger and represent 'worse case' estimates. Table 1 encapsulates the priors used and posteriors obtained in each of the datasets. It is evident that in each case, the threshold value depends on the availability of data near the true (yet unknown) threshold.
Values of s in certain cases of the experimental data may be beyond the set upper bound of 4 [See Supplementary Information] and likely stem from the gradual rise of the sputter yield beyond the threshold. The Q values obtained from the MCMC likelihood estimates, deviate from the Q eff proposed in [19] by at most a factor 2 (molybdenum) and at least 30% (tungsten).

Experimental sputter yields and model parameters: Ne +
Similar experiments using the QCM were carried out for the target elements under neon bombardment. The sputter yields obtained in this work along with references are plotted in Fig. 3. Here, the structural P. Phadke et al.

Fig. 2.
Posterior probabilities for parameters from Eq. (1) for sputter yield data from palladium under Ar + bombardment. The red circle depicts Q eff and s values from Seah [19]. The black cross-hairs show the 'best fit' locations and the contours depict 16th, 50th, and 84th quantiles. The histograms show the marginalized probabilities for each parameter.

Table 1
Parameters used to search the posterior probabilities for the model parameters in Eq. (1). Sputter yields from Ar + bombardment used as input data. The priors are described in the form: U(low, up), denoting uniform probabilities over a range with a lower (low) bound and an upper (up) bound. The expectation values indicating the 'best fit' and 1 uncertainties of the expectations obtained from the Bayesian sampling is compared to the tabulations of the Yamamura model.

Priors
Expectation dependence of the sputter yields is not apparent as the yields are consistent with reports in literature, taken from [36,37] for (bulk) samples immersed in a plasma.
In previously reported studies on the application of the Yamamura model for Ne + ion bombardment, an analysis of the dependence of Q eff on target mass was not performed [20]. The coefficients (a-f in Eq. (6)) being unknown, we rely on the original Q proposed by Yamamura for comparison. The Yamamura model in most cases provides suitable prediction for all elements down to 100 eV. Below 100 eV, the predictions begin to deviate from experiments.
The MCMC sampling was performed much in the same way as for argon sputter yields. A uniform probability density between certain bounds was chosen for the parameters as priors. The results of the ensemble sampling for palladium sputtered by Ne + are shown in Fig. 4. The positive correlation of s with both Q and E th is significantly reduced for this dataset. E th and s, which were previously negatively correlated (for Ar + ) still behave similarly, but the correlation is weaker. The histograms show the marginalized PDFs of the parameters which are Gaussian, but in the case of s, these are asymmetrically skewed. The example here illustrates the difference in values of the maximum likelihood (highest point on the histogram) from the expectation (mean) estimates.
The initial inputs and the expectation values obtained are summarized in Table 2 along with propositions by Yamamura. We notice that the sputter threshold from the MCMC expectation value is higher than as predicted by Eq. (3). A higher experimental threshold is accounted for by the MCMC parameter search. The gradual slope for Ruthenium (s ∼ 4) implies that the threshold is beyond the probed dataset which is also reflected in the relatively lower threshold prediction. We also find from Fig. 4 that the energy scaling given by the exponent s, varies between elements and is correlated to Q (positively) and E th (negatively). The updated Yamamura curves using the expectation values ('best fits') from the MCMC sampling are shown in Fig. 3 along with a shaded region highlighting the 1 variations in the yields from the parameter estimates.
We find that the expectation values for Q represent the experimental data well but in some cases deviate from the predictions of Q eff in Eq. (6) and Yamamura's tabulations as is seen in Figs. 2 and  4, respectively. The elements studied here are heavy (M 2 ∼ 100 or P. Phadke et al. Reference yields are also plotted (common legend out of plots). The Yamamura model using 'best-fit' parameters from the MCMC outputs (solid lines) are compared to Yamamura model (dashed lines). The shaded areas represent the 1 deviation in 'best-fit' parameters' influence on the yield.

Table 2
Parameters used to search the posterior probabilities for the model parameters in Eq. (1). Sputter yields from Ne + bombardment used as input data. The priors are described in the form: U(low, up), denoting uniform probabilities over a range with a lower (low) bound and an upper (up) bound. The expectation values indicating the 'best fit' and 1 uncertainties of the expectations obtained from the Bayesian sampling is compared to the tabulations of the Yamamura model.

Priors
Expectation greater). Eq. (6) predicts for heavy elements, no change in Q eff other than that induced by the inter-atomic spacing. In the following section we shall compare Q obtained from MCMC ensemble sampling to the predictions of Eq. (6) using reference data for a multitude of elements where experimental data is available.

Q MCMC and Q eff
In Sections 5.1 and 5.2, we found that the values of Q and s are correlated using the MCMC ensemble sampling method. The values of s were found to vary from one element to the other, with no distinct pattern observable in relation to target specific parameters. The linear scaling factor Q (Q MCMC in our study), requires some scrutiny as its effect on the sputter yield is relatively clear. Here, we begin by comparing the obtained Q MCMC with Q eff from [19,20]. Eq. (6) predicts that the term Q eff r 3 for argon bombardment saturates to a constant value for target masses > 100 amu. This plateau, while apparent in the evaluation in [19,20], would be unlikely as elemental properties are periodic and any initial fluctuation in Q eff r 3 at low target masses would carry forward to high masses as well. In the evaluation, experimental sputter yield datasets for a variety of target materials were studied, however, there was a lack of populated data for the region between 110 amu and 180 amu. Any fluctuations, if present in this region, would flatten out for the higher target masses.
Here we use experimental data available for a variety of elements obtained from the literature [34,36,[38][39][40][41][42][43]. To supplement the lack of sputter yield data for target masses between 110-180 amu (d and f block elements), we chose to simulate the yields for a selection of lanthanides using TRIDYN [13]. Further, to understand the limitation/validity of this simulation based approach, we randomly selected elements for which literature data was available and used them as inputs for the MCMC ensemble sampler. Fig. 5a depicts the expectation values of Q MCMC for experimental sputter yields from literature data (labelled literature), and for simulated yields for selected targets under Ar + ion bombardment. For comparison, Q values of the Yamamura model [19] are also plotted along with the predictive function for Q eff from Eq. (6). While the scatter in the datasets from MCMC ensemble sampling is quite large, the behaviour of Q MCMC is generally consistent with Eq. (6) for target masses <100 amu. For 2 > 100 amu, the initial oscillation (in the region of 20≤ M 2 ≤ 70 amu), repeats and continues on up to 250 amu. The periodicity is non-uniform and is not captured by Eq. (6), which in this range, approximates the variation to a constant. Q MCMC from Mo-Ar + results from QCM measurements, marked in Fig. 5a, is an outlier.   A significant scatter is seen for data-points between 100 amu and 110 amu, and the cause is yet unknown, however, an overall oscillatory behaviour is observed over the range of target masses.
Data for neon bombardment is more scattered as seen in Fig. 5b. Results not only from our MCMC analysis, but also reports from [20] show elevated scatter in comparison to the argon data. Although a structure for neon datasets similar to argon can be argued, the scatter in the data prevents such a conclusive remark.
In general, considering the oscillatory behaviour of data, the two Gaussian terms in Eq. (6) do not fully encapsulate Q MCMC or Q eff for M 2 > 100 amu. A multitude of Gaussians would be needed under the approach of Eq. (6) which becomes cumbersome with the addition of more terms. We shall discuss in the following section an alternative perspective to the evaluation of the datasets in an effort to encapsulate the behaviour of both ion species.

Dependence of Q on ion and target parameters
The proposition of Eq. (6) was a significant step forward in determining an empirical relation to the scaling factor, Q, used in a multitude of predictive models of sputter yields. Although it was shown [19,20] to be valid for models such as the Matsunami formula [16] (not discussed here), the approach was found to be insufficient for the Yamamura model due to the scatter in the data other than for argon (in addition to correlations between fit parameters). Further, the formulation of Eq. (6) is empirical and contrasts with the original formulation by Sigmund [12], which was also pointed out in [19].
In all previous attempts in qualifying dependence of Q on various parameters, datasets have been heterogeneous. For example, in Yamamura's original formulation [15], sputter yield data was obtained both experimentally as well as using simulation codes. Experimental data in themselves contain scatter due to a multitude of factors: different surface preparation, variability in roughness, crystalline phase and density, poorly characterized ion beams or data obtained at different energy ranges especially for ones far away from the sputter threshold see Fig. 1, and Fig. 3. This results in large spread in available Q and Eth estimates which leads to difficulty in discerning trends in the data.
We attempt to remedy this discrepancy by using data obtained using a consistent, homogeneous method. We source sputter yield data using TRIDYN simulations for 70 common elements (excluding gases, short lived elements and some actinides) in an energy range from the sputter threshold to keV energies. TRIDYN has been shown to reliably simulate sputter yields at low incident energies for a multitude of elements [35]. TRIDYN in this case, represents an ideal behaviour of a target element, restricted by the sublimation energy and assuming binary collision of the ion with a target atom. Data obtained in this manner does not suffer from experimental scatter. Additionally, we can simulate a much larger set of elements than otherwise experimentally available. These combined advantages can assist in finding correlations between parameters in the sputter yield data. We fit the simulated sputter yields using the MCMC algorithm as discussed previously. Using the obtained likelihoods, we examine the dependence of Q on various ion-target dependent parameters. Fig. 6 shows the Pearson correlation coefficients of various parameters with respect to one another. For example, we see that the target mass, M 2 is negatively correlated with the energy transfer parameter ( ) for argon and neon. This is to be expected as the energy transfer decays with increasing target mass after a certain maximum dictated by the ion mass. The actual dependence of on 2 however, is more complex as known from Eq. (4). Similarly, we find that the dependence of Q is stronger (more positive) with respect to E th rather than or . This dependence on th is replicated in the neon data ( Fig. 6(b).
In order to visualize the dependence of Q on three different parameters, we plot the Q values from the argon and neon datasets as a function of in Fig. 7. By varying the size of the datapoints by the corresponding threshold energy (obtained from the MCMC algorithm), P. Phadke et al. we find that the higher the threshold energy, the higher is the value of Q. This dependence is stronger in the argon data than in the neon data as is also seen as the correlation coefficient drops from 0.90 to 0.67. Values of Q obtained from MCMC fits to sputter yields from literature (Q lit MCMC ) are also presented in Fig. 7. Dependence of Q lit MCMC on is qualitatively similar to that from fits to TRIDYN simulated yields (Q TRIDYN MCMC ). The deviations in absolute values of Q MCMC could arise from either experimental scatter (as outlined above), or from physical effects that are unaccounted for in TRIDYN. While these deviations do not allow to reproduce experimental Q, Q TRIDYN MCMC is smoother and covers the full range of which aids in revealing trends in the data. Dependencies of Q on M 2 are similar to those on . The choice between M 2 or is therefore not important for the subsequent analysis.
Sub-dividing the dataset by the location of elements in the periodic table, we find that the peaks in Figs. 7a and 7b are composed of the d-block elements. The p and f-block elements form their own set of peaks while s block elements are mostly confined to valleys. The overall structure of the Q values is similar to the variation observed by Seah (See Supplementary Figure S7) for both Argon and Neon. However, given the correlations, and correspondence of peaks with orbitals of target elements, we find that the description of Q is not limited to the density but rather a combination of multiple parameters. The simplest correlation we find for Q from the pairwise correlation is to the threshold energy. Datasets for the pairwise correlations show a nearly linear dependence of Q on E th , within our best estimate of E th . As a simple formulation exploiting this linear dependence, we find coefficients for a linear fit for Q as a function of the threshold energy, th , calculated using the MCMC method. This is done using a robust linear regression model with a Huber loss function in order to reduce biasing by outliers. The datasets show varying slopes as a function of the electron orbitals as seen in Fig. 8. Elements classified in the s, p and f blocks in the periodic table exhibit Q values that are well approximated by this linear dependence. The s block elements cover a wide range of E th values, with the highest corresponding to Be. Be is also seen as a stark outlier when Q values are plotted against (Fig. 7) in both literature and TRIDYN datasets. However, Q values from the d orbital present significant scatter and the linear dependence is not apparent. Upon further subdividing the d block elements by rows, we find the linear dependence is restored within the subdivisions. Such a subdivision also favours other blocks, but due to fewer number of elements, such a procedure is counter-productive. However, such a stratified behaviour points towards an unaccounted, periodic variable. Understanding such a variable would be needed to gain insights into the physical meaning of Q and its dependencies.
Speculatively, the scatter could also originate from the ensemble clustering around a bound under some cases due to the power law term arising in Eq.  similar linear dependence is also observed in the neon datasets (see Supplementary Information Figure S9) which is in contrast to Seah's approach (see Fig. 5).
As a final step, in an attempt towards approximating Q values without the use of lookup tables, we map the variation of Q to the threshold predicted by the MCMC. In doing so, the number of parameters needed to estimate Q for a given target is significantly reduced. The resulting dependence of Q as a function of the threshold from the MCMC fits is in accordance to: Fits to this simple model to simulated data are summarized for argon and neon in Table 3 [also see Fig. 8 Table 3 Best-fit values approximating Q from the sputter threshold in Eq. (3). The goodnessof-fit (R 2 ) for each ion specie is furnished for the linear approximation as well as the functional form proposed by Seah [19], calculated from Eq. (6) Figure S9 for fit results]. In order to compare the linear approximation with the equation proposed by Seah [19], we also look at the goodness of fit for both functional forms. It is clear from Table 3 that the linear approximation for each orbital type estimates Q just as well, if not better than Eq. (6). Additionally, it reduces the set of parameters from six, down to two for a given orbital. Expanding on the linear approximation, E th depends on and U as per Eq. (3). The Q factor in the Yamamura model (Eq. (1)) can then be approximated as : = ′ ∕ + ′ where ′ and ′ are the slope and intercept relating the sputter threshold (from Table 3), and k follows from Eq. (3) as 6.7 (if 1 < 2 ) or 1 + 5.

and Supplementary information
This approximation is valid for neon as well. Comparing the predictions of Eq. (10) to the overall behaviour of Q, in Fig. 7, shows that such a simplified approach can accurately estimate a wide range of Q values. The validity for both Argon and Neon shows promise in the capacity of such a simplification in estimating Q values. However, the reasoning behind the stratification, as mentioned previously, demands for more investigations into the physical nature of Q and its dependencies.
As is seen from Figs. 6 and 7, the true nature and dependence of Q on target properties is multi-dimensional. Attempts at empirically defining dependence of Q on ion or target parameters have been made in the past. Previous reports described the dependence of Q on density [19], the energy transfer parameter [22] and even a complex dependence on ion-target parameters [44]. We find that Q depends not only on the properties mentioned above, but also on the threshold energy and the filling of electronic orbitals. Bayesian MCMC is shown to be a powerful technique for revealing such dependencies between fit parameters and exploring a simple 3-dimensional parameter space. The flexibility of MCMC also allows for ease of scaling the model to incorporate additional parameters if needed, without increasing complexity.
Using Pearson's pair-wise correlation as we do here to compare dependencies of target parameters not directly included in the model, is limited to observing linear dependencies. Hence, the true nature and dependence of Q on the above parameters could potentially be nonlinear and may require considering additional parameters in the model, similar to the work by García-Rosales [44].
Considering the above, the nature of Q is not merely a scaling of the sputter yields to fit experimental data, but a description of a fundamental interaction of ion species with target materials.

Conclusions and outlook
In this work we have reported the sputter yields of four transition metals relevant to semiconductor, catalysis, photo-lithography and fusion research. We have analysed the description of experimentally obtained yields by the semi-empirical Yamamura model using a Bayesian MCMC parameter search. This allowed us to obtain values of the linear scaling factor Q, considering correlations with the sputter threshold energy, E th and the power-law scaling factor, s.
Further, simulating sputter yields using TRIDYN for ∼70 elemental targets allowed us to populate Q values in a consistent manner without experimental scatter. Using Bayesian MCMC algorithms on simulated data enabled comparison of Q to previous semi-empirical descriptions in literature. It was found that Q has a simple positive correlation with the sputter threshold energy and can be reasonably predicted by a linear relationship, provided that the data is separated according to orbital filling, for example, s, p, d and f block elements. This dependence can be exploited to estimate an unknown Q in a simple manner. However, the nature of Q is deemed to be complex, depending on the density, sputter threshold energy, energy transfer parameters as well as more fundamentally, on the filling of orbitals. This work points towards Q not being a mere scaling parameter, but rather a descriptor of fundamental processes of ion-target interactions.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.