Pseudo-Hydrogen Passivation: A Novel Way to Calculate Absolute Surface Energy of Zinc Blende (111)/(¯1 ¯1 ¯1) Surface

Determining accurate absolute surface energies for polar surfaces of semiconductors has been a great challenge in decades. Here, we propose pseudo-hydrogen passivation to calculate them, using density functional theory approaches. By calculating the energy contribution from pseudo-hydrogen using either a pseudo molecule method or a tetrahedral cluster method, we obtained (111)/ surfaces energies of Si, GaP, GaAs, and ZnS with high self-consistency. This method quantitatively confirms that surface energy is determined by the number and the energy of dangling bonds of surface atoms. Our findings may greatly enhance the basic understandings of different surfaces and lead to novel strategies in the crystal growth.

maintain charge neutrality on the bottom surface, and also stabilize the bottom surface by satisfying electron counting rule (ECR) [27][28][29] . This passivation ensures that states at the bottom surface are localized and have no interactions with top surface. The energy of the top surface can be directly calculated if the pseudo-H passivation energy can be evaluated. Therefore, a natural and intuitive way to calculate the absolute surface energy is to analyze the pseudo-H passivation process. We show that the energy of the passivated surface can be directly calculated from the pseudo chemical potentials of the pseudo-H atoms attached on the bottom surface. Further, our calculations show that simple pseudo-molecules already give reasonably accurate values of the pseudo chemical potentials. Surface energy calculated from this approach shows comparable self-consistency with the wedge structure calculation, while the computation is much simpler. For high accuracy calculations, we construct a tetrahedral cluster with four equivalent (111)/( ) 111 facets to calculate the pseudo-H chemical potentials and the surface energies show improved self-consistency.
Consider a slab of a binary AB compound of zinc blende structure along [111] direction. The bottom ( ) 111 surface with B-termination is passivated with pseudo-H atom carrying fractional charge denote as H B . The absolute surface energy per unit area of the top (111) surface is then given by where E slab is the total energy of the slab with bottom surface passivated, n A (n B ) is the number of A(B) atoms in the slab, μ A (μ B ) is the chemical potential of A(B) atom, µ H B is the chemical potential of pseudo-H H B , α 111 is the area of (111) surface and σ bot pass is the surface energy of the passivated bottom surface. Assuming a thermodynamic equilibrium between the bulk and surface, we can write The pseudo chemical potential describes the energy gain from adding one pseudo-H atom and passivating one dangling bond on the bottom surface with this pseudo-H atom. This pseudo chemical potential can be decomposed into where the former part is the chemical potential of H B atom, and the latter part in bracket is the binding energy between the surface atom and the pseudo-H atom. This binding energy is just the energy of the passivated surface, divided by the number of passivated bond rather than surface area. It can be further decomposed into δ E int due to the intrinsic property of the surface atom, and δ E env due to electronic environment. Since passivated surfaces satisfy charge neutrality and ECR [25][26][27] , contribution from the environment is expected to be localized, and the major contribution comes from the local electronic environment around the pseudo-H atoms. It is difficult to calculate each individual part of the pseudo chemical potential, but the summation of all parts can be estimated under a local electronic environment similar to that of the surface atoms. This transforms the problem of calculating the energy of individual polar surface to a problem of estimating energy of bonds between surface atoms and pseudo-H atoms with a similar electronic environment. Such estimation only requires reproducing a local electronic environment similar to that of the surface atoms and pseudo-H atoms on the surface, but not the overall structure and symmetry of the surface. Therefore, this method is generally applicable to any crystal planes, as long as we can determine the pseudo chemical potential of H B with the similar local environment on the surface. and RHS is from standard slab calculation. Therefore, Eq. (7) can be used to directly determine the difference between the obtained pseudo chemical potentials and the standard slab calculations, which defines the self-consistency of the calculation.
Here, we propose two ways to determine the pseudo chemical potential for the (111)/( ) 111 surface, one from a simple pseudo-molecule and the other from a tetrahedral cluster. For the pseudo-molecule method, we construct a CH 4 -like molecule, with A(B) atom at the center of a tetrahedron bonded to four H A (H B ) atom at the corner of the tetrahedron, as shown in Fig. 1(a). It can be viewed as passivating four dangling bonds of a free-standing atom by pseudo-H. Since there are four bonds between center atom and pseudo-H, we can determine the pseudo chemical potential by Using this method, the chemical potential of HA and intrinsic contribution to the binding energy can be calculated, but it does not reproduce the local electronic environment. This method is straightforward and computationally inexpensive, nevertheless yields a fairly accurate result. Thus, it can be taken as the 0th order approximation for the pseudo chemical potential of the pseudo-H.
The cluster method, in addition, reproduces local electronic environment similar to that on the surface. The structure is shown in Fig. 1(b) (for details of the tetrahedral clusters, please refer to the supporting information). The cluster contains four (111) facets and all the dangling bonds on the surface are passivated by the corresponding pseudo-H atoms. The size of the cluster can be identified by n, the number of atoms on the edge. From the figure, we can identify surface atoms with different local environment. For each A atom on the corner, it is bonded to one B atom and three H A atoms; for each A atom on the edge but not on the corner, it is bonded to two B atoms and two H A atoms; for each A atom on the face of the tetrahedron, it is bonded to three B atoms and one H A atom (similar to (111) surface). We can denote the pseudo chemical potentials under these three conditions as μ H cor A , μ H edge A and μ H face A respectively. Local electronic environment of H A atoms on the face of the clusters is similar to that of H A atoms on (111) surface. Therefore, μ H face A is a good approximation to μ H A on (111) surface. Since number of A atoms, B atoms and pseudo-H atoms can all be expressed by the cluster size n, we write the total energy of the cluster as  23 , we also calculated (111)/(111) surface of GaAs, as well as the wedge structure of GaAs. For (111)/(111) surface, only the absolute energy of Si surfaces can be calculated by constructing symmetric slabs. However, for compound semiconductors, we can construct slabs with both surfaces passivated, and calculate the energy of the fully passivated slabs, where two different kinds of pseudo-H atoms are involved. Then by making use of Eq. (7), we can obtain the sum of the pseudo chemical potentials with the standard slab calculations. The differences per surface area between the sum based on the cluster (or pseudo molecule) method and the slab method were used to check the self-consistency of our method as well as to estimate the errors of the obtained surface energies. The percentage differences were also calculated. Throughout the calculation, the chemical potentials of Ga and that of Zn are taken at the rich limit for GaP (GaAs) and ZnS, respectively. As GGA functional usually gives a smaller band-gap than experimental value, which affects the accuracy of the surface energies 19 , we also performed calculations with screened hybrid functional of Heyd, Scuseria, and Ernzerhof (HSE) 38,39 on slabs and pseudo-molecules of GaP to test the difference between GGA functional and hybrid functional. Results indicate that our proposed method is general and is not functional dependent.
All the slab calculations were performed on (1 × 1) slabs, with (10 × 10 × 1) Monkhorst-Pack 40 k-point mesh for integration over Brillouin zone for GGA calculations and (4 × 4 × 1) for hybrid functional calculations. The slabs and clusters were separated by at least 15 Å vacuum. Pseudo-H atoms with charge q = 0.5e, 0.75e, 1.25e, and 1.5e were used to passivate dangling bonds of S, P (As), Ga, and Zn atoms, respectively. For Si, the passivation is done by true H. All the atoms in the slab and cluster were allowed to relax until forces converged to less than 0.005 eV/Å. Slabs along [111] direction contain 9 bi-layers, with both surfaces passivated by the corresponding pseudo-H atoms. Slabs along [110] direction contain 12 layers, and calculations were done both for slabs with both surfaces un-passivated and slabs with one surface passivated. Convergence tests are performed by increasing the number of layers in the slabs, and the results indicate that the obtained numerical errors of the surface energies are less than 0.5 meV/Å 2 .
After the slab calculations, we calculated pseudo chemical potential by using pseudo-molecule method. We calculated the sum of the pseudo chemical potential of the pseudo H atoms that passivate the anion surface and the potential of the H atoms that passivate the cation surface. A comparison between the sums obtained by our method and that obtained by the slab method is shown in Fig. 2. For 111/(111) surface, difference between the slab calculations and the pseudo chemical potentials calculations are all within 6 meV/Å 2 . Hybrid functional calculations of GaP show a difference of 9.1 meV/Å 2 , slightly larger than that for GGA calculations. Calculations on (110) surfaces of Si, GaP and ZnS also show differences within 6 meV/Å 2 . Previous calculations based on wedge structure have 3 meV/Å 2 3 and 20 meV/Å 2 19 differences for GaAs and zinc blend GaN respectively. Therefore, these results show good accuracies comparable to the wedge structure calculation, whereas the calculations are much simpler.
For the tetrahedral cluster method, although any four clusters can be used to solve for Eq. (8), different selections in fact give different results, as shown in Fig. 3. To make fair comparison between different systems, the percentage differences rather than differences per surface area were used. The variation of the results from clusters of different sizes is because we determine the pseudo chemical potentials from the energy differences between clusters. If energy differences between chosen clusters are large, the errors in the total energies of those clusters will be less significant. Hence the obtained pseudo chemical potentials will be more accurate. Therefore, in all the calculations, two smallest clusters and two largest clusters are chosen in the linear equation set to improve the accuracy. As can be seen from Fig. 3, the last three points show good convergence, with percentage differences less than 0.3% (1 meV/Å 2 in term of energy difference per surface area), and the obtained E tot (AB) calculated from Eq. (9) also shows only a few meV difference with that from bulk calculations. Therefore, we can take theses converged results as the pseudo chemical potentials obtained by cluster method, and the remaining differences as the errors of our method. In addition, the differences between sum of the pseudo H chemical potentials obtained by the cluster method and that obtained by the slab method are summarized in Fig. 2. For both (111)/(111) and (110) surfaces of Si, GaP and ZnS, the differences are within 1 meV/Å 2 with slab calculations, except for ZnS (110) surface. In particular, the difference of GaAs (111)/(111) surface is only 0.24 meV/Å 2 . We also calculated this difference by wedge method with both LDA and GGA functional, to compare the self-consistency of tetrahedral cluster method and wedge method. For both LDA and GGA functional, differences of GaAs (111)/(111) surface based on wedge structure method are around 3 meV/Å 2 , consistent with previous results 3 . Since the difference based on tetrahedral cluster method is one order of magnitude smaller than that based on wedge method, the tetrahedral cluster method shows better self-consistency and a higher accuracy than the wedge method does. To directly compare the absolute surface energies obtained by tetrahedral cluster method and that by wedge method in previous works 3, 23 , we also calculated the absolute surface energies of vacancy-terminated GaAs (111)/(111) surface. By tetrahedral cluster method, we obtained 39.2 meV/Å 2 for V Ga -terminated (111)-2 × 2 surface, and 51.2 meV/Å 2 for V As -terminated (111)-2 × 2 surface, whereas wedge structure method with GGA functional yields 45.2 meV/Å 2 and 52.7 meV/Å 2 respectively. This surface energy based on GGA functional by both tetrahedral cluster method and wedge method agrees reasonably well with recent works 23 by GGA functional (self-consistency check was not provided in ref. 23). Additionally, by using wedge method with LDA functional, we obtained 59.4 meV/Å 2 for V Ga -terminated (111)-2 × 2 surface, and 69.5 meV/Å 2 for V As -terminated (111)-2 × 2 surface, which agrees with values obtained by Zhang and Wei 3 (58 meV/Å 2 and 63 meV/Å 2 respectively). By comparing results of the same wedge structure calculations based on different exchange-correlation functional, we found that LDA and GGA functional yield different absolute surface energies, and large discrepancies between our results and previous work 3 are mainly due to this reason. Additionally, the cluster method also works well for (110) surface as shown in Fig. 2. These results confirm that the number of pseudo-H atoms attached to each surface atom should be the major contribution to the pseudo chemical potential of pseudo-H.
From the estimation of pseudo chemical potential of pseudo-H atoms on (111)/(111) and (110) surfaces, we can see that major contributions of the pseudo chemical potential are from the chemical potential of pseudo-H µ / H A B , and intrinsic contribution to the binding energy δ E int . The contribution from the local electronic environment, δ E env , is not significant. This is because on all the slabs and clusters, pseudo-H atoms have enough space to relax, and the surrounding local electronic environment only serves as a perturbation on the binding energy. Therefore, we have δ µ + δ . Also for δ E env , the contribution from the 1 st nearest neighbors of the surface atoms is much larger than that from the rest. Since the cluster method gives correct 1 st nearest neighbors for surface atoms on both 111/(111) and (110) surfaces, both results are very accurate. This shows general applicability of the pseudo chemical potential to determine the absolute surface energy of polar surfaces, regardless of the overall geometry of the surfaces.
Several advantages can be achieved by using our proposed methods. Since pseudo-molecule method takes both the chemical potential of pseudo-H atoms and the intrinsic contributions to the binding energy between the pseudo H and the surface atoms, it has comparable accuracies as the wedge structure calculations. More importantly, this method is much simpler than the wedge method and can be easily applied to other surfaces, especially to polar surfaces other than (111)/(111) surfaces, where wedges may be difficult to construct 23 . The stability issues of the wedge methods can be avoided by the cluster methods, because of the high symmetry of the tetrahedral structures. The atomic structures are allowed to fully relax without constraints in all the cluster calculations. Therefore, the inaccuracy caused by instability of pseudo-H atoms and the finite size effects in wedge structure calculations can be largely avoided. Even though 4 clusters are essential for determining the surface energies while The smallest clusters with n = 2 are always included and horizontal axis denotes the size of the largest cluster included. The difference between the Si (111) surface energy based on pseudo chemical potential calculation and that based on slab calculation is also included as a reference. Slight increase of differences at n = 9 for GaP and ZnS are mainly due to numerical errors between slab calculations and cluster calculations.
only two wedge structures are needed, two of the clusters are very small and easy to calculate. The large size difference between the large cluster and the small cluster greatly improves the accuracy of the method. Generally speaking, our methods are expected to save computing time and yield high accuracies.
From Table 1, we can conclude the general trends of the surface energies for different compounds. For unreconstructed (111)/(111) surfaces, surface energies for Si, GaP and ZnS follow the trend of their cohesive energies, since the electron redistribution is not significant and surface energies are just directly determined by the dangling bonds on the surface. However, surface energies on (110) surfaces decrease sharply with the increase of the iconicity of the materials, because in those compounds, ECR can probably be better satisfied when electrons in cation dangling bonds are transferred to anion dangling bonds due to the large electronegativity difference between them [27][28][29] . Also, the energy cost from forming the dimer-like structures on the surfaces is smaller for compounds with strong ionicity because the bond strength in such compounds is weaker than that in covalently bonded compounds.
In summary, we have proposed a new method to calculate surface energy of (111)/(111) polar surface of zinc blende structure, based on pseudo-H passivation analysis. Tests on (111)/(111) and (110) surface of Si, GaP and ZnS show very accurate results and good consistency with slab calculations. This method is not restricted to (111)/(111) surfaces, and it is generally applicable to other surfaces of many other types of crystals. The 0 th order approximation of the method yields reasonable accuracy that is comparable with wedge methods, but saves much computing time. The high order approach largely improves the accuracy of the absolute surface energy calculations. Since absolute surface energies are important to determine growth mode, wetting conditions, and crystal quality in hetero-structures of semiconductors, results from our methods are expected to provide very important physical insights in crystal growth techniques, thin film properties controls, and device performance enhancement. In particular, this method can give accurate surface energies of c/-c planes of wurtzite structures, where surface structures are similar to zinc blende (111)/(111) surface, and absolute surface energies of the wurtzite polar surfaces can be determined from the pseudo-molecule or tetrahedral cluster method 41 .

Si
GaP ZnS