A finite element method to simulate dislocation stress: A general numerical solution for inclusion problems

We developed a simple and efficient way to simulate stress from arbitrary shape dislocations within the finite element method (FEM) frame-work. The new method is implemented as a single-step FEM simulation using analytic solutions in an infinite medium as input terms of FEM solver with special internal boundary treatment. It is fundamentally equivalent to the multistep “image method” proposed by Van der Giessen and Needleman [Modell. Simul. Mater. Sci. Eng. 3 , 689 (1995)], but the force equilibration is achieved by a single step FEM method. The singularity at the dislocation core line is removed without mesh modification, which removes many technical difficulties to couple the dislocation stress model to other process/device simulation models.


I. INTRODUCTION
Mechanical stress plays an important role in both detrimental 1 and beneficial 2,3 ways in modern semiconductor device fabrication and, thus, improving stress simulation capability, especially dislocation stress model, has been of great interest in the technology computer-aided design (TCAD) community for decades. However, there has been no practical dislocation stress simulation methodology applicable to routine semiconductor modeling activity.
In principle, solving stress distribution generated by dislocations is Eshelby's inclusion problem, 4 which has a history longer than a half century. Zhou et al. 5 provided a review on the research history of the inclusion problem. Early research on the inclusion problem focused on the analytic solution in special geometries and boundary conditions. In 1993, Lubarda et al. reduced Eshelby's inclusion problem in a finite size domain with traction-free boundary conditions to a boundary value problem, 6 which was extended to the general boundary conditions by Van der Giessen and Needleman 7 and led to several different types of finite element method (FEM) solutions. [8][9][10][11] To our knowledge, four major FEM methodologies have been published so far for dislocation stress simulation in a multidomain structure. The first method by Van der Giessen needs four calculation steps: (1) calculate analytic solution in the region, Ω 1 in Fig. 1, enclosing a dislocation and the accompanied Volterra cutplane, 12 (2) solve FEM for the remaining region Ω 2 in Fig. 1 with the displacement boundary condition at the interface, (3) calculate the sum of traction from the previous steps and solve FEM equations by applying the negative traction to the interface, and (4) sum up solutions from all the previous steps. Although the idea is straightforward, solving FEM stress equations in subdomains of a complex structure separately with the specified boundary conditions and merging the results require complicated correction processes depending on the configuration of the external boundary and internal interfaces. 8 The second method is to insert an "eigenstrain" inclusion plate and let boundary conditions be taken care of by the FEM solver. 9,13 Although this method can successfully calculate the stress distribution in a complex structure, it still has a limitation for practical application because defining a thin region and refining mesh are time-consuming and can cause bad convergence. In addition, it is hard to superpose multiple dislocations. In the third

ARTICLE
scitation.org/journal/adv method used by Gatti et al., 10 analytic solutions were used as the input of FEM solver with user-prepared external function. Their method also requires a special mesh prepared by users to avoid singularity near the dislocation core line. Finally, the extended-FEM (XFEM) uses displacement discontinuity at the Volterra cut-plane, which requires a significant change to the standard FEM. 11 Other variations such as "homogenization method" 14 and "thermoanalog model" 15 also require strain averaging and fine regular mesh preparation, respectively, near the Volterra cut-plane. All four methods produce good results for the tested structures, but they require either special mesh preparation or significant extension to the standard FEM, which can be a hurdle in coupling the dislocation stress model to other general semiconductor process/device simulation models. Therefore, we developed a simple and efficient method compatible with the standard FEM.

A. Governing equations of FEM
The principle of the virtual work equation, also known as the weak form of the force equilibration equation, 16 is used as the governing equation of the FEM solver.
where σ is the stress tensor, B is the differential operator to make Bu the strain tensor, δu is the infinitesimal virtual displacement vector, f is the body force vector, and t is the traction vector at the boundary. In our simulator, the stress, σ, is given by, where D is the elastic stiffness matrix, (u * , σ * ) is the previously equilibrated solution pair, and σ i is the sum of all other stress sources not equilibrated yet, such as thermal expansion and lattice mismatch. In FEM solver, u * is used as the initial guess of the solution u. In case there is no σ i term or change in the boundary condition, the final solution is the same as (u * , σ * ). With the nonzero σ i term or modified boundary conditions, the newly equilibrated (u, σ) deviates from (u * , σ * ). In our method, the analytic stress solution in the infinite medium, (u a , σ a ), is used as the previously equilibrated solution (u * , σ * ) in Eq. (2). A finite size simulation domain guarantees a change in boundary conditions from an infinite size domain. The boundary effect is taken care of by the third term of Eq. (1).
Depending on the structure of the simulation domain, there can be multiple materials or regions, which result in boundaries inside bulk, and treating internal boundary properly is important to get correct solution. Our FEM dislocation stress model requires continuous u at the internal boundary, which is different from continuous Δε at internal boundaries for other stress sources. The details are discussed in Sec. II C.

B. Analytic solutions in an infinite medium
The analytic solutions for curved dislocations can be represented using complicated integral equations, which can be simplified to a set of simple functions for line dislocations, 12 and they are given as, where ⃗ b is the Burgers vector, Ω is the solid angle made by the Volterra cut plane terminated by a closed loop of a dislocation core line, d⃗ s is the line integration element along the dislocation core line, ⃗ R is the vector defined between the dislocation core and the evaluation position, R is the magnitude of ⃗ R, μ is the shear modulus, ν is the Poisson ratio, and I is the unit dyadic.
Equations (3) and (4) can be formulated as a componentwise form to facilitate implementation, where ε ijk is the Levi-Civita symbol and δij is the Kronecker delta.

C. Singularity removal and boundary conditions
The analytic solution σ a diverges near the core line and can cause a convergence problem. In our method, to remove the singularity, a cut-off radius comparable to the interatomic distance is introduced. When the distance R in the denominators of Eqs. (3)-(6) is less than the cut-off distance, R 0 , it is replaced by R 0 . As a result, both u a and σ a are consistently prepared and the stress is smoothly connected across the core point without diverging. The analytic way to avoid the singularity removes the special mesh preparation step required in other methods and guarantees convergence of the simulation. Although Eqs. (5) and (6) are derived assuming isotropic elasticity, our methodology with Eqs. (1) and (2) is still applicable to any consistently prepared (u * , σ * ), regardless of the anisotropy of elasticity because Eqs. (1) and (2) add only corrections from boundary effects to the original input.
At the interface between the analytic solution domain (Ω 1 ) and the remaining region (Ω 2 ), displacement continuation condition is applied and force equilibration is applied after the virtual work equations from both sides of the interface are summed, which takes care of coupling effects between two regions. The displacement continuation at the internal boundary is critical to get the correct stress distribution from dislocations. A different boundary condition, Δε continuation, where ε is the strain, is used for most other cases to keep the process history during multistep process simulations. For example, when a new amorphous layer is deposited on top of stress equilibrated layers containing dislocations, dislocation stress will not affect the new layer. Figure 2 shows details of the simulation step for multistep processes.

III. RESULTS AND DISCUSSION
The new method produces the same result as the eigenstrain method (Fig. 3) and the analytic solution domain can be arbitrarily chosen as long as it includes the dislocation core line and the accompanied Volterra cut plane as shown in Fig. 1. Figure 4 compares two methods in the stacked nanosheet gate-all-around transistor structure 17 with the overlaid mesh, where two dislocation core lines are located at both ends of the Volterra cut-plane [case (c) in Fig. 5]. With the eigenstrain method, defining the thin inclusion plate and remeshing the structure could be a complicated and

ARTICLE
scitation.org/journal/adv time-consuming task especially when it is defined along a curved boundary. With the new method, however, the boundary shape does not add any complexity. Figure 6 shows an example of dislocation termination at the curved boundary. In addition, the new method enables superposition of multiple dislocations. Figure 5 shows some combinations of edge dislocations to make extended forms. Due to the linearity of the equation, the superposed analytic solution pair still works with the governing equations. To our knowledge, a dislocation loop stress simulation result using the analytic method with a proper boundary correction in a finite size domain is reported for the first time (Fig. 7). Even with the advantages of our new methodology described in this report, implementing numerical integrations included in Eqs. (3)-(6) can be time-consuming and difficult for complex shape dislocations. Therefore, we are further extending our method to overcome such difficulties and the result will be reported in a separate paper.

IV. SUMMARY
In summary, we developed a practical FEM stress simulation method applicable to all types of dislocations, which can be also an efficient numerical solution for the inclusion problem. By incorporating the singularity-removed analytic solutions into an FEM stress solver with proper internal interface boundary treatments, a fast converging FEM dislocation stress simulation was accomplished, which removes technical issues of pre-existing methods and thus enables easy coupling to other process/device simulation models.