Numerical Modeling of Electrodeposition Phenomena
Nadia Strusevich, Michael Hughes, Christopher Bailey, Georgi Djambazov
University of Greenwich
Park Row, Greenwich, London SE10 9LS
{n.stroussevitch,m.s.hughes,c.bailey, g.djambazov}@gre.ac.uk
Abstract
Electrodeposition is a widely used technique for the fabrication of high aspect ratio microstructure components.In recent years much research has been focused within this area with an aim tounderstanding the physics behind the filling of high-aspect ratio vias and trenches on PCB’s and in particular how they can be made without the formation of voids in the deposited material. This paper describes some of the fundamental work towardsthe advancement of numerical models that can predict the electrodeposition process and addresses:
i)A novel technique for interface motion based ona variation of a donor-acceptor technique
ii)A methodology for the investigation of stress profiles in deposits
iii)The implementation of acoustic forces to generate replenishing electrolytic flow circulation in recessed features
Introduction
Electrodepositionis a complex process and truly multi-physics in its nature, coupling together the solutions of several physical phenomena such as fluid-flow, temperature, ionic concentrations, electric current and chemical reaction rates. In recent years the drive to produce high quality, tiny yet intricate components has lead research effort in the deposition of micro scale high aspect ratiofeatures.
Key to the deposition quality in these recessed features is the ability to replenish the ionic species. Deposition rates are otherwise limited to the speed of diffusion processes because they are flow dead zones andsubsequently there is a tendency for void formation in the bottom of the trenches as the concentration of ionic species depletes.Towards this aim there has been much interest in using either chemical additives that can produce a ‘superconformal’ accelerationof the deposition from the bottom of the trench, known as CEAC, see [1]; or more recently high frequency sound waves that generate electrolytic flow deep into these recessed features andreplenish the supply of ionic speciesto produce a faster, symmetrical void-free filling of the trench, see [2].
This paper describes some of the fundamental work concerning the development ofnumerical modelsforthe electrodeposition process of copper. This work has been undertakenin collaboration with Herriot-WattUniversity, as part of the MEMSA project in order to explore such scenarios as the deposition of high-aspect ratio features and the use of megasonic agitation to improve the deposition quality.
Key to this research is the development of the underlying numerical models of deposition; this in particular involves the accurate representation of the moving metallic/electrolyte interface. Previous work by Hughes et al. [3], has utilized a level-set technique to represent this interface motion, however in this paperanalternative and more accurate approach has been developedthat is based on a novel variation of the donor-acceptor technique. Results are presented for two types of interface deposition kinetics; driven respectively by Ohm's and Butler-Volmer laws.
Additionally a finite element based approach to examine the intrinsic and residual stresses in the electrodeposits is described; as the strength and durability of these deposits is of crucial importance to the operational expectations of the derived devices. The method described is generic and can be adapted to examine such affects as variation of grain size on stress profiles.
Finally, the effect of acoustic agitation on the electrolytic flow within recessed,high aspect trenches has been considered; as a pre-cursor to their implementation in the full electrodeposition model. CFD simulations have been computed over trenches of varying aspect ratio to examine the influence of the acoustic streaming on the electrolytic transport.
All software development and numerical analysis has been implemented within the unstructured, control volume CFD code PHYSICA [4] and results from the modeling simulations are discussed in the subsequent sections.
Modelling of Electrodeposition: Primary Current Distribution
The formula for the deposition rate that evaluates how much metal is deposited during each time stepis described by the equation below, where v has units(m/s)
Here isthe molar volume of copper (equalto 7.1·10-6m3/mol), i is the current density (A/m2), zis the charge number (equal to 2 for copper),F is Faraday’s constant (equal to 96485.33 Coul/mol), and is the normal to the electrolyte/metal interface pointing into the electrolyte.
Dependent on the electrodeposition regime to be simulated, we distinguish the primary current distribution, governed by Ohm’s law, and the tertiary current distribution, governed by the Butler-Volmer law.
The following assumptions are made: (i) a single cathodic reaction takes place at the cathode; (ii) the physical properties of the electrolyte remain constant; (iii) the activation overpotential at the cathode is constant; (iv) the anode does not change shape during electrodeposition and the cathode is of primary interest, (v) only basic electrodeposition is considered with no use of electrolytic additives or accelerators.
This section discusses the numerical modelling of the electrodeposition process of copper, provided that the current distribution follows Ohm’s Law, also known as the primary current distribution. This distribution is normally applicable if the resistance of the electrolyte is much higher than the resistance of the interface,and takes into account the migratory properties of ions in electrolyte and does not consider diffusion mass transport. The current density in this case can be written as
where σelis the electrolyte’s electric conductivity (A2s3/kg/m3 ) and is the electric potential (V).
We start with the simplest situation, in which copper is deposited on the horizontal surface.This geometry restricts the growth of the deposited metal to 1-D. The governing equation for the potential field in the electrolyte obeys the Laplaceequation. In this case the potential depends linearly on the verticalcoordinate y of the domain, y[h, L], where h is the height of the electrolyte/metal interface that changes as the deposition process develops, and Lis the overall height of the domain. We assume that (L)is kept equal to 200V, while (h) = 0. Thus, as the interface moves higher, the gradient of the potential increases.
The numerical solution is obtained by PHYSICA by applying a variation of the donor-acceptor method. The domain is covered of by a mesh of rectangular elements. In each time step, we determine the layer of elements that contain the interface, and update the domain to become the electrolyte above the new location of the interface. For the interface layer PHYSICA finds the potential, and based on that, the current density i and the deposition rate v are computed to determine a new position of the interface.
Figure 1: Analytically and numerically computed deposition rate for the primary current distribution
The values of the analytically computed and numerically obtained deposition rates are plotted in Figure 1. The graph for the numerical results has a staircase-like shape, different from the linear shape of the analytical solution. This is due to the fact that our numerical method does not update the deposition rate while the interface remains inside the same layer of elements. When such a transition to a new layer takes place, both graphs coincide.Notice that the deposition rate grows as the gradient of the potential grows over time. On the other hand, the graphs for the interface position (analytically and numerically found) are virtually indistinguishable, and therefore are not presented here.
Now we turn to a two-dimensional situation, in which copper is deposited on a surface with a trench. The 2-D geometry is presented in Figure 2; due to symmetry, only the left half is shown.
Figure 2: Two-dimensional geometry (the left half) with the governing equation and the boundary conditions for theprimary current distribution. The interface is marked with a thick line
The governing equation for the potential field in the electrolyte obeys the Laplace equation. The aspect ratio r of the trench is defined as b/(2a).
Figure 3: (a) the interface position after 3.5 min, (b) deposition thickness in the middle part of the side wall (the top line) and at the bottom of the trench (the lower line)
The implementation of the numerical solution by PHYSICA in the 2-D case is more complicated than in the 1-D case. A special code is written for determining the position of the interface in each time step. For each element that contains the interface in the beginning of a time step, the volume of deposited metal by the end of this step is computed. If that volume exceeds the remaining capacity of the element, the excess is recursively redistributed to its neighbouring elements in proportion to the difference of potentials at the elements, until the whole excess is redistributed. Numerical results of electrodeposition for the cell with aspect ratio 1:1, a = 75μm, b = 150μm are shown in Figure 3.
As can be in Figure 3, the deposited metal layer grows faster on the side wall than at the bottom of the trench, thereby leaving a void. A similar behaviour is observed for other aspect ratios, e.g., 2:1 and 3:1.
Modelling of Electrodeposition: Butler-Volmer Current Distribution
The tertiary current distribution model takes into consideration ion transport due to diffusion described by Fick’s second law. Concentration is the variable factor in the diffusion layer and the kinetic expression at the electrode has to take into account the concentration of ions at the surface. For the tertiary current distribution, the current density i is given by the Butler-Volmer equation
where i0is theexchange current density (A/m2), cint and c∞ are the molar concentration of copper at the interface and in the far field (mol/m3), α the transfer coefficient, Ris the gas constant (equal to 8.314 J/mol/ºK), T is the temperature (ºK) and η is the overpotential (V). Unlike the primary current distribution, the current density i is depends on the metal concentration and the overpotential. The governing equation and the boundary conditions are shown in Figure 4.
Figure 4: 2-D geometry with the governing equation and the boundary conditions for the tertiary current distribution. The interface is marked with a thick line.
We performed numerical modelling of electrodeposition for three values of the aspect ratio of the trench equal to 1:1, 2:1 and 3:1. Recall that the aspect ratio ris defined by a/(2b);see Figure 2. In all cases a=75μm, L-b = 500μm. The initial concentration was taken to be 350 mol/m3, α = 0.5, η = -0.3V, i0 = 26 A/m2and the temperature T = 293ºK. The simulation was run until either the trench is filled or the copper ions concentration is depleted.
Figure 5: Results for r=1:1
For each aspect ratio, the numerical results after the simulation was completed are shown in Figures5-7, respectively, where part (a) shows the deposition level and part (b) shows the concentration distribution. The simulation time for r=1:1 was 22 min and, as shown in Figure 5, the trench is completely filled, while the minimum concentration in the neighbourhood of the interface was still fairly large.
Figure 6: Results for r=2:1
Figure 7: Results for r=3:1
The simulation time for r=2:1 was 25 min and, as shown in Figure 6, a void is left in the trench; the minimum concentration in the void is depleted, i.e., too small to get rid of the void.The simulation time for r=3:1 was 23 min and, as shown in Figure 7, the metal is deposited as a layer along the trench, and the minimum concentration is depleted, i.e., too small to get rid of the void.
Prediction of Intrinsic Stress Distribution by Finite Element Analysis
The presence of residual or intrinsic stress in the electrodeposited films is the important factor to consider in devices because such a stress causes the deposited layer depending on its thickness and elastic properties to expand or contract either to a concave or to a convex shape. In some cases intrinsic stress exceeds the strength of the deposited layer, resulting in film cracking, deformation of devices, and interfacial failure. Even moderate tensile or compressive stresses in films may lead to geometric distortions, warpingor shrinkage.Thus, the intrinsic stress in films deposited on substrates should be made as small as possible by controlling the deposition conditions of a film.
Many factors affect internal stress in electrodeposited films, including film composition, nature of the substrate surface and the deposit, current density, the deposit thickness, among others.It is essential to determine effects that process variables may have on stress.This will allow us to find such values of the critical process variables that the deposition process is optimized with respect to a chosen objective, e.g., a desirable stress profileor a high deposition rate.
We are not aware of any established mathematical models that describe the influence of process variables on intrinsic stress. Most of the studies in this area are based on experimental measurements. For example, in [5] Kim et al. study the behaviour of the intrinsic stress of copper film electrodeposited on the Au seed. The authors measure the substrate’s radius of curvature before and after film deposition and compute the stress by Stoney’s equation.They demonstrate how the value of stress depends on the film thickness and on the applied current density.
Engelstad et al. [6] point out that Stoney’s equation can only be reliably used under the assumption on the uniform film stress. To overcome this drawback, they make stress predictions by the Finite Element analysis, taking interpolated measured curvatures as input and deliver intrinsic stress as a function of location and direction.
In this work, we rely on the data presented by Kim et al. [5] for the average values for tensile intrinsic stress obtained for 50-, 100-, and 150-μm thick copper films deposited on an Au seed and with the current density varying within the interval 25mA/cm2-65mA/cm2. For the purpose of predicting the intrinsic stress values we use the analysis software package PHYSICA. The nature of the intrinsic stress does not allow us to compute it directly. Instead, we use the finite element approach to the numerical modelling of another stress-generating process and then convert the obtained solution into the desired values on the intrinsic stress. In particular, we select heat transfer as this alternative process and take either the thermal expansion coefficient or the temperature change as the process parameter. We establish such values of the chosen parameter that generate the values of tensile stress close to the measured intrinsic stress values reported in [5]. Our approach is close to the one used by Chen and Ou [7], who express intrinsic strains in film and substrate in terms of the equivalent thermal expansion coefficients and the equivalent temperature difference.
Having conducted numerical experiments with PHYSICA, we have observed that the average values of the process parameter can be accurately approximated by a function of two variables, the film thickness h and current density i, given by
,
where f1(h)= 0.0279h2- 2.3202h+ 128.45 and f2(h)= 0.0279h2 - 2.3202h+ 128.4; see Figure 8.
Figure 8: Process parameter as
a function of two variables
In order to obtain the values of intrinsic stress reported in [5] for a given pair of values of h and i, we fix the process parameter to be equal to run the heat transfer simulation and output the computed average stress value.
The average value of the process parameter can be expressed as
,
whereis an appropriate distribution function. Differentiating, wefindthis process parameter distribution function for a fixed current densityas
.
Substituting the value of the parameter into this formula we derive the expression for the distribution function
This function is then embedded into the PHYSICA code, which allows us to view the process parameter distribution; e.g., see Figure 9. The figure shows one half of the geometry symmetric with respect to the vertical axis.
Figure 9. Process parameter distribution for the 100-μm thick copper film deposited at i = 35mA/cm2
In turn, based on this distribution function we are able to simulate numerically the values of intrinsic stress in the elements of the deposited copper; see Figures 10 and 11.
Figure 10. Intrinsic stress distribution in copper film deposited on Au seed for h = 100μm:
(a) i = 35mA/cm2; (b) i = 55mA/cm2
Figure 11. Intrinsic stress distribution in copper film deposited on Au seed for h = 150μm:
(a) i = 35mA/cm2; (b) i = 55mA/cm2
As can be seen in Figures10 and 11, the increase in current density during electrodeposition induces larger stress and also that stress of thick copper films (over 50μm) increases with further growth of the deposited copper.The average values of internal stress in all examples shown in Figures 10 and 11 are fairly closed to the data reported by Kim et al. [5].
Acoustic Streaming
Megasonic agitation provides a promising means to solve the problems concerning void formation in high aspect ratio trenches and vias by enhancing a replenishing ionic transport that would otherwise not penetratedeep into these features. The effect of this is to decrease the Nernstdiffusion layer and therefore increase the kinetic rate at the interface. Experiments at Heriot-Watt by Kaufmann et al. [8],have shown an improvement in both the deposition rates and deposit quality for microvias with features of aspect-ratio upto 2:1. being successfully filled; the limitation of trenches with ratios larger than 2:1 being reported as due to the lack of conductive seed layer deep into the microcavity.
The Acoustic Streaming process is generated from the high frequency agitation (~ 1MHz) of a Piezoelectric plate generatingan acoustic intensity of 60kW/m2 with correspondinghigh frequency pressure variations within the feature. These pressure variations generate Reynolds stresses that drive a steady fluid recirculation within the cavity that replenishes the supply of the reacting ionic species from the bulk of the containing bath.
The numerical solution process is in two stages;the first concerns the analytical solution of the time dependent wave equations as developed by Rayleigh and Nyborg for channels between infinite parallel plates and is detailed concisely by Nilson and Griffiths[2]. The second stage involves calculating the time averaged Reynolds stresses that are produced from theseharmonic wave solutions and their numerical implementation into the Navier-Stokes equations.
Simulations have been made to examine the effects of acoustic streaming in extremely narrow features withwidths in the order of tens of microns so that the resulting flow profiles can be compared with those produced by Nilson and Griffiths [2].