Surface Morphology Evolution during Chemical Mechanical Polishing Based on Microscale Material Removal Modeling for Monocrystalline Silicon

Chemical–mechanical polishing (CMP) is widely adopted as a key bridge between fine rotation grinding and ion beam figuring in super-smooth monocrystalline silicon mirror manufacturing. However, controlling mid- to short-spatial-period errors during CMP is a challenge owing to the complex chemical–mechanical material removal process during surface morphology formation. In this study, the nature of chemical and mechanical material removal during CMP is theoretically studied based on a three-system elastic–plastic model and wet chemical etching behavior. The effect of the applied load, material properties, abrasive size distribution, and chemical reaction rate on the polishing surface morphology is evaluated. A microscale material removal model is established to numerically predict the silicon surface morphology and to explain the surface roughness evolution and the source of nanoscale intrinsic polishing scratches. The simulated surface morphology is consistent with the experimental results obtained by using the same polishing parameters tested by employing profilometry and atomic force microscopy. The PSD curve for both simulated surface and experimental results by profilometry and atomic force microscopy follows linear relation with double-logarithmic coordinates. This model can be used to adjust the polishing parameters for surface quality optimization, which facilitates CMP manufacturing.


Introduction
Owing to its excellent mid-to short-spatial-period error control with a root-meansquare (RMS) roughness value in the subnanometer level [1], chemical-mechanical polishing (CMP) has been applied for monocrystalline silicon mirror fabrication in the aerospace industry and high energy beam system domains [2][3][4][5]. Microscale surface morphology is a direct source for silicon mirror evaluations of mid-to short-spatial-period errors during CMP. Optimization of the silicon surface morphology requires precise control of surface features, such as roughness and microscratches, which are determined based on material removal characteristics.
Preston's equation is widely accepted for describing macrolevel material removal, where the material removal rate (MRR) is the product of Preston's constant, the applied load, and the workpiece velocity. However, several aspects of CMP, such as slurry hydrodynamics and microlevel surface features, cannot be explained by using this linear equation. Thus, efforts have been made to expand Preston's equation into theoretical, nonlinear forms. Runnels [6] proposed a hydrodynamic erosion model that considers the fluid layer shear stress by using the steady-state two-dimensional (2D) Navier-Stokes equation. This feature-scale model linked the wafer and abrasive particles together through the layer thickness and erosion law, with good agreement with the experimental erosion

Mechanical Material Removal by Plastic Deformation
During CMP, free abrasives with diameters 10-100 nm in a polishing slurry enter the gap between the polishing pad and workpiece ( Figure 1) under a dynamic load of 0.1-10 µN per particle [11,13]. Most adopted pad materials, such as pitch and polyurethane, exhibit elastic behavior. Contact abrasives generate recoverable elastic indentations within the elastic limit. For brittle workpieces, such as silicon and carbide silicon, both brittle and ductile abrasion exist [22,23]. Sphere-shaped abrasives create nanoscale scratches by plastic deformation with lateral cracks ( Figure 2); this is referred to as the mechanical material removal process in CMP [24]. This microscratching or microcracking process by a high number of free abrasives is similar to a dynamic blunt indent on the workpiece [13,23]. The typical nanoscale scratch depth generated in this process is 0.1-1 nm for brittle materials such as silicon, which is 2-3 orders of magnitude smaller than the silica abrasive size. The final surface morphology of the silicon workpiece is partially determined by this scratching process ( Figure 3).

Micropolishing Model
The surface morphology of the polishing workpiece is determined by the chemical and mechanical material removal processes during CMP. In this section, we discuss the monocrystalline silicon material removal characteristics of these two processes based on the Hertzian contact [14,15] and wet chemical etching models. A micropolishing model combining the two models is proposed to quantitatively explain the roles of the ASD, number density, applied load, and chemical etching rate, together with workpiece, abrasive, and pad characteristics.

Mechanical Material Removal by Plastic Deformation
During CMP, free abrasives with diameters 10-100 nm in a polishing slurry enter the gap between the polishing pad and workpiece ( Figure 1) under a dynamic load of 0.1-10 μN per particle [11,13]. Most adopted pad materials, such as pitch and polyurethane, exhibit elastic behavior. Contact abrasives generate recoverable elastic indentations within the elastic limit. For brittle workpieces, such as silicon and carbide silicon, both brittle and ductile abrasion exist [22,23]. Sphere-shaped abrasives create nanoscale scratches by plastic deformation with lateral cracks ( Figure 2); this is referred to as the mechanical material removal process in CMP [24]. This microscratching or microcracking process by a high number of free abrasives is similar to a dynamic blunt indent on the workpiece [13,23]. The typical nanoscale scratch depth generated in this process is 0.1-1 nm for brittle materials such as silicon, which is 2-3 orders of magnitude smaller than the silica abrasive size. The final surface morphology of the silicon workpiece is partially determined by this scratching process ( Figure 3).  Schematic of plastic deformation and lateral cracks on the silicon workpiece by single abrasive at the pad-workpiece interface; di and di' denote the workpiece plastic and pad elastic deformation depths, respectively, h is the gap between the pad and workpiece, is the contact area radius for the abrasive and workpiece, and is the radius of the abrasive particle.

Figure 2.
Schematic of plastic deformation and lateral cracks on the silicon workpiece by single abrasive at the pad-workpiece interface; di and di' denote the workpiece plastic and pad elastic deformation depths, respectively, h is the gap between the pad and workpiece, a i is the contact area radius for the abrasive and workpiece, and r i is the radius of the abrasive particle. Figure 2. Schematic of plastic deformation and lateral cracks on the silicon workpiece by single abrasive at the pad-workpiece interface; di and di' denote the workpiece plastic and pad elastic deformation depths, respectively, h is the gap between the pad and workpiece, is the contact area radius for the abrasive and workpiece, and is the radius of the abrasive particle. To simulate the evolution of the surface morphology during CMP, a reliable scratch geometry is required. The actual deformation for the plastic-elastic behavior during scratching includes elastic, plastic, and densification deformation. The last two items are the permanent deformation on the silicon surface. In the mechanic material removal model, we focused on the relationship between applied load and the morphology of the scratch. To simplify the problem, scratching was assumed to be the accumulation of Hertzian cracks along the movement direction of abrasive particles [14,15,19,20]. First, we considered single abrasive particles in contact with the pad and the workpiece. To simplify the plastic-elastic behavior for a pad-abrasive-workpiece three-body system, we adopted Suratwala's [15] hypothesis that the formation of applied load is determined by elastic contact mechanics. According to Johnson [14], the combined modulus of the system satisfies [14,15]  To simulate the evolution of the surface morphology during CMP, a reliable scratch geometry is required. The actual deformation for the plastic-elastic behavior during scratching includes elastic, plastic, and densification deformation. The last two items are the permanent deformation on the silicon surface. In the mechanic material removal model, we focused on the relationship between applied load and the morphology of the scratch. To simplify the problem, scratching was assumed to be the accumulation of Hertzian cracks along the movement direction of abrasive particles [14,15,19,20]. First, we considered single abrasive particles in contact with the pad and the workpiece. To simplify the plastic-elastic behavior for a pad-abrasive-workpiece three-body system, we adopted Suratwala's [15] hypothesis that the formation of applied load is determined by elastic contact mechanics. According to Johnson [14], the combined modulus E e f f of the system satisfies [14,15] where d i represents the penetration depth of the abrasive. The plastic deformation depth d i m−p satisfies [14,15] where d i is approximately three orders of magnitude smaller than a i during CMP; therefore, Equation (5) can be rewritten as By substituting Equation (7) into Equation (4), we can obtain the applied load as a function of the abrasive size: For free abrasives, the total applied load is the sum of every abrasive participating in material removal, or more generally, the integration of the product of the applied load and number density within each abrasive size range, as follows: where N i is the number of abrasives of radius r i , f (r i ) is the probability density for the ASD curve of radius r i , and N t is the total number of abrasives at the workpiece-pad interface, which is determined by the slurry concentration. For the polishing slurry, the ASD can be fitted by a logarithmic normal distribution in the following form [9,13]: where r 0 denotes the average abrasive size and σ denotes the half width of the ASD. The abrasives should directly contact the pad and workpiece to form plastic deformation, indicating that the size of the abrasives should be larger than the gap between the pad and workpiece. In addition, larger particles were not permitted to enter the gap in our experiments, indicating a maximum abrasive size restriction. Here, we assumed that the size of the effective abrasives should satisfy the following relationship: where h is the gap between the workpiece and pad, and r 0 and σ denote the mean value and standard deviation from the ASD fitting result, respectively. Abrasives with radii beyond this range cannot create microscratches and do not contribute to the applied load; hence, they can be removed from the integral. Equations (8)- (11) can be combined to rewrite P total as a function of h: From the simple elastic load balance by the pad and workpiece over the abrasive, we obtained the gap value, contact radius, abrasive penetration depth, and plastic deformation depth corresponding to the experiment. This formulation is based on a ductile-regime polishing hypothesis proved by Bifano [23] to be capable of brittle materials such as silicon and silicon carbide, thereby is also adaptable to ductile materials such as fused silica and soda-lime glass.

Chemical Material Removal by Simplified Wet Chemical Etching
Chemical material removal is another basic material removal process that occurs during CMP. In an alkaline atmosphere, which is widely applied for silicon mirror CMP, silicon is considered to be removed through the following reactions [26]: and Silicon atoms at the workpiece surface are oxidized by hydroxyl ions in the aqueous slurry and removed from the bulk material. This process is similar to the wet chemical etching of silicon in KOH [27,28], except that the reaction rate is considerably lower during CMP. Although chemical reaction is not the major part of material removal, this full-plane erosion helps expose the grain boundaries formed during formation, thereby causing orange peel structures on the silicon surface [29][30][31]. The different binding energies of the crystal planes, surfactants with micellization, and chelation make this chemical process almost impossible for precise quantification. To simplify the chemical material removal model, two assumptions were made: (1) monocrystalline silicon has a similar reaction preference as the isotropic substances. The mechanical component of the polishing rate plays a dominant role in determining the polishing rate in the boundary lubrication regime, as is our case [32]. No scratch orientation preference during CMP has yet been reported, proving that the anisotropic characteristic of silicon can be ignored for simplification of the chemical reaction during CMP; (2) the reaction rate at the scratch bottom is the same as that of bulk materials. The simulation error for this assumption can easily be eliminated by linear fitting of the expansion ratio for the RMS roughness or peak-valley value. The trends of the surface morphology evolution and salient features are not influenced.
For hydrofluoric acid-based etching on fused silica, Feit [33] described the morphological evolution of two neighboring cracks by using a surface area-volume model. The width of the parabolic-shaped crack increased with the square root of the etching time before the intersection (Equation (15)). When cracks or scratches intersect with neighboring cracks, the crack widening rate decreases because the width of the coalescing cracks is twice as large. The etching widths of the two neighboring parabolic cracks are shown in Figure 4. The scratch width function for the silicon workpiece could be obtained by adjusting the etching rate r e , the coefficient 2.8, and the initial width w 0 in Equation (15).
where w is the crack width, w 0 is the initial crack width, r e is the etching rate, t is the crack etching time, and the coefficient 2.8 is related to the initial crack shape in their model. We simulated the evolution of 3000 randomly selected scratches on an ideal surface by using Feit's model. The computation period for one iteration over this scratched surface was 1055 ms. Because scratches accumulated over repeated iterations, the time cost rapidly increased over the first few iterations; this is not suitable for MRR simulations during CMP.
To reduce the storage and time cost, a rapid scratch morphology generation method is needed. We proposed a Gaussian blur and smoothing method [35] for scratch widening during CMP. The MATLAB 'fspecial' and 'imfilter' command predefined a rotationally symmetric Gaussian image filter, in which the Gaussian function in 2D satisfies We simulated the evolution of 3000 randomly selected scratches on an ideal s by using Feit's model. The computation period for one iteration over this scratched s was 1055 ms. Because scratches accumulated over repeated iterations, the time co idly increased over the first few iterations; this is not suitable for MRR simulations CMP. To reduce the storage and time cost, a rapid scratch morphology generation m is needed. We proposed a Gaussian blur and smoothing method [35] for scratch wid during CMP. The MATLAB 'fspecial' and 'imfilter' command predefined a rotat symmetric Gaussian image filter, in which the Gaussian function in 2D satisfies The relationship between the scratch width and etching time was establish curve fitting with different Gaussian parameters , , and . The thin lines sho Gaussian fitting results in 2D, which correlate well with the thick lines obtained by the surface area-volume model, as displayed in Figure 4. R-square for fitting is clo for the etching time 1-5 h, which means a Gaussian function is qualified to d scratches during CMP. The eroding profile can easily be obtained by fitting an Equation (16) through time instead of calculating depth ( ) and width ( ) for each s For 3000 randomly selected scratches, the computation period was 22 ms, which times less than that of the surface area-volume model. Using an image filter with a ian blur kernel [35] predefined by the fitting results, we can simulate the scratch evolution in 3D.
To define the simulated surface in numerical values, surface parameters that a correlated to microscale surface morphology should be adopted. The simulated s morphology consists of two parts, scratches and orange-peel texture. RMS rough chosen here because of its adaptability in characterizing orange peel mentioned in giani and Miranda-Medina's works [31,36]; it is also commonly adopted for engin standards. Figure 5 illustrates the evolution of the surface RMS roughness for lo with and without Gaussian blur by 10 or 20 randomly selected abrasive particles pe The relationship between the scratch width and etching time was established by curve fitting with different Gaussian parameters a, b, and c. The thin lines show the Gaussian fitting results in 2D, which correlate well with the thick lines obtained by using the surface area-volume model, as displayed in Figure 4. R-square for fitting is close to 1 for the etching time 1-5 h, which means a Gaussian function is qualified to describe scratches during CMP. The eroding profile can easily be obtained by fitting a and c in Equation (16) through time instead of calculating depth (d) and width (w) for each scratch. For 3000 randomly selected scratches, the computation period was 22 ms, which was 50 times less than that of the surface area-volume model. Using an image filter with a Gaussian blur kernel [35] predefined by the fitting results, we can simulate the scratch width evolution in 3D.
To define the simulated surface in numerical values, surface parameters that are best correlated to microscale surface morphology should be adopted. The simulated surface morphology consists of two parts, scratches and orange-peel texture. RMS roughness is chosen here because of its adaptability in characterizing orange peel mentioned in Rebeggiani and Miranda-Medina's works [31,36]; it is also commonly adopted for engineering standards. Figure 5 illustrates the evolution of the surface RMS roughness for loop 150 with and without Gaussian blur by 10 or 20 randomly selected abrasive particles per loop. Despite slight vibrations due to random particle selection for each loop, the RMS roughness increased over the polishing time on an ideal surface for the four simulation conditions. The RMS roughness increased more slowly with Gaussian blur, which indicated that chemical material removal had an inhibitory effect on the surface during CMP. As the smoothing parameters were adjusted, the chemical effect gradually approached the true level. In addition, the RMS growth curve for 10 abrasives per loop without Gaussian blur had a growth rate similar to that of 20 abrasives per loop with Gaussian blur. This result implies that the abrasive number and chemical etching rate may have similar effects on the RMS roughness, even though they are associated with different material removal processes and surface features. true level. In addition, the RMS growth curve for 10 abrasives per loop without Ga blur had a growth rate similar to that of 20 abrasives per loop with Gaussian blu result implies that the abrasive number and chemical etching rate may have similar on the RMS roughness, even though they are associated with different material re processes and surface features. Simulated RMS roughness evolution with and without chemical material removal loops. The red and black solid lines are for 20 randomly selected particles, which are ra placed for each loop (3000 in total), and the blue and green solid lines are for 10 randomly s particles, which are randomly placed for each loop (1500 in total).

Micropolishing Model during CMP
The microscopic surface morphology was determined by the chemical and me cal material removal processes discussed in the previous section. The simplest ap to combine chemical and mechanical material removal processes in a micropo model is to find two individual parameters in the two processes that have similar on the surface morphology evolution. As shown in Figure 5 in Section 2.2, both the ical etching rate and abrasive number affect the RMS growth rate. The RMS rough linear to the standard deviation chosen for Gaussian blur, which is defined by the cal etching rate. To verify the influence of the number of abrasives, we tested th roughness evolution as a function of the number of abrasive particles. Chemical re prevailed over mechanical removal at an extremely low concentration of 10 partic loop sample (as shown in Figure 6). The RMS roughness changed slightly over t crolevel while a texture similar to an orange peel increased because of the exposu of grain boundaries rose, which worsened the surface quality. For a large concen (100 particles per loop), the mechanical removal surpassed the chemical remov croscratches are the main features of the silicon surface. The RMS roughness inf point appears after 3500 iterations for 100 particles per loop. For 50 particles per loo loops were sufficient to reach the inflection point. It was concluded that for a large concentration, the MRR increased with a decrease in the surface RMS roughness, was consistent with the experimental results. Comparing Figures 5 and 6, the ch etching rate and abrasive particle number changed the RMS roughness grow through time in a similar manner. By normalizing the chemical etching rate and ab

Micropolishing Model during CMP
The microscopic surface morphology was determined by the chemical and mechanical material removal processes discussed in the previous section. The simplest approach to combine chemical and mechanical material removal processes in a micropolishing model is to find two individual parameters in the two processes that have similar effects on the surface morphology evolution. As shown in Figure 5 in Section 2.2, both the chemical etching rate and abrasive number affect the RMS growth rate. The RMS roughness is linear to the standard deviation chosen for Gaussian blur, which is defined by the chemical etching rate. To verify the influence of the number of abrasives, we tested the RMS roughness evolution as a function of the number of abrasive particles. Chemical removal prevailed over mechanical removal at an extremely low concentration of 10 particles per loop sample (as shown in Figure 6). The RMS roughness changed slightly over the microlevel while a texture similar to an orange peel increased because of the exposure rate of grain boundaries rose, which worsened the surface quality. For a large concentration (100 particles per loop), the mechanical removal surpassed the chemical removal. Microscratches are the main features of the silicon surface. The RMS roughness inflection point appears after 3500 iterations for 100 particles per loop. For 50 particles per loop, 2500 loops were sufficient to reach the inflection point. It was concluded that for a large slurry concentration, the MRR increased with a decrease in the surface RMS roughness, which was consistent with the experimental results. Comparing Figures 5 and 6, the chemical etching rate and abrasive particle number changed the RMS roughness growth rate through time in a similar manner. By normalizing the chemical etching rate and abrasive particle number, the two material removal processes can be combined to describe the surface morphology evolution.
Materials 2022, 15, x FOR PEER REVIEW 9 of 18 particle number, the two material removal processes can be combined to describe the surface morphology evolution. First, Gaussian blur over the entire surface was adopted for continuous scratching. The size of the symmetric Gaussian low-pass filter is linear to the contact circle radius of each effective abrasive for the initial scratch. The standard deviation of the filter is determined by its linear relationship with the square root of the deformation depth . The scratch evolution is shown in Figure 7 for 101 loops. The surface morphology evolution initially approached the main problem. However, a nonconvergence RMS roughness evolution with a slight vibration owing to random particle selection, as shown in Figure 8, was obtained after extended iterations. This finding reversed the experimental result, that for specific slurry ASD, a limited RMS was reached after sufficient iterations. Gaussian blur for the entire surface led to a change in the original Gaussian-type height distribution ( ), as follows: where and denote the distribution functions of the original surface height and Gaussian blur, respectively, and the standard deviations, and and the average values. Considering = 0 and = 0 for simplification, we obtain where denotes the standard deviation for the new surface height distribution. and satisfy First, Gaussian blur over the entire surface was adopted for continuous scratching. The size of the symmetric Gaussian low-pass filter is linear to the contact circle radius a i of each effective abrasive for the initial scratch. The standard deviation of the filter is determined by its linear relationship with the square root of the deformation depth d i . The scratch evolution is shown in Figure 7 for 101 loops.
Materials 2022, 15, x FOR PEER REVIEW 9 of 18 particle number, the two material removal processes can be combined to describe the surface morphology evolution. First, Gaussian blur over the entire surface was adopted for continuous scratching. The size of the symmetric Gaussian low-pass filter is linear to the contact circle radius of each effective abrasive for the initial scratch. The standard deviation of the filter is determined by its linear relationship with the square root of the deformation depth . The scratch evolution is shown in Figure 7 for 101 loops. The surface morphology evolution initially approached the main problem. However, a nonconvergence RMS roughness evolution with a slight vibration owing to random particle selection, as shown in Figure 8, was obtained after extended iterations. This finding reversed the experimental result, that for specific slurry ASD, a limited RMS was reached after sufficient iterations. Gaussian blur for the entire surface led to a change in the original Gaussian-type height distribution ( ), as follows: where and denote the distribution functions of the original surface height and Gaussian blur, respectively, and the standard deviations, and and the average values. Considering = 0 and = 0 for simplification, we obtain where denotes the standard deviation for the new surface height distribution. and satisfy The surface morphology evolution initially approached the main problem. However, a nonconvergence RMS roughness evolution with a slight vibration owing to random particle selection, as shown in Figure 8, was obtained after extended iterations. This finding reversed the experimental result, that for specific slurry ASD, a limited RMS was reached after sufficient iterations. Gaussian blur for the entire surface led to a change in the original Gaussian-type height distribution f (z), as follows: where f 1 and f 2 denote the distribution functions of the original surface height and Gaussian blur, respectively, σ 1 and σ 2 the standard deviations, and µ 1 and µ 2 the average values.
Considering µ 1 = 0 and µ 2 = 0 for simplification, we obtain where σ 3 denotes the standard deviation for the new surface height distribution. A and σ 3 satisfy and = √2 + A repeated standard deviation change from to with a stable led to a quasilinear RMS roughness increase after sufficient iterations; hence, a nonconvergence RMS roughness value for the entire surface. To obtain convergent RMS roughness, a Markov chain was applied to compute a Gaussian-type random height distribution series with stationary initial and target surface background RMS. The background height variation series was calculated using the Metropolis-Hastings method [37] for a normal distribution to satisfy the convergent RMS roughness behavior within the scratch-free position during polishing. Silica abrasive particles were randomly selected from the ASD curve and placed on a silicon surface within a 500 × 500 pixels area. The particles moved at the same velocity in a random direction. After the mechanical process, the scratched position was Gaussian blurred with the background morphology offered by the Markov chain to simulate the chemical material removal process.

Experimental Setup
Colloidal SiO2 (JN-30, Qingdao Jiyida, Qingdao, China) was adopted as the polishing slurry. The silica colloid was diluted to 15 wt.% by ultrapure (UP) water (conductivity 1.12 μS/cm). The slurries were magnetically stirred to avoid slurry particle agglomeration during polishing, which could possibly influence the ASD. The ASDs for the polishing slurry were measured by using dynamic light scattering (DLS) techniques (Zetasizer, Malvern Panalytical, Malvern, UK; sizes ranging from 0.3 to 10 4 nm). The results and fitting curves based on Equation (10) are shown in Figure 9. The ASD of the commercial colloidal silica was bimodal, with peaks at ~3 and 107 nm. A repeated standard deviation change from σ 1 to σ 3 with a stable σ 2 led to a quasilinear RMS roughness increase after sufficient iterations; hence, a nonconvergence RMS roughness value for the entire surface. To obtain convergent RMS roughness, a Markov chain was applied to compute a Gaussian-type random height distribution series with stationary initial and target surface background RMS. The background height variation series was calculated using the Metropolis-Hastings method [37] for a normal distribution to satisfy the convergent RMS roughness behavior within the scratch-free position during polishing. Silica abrasive particles were randomly selected from the ASD curve and placed on a silicon surface within a 500 × 500 pixels area. The particles moved at the same velocity in a random direction. After the mechanical process, the scratched position was Gaussian blurred with the background morphology offered by the Markov chain to simulate the chemical material removal process.

Experimental Setup
Colloidal SiO 2 (JN-30, Qingdao Jiyida, Qingdao, China) was adopted as the polishing slurry. The silica colloid was diluted to 15 wt.% by ultrapure (UP) water (conductivity 1.12 µS/cm). The slurries were magnetically stirred to avoid slurry particle agglomeration during polishing, which could possibly influence the ASD. The ASDs for the polishing slurry were measured by using dynamic light scattering (DLS) techniques (Zetasizer, Malvern Panalytical, Malvern, UK; sizes ranging from 0.3 to 10 4 nm). The results and fitting curves based on Equation (10) are shown in Figure 9. The ASD of the commercial colloidal silica was bimodal, with peaks at~3 and 107 nm. Monocrystalline Si<111> wafers (diameter: 30 mm; thickness: 5 mm) were used as the workpiece. Six pieces of silicon wafers were pitch-buttons blocked on an aluminum connector and coarsely ground using diamond abrasive particles (Microgrit 10, 14, and 28 T) on a pig iron lap to obtain a quasiflat coplanar surface ( Figure 10). For pad preparation, a polyurethane pad (LP66, Universal Photonics, New York, NY, USA) was attached to an aluminum polishing lap (diameter: 300 mm; thickness: 20 mm) by using hot melt glue and placed upside down on a flat marble slab to obtain a fully fit no-bubble surface. A 10 mm square pattern with 2 mm wide V-grooves having a depth of 3 mm was engraved on the pad for slurry flow. The workpieces were polished by using pure silica colloid slurry on a polyurethane pad with an applied load of 8 N and a slurry flow rate of 10 mL/min. After sufficient polishing time, the workpieces were rinsed with UP water. The polishing parameters used are listed in Table 1.  The surface morphology was tested using profilometry (50×, scan area of 640 × 480 pixels within 125 μm × 94 μm, narrow-band green light, ContourGT-X3, Bruker™, Billerica, MA, USA) and atomic force microscopy (AFM) (tapping mode, scan area of 256 × 256 pixels within 2 μm × 2 μm, AFM Dimension Icon, Bruker™, USA). The typical surface morphology results are shown in Figure 11, and the RMS roughness over five repetitions of the polishing experiments is listed in Table 2. The size of pixels for profilometry and AFM in the x-y direction is 20 μm and 7.8 nm, respectively. A typical scratch width in Monocrystalline Si<111> wafers (diameter: 30 mm; thickness: 5 mm) were used as the workpiece. Six pieces of silicon wafers were pitch-buttons blocked on an aluminum connector and coarsely ground using diamond abrasive particles (Microgrit 10, 14, and 28 T) on a pig iron lap to obtain a quasiflat coplanar surface ( Figure 10). For pad preparation, a polyurethane pad (LP66, Universal Photonics, New York, NY, USA) was attached to an aluminum polishing lap (diameter: 300 mm; thickness: 20 mm) by using hot melt glue and placed upside down on a flat marble slab to obtain a fully fit no-bubble surface. A 10 mm square pattern with 2 mm wide V-grooves having a depth of 3 mm was engraved on the pad for slurry flow. The workpieces were polished by using pure silica colloid slurry on a polyurethane pad with an applied load of 8 N and a slurry flow rate of 10 mL/min. After sufficient polishing time, the workpieces were rinsed with UP water. The polishing parameters used are listed in Table 1. Monocrystalline Si<111> wafers (diameter: 30 mm; thickness: 5 mm) were used as the workpiece. Six pieces of silicon wafers were pitch-buttons blocked on an aluminum connector and coarsely ground using diamond abrasive particles (Microgrit 10, 14, and 28 T) on a pig iron lap to obtain a quasiflat coplanar surface ( Figure 10). For pad preparation, a polyurethane pad (LP66, Universal Photonics, New York, NY, USA) was attached to an aluminum polishing lap (diameter: 300 mm; thickness: 20 mm) by using hot melt glue and placed upside down on a flat marble slab to obtain a fully fit no-bubble surface. A 10 mm square pattern with 2 mm wide V-grooves having a depth of 3 mm was engraved on the pad for slurry flow. The workpieces were polished by using pure silica colloid slurry on a polyurethane pad with an applied load of 8 N and a slurry flow rate of 10 mL/min. After sufficient polishing time, the workpieces were rinsed with UP water. The polishing parameters used are listed in Table 1.  The surface morphology was tested using profilometry (50×, scan area of 640 × 480 pixels within 125 μm × 94 μm, narrow-band green light, ContourGT-X3, Bruker™, Billerica, MA, USA) and atomic force microscopy (AFM) (tapping mode, scan area of 256 × 256 pixels within 2 μm × 2 μm, AFM Dimension Icon, Bruker™, USA). The typical surface morphology results are shown in Figure 11, and the RMS roughness over five repetitions of the polishing experiments is listed in Table 2. The size of pixels for profilometry and AFM in the x-y direction is 20 μm and 7.8 nm, respectively. A typical scratch width in  The surface morphology was tested using profilometry (50×, scan area of 640 × 480 pixels within 125 µm × 94 µm, narrow-band green light, ContourGT-X3, Bruker™, Billerica, MA, USA) and atomic force microscopy (AFM) (tapping mode, scan area of 256 × 256 pixels within 2 µm × 2 µm, AFM Dimension Icon, Bruker™, USA). The typical surface morphology results are shown in Figure 11, and the RMS roughness over five repetitions of the polishing experiments is listed in Table 2. The size of pixels for profilometry and AFM in the x-y direction is 20 µm and 7.8 nm, respectively. A typical scratch width in CMP is several nanometers, based on Section 2.1. It is much smaller than the pixel size of profilometry, so compared with AFM, the surface topography measured by profilometry is scratch-free. Additionally, both results show orange-peel structures regardless of the field of view, indicating that orange peel has fractal properties. The fractal surface analysis is further analyzed in Section 4.1.
Materials 2022, 15, x FOR PEER REVIEW 12 of 18 CMP is several nanometers, based on Section 2.1. It is much smaller than the pixel size of profilometry, so compared with AFM, the surface topography measured by profilometry is scratch-free. Additionally, both results show orange-peel structures regardless of the field of view, indicating that orange peel has fractal properties. The fractal surface analysis is further analyzed in Section 4.1. Figure 11. Si<111> surface morphology after CMP by performing profilometry within 125 μm × 94 μm (left) and AFM within 2 μm × 2 μm (right).

Verification of the Micropolishing Model and Surface Morphology Evolution
Although the experimental ASD curve is bimodal with two peak centers at 2.8 and 107.2 nm, most of the abrasives are located at the second peak. For simplicity, the first peak was ignored in the simulation. With an of 107.2 nm for the silica slurry adopted in the experimental setup ( Figure 9) and pad, abrasive, and workpiece parameters summarized in Table 3, the calculated applied load as a function of ℎ by using Equation (12) is shown in Figure 12. Before the load by the effective abrasive reached ~12 N, the gap decreased with an increase in the applied load. The maximum applied load was obtained at ~0.06 μm. A load larger than this limit causes direct contact with the workpiece pad; hence, the gap and nondirect contact surface areas covered by abrasives decreased. Because the slurry volume entering the gap is linear to the product of the gap and nondirect contact surface area, the number of abrasive particles in the gap decreased. Therefore, for loads larger than the limit, the load applied by the abrasive decreased with a decrease in the effective abrasive quantity.

Verification of the Micropolishing Model and Surface Morphology Evolution
Although the experimental ASD curve is bimodal with two peak centers at 2.8 and 107.2 nm, most of the abrasives are located at the second peak. For simplicity, the first peak was ignored in the simulation. With an r 0 of 107.2 nm for the silica slurry adopted in the experimental setup ( Figure 9) and pad, abrasive, and workpiece parameters summarized in Table 3, the calculated applied load as a function of h by using Equation (12) is shown in Figure 12. Before the load by the effective abrasive reached~12 N, the gap decreased with an increase in the applied load. The maximum applied load was obtained at~0.06 µm. A load larger than this limit causes direct contact with the workpiece pad; hence, the gap and nondirect contact surface areas covered by abrasives decreased. Because the slurry volume entering the gap is linear to the product of the gap and nondirect contact surface area, the number of abrasive particles in the gap decreased. Therefore, for loads larger than the limit, the load applied by the abrasive decreased with a decrease in the effective abrasive quantity.  Figure 12. Applied load as a function of gap calculated by using Equation (12). As described in Section 3, 8 N was adopted as total applied load, and the gap value was 0.1185 μm (marked by an orange circle), based on the curve.
The gap was 0.1185 μm for an applied load of 8 N. Gap ℎ is two orders of magnitude larger than the first ASD peak center, indicating that the abrasives located at the first ASD peak did not contribute to the mechanical MRR. The mechanical material removal caused by the microscratching process can be obtained by calculating the width and depth distribution of the scratches. The scratched position was then Gaussian blurred, as explained in Section 2.2, with the normalization method presented in Section 2.3, to simulate the total material removal after one pass. The surface morphology evolution and final surface morphology after a sufficient polishing time can be achieved by continuous iterations of the procedure.
To test the validity of the micropolishing model on the RMS roughness behavior, it was necessary to compare the simulated surface roughness with the experimental results under the same polishing conditions. The mechanical parameters adopted in the simulation for Si<111> are listed in Table 3. The simulated surface RMS roughness evolution was analyzed within 10,000 iterations ( Figure 13). An RMS roughness convergence of ~1.4 nm was observed after 5000 loops with a random fluctuation, which was the RMS roughness convergence point for the simulation. The RMS roughness over five repetitions of the polishing experiments was measured as 0.6-0.9 nm by employing profilometry and 0.2-0.4 nm by performing AFM, which has the same dimension as the simulated surface. During the calculation of penetration depth of slurry abrasive, the displacement difference at load/unload due to permanent plastic deformation was ignored, which may explain why the simulated RMS is slightly larger than those obtained by profilometry. Figure 12. Applied load as a function of gap calculated by using Equation (12). As described in Section 3, 8 N was adopted as total applied load, and the gap value was 0.1185 µm (marked by an orange circle), based on the curve.
The gap was 0.1185 µm for an applied load of 8 N. Gap h is two orders of magnitude larger than the first ASD peak center, indicating that the abrasives located at the first ASD peak did not contribute to the mechanical MRR. The mechanical material removal caused by the microscratching process can be obtained by calculating the width and depth distribution of the scratches. The scratched position was then Gaussian blurred, as explained in Section 2.2, with the normalization method presented in Section 2.3, to simulate the total material removal after one pass. The surface morphology evolution and final surface morphology after a sufficient polishing time can be achieved by continuous iterations of the procedure.
To test the validity of the micropolishing model on the RMS roughness behavior, it was necessary to compare the simulated surface roughness with the experimental results under the same polishing conditions. The mechanical parameters adopted in the simulation for Si<111> are listed in Table 3.
The simulated surface RMS roughness evolution was analyzed within 10,000 iterations ( Figure 13). An RMS roughness convergence of~1.4 nm was observed after 5000 loops with a random fluctuation, which was the RMS roughness convergence point for the simulation. The RMS roughness over five repetitions of the polishing experiments was measured as 0.6-0.9 nm by employing profilometry and 0.2-0.4 nm by performing AFM, which has the same dimension as the simulated surface. During the calculation of penetration depth of slurry abrasive, the displacement difference at load/unload due to permanent plastic deformation was ignored, which may explain why the simulated RMS is slightly larger than those obtained by profilometry.
The surface morphology after the RMS roughness convergence point is shown in Figure 14. The simulated surface morphology has features similar to the experimental results obtained by profilometry and AFM, as shown in Figure 11. Scratches and textures similar to orange peel, owing to mechanical scratching and exposure of grain boundaries formed during annealing, can be clearly observed in both experiments.  The simulation and experiment targeted different spatial periods. A direct comparison of the surface height information was not appropriate for morphology evaluation. This was proved by the RMS roughness deviation with profilometry and AFM for the same silicon workpiece ( Figure 11). Additionally, surface uniformity is another fundamental consideration in evaluating the correlation between simulation and experiment. The PSD theory [38] was thereby adopted for morphology uniformity in different spatial periods. The 1D PSD of an isotropic fractal surface obeys the inverse power law by using fractal surface theory: where denotes the spectral density, denotes the PSD value for each spatial frequency , and (1) is the PSD value for = 1. The log-log plot for the PSD is linear with a slope value − .   The simulation and experiment targeted different spatial periods. A direct comparison of the surface height information was not appropriate for morphology evaluation. This was proved by the RMS roughness deviation with profilometry and AFM for the same silicon workpiece ( Figure 11). Additionally, surface uniformity is another fundamental consideration in evaluating the correlation between simulation and experiment. The PSD theory [38] was thereby adopted for morphology uniformity in different spatial periods. The 1D PSD of an isotropic fractal surface obeys the inverse power law by using fractal surface theory: where denotes the spectral density, denotes the PSD value for each spatial frequency , and (1) is the PSD value for = 1. The log-log plot for the PSD is linear with a slope value − . The simulation and experiment targeted different spatial periods. A direct comparison of the surface height information was not appropriate for morphology evaluation. This was proved by the RMS roughness deviation with profilometry and AFM for the same silicon workpiece ( Figure 11). Additionally, surface uniformity is another fundamental consideration in evaluating the correlation between simulation and experiment. The PSD theory [38] was thereby adopted for morphology uniformity in different spatial periods. The 1D PSD of an isotropic fractal surface obeys the inverse power law by using fractal surface theory: where K denotes the spectral density, S 1 denotes the PSD value for each spatial frequency f x , and S 1 (1) is the PSD value for f x = 1. The log-log plot for the PSD is linear with a slope value −n. One-dimensional-surface PSD curves from the profilometer, AFM, and simulation were calculated using PSD theory [38] and plotted in Figure 15. The adopted simulating area was 13.25 µm × 13.25 µm within 500 × 500 sampling points, and the PSD curve covered spatial periods for the profilometer and AFM. The three spatial periods correlated well with the omission of frequency at both ends of the three curves, caused by the field of view, lateral resolution, and white noise [3,39,40]. Additionally, three PSD curves follow the same linear relation with double-logarithmic coordinates, which indicates the three surfaces are isotropic with the same fractal dimension. Linear fitting for the PSD curve of AFM follows: log(S 1 ) = −1.28 log( f x ) − 3.87, with R-square equaling 0.975. For profilometry and simulation, similar linear fitting functions were observed. According to the research by Chen et al. [41], we calculated the fractal dimension, which is~2.6. The PSD curve of the simulated surface provided an adjustable spatial frequency between 10 -1 µm -2 and 10 2 µm -1 by changing the ratio of the contact circle radius a i over the pixels or by increasing the number of pixels. Because of the consistency in the PSD curve, the PSD for other spatial frequency regions can be extrapolated according to the fitting result. Therefore, the pixel resolution can be enlarged for surface morphology predictions on different observation scales at any polishing time.
One-dimensional-surface PSD curves from the profilometer, AFM, and simulation were calculated using PSD theory [38] and plotted in Figure 15. The adopted simulating area was 13.25 μm × 13.25 μm within 500 × 500 sampling points, and the PSD curve covered spatial periods for the profilometer and AFM. The three spatial periods correlated well with the omission of frequency at both ends of the three curves, caused by the field of view, lateral resolution, and white noise [3,39,40]. Additionally, three PSD curves follow the same linear relation with double-logarithmic coordinates, which indicates the three surfaces are isotropic with the same fractal dimension. Linear fitting for the PSD curve of AFM follows: log( ) = −1.28 log( ) − 3.87, with R-square equaling 0.975. For profilometry and simulation, similar linear fitting functions were observed. According to the research by Chen et al. [41], we calculated the fractal dimension, which is ~2.6. The PSD curve of the simulated surface provided an adjustable spatial frequency between 10 -1 μm -2 and 10 2 μm -1 by changing the ratio of the contact circle radius over the pixels or by increasing the number of pixels. Because of the consistency in the PSD curve, the PSD for other spatial frequency regions can be extrapolated according to the fitting result. Therefore, the pixel resolution can be enlarged for surface morphology predictions on different observation scales at any polishing time.

Engineering Application of the Model
The ultrasmooth monocrystalline silicon mirror fabrication is a promising application for CMP. For investigations on CMP, high polishing efficiency and surface quality are often the most concerning issues. Although traditional experiment-based process improvements can achieve good polishing accuracy, a lack of knowledge of the material removal mechanism leads to long-term attempts and is often costly. In addition, the uncertainty of end-process surface morphology limits the application of CMP toward deterministic fabrication. In this paper, the effect of applied load, elastic-plastic behavior of materials, ASD, abrasive concentration, and chemical reaction rate determined by slurry components is evaluated theoretically. The micropolishing model established combines chemical erosion with mechanical scratching to predict morphology evolution on the microscale silicon surface. The simulated surface morphology is adaptable to various surface quality assessment standards because it contains the same information as the actual silicon surface. The model also provides openness with different workpiece and pad materials, provided their elasticity and plasticity match the assumption of the model.

Engineering Application of the Model
The ultrasmooth monocrystalline silicon mirror fabrication is a promising application for CMP. For investigations on CMP, high polishing efficiency and surface quality are often the most concerning issues. Although traditional experiment-based process improvements can achieve good polishing accuracy, a lack of knowledge of the material removal mechanism leads to long-term attempts and is often costly. In addition, the uncertainty of end-process surface morphology limits the application of CMP toward deterministic fabrication. In this paper, the effect of applied load, elastic-plastic behavior of materials, ASD, abrasive concentration, and chemical reaction rate determined by slurry components is evaluated theoretically. The micropolishing model established combines chemical erosion with mechanical scratching to predict morphology evolution on the microscale silicon surface. The simulated surface morphology is adaptable to various surface quality assessment standards because it contains the same information as the actual silicon surface. The model also provides openness with different workpiece and pad materials, provided their elasticity and plasticity match the assumption of the model.
In addition to end-process surface morphology predictions, the model can also be used to balance the polishing accuracy, efficiency, and processing cost. By analyzing Figures 5 and 6, a higher slurry concentration or chemical reaction rate leads to a higher polishing efficiency before abrasive saturation. Increasing the applied load is another approach to achieve high efficiency, but the roughness worsens as scratches are more profound on the workpiece. An even higher applied load brings direct contact with the workpiece pad, which may cause severe damage to silicon mirror and polishing equipment. With the simulated applied load under different gaps, which in our case is the curve shown in Figure 12, maximum load on workpiece can be fixed to avoid this damage.

Conclusions
In this study, chemical and mechanical processes were studied by using elastic-plastic deformation and wet chemical etching based on microscopic material removal characteristics during CMP for a silicon workpiece. A micropolishing model was developed to predict microscale surface morphology during CMP. The predicted silicon mirror surface morphology captured salient features, such as microscratches and textures similar to orange peel, with the PSD curve and RMS roughness results being consistent with the experimental results measured by profilometry and AFM. The results imply valuable insights into the mechanism and prediction of microscale surface morphology evolution. These insights include: (1) The gap between the workpiece and pad was found to depend on the abrasive size distribution and applied load, which in turn determined the size of the abrasive that can generate effective mechanical removal. (2) The chemical reaction during CMP, which was previously considered to assist mechanical removal or to have a similar effect to mechanical material removal, was quantified based on scratch widening from the Gaussian-fitting wet chemical etch model and time-dependent background height transition. (3) Chemical and mechanical material removal processes were linked through the abrasive particle number and chemical etching rate to establish a micropolishing model with convergent RMS roughness evolution. (4) The PSD curves calculated from simulation and experiments indicate that the surface is fractal in mid-to short spatial frequency, which means the material removal mechanism in this region is consistent during CMP. The increased understanding of the CMP mechanism obtained from this model can be used for the impact quantification of different polishing factors and further optimization of the polishing process.