Qualitative Simulation of the Growth of Electrolessly Deposited Cu Thin Films

Electroless deposition for fabricating copper Cu interconnects of integrated circuits has drawn attention due to its low processing temperature, high deposition selectivity, and high coverage. In this paper, three-dimensional computer simulations of the qualitative growth properties of Cu particles and two-dimensional simulations of the trench-filling properties are conducted. The mathematical model employed in the study is a reaction-diffusion equation. An implicit finite difference discretization with a red-black Gauss-Seidel method as a solver is proposed for solving the reaction-diffusion equation. The simulated deposition properties agree with those observed in experimentation. Alternatives to improve the deposition properties are also discussed.


Introduction
Copper is widely used as the interconnecting metal of integrated circuits because of its low resistivity and high resistance to electromigration and stress voiding.Electroless copper deposition has been shown to be a viable technology due to its low processing temperature, deposition selectivity, and high conformity.Generally, the electroless copper deposition solution contains a copper compound, ethylenediaminetetraacetic acid EDTA ligand as a complexing agent for copper, formaldehyde HCHO as a reducing agent capable of reducing the copper compound to metallic copper, and additives as surfactant and stabilizer.Copper is then deposited from the solution onto catalytic seeds on the dielectric layers such as silicon dioxide.The overall reaction on the copper surface is with the overall reaction rate where the constants K, a, b, c, and d have to be determined experimentally.The electroless deposition occurs in stationary solution.Therefore, the convection and drift terms can be ignored.
Under some process conditions, Schacham-Diamand et al. 1 have shown that Cu deposition rate is primarily determined by the concentration of the dissolved Cu II ions.Based on these conditions, Smy et al. 2 proposed a two-dimensional reaction-diffusion equation for simulating the Cu deposition for trench-filling.The reaction-diffusion equation is given as where C is the ionic copper concentration and D is the diffusion coefficient.The boundary conditions for 1.3 are i C C 0 well away from the reaction surface, and ii D dC/dn FC r where F and r are related to the rate constant and reaction order, respectively.Smy et al. 2 used the thin film growth simulator package SIMBAD to simulate the trench-filling process for fabrication of very large scale integration VLSI interconnections.A sequence of grid generation followed by solving the quasi-steady-state of 1.3 was performed.Their numerical results compared well with the experimental films over a range of topography with aspect ratios varying from 1 : 1 to 4 : 1.This paper is concerned with the computer simulation of the qualitative properties of the growth of Cu deposits.Equation 1.3 with the corresponding boundary conditions is extended to a three-dimensional problem.In order to simulate the evolution of the growth of Cu particles, both spatial and temporal discretization are applied.An implicit finite difference scheme with red-black Gauss-Seidel iterative method as a solver is employed in this paper.This numerical method has been successfully applied to reaction diffusion systems of two species with nonlinear reaction terms 3 .It was proven to be unconditionally stable and convergent under some conditions, and the numerical results showed that this numerical method was efficient.In this study, the residual of the linear system arising from discretization is computed to confirm the convergence of the numerical approximations.
The feature size in recent integrated circuits has already reached sub-100 nm scale.Moreover, electroless deposition has been recently introduced in fabrication of nanostructured barrier layers against Cu diffusion.In this application, ultrathin films with thicknesses of tens of nm are required 4, 5 .Although electroless deposition is a promising technique for fabricating Cu-metallization thin films, some problems arise such as low plating rate and void formation.Understanding the mechanism of the film growth helps to minimize the thickness and improve the quality of thin films.Some factors that cause these problems and methods for solving them are studied by numerical simulation.The numerical results are compared with experimental results.
The outline of this paper is as follows.The numerical method is detailed in Section 2. Examples are given with discussion in Section 3. Section 4 provides a brief conclusion.

Numerical Method
Suppose the deposition solution is prepared in a container of dimension 0, α × 0, β × 0, d without agitation, and the substrate is placed horizontally on the bottom of the solution as shown in Figure 1 a .The copper films are required to have a thickness of tens of nm.The reaction takes place only on the surfaces of the Cu particles grown on the substrate.Because the depth of the solution is usually larger than 1 cm, that is, γ > 1 cm in Figure 1 a , and the concentration of the solution is uniform as it is prepared, the computation focuses on a part of the solution below a level c with an appropriate boundary condition imposed to the problem, as proposed in 2 .
Consider a rectangular region E To describe and track the moving boundaries, Γ 1i t ⊂ E, i 1, 2, . . ., N, a level set method 6 is used to formulate the motion of these boundaries.Let v t, x, y, z defined on R × E be the so-called level set function: v t, x, y, z v * for x, y, z inside the interface Γ 1i t , 1 ≤ i ≤ N, and v t, x, y, z 0 otherwise.Here, v * is a prescribed value related to the density of Cu so that the deposition rate is Fu r /v * at any point on the surface 2 .The motion of the interfaces The level set equation 6 is This is equivalent to Consider a uniform rectangular mesh on E with mesh size The increment in t is denoted by Δt and t n nΔt for n ≥ 0. The approximation of u t n , x i , y j , z k is denoted by the standard notation u n ijk .Equation 2.2 is discretized using the backward-time central-space scheme 7 .The backward-time discretization for time derivatives is given by u t t n , x i , y j , z k ≈ u n ijk − u n−1 ijk /Δt, and the central-space discretization for second-order spatial derivatives is given by

2.8
Equation 2.3 is discretized by the central difference method for first derivatives.For example, Therefore, 2.3 gives u n iN 2 −1k .

2.10
Equation 2.4 is also discretized by the central difference method analogous to 2.9 .Let x i , y j , z k ∈ Γ 1p t n with 1 ≤ p ≤ N. Suppose that all of its neighboring grid points are in Ω t n except for x i−1 , y j , z k , which is inside the Cu particle p at t t n .The reaction term in 2.4 is estimated by the backward time approximation Fu r ≈ F u n−1 ijk r .Equation 2.4 is then discretized by

2.12
Suppose two of the neighboring grid points, say x i−1 , y j , z k and x i , y j−1 , z k , are located inside the Cu particle p at t t n .Let Similar approximation is made for u n i 1/2 j 1/2 k .Equation 2.4 for this case is then approximated by

2.14
If three of the neighboring grid points, say x i−1 , y j , z k , x i , y j−1 , z k , and x i , y j , z k−1 , are located inside the Cu particle p.The approximation is analogous

2.15
Hence, the discretization of 2.1 -2.5 is given as follows: and n > 0 with the substitution of 2.10 -2.15 for the grid points outside Ω t .The weight Since the Cu deposits also depend on the size of the surface, the weight w n−1 ijk is determined not only by the coefficients of the reaction terms in 2.12 -2.15 but also the location of x i , y j , z k , which will be explained later.

Mathematical Problems in Engineering
Let v n ijk be the approximation of v t n , x i , y j , z k .From the definition of the level set function ijk , where x i , y j , z k is on ∪ N l 1 Γ 1l t n , has to be computed.For a boundary point x i , y j , z k as described in 2.12 , we have v n−1 i−1jk v * and v n−1 i 1jk 0. Equation 2.7 is discretized using the backwark-time central-space scheme 2.17 For a boundary point as described in 2.14 , we have

2.18
Similarly, the discretization for a boundary point as described in 2.15 is

2.19
From 2.17 -2.19 , the discretization of 2.7 can be writen as

2.20
Again, the weight μ n−1 ijk is determined not only by the coefficients in 2.17 -2.19 but also the location of x i , y j , z k .Rewrite 2.16 as

2.21
Now, w n−1 ijk F u n−1 ijk r /h Δt in 2.21 estimates the Cu ions consumed at x i , y j , z k over t n−1 , t n in the deposition process, while μ n−1 ijk F u n−1 ijk r /h Δt in 2.20 estimates the Cu deposits at x i , y j , z k over t n−1 , t n .In the simulation, we set w n−1 ijk μ n−1 ijk for a boundary grid point x i , y j , z k so that the Cu ions consumed equal the Cu deposits.Let Γ n in,p be the approximation of Cu particle p, Γ n p the approximation of Γ 1p t n , and Ω n the approximation of Ω t n .Γ n p is formed by the grid points satisfying and at least one of x i±ii , y j±jj , z k±jj lies in Γ n in,p , where ii 0, ±1, jj 0, ±1, and kk 0, ±1.
At the beginning of the deposition, let Γ 0 in,p , 1 ≤ p ≤ N, be given for the deposition of catalytic seeds.Then, the set of boundary points ∪ N l 1 Γ 0 p and Ω 0 are determined by the above definition.The initial value v 0 ijk is given as above for n 0 except for v 0 ijk 0 , and add its neighboring grid points which are in Ω n−1 to the set ∪ N l 1 Γ n l .Here, the neighboring points of x i , y j , z k are x i ii , y j jj , z k kk with ii, jj, kk 0, ±1, but are not all zero.In this paper, we assume the rate of Cu deposits at a boundary grid point also depends on the its distance from a nearest grid point in ∪ N l 1 Γ n in,l .For example, suppose x i , y j , z k ∈ Γ n−1 p and x i , y j , z k ∈ Γ n in,p .If x i 1 , y j , z k , x i , y j 1 , z k , and x i 1 , y j 1 , z k are added to Γ n p , then it is assumed that the rate of Cu deposits at x i 1 , y j 1 , z k is smaller than that at x i 1 , y j , z k and x i , y j 1 , z k .Thus, μ n ijk may have three different vales.In this paper, μ n ijk 1, 1/2, or, 1/3 is used in the numerical simulation.A larger μ n ijk value is corresponding to a smaller distance from x i , y j , z k ∈ ∪ N l 1 Γ n l to its nearest grid point in ∪ N l 1 Γ n in,l .The computation procedure is summarized as follows: given the initial conditions Γ 0 in,p for 1 ≤ p ≤ N, and u 0 ijk for 0 ≤ i ≤ N 1 , 0 ≤ j ≤ N 2 , and 0 ≤ k ≤ N 3 , then Γ 0 p , Ω 0 , and v 0 ijk can be determined.At each time step t n , v n ijk is computed by 2.20 at each boundary grid point in ∪ N l 1 Γ n−1 l .Then, v n ijk is used to update the set ∪ N l 1 Γ n in,l , the boundary set ∪ N l 1 Γ n l , and the domain Ω n .It is followed by the computation of u n ijk using 2.21 at the points which are not in ∪ N l 1 Γ n in,l .

Numerical Simulations and Discussion
Since the Cu deposition from an aqueous solution takes place during the plating time, the growth of Cu particles is computed at each time step over the plating time.To simulate a particle size of tens of nm, the mesh step h defined in Section 2 has to be set at about 1 nm and the diffusion coefficient D 10 −5 cm 2 /sec 10 9 nm 2 /sec.Thus, h 2 / DΔt 10 −9 /Δt is very close to zero unless Δt is small enough, for example, Δt < 10 −8 .When h 2 / DΔt is close to zero, the linear system defined by 2.21 may be nearly singular by the Gerschgorin Circle Theorem 8 .On the other hand, using such a small Δt value results in a running time that is too large to be practical.In this paper, qualitative simulations are performed by choosing the proper parameter values in 2.2 -2.4 so that the simulated growth properties agree with those in experimentation and the running time is short.The linear system 2.21 is solved by red-black Gauss-Seidel iterative method 9 .The maximum norm of the residual after thirty iterations at each time step is less than 5 × 10 −9 for all numerical examples presented in this section.Thus, the convergence of the iterative method is assumed.

Simulation of the Growth of Cu Particles
Consider the rectangular domain R 0, 60 × 0, 60 × 0, 44 .Let the parameter values be

3.1
Suppose that four catalytic seeds are deposited on the substrate.Figures 2 a -2 c show the sizes of particles after plating time t 2, 10, and 20, respectively.The in-plane and out-plane growth are plotted in Figure 2 d .At first, conformal deposition is observed.As Cu particles grow larger, the space on the bottom of the film between particles becomes smaller.This situation makes it difficult for ions to diffuse from upper part to the lower part of the particles in the deposition solution.When 70%-80% of the bottom region is covered by the particles, the in-plane growth starts to decrease dramatically.As the deposition time is prolonged, the difference between in-plane growth and out-plane growth becomes significant.At this point, the concentration on the bottom of the film is close to zero and the concentration at the top of the Cu particles is much higher.The out-plane growth remains at a constant rate while the in-plane growth rate is close to zero.The simulated growth properties agree with experimental observations of electroless deposition of Ag films and Co-P barrier layers 10, 11 .The thickness of the film increases due to this nonconformal growth.Moreover, Tong and Wang 10 pointed out that the nonconformal growth may produce voids on the bottom of the film between particles.Let the rate constant F increase to 60. Figure 3 a shows that the nonconformal growth is more significant due to the fast consumption of Cu ions and difficulty of diffusion.Copper ions are deposited on the surface of a Cu particle as they diffuse from the top to the bottom of the particle.Thus, two particles connect to each other at the middle and a void forms at the bottom as shown in Figure 3 b .In contrast, let the rate constant F decrease to 15.The difference between in-plane growth and out-plane growth becomes smaller as shown in Figures 4 a and 4 b .The thickness of the deposits decreases by 17% compared with the case of F 40.
The nonconformal growth occurs when a particle grows faster at the top than at the bottom.The combination of accelerating and inhibiting additives has been proposed to obtain bottom-up growth of the deposits 12 .Thus, the reaction rate can be assumed to be a function of z.Define

3.2
The reaction rate is 2 F at the bottom, F for z ≥ d, and decreases linearly from the bottom to z d.Consider F 20 and d 15.The in-plane growth is faster than out-plane growth as shown in Figures 4 c and 4 d .The thickness of the deposition is 31% less than the case of F 40. Another factor that affects the thickness of the deposits is the population density of the catalytic seeds.Using the same computational domain as above with parameter values in 3.1 , let nine catalytic seeds be deposited.the fabrication of nanostructured barrier layers.This seeding process significantly increased the population density of catalytic seeds, and barrier layers of thicknesses of 10 nm were fabricated.
Further simulations show that the nonconformal growth and voids may also occur as the diffusivity D decreases.For example, let D 300, F 15, and all other parameter values are set as in 3.1 .The growth of the particles is similar to Figure 3, where D 1000 and F 60, with nonconformal growth and void formation.Detailed discussion of the growth problem caused by low diffusivity will not be given here since it is a repeat of the discussion given above.

Simulation of Trench Filling
Electroless copper deposition is a promising method for the fabrication of copper interconnections.In this subsection, simulations of trench filling are performed by directly applying the numerical method in Section 2 to a two-dimensional problem and assuming the seed layer is uniformly deposited on the substrate as shown in Figure 5 a .Consider the rectangular domain R 0, 40 × 0, 40 kw for a k : 1 aspect-ratio feature, and let the parameter values be set as in 3.1 with w 20. Figure 5 b shows that voids occur for the 1 : 1 aspect-ratio feature even at a low reaction rate F 6. To achieve void-free trench filling, the deposition rate has to be very small.Figure 5 c shows F 3 for void-free filling.As aspect ratio increases from 1 : 1 to 2 : 1, F 3 fails to be void-free filling as shown in Figure 6 a .Void-free filling can be achieved if F 0.5 as shown in Figure 6 b .However, the plating time is large when the deposition rate is small.This problem becomes more severe as the aspect ratio increases.Table 1 shows the reaction rate and plating time in void-free trench filling for the k : 1 aspect-ratio feature, k 1, 2, 3, and 4. The simulated trench-filling property agrees with the bottom-up growth of copper by electroless deposition studied in recent years 13 .Inhibiting additives in plating solution have demonstrated successful void-free filling in high aspect-ratio features.However, the resulting deposition rate is too small.Void formation in trench filling is caused by the larger deposition rate at the trench openings than at the trench bottom.Similar to the previous subsection, we may assume a larger reaction rate at the trench bottom than at the trench openings to simulate the effect of accelerating and inhibiting additives.Equation 3.2 is used in this simulation with d kw, where w 20 and k is the aspect ratio.Table 2 shows the reaction rate at trench opening and the plating time in void-free trench filling for the k : 1 aspect ratio feature, k 2, 3, and 4. Larger deposition rate at the bottom than at the trench openings allows bottom-up growth.The overall reaction rate can be set much larger than the constant reaction rate shown in Table 1, and the void-free filling can still be achieved.Thus, the deposition is more efficient.

Conclusion
A numerical method is proposed for simulating the growth properties of electroless Cu deposits.The mathematical model used in this paper is a reaction-diffusion equation of the dominant ionic species.The simulation shows that the copper particles grow conformally at the first stage of the deposition.The out-plane growth becomes faster than the in-plane growth when 70%-80% of the substrate is covered by the particles.The nonconformal growth is more significant if the deposition time is prolonged.This results in large thickness of the films.As the reaction rate increases or the diffusivity decreases, this problem becomes more severe and void formation occurs at the bottom of the film between Cu particles.The growth properties observed in the simulation agree with those observed in experimentation.A lower reaction rate or accelerating and inhibiting additives in the solution can resolve this problem.However, improving the density of seeds is more efficient if ultrathin films are desired.Similarly, the simulation of trench filling shows that void-free trench filling can be achieved with a very small reaction rate or accelerating and inhibiting additives.However, a small reaction rate results in a large plating time.In contrast, the combination of accelerating and inhibiting additives is much more efficient.

Figure 2 :
Figure 2: The growth of Cu particles using parameter values in 3.1 at a t 2, b t 10, and c t 18. d The in-plane growth and out-plane growth over the plating time.

Figure 3 :
Figure 3: The deposition simulated with F 60 and all other parameter values set as in 3.1 .a The nonconformal growth becomes more severe when the reaction rate increases.b Voids form on the bottom between particles.

Figures 4 eFigure 4 :
Figure 4: The growth of Cu particles with a a low reaction rate F 15, c a reaction rate function F z in 3.2 , and e densely populated catalytic seeds.b , d , and f The in-plane growth and out-plane growth for deposition in a , c , and e , respectively.

Figure 5 :Figure 6 :
Figure 5: a The seed layer in the trench.b Void formation occurs at F 6 for a 1 : 1 aspect ratio feature.c Void-free filling is achieved at F 3.
the boundary of E. In the process of electroless Cu deposition, the Cu-complex ions in the solution are reduced to Cu and deposited on the surfaces of Cu particles.The sizes of the Cu particles are dependent on the plating time t.Let Γ 1i t ⊂ E, i 1, 2, ..., N, be the surfaces of the regions occupied by N Cu particles at time t as shown in Figure1b .The set Γ 1N 1 t Γ 1 \ { x, y, 0 | x, y, z 0 ∈ Γ 1i t , i 1, 2, ..., N for somez 0 ∈ 0, c } represents the region on the bottom of the solution that is not covered by Cu particles.The domain Ω t of the problem is formed by the region bounded by∪ N 1 i 1 Γ 1i t ∪Γ 2 ∪Γ 3 .Let u t,x, y, z be the ionic copper concentration at the point x, y, z ∈ Ω t and plating time t.

Table 1 :
Void-free trench filling with constant reaction rate F.

Table 2 :
Void-free trench filling with reaction rate F y defined by 3.2 for d kw and w 20.