Analysis and prediction of sputtering yield using combined hierarchical clustering analysis and artificial neural network algorithms

1.   Introduction

Sputtering is a process whereby target atoms are removed by bombardment with energetic ions. The sputtering behavior of components in various devices, such as electrical thruster engines and fusion plasma chambers, is of significant importance for modeling their long-term performance [1, 2]. Also, plasma processing has been used in semiconductor manufacturing and surface processing [3–5]. The plasma–surface interactions in plasma processing typically involve ion bombardment and surface chemical reactions induced by incident ions, free radicals and photons. The sputtering yield, which represents the average number of atoms removed from the surface per incident ion, is a critical metric for assessing the efficiency of plasma etching processes [6]. Extensive research has shown that the sputtering yield is related to the energy and angle of the incident ions, as well as the surface material and types of ions [7].

The sputtering yield is a critical factor that affects the description of the sputtering process, enhances processing precision and predicts the lifespan of thrusters. Modeling of plasma sputtering phenomena has been developed using various methods, including neural networks and molecular dynamics (MD). Gergs et al [8] used a neural network algorithm to establish an alternative model for plasma–surface interactions. They limited the complexity of the interaction to a minimum by employing an interface model based on an analytical expression. In addition, Gergs et al predicted the time-scale of plasma deposition and sputtering using physics-separating artificial neural networks. This model successfully predicted the phenomenon of sputtering deposition in an Ar/N₂ discharge [9]. Preuss et al proposed a prediction method based on the Bayesian framework, which is mainly used to determine unknown parameters such as surface binding energy (SBE) [10]. This method effectively reduces model errors and parameter orders. MD simulation is considered to be a reliable theoretical tool for analyzing typical plasma processes. Tinacba et al [11] employed MD simulations with a novel S potential model to demonstrate that the etching yields of SF⁵⁺ ions on Si and SiO₂ align well with experimental results, indicating sputtering enhancement by F atom reactions. Cagomoc et al [12] used MD simulations to explore CF_x ion beam etching of SiO₂ with carbon masks, creating 4 nm diameter holes. High ion energy and prolonged etching resulted in tapered holes due to physical sputtering, with redeposited Si-containing species forming sidewalls. Oxygen preferentially removed from SiO₂ surfaces led to Si-rich sidewalls. Simultaneous irradiation with CF_x radicals doubled the flat SiO₂ surface etching yield, but excessive radical flux caused etching to stop. Kawase and Hamaguchi [13] investigated the angular dependence of sputtering yields of SiO₂ substrates under CF₃ beam injection using MD simulations. The obtained data aligned reasonably well with experimental findings. Atomic compositions within the SiO₂–CF mixing layer, as well as kinetic energies and compositions of sputtered species, showed significant dependence on injection angle. Kino et al utilized the hierarchical clustering analysis (HCA) algorithm to investigate the impact of various descriptors on the regression of sputtering yield [14]. However, a critical issue arises from the lack of analysis of physical mechanisms within the HCA algorithm as well as the absence of model validation for the regression process. As mentioned in the literature, most studies have focused on parameter regression problems and the sputtering yield analysis of specific materials. A key issue is that these studies lack correlation analysis between physical variables and high-precision regression models.

The main content of the work is to simplify the physical properties of targets and ions related to sputtering phenomena to obtain a more efficient regression model. Additionally, we discuss the physical significance of the clustering process and the use of neural network algorithms in the regression process to achieve the best model for predicting sputtering yield. The methodologies and experimental setup employed in these analyses are detailed in section 2, with results and discussion presented in section 3, followed by a summary in section 4.

2.   Methods and experimental setup

2.1   Sputtering yield

Figure 1 illustrates the typical relationship between ion incident energy E and sputtering yield (SY). It can be seen that with an increase in ion incident energy the sputtering yield first increases. When the energy reaches a specific value, the sputtering yield reaches its maximum value, that is E_max and SY_max (red line in figure 1). After surpassing E_max, the increase in sputtering yield becomes slower or even decreases [15]. Simplifying the model process by considering only E_th, E_max and SY_max is reasonable, especially in the sputtering process of low-temperature plasma. Based on several works [8, 14, 16], 17 attributes related to sputtering are summarized in table 1. The training dataset used for the research also comes from these literature sources.

Figure  1.  Sputtering yields compared with energy.

DownLoad: Full-Size Img PowerPoint

Table  1.  The features of erosion.

Target	Ion
Atomic mass	Atomic mass
Atomic number	Atomic number
Material density	-
Bulk modulus	-
Melting point	Melting point
Boiling point	Boiling point
Heat of evaporation	Heat of evaporation
Heat of formation	-
Atomic radius	Atomic radius
van der Waals radius	van der Waals radius

 | Show Table

DownLoad: CSV

2.2   HCA algorithm

The HCA algorithm is an unsupervised learning approach that does not necessitate prior knowledge for data analysis and processing. It relies on similarity measurement and achieves target clustering by seeking homogeneity among the targets. The HCA algorithm is currently extensively utilized in various domains such as feature selection [17–19]. A standard HCA pipeline comprises the distinct steps of dimensionality reduction, correlation calculation and merging of similar clusters. Datasets usually encompass many observations distinguished by numerous features. Each feature can be regarded as an additional dimension to the data, rendering it challenging to comprehensively visualize all data in a two-dimensional plot.

In our study, the Pearson correlation coefficient is utilized to compute the correlation matrix. The Pearson correlation coefficient [20] can be used to measure the linear correlation between data, with output values ranging from −1 to 1. A value closer to 1 indicates a stronger positive correlation, while a value closer to −1 indicates a stronger negative correlation; a value of 0 indicates no correlation. The Pearson correlation coefficient, which utilizes the Euclidean distance to represent similarity, effectively eliminates differences in scale between variables during the computation process

where p and q represent distinct physical attributes, ${r_{pq}}$ is Pearson’s correlation function, ${\sigma _{pq}}$ is the covariance matrix between p and q and ${\sigma _p}$ is the standard deviation of p. The distance between physical attributes can be defined as

where P is the number in the same group and i is the value of physical attributes.

2.3   Neural network algorithm

A neural network is an effective tool for nonlinear function fitting, and is particularly suitable for black box fitting with multiple inputs and multiple outputs. Neural networks are widely used in various fitting and regression models. A neural network primarily consists of forward and backward propagation. Forward propagation combines the input data with the weight matrices and produces an output. Backward propagation, typically using a gradient descent algorithm, calculates the output error and propagates it backward to the input layer [21]. The modeling process of an artificial neural network (ANN) involves several key steps: data collection, preprocessing, architecture design, forward propagation, backpropagation, validation, prediction and evaluation, as shown in figure 2. The dataset is divided into training, validation and testing sets in the ratio 0.7, 0.2 and 0.1, respectively. The model learns features and patterns from the training set by adjusting weights and biases to minimize the loss function. The validation set is used to fine-tune hyperparameters (such as learning rate, regularization parameters, network structure, etc.) and to evaluate model performance. The test set is used to assess the final performance of the model, and it is not involved in the training or validation processes. Prior to training, all data undergo standardization. The loss function employed is the mean square error loss function, with regularization utilized to prevent overfitting. The size of the dataset is 267.

Figure  2.  Flowchart of ANN training.

DownLoad: Full-Size Img PowerPoint

2.4   Experimental setup

To validate the developed model, an ion source is employed for data validation. The experimental principle is illustrated in figure 3, and allows for the characterization of sputter erosion rates at different energies and angles. As shown in the diagram in figure 3, the experimental setup consists of a Kaufman ion source (left side of the blue plume) and a sample bracket (right side of the blue plume). The operating conditions for the ion source are listed in table 2. By adjusting the grid voltage of the Kaufman ion source, different energy plumes can be obtained. The sample bracket contains raised supports at different angles for placing metal samples. The samples used in our study are DT4C samples measuring 10 mm × 10 mm × 1 mm. The platform size is 150 mm × 150 mm × 2 mm. The reduction in sample mass is obtained through a high-precision balance and the ion energy and density of the ion source are calibrated using Langmuir probes and retarding potential analyzers [22, 23].

Table  2.  Kaufman ion source parameters.

Parameter	Value
Diameter	16 cm
Flow rate range	5–10 sccm
Beam energy range	100–1000 eV
Gas	Ar
Beam uniformity region	700 mm

 | Show Table

DownLoad: CSV

Figure  3.  Schematic diagram of the setup of the erosion experiment.

DownLoad: Full-Size Img PowerPoint

3.   Results and discussion

3.1   The HCA features

First, the correlation coefficients between various physical quantities are computed using Pearson distance, as illustrated in the figure 4(a). The results indicate that during the initial merging, the enthalpy of vaporization and the heat of vaporization are merged into a similar group. The enthalpy of vaporization and heat of vaporization both reflect the interactions between molecules within a substance and the change in state as it transitions from one phase to another. In the second step, the atomic radius of the substance and the van der Waals force are grouped together. Subsequently, the heat of formation of ions is clustered together with the heat of vaporization, almost simultaneously, while the melting point and boiling point of the target substance are also almost simultaneously grouped together. Indeed, the atomic radius of a substance can influence its van der Waals force through both the distance effect and the polarization effect. In the subsequent clustering, the melting point and boiling point of the target material are regarded as belonging to the same category. Subsequently, the heat of ion evaporation and heat of formation are grouped together, consistent with the clustering of the physical properties of the target material. Similarly, the ion radius and van der Waals radius of the plasma source are grouped together in the following clustering. The melting point and boiling point of ions are grouped together in the next clustering. The properties of ions and targets are clustered using a similar strategy, implying that, according to the HCA algorithm, these properties can be considered the most relevant variables. The atomic number of ions, the mass of ions and the atomic radius and van der Waals radius of ions are grouped together in the next clustering. The mass of the target and the bulk modulus are grouped together as the algorithm continues its calculations. All properties are assigned to at least one category. Next, combinations between categories are formed using the HCA algorithm.

Figure  4.  (a) Heatmap of the first step clustering in the HCA (vdW: van der Waals). (b) Final results of HCA clustering.

DownLoad: Full-Size Img PowerPoint

As revealed in figure 4(b), in the first layer of clustering the atomic mass, atomic number, atomic radius and van der Waals force of ions are grouped together (highlighted by the yellow lines). These properties are closely related to the energy that ions acquire during their motion. Assuming that ions carry one charge each, the mass of ions will affect the speed at which the plasma acquires energy. Changes in properties such as atomic radius and van der Waals force also affect the collisions between ions, leading to variations in ion energy and other related factors. For the green grouping, these include the ion’s heat of evaporation, heat of formation, melting point and boiling point. These four properties share a similar interpretation, qualitatively describing the state transitions of ions. This behavior may be related to partial chemical reactions occurring during collisions between ions and targets. The remaining three groupings are all properties of the target, and are similar to the clustering of ions. The properties connected by red lines are related to collisions involving the target, like the yellow grouping for ions. The remaining three groupings are all properties of the target, and are similar to the clustering of ions. The properties connected by red lines are related to collisions involving the target, like the yellow grouping for ions. Unlike the clustering for ions, the melting point and boiling point of the target are grouped separately. Subsequently, in the next level of clustering, the properties connected by purple and brown lines are grouped together. The brown lines can be considered to have the same clustering basis as the green lines.

Next, we will discuss the various collision processes that occur during sputtering collisions between ions and targets, with a focus on the relevance of atomic number and mass. This will provide a clearer understanding of the mechanisms at play. Various reaction processes have been reported for different masses at the same energy, including inelastic collisions, spikes, surface corrections and cascades [24]. If the ion energy exceeds about 50 keV and the incident ion is heavier than the target atoms, a block of atoms along the initial path of the projectile can be displaced, leading to what is known as the spike regime. In this regime, the transport theory breaks down due to non-linear effects: the next target atom for the incident projectile may be one already in motion from the collision. At very high energies (exceeding 1 MeV) for any projectile–target pair, and at somewhat lower energies when the projectile is lighter than the target, consideration must be given to inelastic energy loss channels. The energy loss equation for the ion sputtering process becomes more complex [24]. At lower energies than those typically associated with the well-understood cascade region (less than 1 keV), energy transfer to the recoil atoms in the cascade is not isotropic. The initial target atom struck by the projectile absorbs most of the lost energy and reacts independently. This atom becomes the primary knock-on and may subsequently produce additional knock-on atoms. The SBE is defined as the energy required to remove an atom from the top surface layer into the vacuum during ion sputtering. It is a crucial parameter for characterizing the ion sputtering process. As the SBE increases, the sputtering yield gradually decreases [25]. This certainly applies to mixtures, where the proportion of two different elements affects the SBE. For a single substance, the SBE has a fixed value. However, analyzing the SBE of mixtures allows us to understand its impact on the sputtering process. The SBE influences the sputtering threshold energy E_th and the extent of this influence can be calculated using specific formulae. Subsequently, the classical sputtering yield fitting formula can be employed to analyze how an increase in SBE affects the sputtering yield [26, 27]

where S_BE is the SBE, E_c(A−B) is the covalent bonding energy, E_i(A−B) is the ionic bonding energy, t_i and t_c are the weighting coefficients, M₁ is the mass of the incident ion, M₂ is the mass of the target atom and A is a constant related to nuclear stopping power.

From the analysis above, the influence of the atomic number of ions on sputtering can be explained by the fact that materials with different atomic numbers have different relative molecular masses (or masses, for materials that are part of a mixture). This leads to different collision mechanisms and motion in electric fields, resulting in different sputtering behaviors. For the target, its molecular weight mainly determines the sputtering results on the target surface, such as cascade sputtering and elastic collisions, thus leading to the clustering behavior of the fifth group (brown in figure 4(b)).

The remaining clustering mainly consists of four categories: melting point, boiling point, heat of vaporization and heat of formation. These four categories of properties are all related to energy transfer during the state change process of the material [28]. These properties are all related to the intermolecular forces within a substance, implicitly reflecting the cohesive energy and ideal enthalpy of the substance’s surface [29, 30]. This theory is considered a classical explanation for the sputtering process. In sputtering, when several components are present, the more volatile substances are lost first, leading to a change in composition [24]. Summarizing the above phenomena, it is highly reasonable to model the sputtering process using ion mass, ion evaporation enthalpy, target mass and target evaporation enthalpy. These four properties encompass all the physical variables that describe the sputtering process, including collisions, energy transfer, target motion modes and yields.

To ensure the accuracy of our predictions of corrosion parameters, we ultimately selected a neural network model construction with a group number of 4. The final selection for the inputs of the neural network prediction model included the atomic mass of ions, the atomic mass of the target material as well as the heat of vaporization for ions and for the target material. This selection is primarily based on considerations of the physical phenomena. The corrosion process can be viewed as an energy exchange process, where the kinetic energy of ions is transferred to the target, causing the target to gain energy and emit substances. Therefore, attributes related to kinetic energy, such as atomic mass, along with properties associated with phase transitions, namely the heat of vaporization, were chosen.

3.2   ANN analysis

First, we compared the results of neural networks with different network architectures. The specific network architectures are outlined in the table 3. The input layer is divided into two categories: for Models 1–3 the input layer consists of four features, while for Models 4–6 the input layer comprises all attributes in the dataset. The output layer consists of three features: E_th, E_max and SY_max.

Table  3.  Neural networks configuration.

Number	Input layer	Hidden layer	Output layer
Model 1	4	10×10	3
Model 2	4	5×5	3
Model 3	4	15×15	3
Model 4	17	10×10	3
Model 5	17	5×5	3
Model 6	17	15×15	3

 | Show Table

DownLoad: CSV

First, a comparison was made between Models 1 and 6 to illustrate the accuracy of feature selection using the HCA algorithm. The coefficient of determination, denoted as R², is used in the model, and its formula is shown in below (figure 5 presents the R² values for different models)

Figure  5.  R² values of the six models.

DownLoad: Full-Size Img PowerPoint

where ${y_i}$ is the true value, ${y_i}$ is the predicted value and $\overline y$ is the average of the true values.

Among the six models, the optimal model, Model 2, is a 5×5 structure, while the lowest performance is observed in Model 5. This indicates that the HCA algorithm effectively increases regression accuracy by eliminating unnecessary interfering descriptors. For models with an input layer of 4, the overall R² values are higher than those with an input layer of 17. Comparing Model 1 with Model 3, the 5×5 structure shows higher accuracy. Generally, more complex network structures exhibit stronger fitting capabilities, but our model shows inconsistency. This inconsistency might be due to the different sources for the datasets. However, overall, the structure with four descriptors outperforms the one with 17 inputs.

3.3   Model accuracy

To validate the accuracy and applicability of the model, a Kaufman ion source was used for sputtering corrosion experiments. The experimental setup was configured as described in section 2.4. Figure 6 illustrates the results of the corrosion experiment. The blue data points represent the experimental test points, while the broken red curve represents the fitted quadratic curve. The zero-crossing and maximum points of the red curve are respectively E₀ = 114.02 eV, and SY_max = 4.3670 atom/ion and E_max = 697 eV.

Figure  6.  Sputtering yield as a function of ion energy, with experimental data points (blue triangles) and polynomial fit (red dashed curve).

DownLoad: Full-Size Img PowerPoint

V₀, SY_max, and V_max were utilized for comparison between actual values and predicted values and the results are depicted in table 4.

Table  4.  Comparison between predicted values and experimental values.

Feature	Experimental value	Predicted value	Relative error
Eth (eV)	114.02	117.8	3.23%
Emax (eV)	697	1995.7	186.32%
SYmax (atoms/ion)	4.3670	4.211	3.577%

 | Show Table

DownLoad: CSV

The experimental and predicted values for E_th and SY_max have a small error, but there is a large error of 186.32% in our prediction for E_max. This discrepancy is attributed to an inability to increase the power of the ion source, leading to inaccuracies in the fitted E_max value. At the same time, DT4C is 99.5% pure iron containing a small amount of carbon that affects its sputtering erosion E_max value.

4.   Conclusions

The research in this paper primarily focuses on establishing a concise yet highly accurate sputtering regression model. This model is based on the hierarchical clustering analysis (HCA) and artificial neural network (ANN) algorithms. Four attributes are selected for parameter regression of the sputtering model based on the HCA algorithm. The effectiveness of the simplified attributes is demonstrated using the ANN algorithm, and the impact of different ANN structures on varying numbers of attributes is also elucidated.

The results demonstrate that a combination of the HCA–ANN algorithms can effectively reduce the regression input parameters of the sputtering process. The physical significance of the HCA algorithm process is also explained. Finally, a 5×5 ANN network is established for the regression of the sputtering process. Validation of the ANN model is conducted based on specially designed sputtering experiments. Apart from E_max, the error of the ANN compared with the experimental results is less than 5%.