Material removal mechanism of copper chemical mechanical polishing with different particle sizes based on quasi-continuum method

In this paper, the material removal mechanism of copper chemical mechanical polishing was studied by the quasicontinuum method that integrated molecular dynamics and the finite element method. By analyzing the abrasive process of different particle sizes on single crystal copper, we investigated the internal material deformation, the formation of chips, the stress distribution, and the change of cutting force. Results showed that shear band deformation was generated along the cutting direction at approximately 45° inside the workpiece material. The deformation was accompanied by dislocations and sliding phenomena in the shear band region. Smaller abrasive particle size led to poor quality of the workpiece, while a larger particle size led to better quality. However, larger particle size resulted in greater plastic deformation and deeper residual stress inside the workpiece. Size change of abrasive particles had little effect on the tangential cutting force.


Introduction
In IC manufacturing, copper chemical mechanical polishing (CMP) is now the most widely used planarization technology for silicon wafer processing because it achieves a super-smooth surface without injury. It is also the ideal multilayer wiring planarization method. Despite this advantage, systematic studies are still lacking, especially, in theoretical research. When comparing CMP with other technologies, CMP technology has been developed from real-world usage. The invention, development, implementation, and application of CMP have mostly occurred in industry rather than in academia. The CMP process involves tribology, chemistry, fluid mechanics, solidstate physics, and many other disciplines. Therefore, its mechanism is very complicated and increases the difficulty of research.
Precise control of the CMP process depends largely on the understanding of its material removal mechanism. When materials are removed by removing several atoms or atomic layers at a time by the CMP process, it is difficult to control and observe the process. Therefore, several difficulties exist in the experiments, calculations, and analysis of the CMP process. Steigerwald et al. provide a survey on the status of current CMP modeling [1]. They investigate different modeling assumptions and discuss the controversial treatment of physical effects. Many scholars use the molecular dynamics method to study the material removal mechanism of CMP. Yokosuka et al. developed a new accelerated quantum chemical molecular dynamics program known as "Colors" to simulate the CMP process. It is more than 5,000 times faster than the original first principles molecular dynamics program, as it is based on their original tight-binding theory [2]. Ye et al. simulate nanoscale polishing of a copper surface with molecular dynamics by utilizing the embedded atom method. Mechanical abrasion produces rough planarized surfaces with a large chip in front of the abrasive particle and with dislocations in the bulk of the crystal [3]. Han et al. have performed molecular dynamics analysis of chemical mechanical polishing of a silicon wafer and the simulation result shows that huge hydrostatic pressure is induced in the local area and leads to silicon atoms transforming from the classical diamond structure (α silicon) to metal structure (β silicon) [4]. Si et al. have used molecular dynamics simulations of a nanoscratching processes to investigate the removal mechanism of single crystalline silicon in the CMP process on the atomic-scale. The simulation results under a scratching depth of 1 nm showed that a thick layer of silicon material was removed by chip formation and an amorphous layer was formed on the silicon surface after nanoscratching [5]. Wu et al. investigated and calculated the wafer topography effects on the contact pressure distribution during CMP with a 3D solid-solid contact model. They proposed a formula that helps to specify the polishing parameters during CMP [6]. Zhou et al. used mixed elastohydrodynamic lubrication model to investigate the CMP process and the results showed that the proposed layer elastic theory is an optimal model to describe the polishing pad deformation behavior in CMP [7]. Chen et al. used molecular dynamics simulation to analyze the material removal mechanism of a silicon substrate under the impact of a large porous silica cluster with different pore diameters. Their findings were inspiring for other researchers who want to optimize the process parameters in order to obtain a lower surface roughness and higher material removal rate during the CMP process [8]. Lee et al. proposed a mathematical model-based evaluation method to determine the environmental burden of the CMP process. They adopt our previously reported material removal rate model for CMP and modified it to incorporate the effect of the slurry flow rate and process temperature. The results showed that slurry consumption strongly impacts the carbon dioxide equivalents of the CMP process [9]. Yang et al. proposed a chip-scale CMP model that is based on the Greenwood-Williamson theory. It combines the effects of both the pad bulk deformation and the feature size, which are mainly responsible for global pressure distribution and the feature-scale removal rate variation, respectively [10].
However, the simulated scale range of the CMP process by molecular dynamics method is limited. The maximum grain radius for simulation generally cannot exceed 10 nm, while the actual average grain radius in CMP process is 50 nm according to the work of Larsen-Basse et al. [11]. The quasicontinuum (QC) method, a multiscale simulation method that integrates the molecular dynamics method and finite element method, is potentially an effective way to solve this problem based on the work of Tadmor [12].
The QC method, proposed by Tadmor et al. [12], is one of the multiscale simulation methods. By mixing a continuum and atomistic approach, this method uses an atomic scale solution near dislocation core regions while using "representative atoms" to conduct a rough description away from the dislocation core regions. By calculating the system's energy and force in the state of a reduced degree of freedom, the mixing of atoms and continuum is achieved. In this paper, a CMP multiscale computer simulation model is used to simulate the CMP process of the microparticle contact interaction with the workpiece, which is the mechanical action process of the CMP material removal mechanism. Figure 1 depicts the basic principle of the QC method. The black filled circles in Fig. 1(a) are the representative atoms, and Fig. 1(b) shows the finite element mesh that is composed of the corresponding representative atoms. As shown in Fig. 1(b), atoms B, C, and D are finite element nodes, and atom A is a nonrepresentative atom. The criteria for judging whether an atom is a representative atom or nonrepresentative atom is based on the formula, it is a representative atom; otherwise, it is a nonrepresentative atom. The displacements of nonrepresentative atoms are derived from the linear interpolation for representative atoms. This purpose can be achieved using the interpolated equations of the finite element method in the QC method. In regions where a fully detailed atom description is required, we can choose all atoms as representative atoms and reduce the density of representative atoms in regions with a smaller deformation gradient. As shown in Fig. 1, we may choose all atoms in the dislocation core region as representative atoms and reduce the density of representation atoms in regions away from the dislocation core region. Currently, the QC method has been used to study a variety of basic phenomenon in crystal deformation such as fracture, grain boundary structure and deformation, nanoindentation, and interaction of threedimenssional dislocations. Miller et al. has considered fracture in two distinct settings. First, cracks in single crystals are investigated, and are then considered as a crack advancing towards a grain boundary in its path. In the investigation of single crystal fracture, the competition between simple cleavage and crack-tip dislocation emission is evaluated [13]. Shenoy et al. have shown how their mixed atomistic and continuum analysis is adapted to the treatment of interfacial deformation. Such calculations demand the generalization of the original quasi-continuum formulation to allow for the existence of more than one grain simultaneously [14]. Tadmor et al. investigate two different crystallographic orientations of the film and two different indenter geometries, a rectangular prism and a cylinder. They obtain both macroscopic load versus indentation depth curves and microscopic quantities, such as Peierls stress and density of geometrically necessary dislocations beneath the indenter [15].

Simulation model
A QC program was adopted in this paper to calculate the core, and the simulation was achieved by secondary development on that basis. The entire program includes parameters to calculate meshing, boundary conditions and outer loads, the QC computation core, and output result modules, as shown in Fig. 2.
In the CMP process, part of the abrasive particles are embedded into the surface of the polishing pad. Each embedded particle is the equivalent of a fixed abrasive, plowing workpiece surface that removes workpiece material, as shown in Fig. 3. This single polishing process was simulated in the paper, and the simulation model was shown in Fig. 4. As shown in Fig. 4, a single crystal copper workpiece was meshed in a finite element grid. The top grid was meshed denser and denser until it became an atomic region, which achieved the integration of the atomistic approach and finite element method. The workpiece dimension was 500 Å × 250 Å. Perpendicular to the paper, the width of a layer of lattice was 3.6 Å. The workpiece contained 40,000 atoms, while the number of representative atoms required for computing was less than 5,000, greatly reducing the computation amount.  Ideally, in order to calculate the total energy, all atoms in the domain need to be visited by: where i E is the energy contribution from site i x . The precise form of i E depended on the potential function used. In the regions where the displacement field was smooth, keeping track of each individual atom was unnecessary. After selecting some representative atoms, the displacement of the remaining atoms could be approximated by linear interpolation.
In the equation below, where the subscript  identified the representative atoms, S  was an appropriate weight function, and rep N was the number of representative atoms involved.
Although this step reduced the degree of freedom needed for calculation, every atom still needed to be visited. To reduce the computational complexity involved in obtaining the total energy, several simplified rules were introduced. Among these rules, the Cauchy-Born rule assumed the deformation gradient, . With these approximations, the evaluation of the total energy was reduced to a summation over the finite elements: where e N is the number of elements. Equation (3) is called the local version of QC. In the presence of defects, the deformation tends to be nonsmooth. Therefore, the approximation made in the local QC is inaccurate. A nonlocal version of the QC has been developed in which the energy is expressed as: where n  is a chosen suitable weight, E  is the energy from each representative atom, and it is computed by visiting the representative atom's neighboring atoms whose positions are generated by the local deformation. Practical implementation usually combines both local and nonlocal versions of the method. A criterion has been proposed to identify the local/nonlocal regions, allowing the whole procedure to be adaptively applied. Figure 5 reviews the entire flow process of the QC method.
The composition of the abrasive particle is usually Al 2 O 3 , SiO 2, or CeO 2 in the actual process. The stiffness of the abrasive particle is much higher than the stiffness of the workpiece. Although the deformation of the particle is extremely small in the simulation process, the particle is considered a rigid body in the paper. Abrasive particles are initially placed at the left position of the workpiece. A displacement load is applied with a load step of 0.3 Å. There are 300 steps in total. The workpiece (001) is cut along the direction of the crystal plane [100] with a cutting depth of 10 Å. Fixed constraints are applied at the bottom and the right end surfaces of the workpiece. For the direction perpendicular to the paper, a periodic boundary condition is adopted. The simulation process is quasi two-dimensional. In fact, due to the static properties of the QC method, energy minimization is executed every time a load step is applied. As a result, the system is in an equilibrium state at each load step. The stiffness of abrasive particles is much higher than of the workpiece. Therefore, the particles are seen as a rigid body due to little deformation in the simulation process. The radiuses of abrasive particles for the simulation are 5 nm, 10 nm, and 20 nm. In comparison with the molecular dynamics method, the particle sizes in this paper are closer to the actual situation. For the copper workpiece, the embedded atom method (EAM) potential function is adopted with a cutoff radius of 4.95 Å. The EAM potential function that describes the interaction between copper atoms can be denoted by: where U represents the system's total energy, ( ) The simulation program runs on a Linux system, and approximately 100 hours were required for a computational process. Slightly more time was consumed when the particle radius increased. The model's parameters were shown in Table 1.  at the contact region of the workpiece and abrasive particles undergo initial deformation, and expand into a shear band along the 45° direction as the particles move forward in the plowing process. The deformation is generated mainly by compression and shear action.
With the particles moving, dislocation and slipping occur in the deformation region. When the strain energy stored in the deformed crystal lattice surpasses a certain value, atoms rearrange to release the unnecessary energy, and then return to the equilibrium state. When the strain energy is insufficient to rearrange the atoms, some atomic region can move in block to present a sidestep shape on the workpiece surface. When the particle radius is 5 nm, as shown in Fig. 6, chips can move forward as the particles move and they accumulate in front of the particles, blocking the forward way of the particles. As the particle radius increases, chips obviously no longer accumulate. That phenomenon results from the decreased angle between  the particles and the workpiece, causing the chips in front of the particles to undergo a squeezing action. As a result, workpiece material is removed by the shearing action of the particles, and the quality of the processed work piece surface is relatively poor.
As the radius increases as shown in Figs. 7 and 8, the workpiece material undergoes the shearing and squeezing actions from the particles, which leads to better surface quality and deeper deformation inside the work piece. Figure 7 and 8 show the same phenomena as described in Fig. 6.
Additionally, we can now know the advantage the QC method has in efficiency. The wafer system in this paper contains approximately 4 × 10 4 atoms, but the number of representative atoms needs to be calculated in an atomic region less than 5,000 for the QC method. Therefore, the calculation quantity can be greatly reduced in the QC method. Figure 9 shows the stress diagrams of the workpiece at the 70 th load step with the particle radiuses of 5 nm, 10 nm, and 20 nm.

Stress analysis
As shown in Fig. 9(b), most of the stress concentrates on the deformation region beneath the particles, and they are distributed along the shear band. The internal stress distribution of the workpiece is complex. In region A beneath the particles, the workpiece material undergoes plastic deformation due to compression action, creating a massive stress concentrate in this region. Only a small amount of stress is distributed in region B. In region C, a certain amount of stress is generated due to stretching action, and the stress increases when the particles move forward. Some stress is distributed in front of the accumulated region, shown in region D.
The phenomena in Figs. 9(a) and 9(c) are similar to those in Fig. 9(b), except as particle radius increases, compression of the particles on the workpiece material becomes stronger and the stress distributes deeper into the internal structure of the workpiece. Unlike the other two situations, with a particle radius of 5 nm, more stress occurs in front of the chip-accumulation region, as shown in Fig. 9(a). This results from the different ways chips are accumulated. Fig. 9 Stress diagrams at the 70 th load step with particle radius of (a) 5 nm, (b) 10 nm, and (c) 20 nm.

Cutting force analysis
In the QC simulation, as the particle is in uniform motion, the interaction force between the abrasive particle and workpiece is equal to the cutting force. This value can be obtained by summing the forces of all atoms that stress the particle. The cutting force is denoted as N is the number of all atoms in a particle, and ( ) f k is the force stressed by the th k atom. Figure 10 is the force-displacement plot when cutting particles of different radiuses, where the abscissa is the displacement of particles, and the ordinate is the component of the cutting force along the cutting direction.
As shown in Fig. 10, the tangential force increased rapidly with the rising cutting displacement at the initial stage of cutting. Then, the tangential force became relatively constant. The twists, turns, and sudden drops resulted from the internal dislocation nucleation of the workpiece. The average tangential cutting forces of 5 nm, 10 nm, and 20 nm particle sizes were 9.78 nN, 11.40 nN, and 8.81 nN, respectively. When the displacement was between 40-100 Å, the calculated average cutting forces of 5 nm, 10 nm, and 20 nm particle sizes were 6.75 nN, 9.36 nN, and 7.87 nN, respectively. Therefore, we could see that the change of the particle radius had little influence on the tangential cutting force.

Conclusions
An abrasive material removal mechanism in the copper CMP process with different particle sizes was studied in this paper. The simulation results are summarized as follows: (1) Internal materials of the workpiece generated a shear band at approximately 45° along the cutting direction, and dislocation and sliding phenomena occurred.
(2) Chips obviously accumulated when the abrasive particles were small, and the quality of the workpiece surface was poor. Higher quality of the workpiece surface was achieved when the particles were relatively large; however, larger plastic deformation was caused inside the workpiece.
(3) The residual stresses inside the workpiece were concentrated mainly on the following regions: chipaccumulating regions along the shear band, in front of the particles, and to the left end of the processed work piece. As the particle size increased, the stress distribution deepened. (4) The change of particle radius had little influence on the tangential cutting force.