# The Elastic Property of Bulk Silicon Nanomaterials through an Atomic Simulation Method

Authors: Jixiao Tao; Yuzhou Sun
Journal: Journal of Nanomaterials
 (2016)

---

## Abstract

This paper reports a systematic study on the elastic property of bulk silicon nanomaterials using the atomic finite element method. The Tersoff-Brenner potential is used to describe the interaction between silicon atoms, and the atomic finite element method is constructed in a computational scheme similar to the continuum finite element method. Young’s modulus and Poisson ratio are calculated for[100], [110], and [111] silicon nanowires that are treated as three-dimensional structures. It is found that the nanowire possesses the lowest Young’s modulus along the [100] direction, while the [110] nanowire has the highest value with the same radius. The bending deformation of [100] silicon nanowire is also modeled, and the bending stiffness is calculated.


---

## Body

## 1. Introduction

Over the past decade, one-dimensional nanostructures such as nanowires have been extensively studied through both theoretical and experimental methods. The most important one is the silicon nanowire (SiNW) that has been successfully applied in the nanoelectromechanical devices such as field effect transistors (FETs) [1–3] due to their unique optoelectronic and mechanical properties [4, 5]. Researchers have made significant contributions in the study of the elastic property of SiNWs, in which they are considered with different diameters and different axial or surface orientations. It is generally accepted that Young’s modulus tends to increase with the increasing diameter [6, 7].The molecular dynamics and density functional theory [6–9] are the common atomic scale simulating methods to study the elastic property of nanostructures. However, their computational cost is very huge and they are valid only for the small size structures. In order to overcome this shortcoming, Liu et al. [10, 11] proposed an order-N atomic finite element method (AFEM) which can achieve the same accuracy with molecular mechanics but is much faster than the latter one. Sun and Tao [12] adopted this method to investigative the elastic properties of carbon, boron nitride, and silicon carbide nanotubes, and Young’s moduli are accurately obtained and the buckling behavior is better displayed. The present work extends its application to three-dimensional nanostructures in order to study the fine elastic property of SiNWs. An appropriate type of Tersoff-Brenner potential is employed to describe the atomic interaction. An atomic element contains 17 atoms for the bulk silicon crystal, and the global stiffness matrix and nonequilibrium force vector are assembled similar to the continuum finite element method (FEM). The nonlinear iteration is implemented to determine the equilibrium.


## 2. Atomic Scale Modeling Method for Bulk Silicon

Bulk silicon has the diamond lattice structure as illustrated in Figure1, and its elastic property has a great dependence on their growth directions. The present work focuses mainly on SiNWs grown along the [100], [110], and [111] crystallographic orientations, and Figure 2 shows the representative cross section and side view of the SiNWs with different growth directions.Figure 1
The diamond lattice structure of bulk silicon.Figure 2
Cross sections and side view for[100], [110], and [111] SiNWs.
(a)
[100] (b)
[110] (c)
[111]In order to describe the Si-Si interaction, the Tersoff-Brenner potential is adopted with the parameter set refined by Erhart and Albe [13]. Tersoff-Brenner potential is a multibody potential which can be used to describe the interaction between atoms including C, B, H, Si, and N: (1)VBrIJ=VRrIJ-B¯IJVArIJ,where VR and VA are the repulsive pair potential and attractive pair potential, respectively; rIJ is the distance from atom I to atom J; B¯IJ is the bond order function. All above functions have analytic forms (2)VRrIJ=D0S-1exp⁡-β2SrIJ-r0fcrIJ,VArIJ=SD0S-1exp⁡-β2SrIJ-r0fcrIJ,B¯IJ=BIJ+BJI2,BIJ=1+∑k≠i,jfcrikexp⁡2μrij-rikgθijk-1/2,where fc(r) is the cutoff function(3)fcr=1,r<R-D,12-12sinπ2r-RD,R-r≤D,0,R+D<r.The angular function g(θ) is given by(4)gθ=γ1+c2d2-c2d2+h+cos⁡θ2.The potential parameters are listed in Table 1.Table 1
Parameter sets for silicon.
D 0 (eV) r 0 (Å) S β (Å−1) γ c d h 2 μ R (Å) D (Å) 3.24 2.222 1.57 1.4760 0.09253 1.13681 0.63397 0.335 0.0 2.90 0.15The basic idea of AFEM is to treat atoms as the continuum FEM nodes, and each element is characterized by a set of discrete atoms [10, 11]. For a system of N atoms, the energy stored in the atomic bonds can be evaluated using Tersoff-Brenner potential, and the atomic positions are determined by minimizing the system energy (5)Etx=Vtx-∑i=1NFi·xi,Vtx=∑i<jNVBxj-xiin which x=(x1,x2,…,xn)T, xi is the position of atom i, and Fi is the external force exerted on atom i.Giving Taylor expansion ofEt(x) and substituting it into(6)∂Et∂x=0we have(7)KΔu=P,where Δu is displacement increment, K and P are, respectively, the stiffness matrix and nonequilibrium force vector given by(8)K=∂2Et∂x∂xx=x0=∂2Vt∂x∂xx=x0,P=-∂Et∂xx=x0=F-∂Vt∂xx=x0.For the present nonlinear system, (7) is solved iteratively until P reaches zero. Here, Newton method is used to solve the problem.In AFEM, the choice of element depends on the atomic structure and nature of atomistic interaction. In the studies for carbon nanotubes by Liu et al. [10, 11], an AFEM element containing 10 atoms is developed in order to capture the multibody interaction. For the current diamond lattice structure and multibody interaction, an AFEM element should contain 17 atoms that include 4 nearest-neighbor atoms (red balls) and 12 second nearest-neighbor atoms (blue balls) as shown in Figure 3. Therefore, the element stiffness matrix and the nonequilibrium force vector are given by(9)Kelement=∂2Vt∂x1∂x13×312∂2Vt∂x1∂xi3×4812∂2Vt∂xi∂x148×3048×48,Pelement=F¯1-∂2Vt∂x13×1048×1.Figure 3
The AFEM element for silicon.In the present study, all work, including assembling the stiffness matrix and force vector and solving the equation system, is performed with our Fortran codes. The following steps are used. Firstly, the initial configuration of SiNWs with a uniform bond length is constructed, and the obtained coordinates and bond information are stored in an array. Using the obtained array, the information about the first and second neighbor atoms for each atom can be found and is stored in another array. Secondly, boundary conditions are applied to SiNWs to get the initial equilibrium coordinates of the system. Here, we restrain one side of SiNWs and make the other side free until the system returns to the equilibrium configuration. Table2 gives the bond lengths of different types of SiNWs when the potential parameters in [13] are employed. These results are similar to the work in [14]. In the third step, displacement field is applied to the initial equilibrium coordinates in order to simulate the axial tension or compression. It is achieved by completely fixing one side of SiNWs and incrementally imposing an axial movement at the other end. The length of the nanowire is changed by 0.01 nm per loading step. In the final step, the total energy of SiNWs is obtained, and it is treated as a function of strain. With the first- and second-order derivatives, Newton iteration method is thus applied to obtain the equilibrium state, and 3–5 iterative steps can generally achieve a good convergence. Therefore, the present method is far faster than molecular dynamic method and is very efficient for the nanostructures with a larger number of atoms.Table 2
Bond lengths of [100], [110], and [111] SiNWs.
Growth directions 100 110 111 Diameter (nm) 1.21 1.83 2.43 1.21 1.47 1.77 0.83 1.20 1.59 Bond length (nm) 0.233 0.234 0.234 0.233 0.233 0.234 0.232 0.233 0.233


## 3. Results and Discussions

In continuum mechanics, the material can be represented by two independent constants, namely, Young’s modulusY and Poisson’s ratio ν. For a material undergoing a uniaxial deformation, Y is defined as [11, 15](10)Y=1Ω∂2Vt∂ε2ε=0,where Ω is the volume corresponding to the initial equilibrium configurations of SiNWs; Vt is the total strain energy under tensile or compressive deformation, and ε is strain.The Poisson ratioν is calculated as (11)ν=1εD-D0D0,where D0 and D are the mean diameters of SiNWs corresponding to the initial and deformed equilibrium configurations, respectively. They can be derived from S=πD2/4, where S is the cross-sectional area.Young’s modulus for the bulk silicon crystal under tension along[100] directions is first calculated in this paper, and it has the bulk structure with the side length 2.16 nm. The total potential energy is plotted versus strain in Figure 4. Using the polynomial curve fitting, the obtained Y is 134.03 GPa. For [110] and [111] directions, the present computational results are, respectively, 190.2 GPa and 178.3 GPa, which drop in the literature reports. Ma et al. [6] used all-electron density functional theory to calculate the binding energy, heat of formation, and Young’s modulus of hydrogen-passivated SiNWs with various diameters and orientations, and Young’s moduli of bulk silicon for the [100], [110], and [111] directions are 133.24 GPa, 157.66 GPa, and 163.58 GPa, respectively. Lee and Rudd [15] found that the bulk value of Young’s modulus is 122.53 GPa for the [100] direction by the density functional theory.Figure 4
Potential versus the length change of[100] bulk Si crystal for axial tensile/compressive deformation.Next, the numerical computation is carried out for Young’s modulus and Poisson ratio of SiNWs with different dimensions. Table3 shows the Poisson ratio for [100], [110], and [111] SiNWs, which have a small dependence on diameter. The calculated Young’s moduli are plotted in Figure 5. Leu et al. [5] have given a systematic study of the mechanical property of strained small diameter silicon nanowires using ab initio density functional theory, and the values of Young’s modulus are 139–153 GPa and 43–131 GPa for [110] and [111] SiNWs, respectively. In the researches by Ma et al. [6] and Lee and Rudd [15], Young’s modulus increases while the wire mean diameter increases, and at certain critical size it approaches the bulk value. The present results agree well with these previous studies. It can also be seen from Figure 5 that, with the same mean diameter, nanowires with different grown directions have distinctly different Young’s moduli, and the wires along [110] direction have the highest value while [100] wires have the lowest value.Table 3
Poisson ratio of [100], [110], and [111] SiNWs.
Growth directions 100 110 111 Wire number 1 2 1 2 1 2 Poisson ratio 0.14 0.14 0.19 0.22 0.22 0.25Figure 5
Young’s modulus versus the mean diameter.Besides the axial strain, the simulation has also been carried out for the bending of SiNWs. The bending stiffness of a SiNW can be calculated as(12)U=1L∂2Vt∂κ2,where L is the length and κ is the curvature of the bent SiNWs and is related to the bending angle θ as κ=θ/L [16].The bending deformation can be achieved by incrementally rotating the two end planes in the opposite directions with respect to a line that is perpendicular to the axis of the undeformed SiNW and passing through its center. During the simulating process, the two end planes are reversely rotated 0.1 rad per loading step. Two[100] SiNWs with different diameters are studied in this paper. The change of the total energy is plotted as a function of the bending angle in Figure 6. Using polynomial curve fitting, we can obtain(13)Vt=a1+a2θ+a3θ2.The bending stiffness can be calculated by replacing θ with κL. The bending stiffness is, respectively, 5432.7 and 1964.2 ev·nm for SiNWs with the diameters 2.43 nm and 1.83 nm. Menon et al. [17] investigated the bending stiffness of T and C type nanowires, and their computed results are 12.2 × 10−6 ev·m and 8.5 × 10−6 ev·m, respectively.Figure 6
The change of the total energy as a function of bending angle for[100] SiNWs.


## 4. Conclusions

The AFEM element is developed for the silicon nanowire, and a basic element contains 17 atoms while it is treated as the bulk structure. Young’s modulus and Poisson ratio are calculated for[100], [110], and [111] SiNWs, and the bending stiffness is also obtained. AFEM is proved to be an accurate atomic simulation method for the bulk nanomaterials, and it is found that the [110] nanowire has the largest Young’s modulus while the [100] nanowire has the smallest value. Newton iteration method is applied to obtain the equilibrium state, in which the first- and second-order derivatives are used and 3–5 iterative steps can achieve a good convergence, and the developed method is much faster than the molecular dynamic. In order to take advantage of the order-N nature of AFEM, further studies can be made to combine AFEM and the classical FEM software to solve problem with a large number of degrees, for example, the nanoindentation of bulk monocrystalline silicon.



---

*Source: 8231592-2016-09-28.xml*