Surface Energy of Diamond Cubic Crystals and Anisotropy Analysis Revealed by Empirical Electron Surface Models

A detailed knowledge of structure and energy of surface contributes to the understanding of many surface phenomena. In this work, the surface energies of 48 surfaces for diamond cubic crystals, including diamond (C), silicon (Si), germanium (Ge), and tin (Sn), have been studied by using the empirical electron surface models (EESM), extended from empirical electron theory (EET). Under the first-order approximation, the calculated results are in agreement with experimental and other theoretical values. It is also found that the surface energies show a strong anisotropy. The surface energy of close-packed plane (111) is the lowest one among all index surfaces. For the low-index planes, the order of the surface energies is γ(111) < γ(110) < γ(001). And surface energy variation of the (hk0) and (hhl) planes with the change of the included angle has also been analyzed. EESM provides a good basis for the surface research, and it also can be extended to more material systems. Such extensive results from the same theoretical model should be useful to understand various surface processes for theorists and experimentalists.


Introduction
Diamond cubic crystals, including diamond (C), silicon (Si), germanium (Ge), and tin (Sn), are important materials for a wide range of applications because of their excellent properties. Diamond has the highest hardness of any bulk materials since it is combined with the strong C covalent bond. And diamond is also a famous decorative material because of its special energy band structure. Si is mainly used as structural compounds, glass, and semiconductor materials. Ge is also used as a commercial semiconductor in transistors. Sn is also a famous alloying element for steel to resist corrosion. The performance and reliability of polycrystalline materials depend strongly on their microstructure, and especially the structure of crystal defects. Crystal surface is an important type of defect, which will affect many physical processes such as fracture, adsorption, corrosion, catalysis, crystal growth, radiation, etc. [1][2][3][4][5][6]. Surface energy, one of the most important quantities in surface science, is still difficult to determine experimentally for clean surface and just few data exist [7][8][9]. For different surfaces with different surface energies, the physical performances are different. For example, the hydrazine etching of Si is anisotropic in that the etching rate for Si (111) surface is slower than the others [10]. Therefore, it is essential to obtain the detailed knowledge of structure and energy of surface for understanding many surface phenomena.
However, first-principle calculations are computationally demanding and have typically been used only for the systems with a few atoms. And most semi-empirical potentials, fitted with the data from experiments and first-principle calculations, usually underestimate or overestimate the value, since the electronic structures of atoms are not taken into account for these semi-empirical models. It is almost impossible to find one type of potentials to reproduce the surface energies of all materials. The empirical electron surface model (EESM) [3,[20][21][22], based on empirical electron theory (EET) [3,20,23], has been found remarkably successful in predicting the surface energies for fcc-, bcc-and hcp-metals [3,[20][21][22]24]. The reliability of EET has been extensively examined in the fields of metals, alloys, metallic compounds and ceramics [25][26][27]. Through analyzing valence electron structure (VES), the contribution of different types of valence shell electrons has been summarized in EESM [22]. The model, without any experiment data input, can easily be performed and completed efficiently.
The main objective of this paper is to calculate and analyze surface energies of diamond cubic crystals revealed by EESM. Therefore, such extensive results from the same theoretical model must be useful to understand various surface processes for theorists and experimentalists. The rest of this paper is organized as follows. In section 2, we will briefly describe EESM and the VES of these diamond cubic crystals will be firstly calculated based on EET. Then the number and the type of dangling bonds on particular crystal surface will be analyzed. In section 3, we will present the results and discuss the anisotropy of surface energies, and the following is the conclusion.

Empirical Electron Surface Model
Surface energy (γ) can be defined as the additional value of the free energy per unit area of increase on a particular crystal surface, γ where ∆E represents the additional value of the free energy with the forming of new surface and ∆S is the added area of the crystal surface. In EET, based on the spatial distribution and the action characteristics of these electrons in the formation of solids or molecules, valence shell electrons are divided into four types, dumb pair electrons, magnetic electrons, covalent bond electrons (bonding electrons) and lattice electrons. The dumb pair electrons either represent the bonding and anti-bonding electrons with their resultant bonding power mutually cancelled by each other and their spins opposite to each other, or represent a pair of non-valence electrons of opposite spins deeply sinking down to the atomic orbit. The magnetic electrons represent the localized electrons from the single-occupied orbit. The covalent electron, also from the single occupied orbit, is paired with another electron of opposite spins from the near atom, contributing to atomic bonding. The lattice electron is present in a space among 3, 4 or 6 atoms. When the s energy band is wide, the electron might reach an energy level above the Fermi energy and is called lattice electron. The number of valence electrons of an atom in solids, including the covalent electrons and the lattice electrons, equals to the chemical valence of the atom in solids. The contribution of these four types of electrons to the surface energy has been analyzed in prior studies [20,22].
Hybridization states of the diamond cubic crystals C, Si Ge and Sn, the number and the distribution of these electrons, are listed in Section 2. 4. The number of the dumb pair electrons, the magnetic electrons, and the lattice electrons of the diamond cubic crystals occupy little proportion of the total electron number for valence shell (see Section 2. 4). Correspondingly, the energy caused by the dumb pair electrons, the magnetic electrons and the lattice electrons take up a very little portion in the crystal cohesive energy for C, Si, Ge and Sn. Therefore, for these diamond cubic crystals, when all covalent bonds between two nearest crystal planes are broken, the two nearest crystal planes will form two new crystal surfaces. So the energy requirement for the procedure is just the additional free energy (∆E).
where E α represents the bond energy of α bond; and Z α is the equivalent dangling bond number of α bond on particular crystal plane, and is obtained by the dangling bond analysis method (DBAM) [20,22]. In EET, the bond energy E α can be obtained from the value of the bonding capability (f) of covalent electron, the screen factor (b) upon the core electron, the bond length (D) and the number of covalent electron pairs (n).
(3) Therefore, the key of this surface energy model is the calculation of VES and the analysis of the equivalent dangling bond number.

Calculation of VES and Bond Energies of the Diamond Cubic Crystals
The diamond cubic crystals are the A4 type crystal structure, and the structure cell is showed in Figure 1  So the experiment bond lengths, which cannot be neglected, can be given as the following:  The equivalent bond number of the bonds can be calculated with the formula. (4) where α (α=A, B,…, F) represents bond name, I M represents the reference atom number in the crystal cell; I S represents the equivalent bond number for a reference atom to form α bond; I K is a parameter, I K equals 1 when the two atoms that form the bond are of the same kind, I K equals 2 equals when the two atoms are of different kinds. Therefore, the equivalent bond numbers of α bond are: According to the Pauling bond length formula in EET [23].
where D α represents the bond length; R (1) represents the single bond radius of the atom which forms the bond; n α represents the number of covalent electron pairs on the bond; and β is a factor with the unit of length. It is determined by Pauling's theory. where 0≤ε<0.050, and # is the largest n α in the structure. Therefore, the following r α′ (α′=B, C,…, F) equation can be obtained: Therefore, The n A equation can be obtained as: where ∑n C represents the total covalent electron number of the calculation system. By selecting a suitable β value and substituting the number of covalent electrons and the single bond radii of the atoms of the corresponding hybrid levels, the number of covalent electron pair n A on the strongest covalent bond can be calculated out. By substituting the obtained n A into the r α′ equations (7), the numbers of covalent electron pairs on all the other covalent bonds n α′ (α′=B, C,…, F) can be calculated out. Substituting the obtained n α into the bond length equations yields the theoretical bond lengths of all the bonds: If the bond length differences of all the bonds Δ @AB = | ? @AB − @AB | < 0.005nm, the applied hybrid levels of the atom are believed being in accordance with the actual state, thus the obtained corresponding parameters are the possible VES of the calculated crystal and the results are shown in Table 2. In the formula (3), f represents the bonding capability of covalent electron.
where the factor g represents the contribution of the coupling effect between spin and orbit of d-electrons to the bonding capability. For 4, 5 and 6 periods, g is equal to 1, 1.35 and 1.70, respectively. αʹ, βʹ and γʹ are contents of s, p and d electrons in covalent bonds, respectively.
where l, m, n and l′, m′, n′ represent the sum of the numbers of s, p, d covalent electron and lattice electron of the h and t states, respectively. The terms τ and τ′ are parameters for the h and t states, respectively, and value 1 when the s electron is covalent electron or 0 when the s electron is lattice electron. The values of these parameters can be taken from Ref [28], and some values of diamond elements are given in Section 2.4. n Tσ , C hσ , and C tσ represent the total numbers of covalent electron and lattice electron, the relative compositions for the h state and for the t state of σ hybrid level.
In the formula (15), b is screen factor upon the core electron.

RS.RTU
AV.RWX (15) The factors n and δ in the denominator reflect the total effect of the screen, coulomb, exchange and interrelated interaction of the inner electrons in solids. The value of n of C, Si, Ge and Sn is 3, the value of δ of C and Si is 2, and the value of δ of Ge and Sn is 1 [28]. So the bond energies can be calculated and are also tabulated in Table 2.

Analysis of Dangling Bond of the Related Crystal Plane in Diamond Cubic Crystal Cells
The equivalent dangling bond number (Z α ) on the particular crystal plane can be deduced out.
where I P represents the reference atom number on particular plane with the area S; I K is a parameter, which equals 1 when the two atoms that form the bond are of the same kind or 2 when the two atoms are different kinds; ʹ represents the equivalent dangling bond number for a reference atom to form one type bond with the atoms of the same neighbor crystal plane; N P (N P =1, 2,…) represents the crystal plane number passed through by the bond. It can be seen in Table 2 that the bonds of E and F can be neglected in analysis of dangling bond because the covalent electron pairs and the bond energies of E E and E F are very small (relatively to the value of E A or E B ). Therefore, only the bond A, B, C, and D are considered in the analysis of dangling bond.
As is shown in Figure. 2, there are two atoms on the (001) crystal plane with the area of a 2 nm 2 in the unit crystal cell of the diamond cubic crystals. Every atom with the atoms of the nearest neighbor crystal plane forms two A bonds and four C bonds, and with the atoms of the second neighbor crystal plane forms four B bond, and with the atoms of the third neighbor crystal plane forms two C bonds, and with the atoms of the fourth neighbor crystal plane forms one D bond. So the equivalent dangling bond numbers on the (001) plane with the area of a 2 nm 2 are Z A =2×1×2×1=4, Z B =2×1×4×2=16, Z C =2×1×(4×1+2×3)=20 and Z D =2×1×1×4=8.  As is shown in Figure.   For the (111) surface, the diamond cubic crystals have some specific structural properties, which differ from bcc-metals and fcc-metals [3,21]. As is shown in Figure 4 (a) (view along the <110> crystal direction), the (111) surface has two possible cleavage planes denoted by s 1 and s 2 , respectively. And the interlayer spacing of s 1 and s 2 , denoted respectively by d 1 and d 2 , has the relationship d 1 =a/ 4 3 <d 2 =3a/4√3. The cleavage plane s 1 and s 2 may become the terminations of the outer surface. However, the calculated results show that the surface energy of s1 is much higher than that of s2. The s2 termination of the outer surface is preferred from the perspectives of surface energy minimization. This may be the reason why the s2 outer surfaces are usually observed in experiment [29,30]. So only the dangling bond analysis of s2 will be given for the (111) surfaces in this paper. As is shown in Figure 4 (b), there are four atoms on the (110) crystal plane with the area of √3a2 nm2 in the diamond cubic crystals. Every atom with the atoms of the nearest neighbor crystal plane forms one A bond and six C bonds, and with the atoms of the second nearest neighbor crystal plane forms three B bonds and three D bonds. So the equivalent dangling bond numbers on the (111) plane are ZA=4×1×1×1=4, ZB=4×1×3×2=24, ZC=4×1×6×1=24 and ZD=4×1×3×2=24.
Based on the similar process on the dangling bond analysis of (001), (110) and (111) crystal planes, a computer program written in Microsoft Visual Basic language has been compiled to analyze and calculate the equivalent dangling bond number of the other 45 high-index crystals, some of which are perpendicular to the (001) plane and belong to the family of site planes (hk0) such as (710), (610), (510) and so on, and some of which are perpendicular to the (110) and belong to the family of site planes (hhl) such as (223), (332), (443) and so on. And the calculation results are listed in Table 3. During the process of analysis, we find that (113),

Hybridization States in EET of the Elements
In EET, dumb pair electrons, magnetic electron, covalent electron, and lattice electron are denoted by ǁ, ↑, •, and Ф, respectively. When the covalent electron and lattice electron are equivalent electrons, they are denoted by •, Ф or s′, p′, respectively. The hybridization states of the related elements (C, Si, Ge, Sn, Pb) are as the follows.   [11-13, 15, 16, 31, 32], listed in the last column for comparison. Based on above analysis, we can suppose that the formula for surface energy (γ) of clean and ideal surface for diamond cubic crystals as follow:

Results and Discussion
where α (α=A, B, C and D) is the bond name, and E α , the bond energy of α bond, and S, which are tabulated in detail in Table 3, the area of the related crystal planes. The physical meaning of the number 2 in this formula is that when these bonds between two nearest crystal planes were broken, there will be two crystal surfaces in the region. So the total area is two times of the area of related crystal plane. The calculation results are tabulated in detail in Table 4 and Table 5. The calculation values of the surface energy of various low-index planes (100), (110) and (111) are shown in Table 4 in eV·nm -2 . For comparison we also present some theoretical results [11-13, 15, 16, 31] and experimentally derived values [32]. Under the first-order approximation, the surface energy values calculated in this paper are close to the values which come from the other references. But for Sn, there is a lack of surface energy data in the literature. A further approximation in the present calculation is the neglect of the relaxation and reconstruction of the atomic positions which may lead to errors of up to a few per cent.

crystals (100)
The calculation values of the surface energy of various high-index planes, belonging to the family of site planes (hk0) or the family of site planes (hhl), are shown in Table 5 in eV·nm -2 together with included angles θ (hk0) between the (hk0) planes and (100) plane, and θ (hhl) between the (hhl) planes and (001) plane. The included angle θ (hk0) and θ (hhl) is calculated with the below formula:  For the planes of (hk0) family, surface energy variation of the diamond cubic crystals with included angle θ (hk0) between the (hk0) planes and (100) plane is shown in Figure. 5. It can be seen that the surface energy of the (hk0) surface decreases gradually with the increase of included angle in the interval 0°<θ (hk0) <45°. And the surface energy values of the rest angle range can be deduced according to the symmetry of crystal lattice. So the surface energy reaches the maximum values for the (100) plane at θ=0°, and the minimal values for the (110) plane at θ=45° among all these (hk0) planes which are normal to (001) plane.
For the planes of (hhl) family, Figure. 6 shows surface energy variation of the diamond cubic crystals with included angle θ (hhl) between the (hhl) planes and (001) plane. The calculated values of surface energy for the (hhl) planes are found to differ in character from those of the (hk0) planes. We can see that variation relationship of surface energy with included angle is not found to be smooth but has revealed an abrupt change when passing from one plane to another. This leads to the conclusion that the formation of a free surface (e. g. (111), (113), (553)) will take place strictly of these planes which are normal to (110), because even a small deviation of included angle yields a significant energetic disadvantage, i. e. an abrupt rise in surface energy. These planes have a similar property, which they have two possible cleavage planes and their Miller indexes are all odd. According to the symmetry of the lattice, we can also deduce the surface energy values of the rest angle range.

Conclusions
The surface energies of 48 surfaces for diamond cubic crystals (C, Si, Ge, and Sn) have been calculated by using the empirical electron surface models, and the spatial distribution of dangling bonds of various index surfaces has been analyzed with the dangling bond analysis method. The dangling bond electron density and the spatial distribution of covalent bonds have a great influence on surface energy of various index surfaces. Under the first-order approximation, the results are close to the other theoretical values and the experimental values. These surface energy values may be used as a consistent starting point for study of surface science phenomena. The calculated surface energy shows a strong anisotropy. And the close-packed (111) surface energy is the lowest of these index surfaces as predicted. For the low-index planes, the order of the surface energy result is γ (111) < γ (110) < γ (001) . So according to the surface energy minimization, the (111) texture should be favorable in the diamond cubic crystals. And surface energy variation of the planes of the (hk0) and (hhl) families with the change of included angle has been analyzed. This research model can be further extended to the surface energy estimation of more crystals, which will be the next research objects.