Multi-Physics Modeling of Electrochemical Deposition

Abstract

Electrochemical deposition (ECD) is a common method used in the field of microelectronics to grow metallic coatings on an electrode. The deposition process occurs in an electrolyte bath where dissolved ions of the depositing material are suspended in an acid while an electric current is applied to the electrodes. The proposed computational model uses the finite volume method and the finite area method to predict copper growth on the plating surface without the use of a level set method or deforming mesh because the amount of copper layer growth is not expected to impact the fluid motion. The finite area method enables the solver to track the growth of the copper layer and uses the current density as a forcing function for an electric potential field on the plating surface. The current density at the electrolyte-plating surface interface is converged within each PISO (Pressure Implicit with Splitting Operator) loop iteration and incorporates the variance of the electrical resistance that occurs via the growth of the copper layer. This paper demonstrates the application of the finite area method for an ECD problem and additionally incorporates coupling between fluid mechanics, ionic diffusion, and electrochemistry.

1. Introduction

Electrochemical deposition (ECD) is an integral process in the fabrication of semiconductor devices. In addition, the applications spread through a variety of other fields such as the deposition of coatings for corrosion resistance [1] and water treatments [2]. ECD is not the only deposition technique used for growing copper layers in semiconductor devices, but it is a popular and inexpensive method, since it does not require a vacuum or heating of the system [3]. Additionally, this method is able to deposit high-aspect ratio features with a high degree of precision.

The process of ECD is shown in a simplified schematic in Figure 1.

The process applies a current I that flows from the anode to the cathode. Positive ions are shed from the anode and added to the electrolyte. Diffusion and migration of ions occurs in the electrolyte. The positive ions are deposited on the cathode surface and depletes the ionic concentration near the cathode-electrolyte interface. As time evolves, the copper layer continues to grow, and the rate of this growth is directly proportional to the applied current density. The copper layer is deposited edge high, center low. So it is imperative to control how the copper layer is deposited and one way to control the uniformity of the copper layer is to induce mixing of the electrolyte. In this paper the anode and cathode are copper and the electrolyte is a copper sulfate-acid solution.

Mathematical models of electrochemical phenomena in the literature showcase a spectrum of weak to strong coupling between electric potential, ionic concentration, and fluid motion. The work by Karcz [4], which investigates a tubular fuel cell, chose to account for the current density with a zero-dimensional or one-dimensional model, while maintaining a three-dimensional model for the fluid motion and ionic transport. Karcz concluded that the characteristics of the current density were similar in both models and the one-dimensional model compared favorably to experimental data. Introducing relationships such as the Butler-Volmer equation, which computes a current density based on a nonlinear dependence of the voltage difference between the anode and cathode and ionic concentration of the species that participates in chemical kinetics [5], strengthens the coupling of the electrochemical model and can be observed in work performed by Ritter et al. [3], Drese [6], Hughes, Baily and McManus [7], and Hughes et al. [8], where four different current distributions, at the cathode-electrolyte interface, are presented. These four regimes are tertiary, secondary, primary, and diffusion limited current distributions. This work uses the tertiary current distribution. The strongest coupling presented in the literature shows dependence of ionic concentration and electric potential in the equations that govern ionic transport and electrochemistry, see [2,9,10,11].

The overall objective of this work is to determine how turbulent mixing impacts copper layer growth. With that being said, this paper is the first building block to achieve that higher goal, where the focus will be on the electrochemistry and ionic transport. There have been others that have looked into mixing via a reciprocating paddle [10]. Motion of a paddle within the fluid or a deforming boundary due to depletion or plating of copper require a numerical approach that will account for these changes. One such approach is to use the Arbitrary Lagrangian-Eulerian (ALE) [10] or a dynamic mesh method [2] that moves the mesh based on the motion of a boundary. Other approaches include a level set method [7,8,9], an explicit interface tracking method (EITM) [12,13], and the finite area method (FAM) [1]. Here the FAM is used because the growth of the copper layer is small in comparison to the domain size, and, therefore, would have negligible effects on the motion of the fluid.

The paper is organized as follows: Section 2 provides a description of the mathematical models and key assumptions used in this work. Section 3 summarizes the numerical methods employed to solve the physics model. Section 4 presents the results of a validation study performed to assess the accuracy of the proposed approach. Section 5 extends the validation results by considering concentration dependence. Finally, Section 6 provides a summary of the research and some possible future extensions.

2. Problem Formulation

The physics that governs ECD couples fluid mechanics, ionic diffusion, and electrochemistry. To computationally model a complex system such as ECD, several mathematical, physical, and geometric assumptions must be made that dictate the level of coupling. It is necessary to have a strong enough level of coupling between the different physical models to demonstrate and predict the required physical behavior observed in experiments and practice.

2.1. Governing Equations

The equations that govern the fluid mechanics are conservation of mass and conservation of momentum. The conservation of mass for an incompressible fluid is

$$\nabla ·\mathbf{u}=0,$$

(1)

and the conservation of momentum, the Navier-Stokes equations, are

$$\frac{\partial \mathbf{u}}{\partial t}+\nabla ·(\mathbf{u}\otimes \mathbf{u})-\nabla ·(\nu \nabla \mathbf{u})=-\frac{1}{\rho }\nabla p+\mathbf{F},$$

(2)

where $\mathbf{u}$ is the fluid velocity, $\nu$ is the kinematic viscosity, $\rho$ is the fluid density, p is the fluid pressure, and $\mathbf{F}$ is any body force acting on the fluid.

To mathematically model the transport of copper ions within the electrolyte, the concentration of copper ions is tracked via diffusion, migration, and convection (if fluid motion is present). The governing equation for the transport of Cu²⁺ ions is

$$\frac{\partial C}{\partial t}=-\nabla ·\left[-D\nabla C-\left(\frac{FnD}{RT}C\right)\nabla \phi -C\mathbf{u}\right],$$

(3)

where C is the mass fraction of copper ion concentration, D is the diffusivity of copper ions, F is Faraday’s constant, n is the valence electrons for copper ions, R is the ideal gas constant, T is chamber temperature, and $\phi$ is the electric potential of the electrolyte. This equation is known as the Nernst-Planck equation. The first term on the right-hand side of Equation (3) represents diffusion, the middle term is the migration of the ionic concentration (this is a coupling term between the electric potential and ionic concentration), and the last term is convection due to fluid motion.

The electrochemistry is governed by ensuring that the electric flux in the system is divergence-free, where the flux is defined as

$$\mathbf{N}=-Fn\left[D\nabla C-\left(\frac{FnD}{RT}C\right)\nabla \phi +C\mathbf{u}\right].$$

(4)

Therefore, the governing equation is $\nabla ·\mathbf{N}=0$ or

$${\nabla }^2(FnDC)+\nabla ·\left[\left(\frac{F^2{n}^{2}D}{RT}C\right)\nabla \phi \right]=0,$$

(5)

since the fluid motion is quiescent and Equation (1) still holds $\nabla ·(C\mathbf{u})$ vanishes. If fluid motion was present the term $\nabla C·\mathbf{u}$ would still be included.

The equations that have been presented to this point only account for the electrolyte and do not consider any electrodynamics that occur from the growth of the copper layer. To account for the changes in electric potential on the plating surface the following equation must be solved:

$$\nabla ·({\sigma }_{eff}\nabla {\phi }_s)=-i,$$

(6)

where ${\sigma }_{eff}$ is the effective surface conductivity at the interface between the seed layer and electrolyte and i is the current density on the plating surface. The effective conductivity is calculated by considering the initial seed layer in addition to the growing copper layer; it is defined as:

$${\sigma }_{eff}=\sigma \delta +{\sigma }_0{\delta }_0,$$

(7)

where $\sigma$ is the conductivity of solid copper, ${\sigma }_0$ is the conductivity of the initial seed layer material, and ${\delta }_0$ is the initial seed layer thickness. The current density at the electrolyte-cathode interface is computed using the following:

$$i=\frac{FD{c}_b}{M_{Cu}}\left(\nabla C·\mathbf{n}-\frac{nF}{RT}C\nabla \phi ·\mathbf{n}\right),$$

(8)

where $c_b$ is the bulk concentration of copper ions in the electrolyte and $M_{Cu}$ is the molar mass of copper.

The Butler-Volmer equation Equation (9) is another expression for the current density, but here it used to calculate the overpotential $\eta$, which requires using a Newton-Raphson [14] method to solve the nonlinear equation. The Bulter-Volmer equation is

$$i=i_0\left[exp\frac{{\alpha }_{A}nF\eta }{RT}-Cexp\frac{-{\alpha }_{C}nF\eta }{RT}\right],$$

(9)

where $i_0$ is the exchange current density (cathodic current at equilibrium [15]), and ${\alpha }_A$ and ${\alpha }_C$ are the charge transfer coefficient for the anode and cathode, respectively. $\eta$ is then used to update the interface potential,

$${\phi }_I={\phi }_s-\eta ,$$

(10)

which is used as a boundary condition for the current algorithm. Equation (9) also enforces convergence of the current density inside the PISO loop, where the full algorithm will be discussed in more detail in Section 3.2.

Finally, the equation that governs the deposited copper is

$$\frac{\partial \delta }{\partial t}=\frac{i{M}_{Cu}}{{\rho }_{Cu}nF},$$

(11)

where $\delta$ is the deposited copper thickness and ${\rho }_{Cu}$ is the density of copper. All of the parameters that need to be set in the simulations are defined in

Table 1. Material and Fluid properties that will be calculated through the course of each simulation. Symbols, descriptions, and dimensions are provided here.

Constants	Description	Dimensions (SI)
$\rho$	Fluid density	$\left[\mathrm{kg}/{\mathrm{m}}^3\right]$
$\nu$	Kinematic viscosity of the fluid	$\left[{\mathrm{m}}^2/\mathrm{s}\right]$
D	Cu ion diffusivity	$\left[{\mathrm{m}}^2/\mathrm{s}\right]$
F	Faraday’s constant	$\left[\mathrm{A}\mathrm{s}/\mathrm{mol}\right]$
n	Ion valence	$\left[-\right]$
R	Universal gas constant	$\left[\mathrm{kg}{\mathrm{m}}^2/(\mathrm{K}\mathrm{mol}{\mathrm{s}}^2)\right]$
T	Operating temperature	$\left[\mathrm{K}\right]$
$c_b$	Bulk concentration of copper ions in the electrolyte	$\left[\mathrm{kg}/{\mathrm{m}}^3\right]$
$M_{Cu}$	Molar mass of copper	$\left[\mathrm{kg}/\mathrm{mol}\right]$
${\rho }_{Cu}$	Density of copper	$\left[\mathrm{kg}/{\mathrm{m}}^3\right]$
${\delta }_0$	Initial seed layer thickness	$\left[\mathrm{m}\right]$
${\sigma }_0$	Initial seed layer conductivity	$\left[{\mathrm{s}}^3{\mathrm{A}}^2/(\mathrm{kg}{\mathrm{m}}^3)\right]$
$\sigma$	Electric conductivity of solid copper	$\left[{\mathrm{s}}^3{\mathrm{A}}^2/(\mathrm{kg}{\mathrm{m}}^3)\right]$
$i_0$	Exchange current density	$\left[\mathrm{A}/{\mathrm{m}}^2\right]$
${\alpha }_{A;C}$	Charge transfer coefficients for the anode and cathode	$\left[-\right]$

.

2.2. Material and Fluid Properties

Table 1 and

Table 2. Material and Fluid fields that will be calculated through the course of each simulation. Symbols, descriptions, and dimensions are provided here.

Field	Description	Dimensions (SI)
$\mathbf{u}$	Fluid velocity	$\left[\mathrm{m}/\mathrm{s}\right]$
p	Fluid pressure	$\left[\mathrm{kg}/(\mathrm{m}{\mathrm{s}}^2)\right]$
$\mathbf{F}$	Body force acting on the fluid	$\left[\mathrm{m}/{\mathrm{s}}^2\right]$
C	Mass fraction of copper ion concentration	$\left[-\right]$
$\phi$	Electric potential of the electrolyte	$\left[\mathrm{kg}{\mathrm{m}}^2/(\mathrm{A}{\mathrm{s}}^3)\right]$
${\phi }_s$	Electric potential copper layer	$\left[\mathrm{kg}{\mathrm{m}}^2/(\mathrm{A}{\mathrm{s}}^3)\right]$
${\phi }_I$	Interface potential	$\left[\mathrm{kg}{\mathrm{m}}^2/(\mathrm{A}{\mathrm{s}}^3)\right]$
$\mathbf{N}$	Electric flux	$\left[\mathrm{A}\mathrm{m}/\mathrm{mol}\right]$
$\delta$	Deposited copper layer thickness	$\left[\mathrm{m}\right]$
i	Current density	$\left[\mathrm{A}/{\mathrm{m}}^2\right]$
${\sigma }_{eff}$	Effective conductivity to account for ${\delta }_0$ and $\delta$	$\left[{\mathrm{s}}^3{\mathrm{A}}^2/(\mathrm{kg}{\mathrm{m}}^2)\right]$
$\eta$	Overpotential	$\left[\mathrm{kg}{\mathrm{m}}^2/(\mathrm{A}{\mathrm{s}}^3)\right]$

list the fluid and material properties and the fields being computed throughout the ECD simulations, respectively.

3. Numerical Methods

3.1. Discretization

In this work, a cell-centered finite volume method (FVM) is used to discretize the bulk fluid (Equations (1)–(3)) and electric potential (Equation (5)) equations and a face-centered finite area method (FAM) is used to discretize the electric potential equation along the plating surface (Equation (6)). Both the FVM and FAM approaches are based on integral forms of the governing equations. Unless otherwise noted, all time derivatives are discretized using the second-order accurate, implicit scheme referred to as the backward scheme [16]. The following sections will briefly describe the FVM and FAM spatial discretizations.

3.1.1. Finite Volume Method (FVM)

In a finite volume approach, the computational space is divided into non-overlapping control volumes (CV) that completely fill the space. Figure 2a shows an example of a hexahedral control volume, $V_p$, centered around a point, P. The face f with surface area $S_f$ and unit normal ${\mathbf{n}}_f$ is shared by a neighboring CV with centroid N.

Integral forms of the governing equations are discretized in the FVM by transforming continuous surface integrals into discrete summations over CV faces. For example, the FVM form of the Laplacian term in Equation (16) is given by

$${\int }_V\nabla ·\left(\frac{F^2{n}^{2}D}{RT}\nabla \phi \right)dV={\int }_S\frac{F^2{n}^{2}D}{RT}\nabla \phi ·d\mathbf{S}=\sum _f\frac{F^2{n}^{2}D}{RT}{\mathbf{S}}_f·{(\nabla \phi )}_f.$$

(12)

The face normal gradient, ${\mathbf{S}}_f·{(\nabla \phi )}_f$, is evaluated according to

$${\mathbf{S}}_f·{(\nabla \phi )}_f=|S_f|\frac{{\phi }_N-{\phi }_P}{|\mathbf{d}|},$$

(13)

where $\mathbf{d}$ is the length vector between the center of cell P and neighboring cell N and ${\phi }_P$ and ${\phi }_N$ are the value of $\phi$ at the points P and N, respectively. In the case of non-orthogonal meshes, an additional correction term is introduced (see [16]). Similar discretizations can be introduced for the remaining terms in the governing equations. The interested reader is referred to [18] for additional information on the FVM discretizations used within the OpenFOAM^® framework [19]. For more detailed information on the FVM the interested reader is referred to [20].

3.1.2. Finite Area Method (FAM)

The finite area method can be thought of as a two-dimensional version of the FVM over curved surfaces, and is well-suited to problems for which the thickness of the region of interest is much less than the lateral dimensions of the domain and through-thickness gradients can be neglected. For example, in [17], the FAM is employed to solve the transport of surfactants along the surface of a free-rising bubble. In this work, a face-centered FAM is used to solve for the electric potential field within the thin deposited copper layer along the plating surface, see Equation (6).

In the FAM, similar to the FVM, the surface is decomposed into non-overlapping polyhedral control areas. Figure 2b provides an illustration of two adjacent quadrilateral control areas with centroids P and N, areas $S_P$ and $S_N$, and outward facing normals ${\mathbf{n}}_P$ and ${\mathbf{n}}_N$. Integral equations are then discretized as summations over the control area edges, analogous to the FVM procedure. For example, the FAM discretization of the Laplacian term in Equation (17) is expressed as

$${\int }_{\partial S}\nabla ·({\sigma }_{eff}\nabla {\phi }_s)=\sum _e{\sigma }_{eff}{L}_e{\mathbf{m}}_e·{(\nabla {\phi }_s)}_e,$$

(14)

where the subscript e denotes edge-centered values, $L_e$ is the length of the edge, and ${\mathbf{m}}_e$ is the outward pointing unit bi-normal of the edge. For more detailed information on the FAM the interested reader is referred to [21].

3.2. Solver Algorithm

Figure 3 describes the full solution algorithm that the ECD simulation follows. It should be noted that there is logic that allows the user to turn on or off the ionic transport of copper. The inputs into the algorithm are the initial concentration, seed layer thickness, and electric current. Here the overpotential is initially chosen to be zero and iterated to convergence via a Newton-Raphson [14] scheme.

3.3. Solver Parameters

The numerical solvers and schemes used in the ECD solver are noted in

Table 3. Numerical solvers used in the ECD simulations. LU = Lower-Upper; GAMG = Geomtric-Alegbraic Multi-Grid; FDIC = Faster Diagonal Incomplete-Cholesky; DICGaussSeidel = Diagonal Incomplete-Cholesky with Gauss Seidel; PBiCGStab = Stabilized Preconditioned Bi-Conjugate Gradient; DILU = Diagonal Incomplete-LU; DILUGaussSeidel = Diagonal Incomplete-LU with Gauss Seidel; PCG = Preconditioned Conjugate Gradient; DIC = Diagonal Imcomplete-Cholesky.

Finite Volume Method
Field	Equation	Solver	Preconditioner	Smoother
p	(1), (2)	GAMG	FDIC	DICGaussSeidel
$\mathbf{u}$	(1), (2)	PBiCGStab	DILU	DILUGaussSeidel
C	(3)	PBiCGStab	DILU	DILUGaussSeidel
$\phi$	(5)	PCG	DIC	-
Finite Area Method
${\phi }_s$	(6)	PCG	DIC	-
$\delta$	(11)	PBiCGStab	DILU	DILUGaussSeidel

and

Table 4. Numerical schemes to approximate the differential operators. For the divergence operators, $\nabla ·$ the choice is either Gauss upwind or Gauss linear, which depends on the specific form of the divergence.

Operator	Scheme	Description
$\partial /\partial t$	backward	Transient, $2\mathrm{nd}$ order, potentially unbounded, implicit
$\nabla (·)$	Gauss linear 1.0	Central differencing, bounded
$\nabla ·$	Gauss upwind/Gauss linear	$1\mathrm{st}$ order, bounded/ $2\mathrm{nd}$ order, unbounded
$\nabla (·)·\mathbf{n}$	corrected	Explicit non-orthogonal correction at cell faces, $2\mathrm{nd}$ order, conservative
${\nabla }^2$	Gauss linear corrected	$2\mathrm{nd}$ order, unbounded, non-orthogonal correction, conservative

. Numerical solvers are the linear algebraic solvers used to solve the discretized equations and numerical schemes are the methods used to approximate the differential operators in those equations. The relaxation factors for the equations that govern the fields $\phi$ and ${\phi }_s$ are set to $0.7$. These relaxation factors were chosen to reduce oscillations that occur in the electric potential fields during the loop that ensures the current density converges, but a more formal investigation into the optimal relaxation factors is necessary and is part of future work.

4. Validation Test Case

The geometric assumptions for this problem will limit the computational domain to two-dimensions (2D), as a simplification. Another simplification for the present state of the computational model is to have quiescent fluid motion. In other words there will be no mixing of the fluid, and therefore, it has no motion. The focus on this effort was to describe the coupling between ionic transport and electrochemistry. First, it is demonstrated that the approach agrees with previously published resultsin [1], but without ionic transport.

4.1. Description

The performance of the proposed model described in Section 2 and Section 3 is validated against the published results of [1] for copper electroplating of an L-shaped surface at the bottom of an electrolyte bath. In this simple test case, a constant electric current is supplied along the top surface of the bath and copper is deposited on the bottom L-shaped surface, highlighted in orange as shown in Figure 4. Due to the configuration of the plating surface, the thickness of the Cu layer is expected to be non-uniformly distributed along the plating surface, with the highest deposition rate occurring at the exposed corner of the L-shaped surface. The sampling line for used in the figures below starts in the center of the cyan line (this is the zero distance in the line plots below) then extends up and around the corner of the orange L-shaped plating surface and ends at the left most point on the orange plating surface.

4.2. Model Formulation and Assumptions

For this test case, the bulk fluid in the electrolyte bath is assumed to be quiescent ($\mathbf{u}=0$), and the electrolyte solution is sufficiently mixed such that ionic transport within the bulk fluid can be safely neglected. Under these assumptions, the governing equations simplify to

$$\nabla p=-\rho \mathbf{g},$$

(15)

$$\nabla ·\left(\frac{F^2{n}^{2}D}{RT}\nabla \phi \right)=0,$$

(16)

$$\nabla ·({\sigma }_{eff}\nabla {\phi }_s)=-i,$$

(17)

where $F^2{n}^{2}D/RT=\kappa$ when comparing to [1] which is representative of the conductivity of the electrolyte and the Butler-Volmer equation reduces to

$$i=i_o\left[exp\left(\frac{{\alpha }_{A}nF\eta }{RT}\right)-exp\left(\frac{-{\alpha }_{C}nF\eta }{RT}\right)\right].$$

(18)

Equation (15) expresses the hydrostatic equilibrium existing within the bulk fluid with the gravitational force vector defined as $\mathbf{g}=(0,0,g)$. Equations (16) and (17) describe the electric potential in the bulk fluid and deposited copper layer, respectively. The conductivity of the growing copper layer is given by

$${\sigma }_{eff}=\sigma \delta +{\sigma }_{s0},$$

(19)

where $\sigma$ is the electric conductivity of the deposited metal, $\delta$ is the thickness of the deposited layer, and ${\sigma }_{s0}$ is the electric conductivity of the initial deposited surface. This formulation allows for the initial seed layer and deposited metal to be of different materials with differing electric conductivities. Equation (18) is simply the Butler-Volmer Equation without the dependence on Cu²⁺ concentration. The solution algorithm is still given by Figure 3 and the expressions for interface potential and deposition rate, given by Equations (10) and (11), remain unchanged.

4.3. Domain and Computational Mesh

Referring back to Figure 4, the overall domain is a 1.0 m × 0.5 m × 0.1 m box with a 0.25 m × 0.25 m × 0.1 m section removed. A fixed, uniform current density boundary condition ($\nabla \phi ·\mathbf{n}=i/\kappa$) is applied along the top surface, colored in blue. Here, $\mathbf{n}=(n_x,n_y,n_z)$ is the boundary surface normal. The orange surface indicates the area over which electroplating occurs, and the bottom edge of the plated surface, colored in cyan, is grounded (${\phi }_s=0$). All other surfaces are considered to be perfectly insulated ($\nabla \phi ·\mathbf{n}=0$). The computational domain is discretized with uniform hexahedral cells with a cell size of 1.25 cm. The total number of cells is approximately 22,400.

4.4. Parameters

The results presented in [1] included two test cases, case-1 and case-2, to explore model parameter sensitivity. In this work, the proposed model is parameterized in accordance with values from case-1 only, as the results from case-2 do not exercise any additional physics and are qualitatively similar to the results of case-1. The relevant model parameters are given in

Table 5. Model parameters used in the L-shaped validation case [1]. Additional parameters are presented to account for the modifications to the model.

Parameter	Symbol	Value
Electric conductivity of electrolyte	$\kappa$	1.0 S/m
Electric conductivity of deposited metal	$\sigma$	5.95 × 10${}^7$ S/m
Surface electric conductivity of initial deposited surface	${\sigma }_{s0}$	59.5 S
Molar weight of deposited metal	m	63.546 g/mol
Valence of metallic ion in electrolyte	n	2
Mass density of deposited metal	$\rho$	8940 kg/m${}^3$
Temperature	T	300 K
Total current	J	5 A
Exchange current density	$i_0$	150 A/m${}^2$
Charge transfer coefficients	${\alpha }_A$ and ${\alpha }_C$	0.5, 0.5
Time step	$\Delta t$	10 s
Tolerance of potential relaxation subcycle	$ϵ$	1.0 × 10${}^{-6}$
Copper diffusivity	D	0.67 × 10${}^{-9}$ m${}^2$/s
Bulk copper concentration	$c_b$	6.3536 kg/m${}^3$

. In order to account for the difference between the model presented here and the model presented in [1], namely that $F^2{n}^{2}D/RT=\kappa$, the inclusion of the diffusivity and bulk copper concentration of the electrolyte are included in the model parameters.

4.5. Solution

Model predictions for the validation case described above are shown in Figure 5, Figure 6 and Figure 7. The model is run for a total time of 1500 s. Qualitatively a comparison between the results of the electric potential in the bulk electrolyte can be made between the predictions from [1] and to the predictions from the proposed model, the results of the electric potential from the proposed model are shown in Figure 5, a time $t=$ 1500 s. From this qualitative analysis the proposed model is in agreement with the results presented in [1].

Quantitatively comparing the time distributions of the current density and the copper layer thickness from [1] the proposed model time-dependent behavior of the interface current and deposition thickness distributions are provided in Figure 6a,b. The time-dependent behaviors correspond well with the results presented in [1].

Figure 7a,b show distributions of interface current density and deposition layer thickness on the deposition surface, also at time $t=$ 1500 s. A qualitative comparison between the plots produced by this model and those presented in [1] again demonstrate good agreement.

4.6. Analysis and Discussion

Model predictions of electric potential in the bulk electrolyte, shown in Figure 5, agree qualitatively with the published results of [1]. Unfortunately, a more quantitative comparison of the electrolyte potential is not possible; however, based on the range of scales published in [1], maximum and minimum values between the two model predictions for electrolyte potential (interface potential) agree to within 0.8% (6%) and 5.5% (187%), respectively. It is important to note that the large percent difference between the minimum value of interface potential between the two models is relative to a value $<<1$, and that small differences in meshing or solver parameters can lead to small numerical differences that will accumulate over time.

Figure 6a,b shows the interface current distribution and deposition thickness evolution over time, respectively, and provide a quantitative comparison between the two model predictions. Both models predict an interface current distribution that is constant over time, with the highest values occurring at the exposed corner of the plating surface, around a distance of 0.25 m. The proposed model predicts a peak current density of approximately 251 A/m${}^2$, while the model from [1] predicts a value of 250 A/m${}^2$.

Qualitative comparisons are shown in Figure 7a,b for the current density and Cu layer thickness distributions along the plating surface. Maximum and minimum values for both interface current density and Cu thickness agree with the results in [1] to within 0.2%.

Qualitatively similar results are present in the prediction of deposition thickness, as shown in Figure 6b. Both models predict increasing Cu layer thicknesses over time, with the thickest deposition occurring at the corner of the plating surface. These results are consistent with the deposition rate’s linear dependence on interface current density (see Equation (11)).

The results of this code-to-code comparison provide confidence that the proposed approach can successfully model the electrochemical deposition process in cases where fluid motion is quiescent and ionic transport can be neglected. The performance of the proposed model which includes ionic transport will be discussed in the following section.

5. Ionic Transport Case

An extension case from Section 4 is presented to demonstrate the coupling of copper ionic transport to the electric potential, current density, and copper layer growth. This case uses the same problem parameters presented in Table 5. The difference here is that the branch shown in Figure 3 for ionic transport is turned on so the transport of copper ions, and depletion near the plating surface, impacts the current density and electric potential fields near the electrolyte-copper layer interface. This depletion of copper ions at the plating surface also indicates that the electroltye is no longer well-mixed.

Analysis and Discussion

Figure 8 shows the electric potential field contours at $t=$ 1500 s, when copper ionic transport is included in the physical system. From this, it can be observed that there are slight changes in the surface electric potential at the corner of the plating surface, and the bulk electric potential is decreased.

This change in electric potential at the interface of the plating surface leads to a decrease of the current density over time, but overall the current density at this interface is larger (because of the copper ion concentration coupling), see Figure 9a. The main change in current density occurs where the majority of copper is deposited (see Figure 9b), which is at the corner of the plating surface.

Figure 9a,b show the time dependent evolution of the current density and copper layer growth, respectively. It can be observed from these plots that the peak value of the current density, at the corner of the plating surface, decreases because of the increased resistivity due the growth of the copper layer. It is much less obvious, but the rate at which the copper layer is growing is starting to decrease, meaning that less copper is being deposited in the same amount of time. This occurs because there is less copper ions in the vicinity and because the current density has decreased. The overall current density is greater than in the case presented in Section 4, but it can be observed in Figure 10a that the current density has decreased as the simulation progressed. Figure 10b shows a similar trend to that presented in Section 4.5.

Additionally, Figure 11a,b shows the concentration surface contour along the plating surface at $t=1500$ s and the copper ion deposition at the plating surface interface as time evolves. It can be observed in Figure 11b, unsurprisingly, that the regions of the plating surface that correspond to large deposition result in less copper concentration near the interface of the electrolyte-copper layer. It can also be observed that the growth on the corner of the plating surface results in a shadowing effect along the vertical portion of the plating surface which can be observed in the left side of Figure 9a,b, and Figure 11b, based on the start of the sampling line.

6. Summary

This paper details the development, testing, and application of a multi-physics simulation code capable of predicting the spatio-temporal deposition of copper during an electrochemical deposition process. The details of the underlying mathematical models, several implementation details, and modeling considerations have been presented. The results of a validation study are also included, in which the proposed model is used to predict copper deposition along an L-shaped plating surface. The proposed model results presented for validation are in agreement with the results presented in [1], but the additional assumption that the electrolyte remained well mixed was too restrictive. To address this, the proposed model includes the addition of copper ion transport, which demonstrates how deposition changes as the electrolyte is depleted of copper ions near the plating surface. An extension of the validation problem, which incorporates ionic transport, has also been studied. It can be observed that the depletion of copper ions results in a larger resistivity on the plating surface (the deposited copper layer), which, in turn, decreases the current density on the plating surface, specifically at the corner (where the copper concentration is being depleted). This also slows the rate at which copper is deposited as time evolves.

As mentioned in Section 1 the overall objective is to study how turbulent mixing affects the deposition of copper, and is future work to build upon this initial effort. Mixing the electrolyte will ensure that copper ions are not depleted at the plating surface-electrolyte interface. Mixing is also believed to help with the uniformity of deposition. Incorporating either a forcing function within the fluid, or a moving mesh approach, such as the ALE method, introduces additional challenges in the numerical method.