Atomic-Scale Analysis of Self-Positioning Nanostructures∗

We performed an atomic-scale analysis of self-positioning nanostructures using the atomic-scale finite element method (AFEM) and Kriging interpolation technique. Equilibrium atomic configurations are determined for varying structure thicknesses and results are compared with the continuum mechanics solution under plane strain conditions. We observed significant decrease of the equilibrium curvature radius when the structure thickness is less than 40 nm. It is found that large compressive strains near the free surface affect the distribution of strains as the structure size decreases. [DOI: 10.1380/ejssnt.2008.301]


I. INTRODUCTION
Fabrication of nanoscale structures attracts substantial attention due to possibilities for development of new nanodevices. However, fabrication and manipulation of nanoscale structures are usually difficult to control. Simple and robust procedures for formation of three-dimensional nanostructures are desired for designing smart integrated micro-/nano electro mechanical systems (MEMS/NEMS). One promising approach is the method utilizing multilayer structures with crystal lattice mis-match [1][2][3]. The self-rolling phenomenon of the multilayer structure is referred to as the self-positioning [3]. The self-positioning occurs when a structure is subjected to strain/stress imbalance. For epitaxially grown samples, typical source of the deformation is a crystal lattice mis-matching strain in different layers. A typical structure consists of a substrate, a sacrificial layer, and several crystal lattice mis-matching layers (Fig. 1).
After selective etching of the sacrificial layer, a rolledup nanostructure is fabricated due to the self-positioning. This approach is simple and robust to create threedimensional structures from two-dimensional membranes, and deformation is controllable by selecting layers thickness, material properties, and crystal orientation. There are many potential applications of self-positioning structures [4]. While this technology is at early stage, several applications have been already developed [5][6][7][8][9].
For the self-positioning structure shown in Figure 1, the main structural parameter of interest is the local curvature radius in the equilibrium configuration. Analytical solutions have been derived based on continuum mechanics theory to predict the curvature radius [10][11][12]. Computational modeling has been also performed to understand characteristics of self-positioning nanostructures [13][14][15][16]. However, just few publications investigate selfpositioning structures taking into account atomic-scale effects. In this article, we present atomic-scale analysis of self-positioning nanostructures. We employ the atomicscale finite element method (AFEM) [16][17][18] for investi-gating dependence of self-positioning structure curvature radius on the structure thickness. In addition to previous results [16], we investigate strain distributions using mesh free finite element interpolation technique. AFEM results are compared with a continuum mechanics solution under plane strain condition.

II. THE ATOMIC-SCALE FINITE ELEMENT METHOD
The atomic-scale finite element method (AFEM) has been proposed for analysis of carbon nanotubes [17,18]. The AFEM has similarities to the conventional FEM, so that several computational algorithms for the FEM can also be applied to the AFEM. In the AFEM models, nodes and its neighbors reproduce atomic structures. The AFEM equation system is derived from the total energy minimization. The total energy E is approximated around current configuration x i by the following expansion: It is minimized by equating its vector derivative to zero: Using equations (1) and (2), the AFEM global equation system can be expressed in a form similar to conventional finite element equation system: where K is a global stiffness matrix, u is a displacement vector, and f is a load (force) vector. In the AFEM, the global stiffness matrix K and the load vector f are composed of second and first derivatives of the system energy with respect to atom positions: where E is the total energy of the atomic system, n is a number of atoms (AFEM nodes), and x i is the coordinate vector of an i-th atom. The energy E coincides with a potential energy for problems of finding static equilibrium configuration of atomic structures without external loads. Expressions of energy derivatives should be obtained to formulate the AFEM global equation system, so an energy function E should be at least twice differentiable with respect to positions.

III. MESH FREE FINITE ELEMENT INTERPOLATION
In the conventional finite element method (FEM), finite elements and local continuous functions are constructed at nodes for interpolation inside elements. In the AFEM, atomic bonds are used to construct element stiffness matrices, and there are no solid elements and interpolation functions as in the conventional FEM. It is not possible to perform displacement differentiation without continuous function approximating discrete nodal displacements, so we employ an interpolation method used in mesh free FEMs [19]. Mesh free FEMs are developed to overcome difficulties for construction of finite elements from nodes which are distributed arbitrarily. After solution of the AFEM equation system (3), displacement interpolation functions are constructed at nodes. In the mesh free FEM, shape functions are constructed by special interpolation techniques like moving least squares (MLS), partition of unity (PU), radial basis functions (RBF), and Kriging interpolations [19,20].
The MLS approximation is a popular method which is incorporated in several mesh free FEMs, but there are fundamental difficulties to specify essential boundary conditions. Recently, the Kriging interpolation technique has been introduced to formulate mesh free FEM. Since the Kriging interpolation exactly go through nodal displacements, it is easy to assign displacement boundary conditions.
The Kriging approximation finds unknown fitting coefficients for a given polynomial. In the moving Kriging interpolation, displacement u h is expressed as: where p is the vector consisting of terms of a given polynomial, r is the vector consisting of correlation functions, u is the vector consisting of one of a displacement components at all nodes in a partial region, and A and B are matrices given as:.
Here matrix P contains the term vector p at nodes, and matrix R is composed of correlation function values for a pair of nodal coordinates: where s i is the coordinate of i-th node, n is the number of nodes, and R is a correlation function. Polynomial term vector p and correlation vector r have the following appearance: where p i is the i-th term of a given polynomial and m is the number of terms. There are numerous possibilities for selection of polynomial and correlation function. The polynomial controls order of interpolation while correlation function determines relation between nodes. Quadratic polynomial basis is selected in this study. The term vector p for quadratic basis is given by: We select a quartic spline correlation function: and its parameters developed in recent study [21]: where d is the largest distance of possible pair of nodes in a partial region, θ is the positive correlation parameter. We employ the adaptive θ value (15) in this study. Strain components are calculated by differentiation of Kriging interpolation (6) with respect to coordinates x: FIG. 2: Bi-layer self-positioning structure composed of GaAs top and InAs bottom layers. We model only lattice mismatching layers (Layers 1 and 2).

IV. MODELING OF GAAS-INAS BI-LAYER STRUCTURES
We model bi-layer self-positioning structures composed of GaAs top and InAs bottom layers (Fig. 2).
Corresponding AFEM meshes are constructed in accordance with the real GaAs and InAs crystalline structure called zinc-blende crystal structure. The zinc-blende GaAs and InAs crystals are 3D cubic structures with arsenide atoms placed at eight corners and six face centers, and indium or gallium atoms occupying coordinates (0.  (Fig. 3).
We employ the Tersoff empirical inter-atomic potential function [22] and its parameters developed by Nordlund [23] for modeling In-Ga-As atomic interactions. In the Tersoff model, total potential energy E is given by the following function: where E i is the potential energy of atom i, V ij the potential energy of a bond i-j, r ij the distance from atom i to atom j, f C the cut-off function to disregard effects from distant atoms, f R a reactive component, f A an attractive component, and b ij a bond order term to represent multiatom interaction effects characterized by bonding angles.
Appearance of the potential function (17) is slightly different from original one after parametrization by Nordlund. The parameters fit by Nordlund correspond to basic elastic and melting crystal properties. We tested suitability of the Nordlund parameters for AFEM modeling of GaAs and InAs self-positioning structures [16]. Comparison of the numerical results is shown in Table I.

A. The plane strain solution
We analyze atomic-scale effects in bi-layer selfpositioning structures using the AFEM. Results are compared with those obtained by the analytical continuum mechanics solution under plane strain conditions [10]. In the plane strain solution, curvature K, and strain components ε x , ε y , and ε z are written as follows: where is the Young's modulus, ν i is the Poisson's ratio, t i is a thickness, ε 0 i is a lattice mis-matching strain, R is the curvature radius which is inverse of the curvature K, and y b is the neutral local y coordinate where bending component of strain is zero. The subscript i indicates which layer the parameter is assigned to. If the bottom surface corresponds to the zero y level (y = 0), then y b can be calculated as: We use structural properties found by the AFEM test of Nordlund parameter suitability (Table I) to calculate plane strain solutions. Definition of thickness for structures consisting of just a few crystal layers in the thickness direction should be done with care when calculating curvature radii with an analytical technique. We add some offset equal to the 'radius' of an atom at each free surface. While adding such an offset is not significant for thick structures, it can be important for problems with small number of unit crystals in y direction. If we adopt half of the atom connectivity length as an offset, then for zinc-blende crystal structures the offset is equal to √ 3a/8, where a is a lattice period.
Corresponding offsets are 0.1224 nm for GaAs and 0.1425 nm for InAs.

B. Comparison of curvature radius
AFEM calculations have been performed to obtain equilibrium atomic configuration of bi-layer GaAs/InAs self-positioning structures. Structures consist of 3c and c number of GaAs and InAs unit crystals in the thickness direction, so totally 4c unit crystals where c is the problem size parameter (Fig. 4). We determine the curvature radius of bi-layer structures with problem sizes c = 1, 2, 4, 8, 12, 16, 24, and 36, which correspond to thicknesses 2.56, 4.86, 9.46, 18.65, 27.84, 37.03, 55.41 and 82.98 nm. Number of unit crystals in the length direction is 16c. Periodic boundary condition is employed to imitate infinite width corresponding to plane strain solution. One of the structure ends is fixed to restrict displacements in the x direction. The largest AFEM model consists of 1, 329, 986 atoms that corresponds to almost 4 million degrees of freedom. We employed the preconditioned conjugate gradient (PCG) algorithm with sparse row matrix format for solution of the AFEM equation system (3) at each computational step. Figure 5 shows the ratio of the curvature radius determined by the AFEM and by the continuum mechanics solution for varying structure thickness. The asymptotic value of the ratio is found by least square fitting of a power function R(c) = α 1 (α 2 c + α 3 ) −β + γ, where c indicates the size of an atomic system, and α 1 , α 2 , α 3 , β, and γ are fitting parameters. The parameter γ corresponds to the converged value of curvature radius for infinite number of crystal layers.
According to the fitting, ratio of curvature radius is converged to 1.0037 for infinite number of crystals. Relative difference is −0.36% for the biggest problem (c = 36) we investigated. Therefore, the AFEM result converges to continuum mechanics solution as the structure thick- ness increases. In the smallest problem (c = 1), thickness is 2.56 nm, and relative difference is −18.4%. Difference between atomic-scale and continuum mechanics solutions for structures of small thickness can be attributed to inability of continuum mechanics methods to model surface effects of neighboring atom absence.

C. Comparison of strain components
The Kriging interpolation is incorporated in the AFEM procedure to calculate strains. Strain distributions are obtained as contour plots for atoms near the mid plane normal to z. The Delaunay triangulation [25] is used for rendering of contour plots. We used a color function and a texture mapping technique for generating contour plots [26]. Figure 6 shows strain distributions and contour plots in problems c = 1 and 8. Strain in the z (width) direc- tion is zero, so strains ε x and ε y are compared with the plane strain solution. According to the plane strain solution, strains are linear function of local y coordinate, and the strain ε y is discontinuous at the interface of bi-layer due to difference in material properties and lattice mismatching strains. Figures 6(a) Figures 6(b) and (e) are contour plots of the ε x distribution for problems c = 1 and 8. Distribution for the problem c = 1 looks a little rough since there are less atoms in the whole structure, but it can be recognized that fixed and free end effects are more significant in the smaller structure. In Figs. 6(c) and (f), contour plots of the ε y distribution are shown. At the bottom and top surface, the contour plot is enlarged to show details of compressive strain distribution. The contour plot shows that the compressive strain takes place at the bottom and top surfaces. At the bottom and top surfaces, atoms outside of the structure is absent. This leads to attractive forces in the structure thickness direction. Figure 7 show atomic configurations at the bottom and top surfaces. Smaller atomic spacing are observed near the top and bottom surfaces which are consistent with strain distributions.

VI. CONCLUSION
The atomic-scale finite element analysis is performed for modeling bi-layer self-positioning structures consisting of GaAs and InAs layers. We employed Tersoff-Nordlund potential for modeling atomic interactions of GaAs and InAs materials. Equilibrium atomic positions are calculated for varying thickness using the AFEM. We found that the GaAs/InAs bi-layer film radii of curvature are significantly different from the continuum mechanics solution when total thickness becomes smaller than 40 nm. On the other hand, the AFEM results converge to the continuum mechanics solution when thickness increases over 100 nm. We also calculated strains in the self-positioning structures using Kriging interpolation technique. We found that strain distributions obtained by the atomic-scale modeling are similar to the results of plane strain solution. However, AFEM results show that surface effects are significant for the distribution of strains in the smaller structures. Attractive forces due to absence of atoms lead to significant difference from strains predicted by continuum mechanics theory. In order to calculate strains near the free surfaces, surface effects should be taken into account.