Multiscale Simulation of Surface Defect Influence in Nanoindentation by a Quasi-Continuum Method

Microscopic properties of crystal aluminum thin film have been investigated using the quasi-continuum method in order to study the influence of surface defects in nanoindentation. Various distances between the surface pit defect and indenter and various sizes of the pit have been calculated. In this simulation, as the distance between the pit and indenter increases, the nanohardness increases in a wave that goes up in a period of three atoms, and it is found closely related to the crystal structure of periodic atom arrangement on {1 1 1} atomic close-packed planes of FCC metal; there is almost no influence on the nanohardness when the adjacent distance between the pit and indenter is more than 16 atomic spacing. We have modified the theoretical equation of the necessary load for elastic-to-plastic transition of Al film with the initial surface defect size. Furthermore, when the size coefficient of width (of height) equals about one unit (half unit), the yield load experiences an obvious drop. When it reaches about two units (one unit), the yield load is nearly close to that of the nanoindentation on a stepped surface. Addition-ally, compared to the width, the height of surface pit defect displays a greater influence on the yield load of thin film.


Introduction
Nanoindentation [1], relatively simple and effective, has already been used as a standard technique for evaluating mechanical properties of thin films [2][3][4][5][6][7]. In the recent years, a number of scientists have focused on thin films with defects through simulations and experiment [8][9][10][11]. Yu has analyzed the effects of nanocavity on nanoindentation, which is one kind of defect [12]. Surface roughness, grouped by the pits and steps, has already been a popular topic, and lots of nanoindentation simulations on a step have been studied [13][14][15][16][17][18][19][20]. It is well-known that the surface pits are very common in polycrystalline surfaces on microchips or microelectromechanical systems (MEMS). Therefore, it is necessary and significant to observe the nanoindentation on the pitted surface. Ni et al. [21] have compared nanoindentation of Al thin film with and without surface defects by multiscale simulation, and Zhang et al. [22] have probed the delay effect of dislocation nucleation in nanoindentation due to the surface pit defect. However, the distance effect of the pit and the size effect on elastic-plastic transition, which is especially important to thin film performance in microchips and MEMS, have not been taken into account yet. Now, this chapter is further to study and represent the distance effect [23] and the size effect of the pit [24] on nanohardness by quasi-continuum (QC) method [25], which is a relatively effective way to investigate large-scale model, where the molecular dynamics (MD) is mostly limited by the model scale.

Distance effect of surface pit defect 2.1 Methodology
The quasi-continuum (QC) method [26] is an effective mixed continuum and atomistic approach for simulating the mechanical response, especially in large-scale materials. The Ercolessi-Adams potential, which is one of the EAM potentials [27][28][29][30], is used in this QC method, in order to describe the atomistic behavior. The parameters for Al in this potential are as shown in Table 1.
The schematic illustration of nanoindentation model with a pit defect and its unit cell model of Al in the selected directions is as shown in Figure 1, where the width of the rigid rectangular indenter is 0.932 nm (four times the lattice constant of Al in [1 1 1] direction (d 0 )). Such indenter width is chosen based on the simulation example "Nano-Indentation by a Square Punch" in QC tutorial document [32] and others' works [21,22]. It is convenient to simulate and analyze that if the indenter shape is set rectangular, result from that the boundary of energy field (displacement field) and the distance between the pit and the indenter will not changed during the indenter being driven down into the 110 À Á surface in this simulation, and that is exactly necessary to probe the distance effect of the pit. The width D and depth H of the pit in the model are 0.688 and 0.730 nm, respectively. Such pit size is chosen relatively moderate and proper, more sensitive to the distance effect The experimental values extrapolated to T = 0 K are C 11 = 118.0 GPa, C 12 = 62.4 GPa, and C 44 = 32.5 GPa [27]. based on our previous works [21,22] and some pre-simulations. The thickness of this model in the out-of-plane direction is 0.4938 nm, which equals to the minimal repeat size applying the periodic boundary condition. Fifteen different d (defined as the adjacent distance between the pit and indenter) have been calculated in this chapter, respectively, 1d 0 , 2d 0 , 3d 0 , 4d 0 , 5d 0 , 6d 0 , 7d 0 , 8d 0 , 9d 0 , 10d 0 , 11d 0 , 12d0, 13d 0 , 17d 0 , and 21d 0 , in order to make a more comprehensive investigation. The parameter "PROXFACT" in QC method is applied to judge the range of nonlocal effects. The repatom is made nonlocal when any element is in its range of PROXFACT*rcut, where rcut is defined as the atomistic potential cutoff radius. Nonlocality is triggered if max a, b; k λ a k À λ b k >epscr (k = 1, 2, 3), where λ a k and λ b k are the eigenvalues of right Cauchy-Green deformation tensor C in two elements a and b, and the factor "epscr" is applied to judge whether a repatom has to be made nonlocal because of the significant variations in the deformation gradients around the repatom.
For normal atomistic modeling standards, the dimensions of this simulation thin film is quite large with approximately 1.3 million atoms or 4 million degrees of freedom (as shown 0.1 μm in height and 0.2 μm in width in Figure 1). Fortunately, QC method applies the molecular dynamics model at the intense deformation region and a finite element model elsewhere in order to reduce the degrees of freedom without losing atomistic details, where only 4000 atoms or 12,000 degrees of freedom have to be treated in this model by comparison, and can be easily finished in a few days through personal computer.

Nanoindentation without defect
Nanoindentation without defect is necessary to be studied for comparison, and the load-displacement curve of the nanoindentation on a defect-free surface is shown in Figure 2, where load (N/m) is presented as per unit length of indenter in the out-of-plane direction. It can be easily found out that the load curve gradually To find out the reason of such abrupt load decline in AB segment, the atoms structure and out-of-plane displacements have been probed and shown in Figure 3, where the step of 0.48 and 0.50 nm is, respectively, corresponding to the point A and point B in Figure 2.
Through Figure 3, we can make a conclusion here that the load reaches the critical value for dislocation emission at point A, which indicates the onset of the plastic stage. After that two Shockley partial dislocations are emitted at point B. Therefore the nanohardness of Al thin film without defect is 16.24 GPa, calculated by equation [33]: H ¼ P max A , where P max is the maximum value of load and A is the contact area of the indenter.  Due to the indenter width 0.932 nm in this simulation, the yield load is obtained approximately 15.14 N/m, which is smaller than 24.7 N/m acquired by Tadmor and Miller [31] with the indenter width 2.5 nm. It is reasonable that the paper [34] shows that as the indenter width decreases, the yield load would decline because of the requirement decline of the necessary strain energy. Figure 4 shows the nanohardness with various distances cases of the surface pit defects. Compared to the defect-free situation shown as a red horizontal, it indicates that the nanohardness of pitted surface has been declined. That is because a discontinuity at the boundary and the structure may lead to the reduction of stain energy storage when indenting. Additionally, when the adjacent distance between the pit and indenter (d) increases, the nanohardness increases in a wave that goes up in a period of three atoms (donated by circle in Figure 4) and finally tends to the case of nanoindentation on a defect-free surface.

Distance effect of the surface pit
To make a further probe, it is well-known that many physical properties mostly depend on the stacking patterns of atoms, such as cleavage, electronic band structure, and optical transparency [35]. Based on this simulation, the periodic arrangement of atoms "ABCABC" on {1 1 1} atomic close-packed planes of face-centered cubic metal (the illustration as shown in Figure 4) is exactly corresponding to increasing distances d on [1 1 1] direction. That is to say, when the pit moves each atom in the [1 1 1] direction away from the indenter, the strain energy at the pit surface on the {1 1 1} stacking fault energies (SFE) changes because of the proximity of the pit [36]. Consequently, such wave pattern associated with a cycle of three atoms is closely related to the crystal structure of periodic atom arrangement on {1 1 1} atomic close-packed planes of FCC metal.
A further discussion has been made in order to figure out the spatial extent of surface pit influence on nanohardness. It can be found out from Figure 4 that when the distance between the pit and indenter increases, the nanohardness gradually close to the nanohardness of defect-free case (16.24 GPa). If we set 1.5% to determine whether the nanohardness influence exists, it can be found out that when the adjacent distance (d) goes 16 atomic spacing far away from the indenter, there is almost no effect on nanohardness (as shown in Figure 4). Moreover, it can be predicted that each material has its critical value of such spatial extent of surface pit influence on nanohardness, which might have great significance to the size design of thin film in microchips without obvious reduction of the hardness.
However, it can be easily found out that the first three distances cases, respectively, d = 1d 0 , d = 2d 0 , and d = 3d 0 , do not match such wave pattern well. Atomic structure and corresponding strain distribution of these three cases have been further carried out to explain such unusual phenomenon.
Von Mises strain distribution and a strain comparison before and after the notch propagation in the distance cases of, respectively, 1d 0 , 2d 0 , and 3d 0 , are shown in Figure 5. It can be easily found out that a notch formed at the left side of surface pit in the distance cases of 1d 0 and 2d 0 , which actually induces serious damage to the structure of materials and severe strain concentration (as shown in Figure 5A-D), while it does not if the distance d equals 3d 0 (as shown in Figure 5E and F). That is to say, due to the great reduction of the nanohardness in the cases of 1d 0 and 2d 0 , the first three distance cases in nanohardness curve as shown in Figure 4 will not match the wave pattern that goes up in a period of three atoms.

Critical load for dislocation emission with initial surface pit
As is known to all that the pit influences the nanohardness is actually through the way of affecting nucleation and emission of the dislocation. It is necessary and significant to make a further probe on the critical load for elastic-to-plastic transition in the case of nanoindentation on the pitted surface according to the formula of defect-free model, where the formula to calculate the critical load for dislocation emission is carried out by Tadmor [31], shown as Eq. (1): where P cr is the critical load of dislocation emission, h is the depth of dislocation dipole emitted down beneath the indenter, a is the half width of indenter, γ 111 is the (1 1 1) surface energy of Al single crystal, and k is the slope of the load-displacement curve during the elastic section.
The critical load calculated by the simulation results of QC method and by Eq. (1) of dislocation theory has been displayed, respectively, as "QC data" and "theory load" in Table 2. It is necessary to note that the simulation data in the case of 1d 0 and 2d 0 is not suitable to be taken into account because of the notch propagation. The differential of the critical data between QC method and dislocation theory is fluctuant as the distance (d) changes. Based on the result that the nanohardness goes up in a period of three atoms, the critical load for elastic-toplastic is also in such periodicity. Thus, the correction form (set as Δ) might be reasonably defined as the following: where A(d) is the hardness reduction due to the surface pit and BÁSin(d) is specially set for the periodic arrangement of atoms. It is well recognized that the critical load of elastic-to-plastic transition will decrease [37], when the pit size (D, H as shown in Figure 1) increases. We use a dimensionless factor D a 1 Á H a 1 (dividing by the crystallographic lattice constant) to express the size influence of the pit, which has already been demonstrated reasonable in published article [24]. If the distance between the indenter and the pit is infinitely large, the influence on nanohardness can be almost ignored, and if the pit size increases, the rate of nanohardness change decreases with the distance variation. According to the function property and calculation formula designed by Tadmor [31], we take the form of ln 1 to express the distance effect of surface pit, considering that it is relevantly reasonable. Moreover, the affection of surface pit is closely due to the material property such as Burgers vector b ! , shear modulus μ, and Poisson ν. It is reasonable to apply μb 4π 1Àν ð Þ to express the influence of material property based on Eq. (1). In addition, the atomic periodical arrangement is actually three atoms "ABCABC" on {1 1 1} atomic close-packed planes of FCC metal.
Namely, a form of Sin 2π 3d 0 Á d þ φ is proper to express the periodicity of atom arrangement. Besides, we know that the unit of correction term (Δ) is exactly N/m. So, the correction can be defined as the following based on the discussion above: where α, β, and φ are three constants that need to be matched and fitted. Based on the simulation data in Table 2, these three constants α, β, and φ can be acquired by calculation approximately 3 2 , 2 15 , and À π 3 , respectively. So the critical load for the first dislocation emission of Al film has been revised with initial surface pit as follows:  Table 2.
The comparison of critical load between QC method and dislocation theory. Figure 6 shows the comparison of the critical load for dislocation emission of Al thin film in different distance cases calculated by the theoretical formula before and after modification. Though there is no parameter d in Eq. (1), the curves with blocks are calculated by depth h corresponding each distance case in this simulation. The simulation QC data is closer to the theoretical results which are calculated by Eq. (4) after modification. That is to say, such modification to the theoretical formula is justified as the pit size and the distance between the pit and indenter have both been taken into account.
The modified formula displays the decreasing trend of nanohardness as the distance between the pit and indenter increases, which quite agrees with the experimental results of nanoindentation on the surface step with different distances [14]. Further, such study might be referential to the research of material properties with defects, especially in microchips and MEMS.
3. Size effect of surface pit defect

Method and model
The quasi-continuum method is adopted in this simulation, which is one of the multiscale approaches that keeps an atomistic description at highly deformed regions, whereas a linear elastic continuum method is implemented far away from this dislocation core. In this simulation, the Ercolessi-Adams potential (EAM) [38] is also applied to describe the atomistic behavior of the system.
The nanoindentation model used in the simulation has been shown in Figure 7, where the x-axis direction is [1 1 1] direction, the y-axis direction is [1 1 0], and the outer-of-plane z direction is [1 1 2] direction. And its schematic of local and nonlocal representative atoms with initial surface pit defect has been shown in    and the pit in this model is chosen to be 4d 0 (as shown in Figure 7), which is proved to be reasonable referring to the research on nanoindentation on a stepped surface [14], where the spatial extent of the step's influence has been figured out to be approximately three times the contact radius. Therefore, the contact radius of this simulation is 2d 0 (half of the indenter width), and the distance between the left side of surface pit defect and the center of the indenter is exactly 6d 0 (4d 0 +4d 0 /2 = 6d 0 ), which quite agrees with this reference [14] (6d 0 /2d 0 = 3). Although the indentation tip shape is different between this simulation and the nanoindentation experiment, where it is square and round-like, respectively, actually, the nanohardness are both calculated by the real contact radius. That is to say, the change of the nanohardness in these two cases is quite similar if considering the ratio of defect distance to the real contact radius (as well as the experiment discussion). The simulation models of the width effect and height effect of surface pit defect are shown as Figure 7, and relevant parameter values are shown in Table 3 in order to make a more comprehensive investigation. Besides, the comparison model of nanoindentation on the stepped surface step with the height of the step 10h 0 has also been carried out as shown in Figure 7(d). Further, this model keeps the boundary condition rigid at the bottom and free at the sides, and the thickness is equal to the minimal repeat distance. The displacement-imposed boundary condition forces the atoms under the indenter to move into the material gradually. Each load step of indentation has set 0.02 nm in order to be more proper and effective to catch the dislocation nucleation and mission. Moreover, the final depth is 1.2 nm, which ensures that the behavior in the vicinity of the indenter will not be affected by the far-field boundary conditions.

Width effect on the yield load due to the pit defect
As is known to all that the yield load of materials is one of the most important factors of the material properties, however it can be obviously affected by defects such as surface pit defect. Normally, the yield load can be easily obtained from the first peak load in the load-displacement curve, which suggests onset of the elasticto-plastic transition. In this chapter, we have taken 10 different widths of the pit from D = 1d 0 to 10d 0 with a fixed height H = 5h 0 , in order to investigate the width  effect of surface pit defect on the yield load. The change of the yield load curve as pit width has been revealed in Figure 9. Generally, a reduction tendency of the yield load of Al thin film with the pit defect displays result from more and more serious destruction to the atomic structure by the increase of the pit width. Further, the yield load experiences an extremely slow reduction when the pit width increases from D = 1d 0 to 7d 0 ; after that it obviously drops from 14.8 to 14.24 N/m when the pit width reaches 7d 0 . Then, the yield load decreases slowly again. It is necessary and significant to compare the nanoindentation on a stepped surface with H = 5h 0 as shown in Figure 7(b), where the pit width can be treated as infinitely large. The result is approximately 14.23 N/m, very close to the yield load in the case of D = 10d 0 (the red point in Figure 9). This implies that the yield load of Al thin film nearly equals the yield load value in the case of stepped surface when the pit width increases to 10d 0 .

Height effect on the yield load due to the pit defect
We have taken 10 different heights of the pit from H = 1h 0 to 10h 0 , with a fixed width D = 5d 0 , in order to investigate the height effect of surface pit defect on the yield load as shown in Figure 7(c). The change of the yield load curve as pit height has been revealed in Figure 10. Similarly, the yield load experiences an extremely slow reduction when the pit height increases from H = 1h 0 to 5h 0 ; after that it obviously drops from 14.79 to 14.14 N/m when the pit width reaches 6d 0 . Then, the yield load decreases slowly again.
In the same way, it is also necessary and significant to compare the nanoindentation on a stepped surface with H = 10h 0 as shown in Figure 7(d). The result is approximately 13.75 N/m (the red point in Figure 10), which is quite near the yield load 13.93 N/m in the case of H = 10h 0 .

The investigation of dislocation nucleation and the estimation of Peierls stress
We have carried out a further probe of atomic snapshot and corresponding outof-plane displacement plot, in order to explain the obvious drops of yield load (D = 7d 0 to 8d 0 segment in Figure 9, H = 5h 0 to 6h 0 segment in Figure 10). Taking the cases of D = 1d 0 in width effect simulation and H = 1h 0 in height effect simulation, for example, it can be easily found out from nucleated dislocations and UZ contours displayed in Figure 11 that two dissociated <1 1 0> edge dislocations are emitted beneath the indenter after nucleation during the thin film yielding. Moreover, the dislocated structure along with the out-of-plane displacements experienced by the atoms has also been displayed in Figure 11, with the dimension 0.1 nm, where a fingerprint of the dislocations has been clearly shown between the partials in the stacking fault regions. According to the structure of FCC metal, it can be easily found out that the dislocations are composed of 1/6 <1 1 2> Shockley partials. On the left, and on the right,   Comparing all these width and height effect cases in this simulation, we find that emission depth of dislocations changes due to the size of the pit. Take Figure 11(a) and (b), for example, the dislocation dipole similarly travels into bulk after nucleation at the load step of 0.5 nm; however, its center of the emission depth settles, respectively, at the depth of 5.2 and 6.08 nm. It might be predicted that the yield load of thin film in macroscopy corresponds to the emission depth of dislocation in microscopy.
Considering that Peierls stress is exactly the resisting force during the dislocation movement due to the lattice structure, all these emission depths of dislocations have been adopted as an equilibrium distance to further calculate the Peierls stress predicted by the EAM potential [38]. Except the lattice friction, there are two forces acting on the dislocation: (i) the Peach-Koehler force (F PK ) due to the indenter stress field driving the dislocation into bulk and (ii) the image force (F I ) pulling the dislocation up to the surface. The dislocation, which is forced by the sum of these two forces, escapes the attractive region and propagates into the bulk and is finally stopped by lattice friction. Consequently, the force on the dislocation will be balanced at the equilibrium depth by the lattice friction force that is due to the Peierls stress (σ p ) [31].
Shear stress beneath the indenter is necessary to be further obtained to calculate the Peach-Koehler force. In this simulation, the rectangular indenter is frictionless, applying to an elastic thin film occupying the lower half-plane. When y < 0, the shear stress in bipolar coordinates is [39] σ xy ¼ À Pr 2 sin θ π r 1 r 2 ð Þ 3=2 sin θ À where P is the indentation load. According to the coordinate system of 2a indentation contact (the width of indenter is 2a), as shown in Figure 12, at a depth h beneath the right indenter tip, there is r ¼ The Peach-Koehler force is where σis the applied stress tensor, b is the Burgers vector, and ℓ is the dislocation line vector.
The image force applying to the dislocation dipole with width d = 2a at the emission depth h can be displayed as follows: Peierls stress of all these cases of pit size has been calculated and plotted based on the discussion above, where Figures 13 and 14, respectively, show Peierls stress of width effect and height effect. In the case of width effect, Peierls stress maintains  around the value of 100 MPa with tiny fluctuation from D = 1d 0 to 7d 0 ; after that the Peierls stress displays a sudden obvious dropdown to about 70 MPa, which is quite similar compared with yield load curve in Figure 9. In the case of height effect, Peierls stress maintains around the value of 70 MPa with tiny fluctuation from H = 1h 0 to 5h 0 ; after that the Peierls stress displays a sudden obvious dropdown to about 50 MPa, which is also greatly in accordance with the yield load curve in Figure 10. The conclusion can be drawn that such obvious decline of yield load (D = 7d 0 to 8d 0 segment in Figure 9, H = 5h 0 to 6h 0 segment in Figure 10) is closely related to the severe reduction of the Peierls stress, suggesting that it is reasonable and effective to explain the variation of yield load through the Peierls stress.

Size coefficient
We make a further probe on the difference of turning point between width effect and height effect, corresponding to D = 7d 0 in the width effect simulation and H = 5h 0 in the height effect simulation. It can be predicted that the influence degree of width factor is different from the height factor. It is necessary to quantify the size effect of surface pit defect to explain the reason of these differences. It can also be easily recognized that if the distance between the pit and the indenter decreases, the influence would be much more severe on the hardness and yield load. Namely, controlling the same influence of the pit on the nanohardness, the larger size of the pit is required when the pit goes far away from the indenter. Therefore, we defined a size coefficient α as the following, which is dimensionless in order to explain the size effect of surface pit defect: where "L*" is the characteristic length of the pit, such as D in the width effect simulation or H in the height effect simulation and "d*" is the distance between the center of the indenter and the left boundary of the pit, namely, 6d 0 in this simulation.
According to the simulation result of width effect, the critical width value to make a sudden obvious drop of yield load is 7d 0 (as shown the point D = 7d 0 in Figure 9). Consequently, the size coefficient α is approximately 1.17 , as shown in the point D = 10d 0 in Figure 9, the yield load of thin film with surface pit defect nearly equals the one of nanoindentation with surface step as shown in the red point in Figure 9.
According to the simulation result of height effect, the critical height value to make a sudden obvious drop of yield load is 5h 0 (as shown the point H = 5h 0 in Figure 10). Consequently, the size coefficient α is approximately 0.51 , as shown in the point H = 10h 0 in Figure 10, the yield load of thin film with surface pit defect nearly equals the one of nanoindentation with surface step as shown in the red point in Figure 10.
By contrast, the size coefficient of height is approximately half of the one of width to boost the sudden decline of yield load, implying that the height of the pit has a greater influence on the yield load than the width.
Moreover, the change of yield load of thin film as the pit area has been plotted and shown in Figure 15, where we can easily find out that the slope of yield load curve by height increasing is larger than the one by width increasing. That is to say, the height increasing makes the yield load decrease faster. Besides, the yield load by width increasing is smaller than the one by height increasing during 5h 0 d 0 to 25h 0 d 0 segment of the pit area, which results from the height of the pit in width increase curve (as shown red curve in Figure 15) is larger than the other one (black curve). If the area increases over than 25h 0 d 0 , the yield load by height increase is smaller than the one by the width increase. It can be well explained that the height of the pit in the curve of height increase goes up over 6h 0 , while the one in the curve of width increase still maintains 5h 0 , suggesting that the height of the pit has played a more important role on yield load.

Conclusions
Through this whole chapter, we apply the QC method to study the influence of surface pit defect in nanoindentation. Fifteen distances of adjacent boundaries between the pit and indenter have been taken into account to study the distance effect, compared with the nanoindentation on defect-free surfaces, while various sizes of the pit have been taken into account to study the size effect, compared with the nanoindentation on the stepped surfaces. Based on the discussion, we can make some conclusions as the following: i. Compared with the nanoindentation on the defect-free surface, the pit defect plays a significant role in the elastic-to-plastic transition that leads to the decrease of the nanohardness.
ii. The nanohardness increases in a wave that goes up in a period of three atoms when the pit moves far away from the indenter atom by atom, which is strongly linked to the crystal structure of periodic atom arrangement on {1 1 1} atomic close-packed planes of FCC metal. Moreover, it can be nearly considered as no influence on the nanohardness if the adjacent distance between the pit and indenter increases over than 16 atomic spacing.
iii. The formula for critical load of dislocation emission of Al thin film has been effectively revised with initial pit defect, where the reduction trend of nanohardness, as the increasing distance between the pit and the indenter, is in great accordance with the experimental results of nanoindentation on the stepped surface with various distances. Such study might be referential to the research of material properties with defects, especially in microchips and MEMS.
iv. When the size coefficient of width (of height) equals about one unit (half unit), the yield load experiences an obvious drop. When it reaches about two units (one unit), the yield load is nearly close to that of the nanoindentation on a stepped surface.
v. Compared to the width, the height of surface pit defect shows a greater influence on the yield load of thin film, implying that the height of the pit is a leading factor on the influence of yield load, which might have great significance to the defect design and applications in artificial materials.