Uncertainty Quantification Analysis on Silicon Electrodeposition Process Via Numerical Simulation Methods

Silicon is one of the most important and commonly used semiconductors for the applications of computer chips, photovoltaic cells, and optoelectronic devices, due to its chemical stability and unique optical and electronic properties [1–3]. Various methods have been used to synthesize the Si-based nanostructures, such as sputtering, plasma deposition from silane, laser ablation, and electron-beam evaporation. [4–7]. However, these methods are usually expansive, complicated, and time-consuming procedures, which are based on the use of high-temperature vapors and melts, and cannot uniformly deposit silicon onto complex three-dimensional substrates [8]. With the continuously growing demand for high-purity Si, the development of advanced deposition techniques to produce low-cost Si devices has become a research interest.

Nevertheless, even though electrodeposition is a relatively simple method, various sources of uncertainty may inevitably be introduced during the operating process, such as the variation of the experimental setting, deviation of material properties from the nominal values, and the edge effect, which could result in electrochemical heterogeneity and uneven deposits [16]. Therefore, to improve the quality of the electroplated products, i.e., the uniformity of deposits, it is crucial to better understand the influences of different sources of uncertainty on the uniformity of the resulting products through uncertainty quantification (UQ) analysis and further obtain the optimal experimental conditions. Previously, different experimental and numerical methods have been applied to investigate the influences of various experimental operation parameters on the nonuniformity of the electrodeposition process [17–19]. However, these methods are often expensive and time-consuming, which are not suitable to be used for UQ study. In this study, five experimental operation parameters are selected as the primary sources of uncertainty [20], including the distance between the anode (counter electrode) and cathode (working electrode) d, the area ratio of cathode over anode r, the potential of the anode V, the temperature of the electrolyte bath T, and initial concentration of Si⁴⁺ in the electrolyte c. The UQ analysis is then performed on the silicon electrodeposition process, which includes two steps. Firstly, a finite element (FE) model is constructed to simulate the silicon electrodeposition process under various experimental conditions and to predict the thickness and uniformity of the deposited layer; second, a Gaussian process (GP) based surrogate model is developed using the simulation results as training data to reproduce the outputs of the FE model efficiently and accurately and to conduct the UQ on the silicon electrodeposition process.

Numerical simulation models like FE analysis have been widely adopted in the electrochemistry field to avoid conducting a large number of experiments that are time-consuming and expensive [21–25]. These models can help understand various complex coupled phenomena and the underlying mechanisms, such as electrochemical metal-reduction reactions on moving deposition surfaces, and transport of charged and neutral species via diffusion, electromigration, and convection in the electrolyte, etc. More importantly, through precise prediction of the electrodeposition process, the size and microstructure of deposits can be directly mapped to the operation parameters, e.g., geometries of the electrodes, applied potential, bath composition, pH, and temperature, etc. Chen et al. [26] reported the development of FE models to simulate the two-dimensional (2D) processes of nickel electrodeposition in the lithography, electroforming and molding microfabrication. Species concentrations, electrolyte potential, flow field, and positions of the moving deposition surfaces were computed by coupling the arbitrary Lagrangian–Eulerian approach to track the entire transient deposition processes and to calculate local current densities that influenced the microstructure and functional/mechanical properties of the deposit. Lupo et al. [27] investigated the dendrite formation phenomena during the electrodeposition process using FE methods and validated the model with experimental results. It was found that the FE models could precisely estimate the growth of dendrites and their detailed size distribution over the whole electrode surface. Wheeler et al. [28] presented the modeling of superconformal electrodeposition/superfill using the level set method in the Eulerian framework to simulate the additive accumulation and conservation on a deforming copper/electrolyte interface. However, one critical disadvantage of the FE models is that they are often composed of multiple components, partial differential equations, and physical modules, with the degrees-of-freedom of million level; as a result, the computational cost is rather high, and it is not feasible to implement the UQ analysis via this technique.

Machine learning-based surrogate models have drawn great attention in the material science field recently [29], which are constructed through data-driven, bottom-up techniques. Different from FE models, Machine learning surrogate models can directly capture the general underlying trends of the performances of the systems and efficiently approximate the outputs of large-scale complex systems given any desirable input variables through the realization of a regression model and a random process [30,31]. The hyperparameters and random process in the regression model are trained using the FE simulation results and experimental data. Surrogate models with different predictors, e.g., the response surface method [32], neural networks [33], and GP (as known as Kriging) [34–36], have been applied for the structural dynamic analysis of various engineering systems. In the previous study, we adopted GP-based surrogate models to help understand the lithiation-induced stresses [22] and crack initiation phenomenon [25,37] in the Si-based lithium-ion battery anodes over the design space. It was found that the GP model could effectively estimate the lithiation-induced stress concentration and crack formation in the Si layers without losing fidelity.

In this study, the UQ analysis is performed on the silicon electrodeposition process via numerical simulation methods. The overall modeling procedures contain two steps. In the first step, a 2D multiphysics-based FE model is developed to simulate the electrodeposition process; the thickness variations of electroplated silicon layers are calculated with different combinations of operation parameters as input variables, including the distance d and area ratio r between the counter and working electrodes, the potential of the counter electrode V, the temperature of the electrolyte bath T, and initial concentration of Si⁴⁺c; the design sensitivity analyses of these parameters are conducted. In the second step, the GP-based surrogate model is trained with the FE simulation results; the propagation of uncertainty from inputs over the output space is evaluated; and the improvement strategy of the experimental setup for silicon electrodeposition is provided to minimize the thickness variations of the coated silicon layer.

In this section, the multiphysics-based FE model is established in the commercially available FE simulation tool (comsolmultiphysics) to simulate the silicon electrodeposition process. The electrochemical kinetic on the electrode surfaces, the deformation of the electrode boundary, the transportation of Si⁴⁺ and Cl⁻ and potential change in the electrolyte are calculated. In the following, governing equations that control the reactions and deformation are discussed; then, the model implementation in the fe software is introduced.

where i₀ is exchanged current density in A/m², T is the temperature in K, R = 8.314 J/(K·mol) is the universal gas constant, α_c and α_a are the cathodic and anodic charge transfer coefficient, respectively.

where d_loc is the coating thickness, t is the plating time in s, M = 28 g/mol, and ρ = 2300 kg/m³ are the atomic weight and density of silicon, respectively.

where d_min and d_avg are the minimum and average coating thickness, respectively.

The silicon electrodeposition process is solved in 2D continuum space using a tertiary current distribution module. Figure 1 illustrates the schematic of the model geometry, containing the finite element mesh. The deposition at the working electrode (cathode) and the electricity supply at the counter electrode (anode) are set to take place at 100% current yield. The process is assumed to be time dependent, where the silicon electrodeposition is simulated as the moving boundary of the cathode. The model is governed by the mass conservation for Si⁴⁺ and Cl⁻ ions and the electroneutrality condition. The other boundaries are assumed to be insulated.

where D₀ = 4.13 × 10⁻¹⁰ is a pre-exponential factor, and E_D = 19,200 J/mol is the activation energy for diffusion.

The mesh size of the cathode boundary is set to be 0.05 mm to precisely predict the growth of the deposited silicon layer. The electrolyte domain is meshed using free triangular elements. Complete geometry consists of 7041 domain elements and 373 boundary elements.

To further improve the efficiency to explore the operation parameter space without losing accuracy, a GP-based surrogate model is constructed. GP surrogate model is a statistical procedure that generates an estimated surface from a scattered set of points via a covariance governed Gaussian process interpolation method. It requires training samples for model construction and predicting responses of new sample points [42].

where r is the correlation vector between $x′$ and the input parameters x, A is an n × 1 unit vector.

where CL(x_m,i) represents the confidence level of the classification for sample x_m,i, Φ is a standard normal cumulative distribution function, and |.| is the absolute operator. Note that CCL(M, X_m) is always a positive value within [0.5, 1]. Eventually, the updated GP model acquires enhanced fidelity, which can be used for the thickness variation prediction and optimization design.

In this section, the simulation results from the FE model are first discussed. Then, the design sensitivity analysis is applied to the five operation parameters to evaluate their influences. Afterward, the GP-based surrogate model results are introduced; and the uncertainty quantification study results of the silicon electrodeposition process are discussed. Finally, a design improvement strategy for the experimental condition is provided based upon the surrogate model.

Figure 2 illustrates the thickness change of the electrodeposited silicon layer along the cathode surface from the developed FE model when the simulation runs from 0 to 600 s. As can be seen, a layer of silicon is steadily growing on the cathode substrate, whereas a rather uneven coating is shown: the deposits at the left (<2 mm) and right (>18 mm) edges of the cathode substrate are much thicker than those in the center. For instance, after 600 s deposition, the center area is about 12.7 mm thick, and remains relatively consistent and uniform; however, the two edge regions can reach up to 56.6 mm. Consequently, the thickness variation of the deposition process is calculated to be 11.1% via Eq. (7). The nonuniformity phenomenon can be explained by the edge effect of electrode current density [47]. As illustrated in the inset of Fig. 2, there are concentrated electrode and electrolyte current densities near the edges, leading to high electrochemical reaction rates and thick coatings; on the contrary, the electric field and current density flow in the center area keep consistent, resulting in the uniform coating. It is worth noting that, the nonuniformity is not getting exaggerated with the increase of simulation time. This may be because the electrical conductivity of silicon is low compared to the metal materials; as the silicon layer becomes thicker near the edges, it becomes more and more difficult to transfer charges from the silicon–electrolyte interface to the cathode substrate, which gradually reduces the electrode current density and limits the electrodeposition reaction rates. In the following discussions, the thickness variation of the coating layer is treated as the target property to investigate the influences of different operation parameters and to conduct the uncertainty quantification analysis on the silicon deposition process.

Before the construction of the surrogate model, a sensitivity analysis is conducted on the silicon deposition process to evaluate the influences of different operation parameters on the thickness variation of coated Si layer, to identify the critical parameters with high sensitivity, and to minimize the discrepancy between simulation results and surrogate model.

Figure 3 shows the effects of five operation parameters, including the distance between the two electrodes d, area ratio of cathode over anode r, the potential of the anode (counter electrode) V, temperature T, and initial concentration of Si⁴⁺c, on the silicon layer thickness variation k after 600 s simulation. It can be observed that five different operation parameters show distinct impacts.

First, the thickness variation k has a nonlinear relation with the distance between two electrodes d, as shown in Fig. 3(a). When d = 5 mm, k is rather small, only around 11.1%; with the increase of d, k increases simultaneously; when d ≥ 20 mm, however, further extending d will no longer affect k, which becomes a constant, 23.9%. This can be explained as follows. While the distance between the working and counter electrodes is small, the electric field in the electrolyte is strong so that the charge transfer and ion transportation between two electrodes are constrained and keep uniform, specifically near the cathode substrate; therefore, the current density concentration and sharp spikes of deposited silicon near the electrode edges are limited. On the contrary, extending the distance between electrodes could lower the electric field strength, release the constraints on ion movements, and enlarge the gradient of electrolyte potential field near the electrode edges, leading to an enhanced edge effect and a high thickness variation. Nevertheless, when the electrodes are separated more than 20 mm, the effects of the reduction of the electric field strength become minimal, thus the uneven distribution of the electrode current density will no longer be exaggerated; as a result, k becomes independent of d.

Second, the area ratio between two electrodes r shows a significant impact on the thickness variation k, as illustrated in Fig. 3(b). On one hand, when the area of cathode is smaller than that of the anode, i.e., r < 1, the edge effect is severe during the electrodeposition process, causing large electrolyte potential field gradients (inset of Fig. 3(b)) and thicker silicon coatings around the edges of cathode substrate. The edge effect can be greatly mitigated when the cathode has a slightly larger area (i.e., r ≈ 1.25); consequently, k is minimized, approaching 9.8%. On the other hand, when the cathode is oversized compared to the anode (r > 1.5), the electrolyte potential gradient outside the anode range (circled in the inset of Fig. 3(b)) is rather small, suggesting that the driven force for the silicon electrodeposition reaction is largely reduced; whereas, the electrolyte potential gradient in the center region remains high. Eventually, different from the previous scenarios presented in Fig. 2, the coating layer has a thicker center but thinner edges; furthermore, due to the large variation of the potential gradient between the center and edge regions, the thickness variation k of the coating layer is significantly increased.

Third, the potential of the anode V and initial concentration of Si⁴⁺c have little influences on the thickness variation k, as shown in Figs. 3(c) and 3(e), respectively, where k is fluctuating within the examined ranges and shows nonlinear relationships with V and c. This is because, even though higher V or c could lead to the improved driven force for the electrochemical reaction and correspondingly higher silicon coating rate, they will barely change the electrolyte potential gradient pattern, and thus the thickness variation is only merely affected.

Lastly, as presented in Fig. 3(d), the thickness variation k increases along with the rising temperature T in a nearly linear manner. This is caused by two reasons: first, as introduced in Eq. (8), the diffusion coefficients of different ions will increase at elevated temperature, indicating easier transportation of these species in the electrolyte; in addition, higher T could enhance the electrochemical reaction rate and increase the current density at the electrode–electrolyte interface (Eq. (5)). Both these factors will exaggerate the edge effect, resulting in the nonuniformity of the electrodeposited products.

where $q¯i$ is the value of the ith operation parameter in the baseline design, $k¯i$ is the calculated thickness variation with respect to the ith operation parameter, and $Δqi$ is the difference of the parameter. The sensitivity results are listed in Table 3. The results reveal that the silicon electrodeposition process is sensitive to the geometric setting of the experiments, including the distance and area ratio between two electrodes. Meanwhile, the other parameters, i.e., the potential of the anode, initial Si⁴⁺ concentration, and the temperature of the electrolyte bath, show relatively low influences. In the following, these five operation parameters are used to construct the GP-based surrogate model.

The GP-based surrogate model is constructed based upon the FE simulation and the sensitivity analysis results using the uniform grid sampling within the design space of the operation parameters. Particularly, from the sensitivity analysis, two critical variables are identified that have big impacts on the uncertainty of the system performance. Afterward, more densely arranged samples for them within the examined ranges are chosen to calculate thickness variation k through the FE simulations, to model the complex nonlinear relationship between the inputs and outputs. Therefore, the sensitivity analysis can support the construction of the surrogate model and reduce the error. To have better visualization and explanation of the surrogate model, two predicted performance surface figures are presented as examples. Figures 4(a) and 4(b) illustrate the performance surfaces of silicon electrodeposition thickness variation k with respect to two selected combinations of operation parameters, i.e., distance d and area ratio r between two electrodes, as well as the potential of anode V and temperature of the electrolyte bath T, respectively. The scatters represent the FE simulation results. A highly nonlinear performance surface of k in relation to the geometric setting of the experiments, i.e., d and r, is observed in Fig. 4(a). First, when r is less than 0.75, k remains almost constant at a high level, approximately 25%, and becomes independent of d. Besides, with the variation of d, r has distinct influences: when the electrodes are close to each other, the valley of the performance surface lies at r = ∼1.25, and a sharp peak is shown at the corner of d = 0 mm and r = 2; however, when two electrodes are far apart, the performance surface valley is found at r = ∼1.75, and k will monotonically increase with the decrease of r. In addition, as shown in Fig. 4(b), k is fluctuating within the range of 11.5%–12% in the examined parameter space and is insensitive to the change of T and V. By reducing T, k can be slightly lowered.

The results from the GP surrogate model are in consistent with FE model predictions and sensitivity analysis results, suggesting the high fidelity of the surrogate model. In the next step, the GP surrogate model is used for the uncertainty quantification analysis of the silicon electrodeposition process.

In this study, the Monte Carlo simulation technique is adopted to perform the UQ study using the GP surrogate model. This approach provides several merits [49]. First, the Monte Carlo technique could directly map the input parameter dataset to the output performance space, and acquire the statistical distribution of the output given the distribution of the inputs, which is the most reliable approach to handle the complex nonlinear relationship between input and output [50]. Second, the GP surrogate model is used as a mapping function instead of the computationally expensive FE simulations, which could greatly reduce the computational burden without losing much accuracy when dealing with a large amount of data points.

The five input operation parameters in Table 1 are assumed to follow the truncated normal distribution [51]. Figure 5 shows the histograms of two input parameters, i.e., distance d and area ratio r between two electrodes, as examples. Correspondingly, 10⁵ data sets are formed using the random number generator, each consisting of the aforementioned five operation parameters, where the combination of the value for each parameter has no preference and is totally random. Afterward, the distribution of the output, i.e., thickness variation k, is obtained by feeding the input data sets into the GP surrogate model. Figure 6 illustrates a histogram of the estimated output from the GP surrogate model given the generated input data sets. The output distribution is not Gaussian due to the nonlinear relationship between the input parameters and output performance. The average value and standard deviation are 0.185 and 0.030, respectively.

The probability distribution of the silicon electrodeposition thickness variation k is helpful to understand the uncertainty and variability of the electrodeposition process and could provide a guideline for the design improvement strategies to reduce the noise and enhance the stability. Figures 7(a) and 7(b) present the histogram of the silicon electrodeposition thickness variation k with respect to the change of coefficient of variations of the distance between the electrodes d and initial concentration of Si⁴⁺c, respectively, as examples. As can be seen, these two variables show distinct influences on the uncertainty of the electrodeposition process. On one hand, with the decrease of coefficient of variation of d from 14.4% to 3.6%, the average, and deviation of thickness variation can be decreased greatly. On the other hand, by reducing the coefficient of variation of c, the uncertainty of the system is not much affected. These indicate that different input variables have various contributions to the uncertainty of the Si anode; as a result, in order to acquire a better design of the silicon electrodeposition process with more stable performances and less uncertainty, the most effective approaches are to control the variations of distance d and area ratio r between two electrodes.

In addition to the uncertainty quantification study, which can quantify the uncertainty of the silicon electrodeposition process, the surrogate model could offer another critical advantage: to guide the design optimization of the process to minimize the thickness variation k and enhance the reliability performances. By scanning the surrogate model results within the examined parameter space, the optimal operation parameters to electroplate silicon can be identified, which could lead to the smallest k. Table 4 compares k between the predicted optimal operation parameters and the baseline condition. It can be observed that, compared to the baseline condition, the optimal condition is updated, where r, V, and c are hardly changed, d is slightly extended and T is modified to be 322.3 K; as a result, the k is decreased from 11.1% to 9.97%, suggesting that the nonuniformity phenomenon can be effectively mitigated so that the silicon deposition process is upgraded. Moreover, the optimal condition is imported into the FE model to validate the GP surrogate model estimation. The thickness variation k from the FE model is almost identical to the predicted value from the surrogate model, further demonstrating the accuracy of the GP-based surrogate model.

In this paper, we conduct the UQ analysis on the silicon electrodeposition process to explore the influences of various sources of uncertainty on the quality of the electrodeposited products and to find the optimal experimental conditions providing minimized thickness variation of the coated silicon layer. Five experimental operation parameters are chosen as the primary sources of uncertainty and investigated, including the distance d and area ratio r between the working (cathode) and counter (anode) electrodes, the potential of the counter electrode V, temperature of the electrolyte bath T, and initial concentration of Si⁴⁺ in the electrolyte c. The overall modeling procedures contain two steps. In the first step, a 2D multiphysics-based FE model is developed to simulate the electrodeposition process with different combinations of operation parameters as input; the thickness variation of the electroplated silicon layer is calculated; the design sensitivity analyses are then conducted to identify the critical parameters with high sensitivity. In the second step, the GP-based surrogate model is constructed with the FE simulation results as training data; the propagation of uncertainty from inputs over the output space is analyzed, and the optimal experimental condition is proposed to minimize the thickness variations of the coated silicon layer.

It is found that a rather uneven coating is obtained, where the deposits near the left and right edges of the working electrode are much thicker than those in the center. This is mainly due to the edge effect, i.e., the concentrated electrode and electrolyte current densities near the electrode edges during the reaction. Besides, the electrodeposition process is rather sensitive to the geometric settings of the experiments, i.e., distance and area ratio between the two electrodes; both extending the distance between two electrodes and using oversized or undersized working electrode compared to the counter electrode can significantly exaggerate the thickness variation of the electrodeposited product. On the contrary, other conditions, such as the potential of the counter electrode, temperature, and ion concentration in the electrolyte bath, show less impacts. In addition, the developed GP surrogate model can efficiently and accurately estimate the performance of the electrodeposition given the experimental operation parameters. Different operation parameters are found to have distinct contributions to the uncertainty of the system; by assuming that the parameters follow a normal distribution, the average thickness variation and standard deviation are 0.185 and 0.030, respectively. Furthermore, the optimal operating condition is proposed by the GP model and validated using FE simulation to minimize the thickness variation of the coated silicon layer and to enhance the reliability of the electrodeposition experiment. To narrow the distance between two electrodes, to use a slightly larger working electrode, and elevating the temperature of the electrolyte bath are key to reduce the thickness variation.

• Office of Naval Research (ONR) through the Defense University Research-to-Adoption (DURA) Initiative (Grant No. N00014-18-S-F004; Funder ID: 10.13039/100000006).

• National Science Foundation (NSF) Engineering Research Center for Power Optimization of Electro-Thermal Systems (POETS) with Cooperative Agreement (Grant No. EEC-1449548; Funder ID: 10.13039/100000001).