Qualitative Simulation of the Growth of Electrolessly Deposited Cu Thin Films

Abstract

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.

1. 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 isCuEDTA²⁻+2HCHO+4OH⁻⟶Cu+H₂+2H₂O+2HCOO⁻+EDTA⁴⁻,(1.1) with the overall reaction rate𝑅=𝐾Cu²⁺^𝑎[HCHO]^𝑏[OH⁻]^𝑐Lig^𝑑,(1.2) where the constants 𝐾, 𝑎, 𝑏, 𝑐, and 𝑑 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𝐶_𝑡=𝐷𝐶_𝑥𝑥+𝐶_𝑦𝑦,(1.3) where C is the ionic copper concentration and 𝐷 is the diffusion coefficient. The boundary conditions for (1.3) are (i) 𝐶=𝐶₀ well away from the reaction surface, and (ii) 𝐷(𝑑𝐶/𝑑𝑛)=𝐹𝐶^𝑟 where 𝐹 and 𝑟 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.

2. Numerical Method

Suppose the deposition solution is prepared in a container of dimension [0,𝛼]×[0,𝛽]×[0,𝑑] 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 𝑐 with an appropriate boundary condition imposed to the problem, as proposed in [2].

(a)

(b)

Figure 1 

The computational domain.

Consider a rectangular region 𝐸=[0,𝑎]×[0,𝑏]×[0,𝑐]. Let Γ₁={(𝑥,𝑦,𝑧)∈𝐸∣𝑧=0}, Γ₂={(𝑥,𝑦,𝑧)∈𝐸∣𝑥=0}∪{(𝑥,𝑦,𝑧)∈𝐸∣𝑥=𝑎}∪{(𝑥,𝑦,𝑧)∈𝐸∣𝑦=0}∪{(𝑥,𝑦,𝑧)∈𝐸∣𝑦=𝑏} and Γ₃={(𝑥,𝑦,𝑧)∈𝐸∣𝑧=𝑐} the boundary of 𝐸. 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 𝑡. Let Γ_1𝑖(𝑡)⊂𝐸, 𝑖=1,2,…,𝑁, be the surfaces of the regions occupied by 𝑁 Cu particles at time 𝑡 as shown in Figure 1(b). The set Γ_1𝑁+1(𝑡)=Γ₁⧵{(𝑥,𝑦,0)∣(𝑥,𝑦,𝑧₀)∈Γ_1𝑖(𝑡),𝑖=1,2,…,𝑁forsome𝑧₀∈[0,𝑐]} represents the region on the bottom of the solution that is not covered by Cu particles. The domain Ω(𝑡) of the problem is formed by the region bounded by ∪^𝑁+1𝑖=1Γ_1𝑖(𝑡)∪Γ₂∪Γ₃. Let 𝑢(𝑡,𝑥,𝑦,𝑧) be the ionic copper concentration at the point (𝑥,𝑦,𝑧)∈Ω(𝑡) and plating time 𝑡. At 𝑡=0, Γ_1𝑖(0), 𝑖=1,2,…,𝑁, are the surfaces of the catalytic seeds deposited on the substrate to initiate the plating process. Let𝑢(0,𝑥,𝑦,𝑧)=𝑢₀,(𝑥,𝑦,𝑧)∈Ω(0),(2.1) be the initial concentration of the solution. Equation (1.3) is extended to the following reaction-diffusion problem with a nonlinear reaction acting on ∪^𝑁𝑖=1Γ_1𝑖:𝑢_𝑡=𝐷𝜕𝑢𝑢_𝑥𝑥+𝑢_𝑦𝑦+𝑢_𝑧𝑧,inΩ(𝑡),𝑡>0,(2.2)𝜕𝑛=0,onΓ^1𝑁+1𝐷𝜕𝑢∪Γ₂,𝑡>0,(2.3)𝜕𝑛=𝐹𝑢^𝑟,on^𝑁_𝑖=1𝑢=𝑢₀Γ_1𝑖,𝑡>0,(2.4),onΓ₃,𝑡>0.(2.5) To describe and track the moving boundaries, Γ_1𝑖(𝑡)⊂𝐸, 𝑖=1,2,…,𝑁, a level set method [6] is used to formulate the motion of these boundaries. Let 𝑣(𝑡,𝑥,𝑦,𝑧) defined on ℜ×𝐸 be the so-called level set function: 𝑣(𝑡,𝑥,𝑦,𝑧)=𝑣^∗ for (𝑥,𝑦,𝑧) inside the interface Γ_1𝑖(𝑡), 1≤𝑖≤𝑁, and 𝑣(𝑡,𝑥,𝑦,𝑧)=0 otherwise. Here, 𝑣^∗ is a prescribed value related to the density of Cu so that the deposition rate is 𝐹𝑢^𝑟/𝑣^∗ at any point on the surface [2]. The motion of the interfaces Γ_1𝑖(𝑡), 1≤𝑖≤𝑁 is transported under the velocity field 𝑣⃗=−𝐹𝑢^𝑟/𝑣^∗𝑛,⃗ where 𝑛⃗=∇𝑣/|∇𝑣|. The level set equation [6] is𝑣_𝑡+𝑣⃗⋅∇𝑣=0.(2.6) This is equivalent to𝑣_𝑡=𝐹𝑢^𝑟𝑣^∗||∇𝑣||.(2.7)

Consider a uniform rectangular mesh on 𝐸 with mesh size ℎ=𝑎/𝑁₁=𝑏/𝑁₂=𝑐/𝑁₃. Let (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘)=(𝑖ℎ,𝑗ℎ,𝑘ℎ)∈𝐸 be a grid point where 0≤𝑖≤𝑁₁, 0≤𝑗≤𝑁₂, and 0≤𝑘≤𝑁₃. The increment in 𝑡 is denoted by Δ𝑡 and 𝑡^𝑛=𝑛Δ𝑡 for 𝑛≥0. The approximation of 𝑢(𝑡^𝑛,𝑥_𝑖,𝑦_𝑗,𝑧_𝑘) is denoted by the standard notation 𝑢^{𝑛𝑖𝑗𝑘}. Equation (2.2) is discretized using the backward-time central-space scheme [7]. The backward-time discretization for time derivatives is given by 𝑢_𝑡(𝑡^𝑛,𝑥_𝑖,𝑦_𝑗,𝑧_𝑘)≈(𝑢^{𝑛𝑖𝑗𝑘}−𝑢^{𝑛−1𝑖𝑗𝑘})/Δ𝑡, and the central-space discretization for second-order spatial derivatives is given by𝑢_𝑥𝑥=𝑡^𝑛,𝑢𝑥_𝑖,𝑦_{𝑗𝑛𝑖−1𝑗𝑘},𝑧_𝑘+𝑢_𝑦𝑦𝑡^𝑛,𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+𝑢_𝑧𝑧𝑡^𝑛,𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−2𝑢^{𝑛𝑖𝑗𝑘}+𝑢^{𝑛𝑖+1𝑗𝑘}ℎ²+𝑢^{𝑛𝑖𝑗−1𝑘}−2𝑢^{𝑛𝑖𝑗𝑘}+𝑢^{𝑛𝑖𝑗+1𝑘}ℎ²+𝑢^{𝑛𝑖𝑗𝑘−1}−2𝑢^{𝑛𝑖𝑗𝑘}+𝑢^{𝑛𝑖𝑗𝑘+1}ℎ²+𝑂ℎ².(2.8) Equation (2.3) is discretized by the central difference method for first derivatives. For example,𝜕𝑢𝜕𝑛𝑡^𝑛,0,𝑦_𝑗,𝑧_𝑘=𝜕𝑢𝜕𝑥𝑡^𝑛,0,𝑦_𝑗,𝑧_𝑘=𝑢^{𝑛1𝑗𝑘}−𝑢^{𝑛−1𝑗𝑘}2ℎ+𝑂ℎ².(2.9) Therefore, (2.3) gives𝑢𝑢^{𝑛−1𝑗𝑘}_{𝑛𝑁1+1𝑗𝑘}=𝑢^{𝑛1𝑗𝑘},𝑢^{𝑛𝑖−1𝑘}=𝑢^{𝑛𝑖1𝑘},𝑢^{𝑛𝑖𝑗−1}=𝑢^{𝑛𝑖𝑗1},=𝑢^𝑛𝑁_{1−1𝑗𝑘},𝑢^𝑛𝑖𝑁_2+1𝑘=𝑢^𝑛𝑖𝑁_2−1𝑘.(2.10) Equation (2.4) is also discretized by the central difference method analogous to (2.9). Let (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘)∈Γ_1𝑝(𝑡^𝑛) with 1≤𝑝≤𝑁. Suppose that all of its neighboring grid points are in Ω(𝑡^𝑛) except for (𝑥_𝑖−1,𝑦_𝑗,𝑧_𝑘), which is inside the Cu particle 𝑝 at 𝑡=𝑡^𝑛. The reaction term in (2.4) is estimated by the backward time approximation 𝐹𝑢^𝑟≈𝐹(𝑢^{𝑛−1𝑖𝑗𝑘})^𝑟. Equation (2.4) is then discretized by𝑢^{𝑛𝑖+1𝑗𝑘}−𝑢^{𝑛𝑖−1𝑗𝑘}2ℎ=𝐷𝐹𝑢^{𝑛−1𝑖𝑗𝑘}^𝑟.(2.11) So,𝑢^{𝑛𝑖−1𝑗𝑘}=𝑢^{𝑛𝑖+1𝑗𝑘}−2ℎ𝐹𝐷𝑢^{𝑛−1𝑖𝑗𝑘}^𝑟.(2.12) Suppose two of the neighboring grid points, say (𝑥_𝑖−1,𝑦_𝑗,𝑧_𝑘) and (𝑥_𝑖,𝑦_𝑗−1,𝑧_𝑘), are located inside the Cu particle 𝑝 at 𝑡=𝑡^𝑛. Let 𝑢^{𝑛(𝑖−1/2)(𝑗−1/2)𝑘}≈𝑢(𝑡^𝑛,𝑥_𝑖−ℎ/2,𝑦_𝑗−ℎ/2,𝑧_𝑘), and 𝑢^{𝑛(𝑖−1/2)(𝑗−1/2)𝑘}=(𝑢^{𝑛𝑖−1𝑗𝑘}+𝑢^{𝑛𝑖𝑗−1𝑘})/2. Similar approximation is made for 𝑢^{𝑛(𝑖+1/2)(𝑗+1/2)𝑘}. Equation (2.4) for this case is then approximated by𝑢^{𝑛(𝑖+1/2)(𝑗+1/2)𝑘}−𝑢^{𝑛(𝑖−1/2)(𝑗−1/2)𝑘}√2ℎ=𝐷𝐹𝑢^{𝑛−1𝑖𝑗𝑘}^𝑟.(2.13) So,𝑢^{𝑛𝑖−1𝑗𝑘}+𝑢^{𝑛𝑖𝑗−1𝑘}=𝑢^{𝑛𝑖+1𝑗𝑘}+𝑢^{𝑛𝑖𝑗+1𝑘}−2√2ℎ𝐹𝐷𝑢^{𝑛−1𝑖𝑗𝑘}^𝑟.(2.14) If three of the neighboring grid points, say (𝑥_𝑖−1,𝑦_𝑗,𝑧_𝑘), (𝑥_𝑖,𝑦_𝑗−1,𝑧_𝑘), and (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−1), are located inside the Cu particle 𝑝. The approximation is analogous𝑢^{𝑛𝑖−1𝑗𝑘}+𝑢^{𝑛𝑖𝑗−1𝑘}+𝑢^{𝑛𝑖𝑗𝑘−1}=𝑢^{𝑛𝑖+1𝑗𝑘}+𝑢^{𝑛𝑖𝑗+1𝑘}+𝑢^{𝑛𝑖𝑗𝑘+1}−2√3ℎ𝐹𝐷𝑢^{𝑛−1𝑖𝑗𝑘}^𝑟.(2.15) Hence, the discretization of (2.1)–(2.5) is given as follows: 𝑢^{0𝑖𝑗𝑘}=𝑢₀ for 0≤𝑖≤𝑁₁, 0≤𝑗≤𝑁₂, and 0≤𝑘≤𝑁₃; 𝑢^{𝑛𝑖𝑗𝑁}₃=𝑢₀ for 0≤𝑖≤𝑁₁, 0≤𝑗≤𝑁₂, and 𝑛>0𝑢^{𝑛𝑖𝑗𝑘}−𝑢^{𝑛−1𝑖𝑗𝑘}Δ𝑡=ℎ𝐷²𝑢^{𝑛𝑖−1𝑗𝑘}+𝑢^{𝑛𝑖+1𝑗𝑘}+𝑢^{𝑛𝑖𝑗+1𝑘}+𝑢^{𝑛𝑖𝑗−1𝑘}+𝑢^{𝑛𝑖𝑗𝑘−1}+𝑢^{𝑛𝑖𝑗𝑘+1}−6𝑢^{𝑛𝑖𝑗𝑘}−𝑤^{𝑛−1𝑖𝑗𝑘}𝐹𝑢^{𝑛−1𝑖𝑗𝑘}^𝑟ℎ,(2.16) for 0≤𝑖≤𝑁₁,0≤𝑗≤𝑁₂,0≤𝑘<𝑁₃, and 𝑛>0 with the substitution of (2.10)–(2.15) for the grid points outside Ω(𝑡). The weight 𝑤^{𝑛−1𝑖𝑗𝑘}=0 if (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘) is in Ω(𝑡^𝑛), and 𝑤^{𝑛−1𝑖𝑗𝑘}>0 if (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘) is on ∪^𝑁𝑙=1Γ_1𝑙(𝑡^𝑛). Since the Cu deposits also depend on the size of the surface, the weight 𝑤^{𝑛−1𝑖𝑗𝑘} is determined not only by the coefficients of the reaction terms in (2.12)–(2.15) but also the location of (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘), which will be explained later.

Let 𝑣^{𝑛𝑖𝑗𝑘} be the approximation of 𝑣(𝑡^𝑛,𝑥_𝑖,𝑦_𝑗,𝑧_𝑘). From the definition of the level set function 𝑣, let 𝑣^{𝑛𝑖𝑗𝑘}=𝑣^∗ if (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘) is inside the Cu particle 𝑝, 1≤𝑝≤𝑁, and 𝑣^{𝑛𝑖𝑗𝑘}=0 if (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘) is in Ω(𝑡^𝑛). Here, let 0≤𝑣^{𝑛𝑖𝑗𝑘}<𝑣^∗ if (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘) is on ∪^𝑁𝑙=1Γ_1𝑙(𝑡^𝑛) for computational purpose. Only 𝑣^{𝑛𝑖𝑗𝑘}, where (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘) is on ∪^𝑁𝑙=1Γ_1𝑙(𝑡^𝑛), has to be computed. For a boundary point (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘) as described in (2.12), we have 𝑣^{𝑛−1𝑖−1𝑗𝑘}=𝑣^∗ and 𝑣^{𝑛−1𝑖+1𝑗𝑘}=0. Equation (2.7) is discretized using the backwark-time central-space scheme𝑣^{𝑛𝑖𝑗𝑘}−𝑣^{𝑛−1𝑖𝑗𝑘}Δ𝑡=𝐹𝑢𝑣^{𝑛−1𝑖𝑗𝑘∗}𝑣^{𝑛−1𝑖−1𝑗𝑘}−𝑣^{𝑛−1𝑖+1𝑗𝑘}=𝐹𝑢_{𝑛−1𝑖𝑗𝑘}2ℎ2ℎ.(2.17) For a boundary point as described in (2.14), we have 𝑣^{𝑛−1(𝑖−1/2)(𝑗−1/2)𝑘}=(𝑣^{𝑛−1𝑖−1𝑗𝑘}+𝑣^{𝑛−1𝑖𝑗−1𝑘})/2=𝑣^∗ and 𝑣^{𝑛−1(𝑖+1/2)(𝑗+1/2)𝑘}=(𝑣^{𝑛−1𝑖+1𝑗𝑘}+𝑣^{𝑛−1𝑖𝑗+1𝑘})/2=0. The discretization of (2.7) becomes𝑣^{𝑛𝑖𝑗𝑘}−𝑣^{𝑛−1𝑖𝑗𝑘}Δ𝑡=𝐹𝑢𝑣^{𝑛−1𝑖𝑗𝑘∗}𝑣^{𝑛−1(𝑖−1/2)(𝑗−1/2)𝑘}−𝑣^{𝑛−1(𝑖+1/2)(𝑗+1/2)𝑘}=𝐹𝑢_{𝑛−1𝑖𝑗𝑘}√2ℎ√2ℎ.(2.18) Similarly, the discretization for a boundary point as described in (2.15) is𝑣^{𝑛𝑖𝑗𝑘}−𝑣^{𝑛−1𝑖𝑗𝑘}Δ𝑡=√3𝐹𝑢^{𝑛−1𝑖𝑗𝑘}2ℎ.(2.19) From (2.17)–(2.19), the discretization of (2.7) can be writen as𝑣^{𝑛𝑖𝑗𝑘}=𝑣^{𝑛−1𝑖𝑗𝑘}+𝜇^{𝑛−1𝑖𝑗𝑘}𝐹𝑢^{𝑛−1𝑖𝑗𝑘}^𝑟ℎΔ𝑡.(2.20) Again, the weight 𝜇^{𝑛−1𝑖𝑗𝑘} is determined not only by the coefficients in (2.17)–(2.19) but also the location of (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘). Rewrite (2.16) as 1+6𝐷ℎΔ²𝑡𝑢^{𝑛𝑖𝑗𝑘}−𝐷Δ𝑡ℎ²𝑢^{𝑛𝑖−1𝑗𝑘}=𝑢_{𝑛−1𝑖𝑗𝑘}+𝑢^{𝑛𝑖+1𝑗𝑘}+𝑢^{𝑛𝑖𝑗−1𝑘}+𝑢^{𝑛𝑖𝑗+1𝑘}+𝑢^{𝑛𝑖𝑗𝑘−1}+𝑢^{𝑛𝑖𝑗𝑘+1}−𝑤^{𝑛−1𝑖𝑗𝑘}𝐹𝑢^{𝑛−1𝑖𝑗𝑘}^𝑟ℎΔ𝑡.(2.21) Now, 𝑤^{𝑛−1𝑖𝑗𝑘}(𝐹(𝑢^{𝑛−1𝑖𝑗𝑘})^𝑟/ℎ)Δ𝑡 in (2.21) estimates the Cu ions consumed at (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘) over [𝑡^𝑛−1,𝑡^𝑛] in the deposition process, while 𝜇^{𝑛−1𝑖𝑗𝑘}(𝐹(𝑢^{𝑛−1𝑖𝑗𝑘})^𝑟/ℎ)Δ𝑡 in (2.20) estimates the Cu deposits at (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘) over [𝑡^𝑛−1,𝑡^𝑛]. In the simulation, we set 𝑤^{𝑛−1𝑖𝑗𝑘}=𝜇^{𝑛−1𝑖𝑗𝑘} for a boundary grid point (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘) so that the Cu ions consumed equal the Cu deposits.

Let Γ^{𝑛𝑖𝑛,𝑝} be the approximation of Cu particle 𝑝, Γ^𝑛𝑝 the approximation of Γ_1𝑝(𝑡^𝑛), and Ω^𝑛 the approximation of Ω(𝑡^𝑛). Γ^𝑛𝑝 is formed by the grid points satisfying (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘)∉∪^𝑁𝑙=1Γ^{𝑛𝑖𝑛,𝑙} and at least one of (𝑥_{𝑖±𝑖𝑖},𝑦_{𝑗±𝑗𝑗},𝑧_{𝑘±𝑗𝑗}) lies in Γ^{𝑛𝑖𝑛,𝑝}, where 𝑖𝑖=0,±1, 𝑗𝑗=0,±1, and 𝑘𝑘=0,±1. Then, 𝑣^{𝑛𝑖𝑗𝑘}=0 if (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘)∈Ω^𝑛, 0≤𝑣^{𝑛𝑖𝑗𝑘}<𝑣^∗ if (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘)∈∪^𝑁𝑙=1Γ^𝑛𝑙, and 𝑣^{𝑛𝑖𝑗𝑘}=𝑣^∗ if (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘)∈∪^𝑁𝑙=1Γ^{𝑛𝑖𝑛,𝑙}. At the beginning of the deposition, let Γ^{0𝑖𝑛,𝑝}, 1≤𝑝≤𝑁, be given for the deposition of catalytic seeds. Then, the set of boundary points ∪^𝑁𝑙=1Γ^0𝑝 and Ω⁰ are determined by the above definition. The initial value 𝑣^{0𝑖𝑗𝑘} is given as above for 𝑛=0 except for 𝑣^{0𝑖𝑗𝑘}=0 when (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘)∈∪^𝑁𝑙=1Γ^0𝑙. At 𝑡=𝑡^𝑛, let (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘)∈∪^𝑁𝑙=1Γ^𝑙𝑛−1. Then, 𝑣^{𝑛𝑖𝑗𝑘} is computed by (2.20). Now, suppose 𝑣^{𝑛𝑖𝑗𝑘}≥𝑣^∗. Set 𝑣^{𝑛𝑖𝑗𝑘}=𝑣^∗, add (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘) to ∪^𝑁𝑙=1Γ^{𝑛𝑖𝑛,𝑙}, delete it from ∪^𝑁𝑙=1Γ^𝑛𝑙, and add its neighboring grid points which are in Ω^𝑛−1 to the set ∪^𝑁𝑙=1Γ^𝑛𝑙. Here, the neighboring points of (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘) are (𝑥_{𝑖+𝑖𝑖},𝑦_{𝑗+𝑗𝑗},𝑧_{𝑘+𝑘𝑘}) with 𝑖𝑖, 𝑗𝑗, 𝑘𝑘=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 ∪^𝑁𝑙=1Γ^{𝑛𝑖𝑛,𝑙}. For example, suppose (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘)∈Γ^𝑝𝑛−1 and (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘)∈Γ^{𝑛𝑖𝑛,𝑝}. If (𝑥_𝑖+1,𝑦_𝑗,𝑧_𝑘), (𝑥_𝑖,𝑦_𝑗+1,𝑧_𝑘), and (𝑥_𝑖+1,𝑦_𝑗+1,𝑧_𝑘) are added to Γ^𝑛𝑝, then it is assumed that the rate of Cu deposits at (𝑥_𝑖+1,𝑦_𝑗+1,𝑧_𝑘) is smaller than that at (𝑥_𝑖+1,𝑦_𝑗,𝑧_𝑘) and (𝑥_𝑖,𝑦_𝑗+1,𝑧_𝑘). Thus, 𝜇^{𝑛𝑖𝑗𝑘} may have three different vales. In this paper, 𝜇^{𝑛𝑖𝑗𝑘}=1,1/2,or,1/3 is used in the numerical simulation. A larger 𝜇^{𝑛𝑖𝑗𝑘} value is corresponding to a smaller distance from (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘)∈∪^𝑁𝑙=1Γ^𝑛𝑙 to its nearest grid point in ∪^𝑁𝑙=1Γ^{𝑛𝑖𝑛,𝑙}. The computation procedure is summarized as follows: given the initial conditions Γ^{0𝑖𝑛,𝑝} for 1≤𝑝≤𝑁, and 𝑢^{0𝑖𝑗𝑘} for 0≤𝑖≤𝑁₁, 0≤𝑗≤𝑁₂, and 0≤𝑘≤𝑁₃, then Γ^0𝑝, Ω⁰, and 𝑣^{0𝑖𝑗𝑘} can be determined. At each time step 𝑡^𝑛, 𝑣^{𝑛𝑖𝑗𝑘} is computed by (2.20) at each boundary grid point in ∪^𝑁𝑙=1Γ^𝑙𝑛−1. Then, 𝑣^{𝑛𝑖𝑗𝑘} is used to update the set ∪^𝑁𝑙=1Γ^{𝑛𝑖𝑛,𝑙}, the boundary set ∪^𝑁𝑙=1Γ^𝑛𝑙, and the domain Ω^𝑛. It is followed by the computation of 𝑢^{𝑛𝑖𝑗𝑘} using (2.21) at the points which are not in ∪^𝑁𝑙=1Γ^{𝑛𝑖𝑛,𝑙}.

3. 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 ℎ defined in Section 2 has to be set at about 1 nm and the diffusion coefficient 𝐷=10⁻⁵cm²/sec=10⁹nm²/sec. Thus, ℎ²/(𝐷Δ𝑡)=10⁻⁹/Δ𝑡 is very close to zero unless Δ𝑡 is small enough, for example, Δ𝑡<10⁻⁸. When ℎ²/(𝐷Δ𝑡) 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 Δ𝑡 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⁻⁹ for all numerical examples presented in this section. Thus, the convergence of the iterative method is assumed.

3.1. Simulation of the Growth of Cu Particles

Consider the rectangular domain 𝑅=[0,60]×[0,60]×[0,44]. Let the parameter values be𝐹Δ𝑡==400.,001,𝐷=𝑢1₀000,𝑟=0.8,ℎ=1,=0.03,𝑢^∗=1.4.(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 𝑡=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 𝐹 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 𝐹 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 𝐹=40.

(a)

(b)

(c)

(d)

Figure 2 

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

(a)

(b)

Figure 3 

The deposition simulated with 𝐹=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.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4 

The growth of Cu particles with (a) a low reaction rate 𝐹=15, (c) a reaction rate function 𝐹(𝑧) 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.

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 𝑧. Define𝐹(𝑧)=𝐹𝐹,2−𝑑𝑧,iiff𝑧𝑧<≥𝑑𝑑,.(3.2) The reaction rate is 2𝐹 at the bottom, 𝐹 for 𝑧≥𝑑, and decreases linearly from the bottom to 𝑧=𝑑. Consider 𝐹=20 and 𝑑=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 𝐹=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. Figures 4(e) and 4(f) show the growth of the particles are more conformal, the particle size is smaller, and the plating time is shorter. The thickness and plating time of the deposits are reduced by nearly half compared with Figures 2(c) and 2(d). Thus, improving the density of the catalytic seeds is the most effective method for reducing the film thickness. Chen et al. [5] developed a new seeding process in 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 𝐷 decreases. For example, let 𝐷=300, 𝐹=15, and all other parameter values are set as in (3.1). The growth of the particles is similar to Figure 3, where 𝐷=1000 and 𝐹=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.

3.2. 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 𝑅=[0,40]×[0,40+𝑘𝑤] for a 𝑘 : 1 aspect-ratio feature, and let the parameter values be set as in (3.1) with 𝑤=20. Figure 5(b) shows that voids occur for the 1 :  1 aspect-ratio feature even at a low reaction rate 𝐹=6. To achieve void-free trench filling, the deposition rate has to be very small. Figure 5(c) shows 𝐹=3 for void-free filling. As aspect ratio increases from 1 : 1 to 2 : 1, 𝐹=3 fails to be void-free filling as shown in Figure 6(a). Void-free filling can be achieved if 𝐹=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 𝑘 : 1 aspect-ratio feature, 𝑘=1,2,3,and4. 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.

Table 1 

Void-free trench filling with constant reaction rate 𝐹.

(a)

(b)

(c)

Figure 5 

(a) The seed layer in the trench. (b) Void formation occurs at 𝐹=6 for a 1 : 1 aspect ratio feature. (c) Void-free filling is achieved at 𝐹=3.

(a)

(b)

Figure 6 

(a) Void formation occurs at 𝐹=3 for a 2 : 1 aspect ratio feature. (b) Void-free filling is achieved at 𝐹=0.5.

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 𝑑=𝑘𝑤, where 𝑤=20 and 𝑘 is the aspect ratio. Table 2 shows the reaction rate at trench opening and the plating time in void-free trench filling for the 𝑘 : 1 aspect ratio feature, 𝑘=2,3,and4. 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.

Table 2 

Void-free trench filling with reaction rate 𝐹(𝑦) defined by (3.2) for 𝑑=𝑘𝑤 and 𝑤=20.

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

Acknowledgment

The author would like to thank Ms. Cheryl Robbins for her help in editing this paper.