## Abstract

Silicon is one of the commonly used semiconductors for various industrial applications. Traditional silicon synthesis methods are often expensive and cannot meet the continuously growing demands for high-purity Si; electrodeposition is a promising and simple alternative. However, the electrodeposited products often possess nonuniform thicknesses due to various sources of uncertainty inherited from the fabrication process; to improve the quality of the coating products, it is crucial to better understand the influences of the sources of uncertainty. In this paper, uncertainty quantification (UQ) analysis is performed on the silicon electrodeposition process to evaluate the impacts of various experimental operation parameters on the thickness variation of the coated silicon layer and to find the optimal experimental conditions. To mitigate the high experimental and computational cost issues, a Gaussian process (GP) based surrogate model is constructed to conduct the UQ study with finite element (FE) simulation results as training data. It is found that the GP surrogate model can efficiently and accurately estimate the performance of the electrodeposition given certain experimental operation parameters. The results show that the electrodeposition process is sensitive to the geometric settings of the experiments, i.e., distance and area ratio between the counter and working electrodes; whereas other conditions, such as the potential of the counter electrode, temperature, and ion concentration in the electrolyte bath are less important. Furthermore, the optimal operating condition to deposit silicon is proposed to minimize the thickness variation of the coated silicon layer and to enhance the reliability of the electrodeposition experiment.

## Introduction

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.

_{4}as solute [8]. During the reactions, SiCl

_{4}is electrochemically reduced to form silicon through the following equations:

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^{4+} 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, *p*H, 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^{4+}*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.

## Finite Element Model Development

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^{4+} 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.

### Governing Equations.

^{4+}and Cl

^{−}ions in the electrolyte

**N**are described by the Nernst-Planck Eq. (3), including the diffusion term in concentration gradient and the migration term in an electric field. The convection of electrolyte in a flow field is neglected

_{i}**N**is the flux of ion

_{i}*i*in mol/(m

^{2}·s),

*D*is the diffusion coefficient in m

_{i}^{2}/s,

*c*is the concentration of the ion

_{i}*i*, in mol/m

^{3},

*z*is the valence,

_{i}*u*is the

_{i}*i*th ion mobility in s·mol/kg,

*ϕ*is the electrolyte potential in V, and

_{l}*F*= 96487 C/mol is the Faraday's constant. Furthermore, the assumption of electroneutrality is held for species concentrations [39]. The current density

*i*and current balance in the electrolyte are described in the following equation:

_{l}*i*

_{loc}is expressed as a function of the overpotential

*η*using the Butler–Volmer equation

where *i*_{0} is exchanged current density in A/m^{2}, *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^{3} are the atomic weight and density of silicon, respectively.

*k*is defined as

where *d*_{min} and *d*_{avg} are the minimum and average coating thickness, respectively.

### Model Implementation.

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^{4+} and Cl^{−} ions and the electroneutrality condition. The other boundaries are assumed to be insulated.

*d*is set to be 5 mm. The area of the anode is fixed as 12 mm; and the area ratio of cathode over anode

*r*is 1 as the baseline condition. The potential of the anode

*V*is 0 V, whereas the potential of the cathode is fixed at −1.7 V. The temperature

*T*is 293 K. The initial concentration of Si

^{4+}in the electrolyte

*c*is 200 mol/m

^{3}. The diffusion coefficients

*D*of Si

_{i}^{4+}and Cl

^{−}ions are assumed to obey the Arrhenius type relationship (Eq. (8)) [41]. Table 1 summarizes the baseline condition, as well as the design space of five operation parameters as the primary sources of uncertainty. Table 2 lists the other electrochemical parameters used in the model

where *D*_{0} = 4.13 × 10^{−10} is a pre-exponential factor, and *E _{D}* = 19,200 J/mol is the activation energy for diffusion.

Operation parameters | Baseline | Range |
---|---|---|

Distance between the electrodes d (mm) | 5.0 | 5.0–15.0 |

Area ratio of cathode over anode r | 1.0 | 0.2–1.8 |

Potential of the anode V (V) | 0.0 | 0.0–2.0 |

Temperature T (K) | 293.0 | 293.0–333.0 |

Initial concentration of Si^{4+}c (mol/m^{3}) | 200.0 | 120.0–240.0 |

Operation parameters | Baseline | Range |
---|---|---|

Distance between the electrodes d (mm) | 5.0 | 5.0–15.0 |

Area ratio of cathode over anode r | 1.0 | 0.2–1.8 |

Potential of the anode V (V) | 0.0 | 0.0–2.0 |

Temperature T (K) | 293.0 | 293.0–333.0 |

Initial concentration of Si^{4+}c (mol/m^{3}) | 200.0 | 120.0–240.0 |

Parameter | Value |
---|---|

Potential of the cathode (V) | −1.7 |

Equilibrium potential at the silicon surface (V) | 0.0 |

Cathodic charge transfer coefficient α_{c} | 0.5 |

Anodic charge transfer coefficient α_{a} | 0.5 |

Valence of Si^{4+} | 4.0 |

Valence of Cl^{−} | −1.0 |

Electrical conductivity of silicon | 1000 S/m |

Parameter | Value |
---|---|

Potential of the cathode (V) | −1.7 |

Equilibrium potential at the silicon surface (V) | 0.0 |

Cathodic charge transfer coefficient α_{c} | 0.5 |

Anodic charge transfer coefficient α_{a} | 0.5 |

Valence of Si^{4+} | 4.0 |

Valence of Cl^{−} | −1.0 |

Electrical conductivity of silicon | 1000 S/m |

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.

## Gaussian Process-Based Surrogate Model

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].

*G*(

*x*) has a form of

*S*(x) is a Gaussian stochastic process with zero mean and certain covariance matrix. The covariance function between two inputs

*x*and

_{i}*x*is expressed as the squared-exponential form

_{j}*a*and

_{p}*b*are parameters of the GP model, and

_{p}*k*is the number of the input parameters

**x**.

**R**(

**x**,

_{i}**x**) is an

_{j}*n × n*symmetric correlation matrix and is defined based on the distance between two sample points

**x**and

_{i}**x**. It is worth noting that, other forms of the covariation function, like the Mercer kernel, can be utilized to address more complex problems [43]. The GP model is then trained with the imported simulation results (

_{j}**X**,

_{S}**G**) to generate the initial thickness variation function

_{S}*G*(

*x*). For any new point $x\u2032$, its performance can be estimated via maximizing the likelihood function [44] as a normal distributed variable with mean $G\u0302(x\u2032)$ and variance $e\u2322(x\u2032)$

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

*X*

_{S},

*G*

_{S}) are generated by the FE model to develop the initial GP model, where

*X*

_{S}is the model input parameters, i.e., the five operation parameters in the study, and

*G*

_{S}is the corresponding output performance, i.e., the thickness variation. In this study, a uniform grid sampling strategy is used to select the input

*X*

_{S}. (2) In the next step, the simulation results (

*X*

_{S},

*G*

_{S}) are stored in a sample data pool and used to construct the predictive models

**employing the GP-based surrogate modeling technique [44–46]. (3) The most important sample points**

*M**x** are then selected using the maximum expected information value-based adaptive sampling technique. The performances of these points

*G** are evaluated with the FE model, and the new sample points (

*x**,

*G**) are imported into the sample data pool to upgrade the surrogate model. (4) Cumulative confidence level (CCL) is used to qualify the accuracy of reliability of the surrogate model via Eqs. (13) and (14) [44]

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.

## Results and Discussion

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.

### Finite Element Model Results.

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.

### Sensitivity Analysis.

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^{4+}*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^{4+}*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.

*k*with respect to the five operation parameters. To avoid scaling issues, the differences are normalized to the values used in the baseline design [48]. Equation (15) is used to calculate the normalized sensitivity

*σ*of the thickness variation

*k*

where $q\xafi$ is the value of the *i*th operation parameter in the baseline design, $k\xafi$ is the calculated thickness variation with respect to the *i*th operation parameter, and $\Delta 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^{4+} 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.

### Surrogate Model Prediction Results.

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.

### Uncertainty Quantification Analysis.

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^{5} 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^{4+}*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.

### Design Improvement Strategy.

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.

Optimal condition | Baseline | ||
---|---|---|---|

Input variables | Distance between the electrodes d (mm) | 5.5 | 5.0 |

Area ratio of cathode over anode r | 1.1 | 1.0 | |

Potential of the anode V (V) | 0.2 | 0.0 | |

Temperature T (K) | 322.3 | 293.0 | |

Initial concentration of Si^{4+}c (mol/m^{3}) | 200.7 | 200.0 | |

Thickness variation k from surrogate model | 9.97% | 11.10% | |

Thickness variation k from FE simulation | 9.96% | 11.10% |

Optimal condition | Baseline | ||
---|---|---|---|

Input variables | Distance between the electrodes d (mm) | 5.5 | 5.0 |

Area ratio of cathode over anode r | 1.1 | 1.0 | |

Potential of the anode V (V) | 0.2 | 0.0 | |

Temperature T (K) | 322.3 | 293.0 | |

Initial concentration of Si^{4+}c (mol/m^{3}) | 200.7 | 200.0 | |

Thickness variation k from surrogate model | 9.97% | 11.10% | |

Thickness variation k from FE simulation | 9.96% | 11.10% |

## Conclusion

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^{4+} 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.

## Funding Data

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).

## References

**10533**, p. 105332S.10.1117/12.2294944