A brief review of formation energies calculation of surfaces and edges in semiconductors

Chuen-Keung Sin, Jingzhao Zhang, Kinfai Tse, and Junyi Zhu

The Chinese University of Hong Kong, Hong Kong, China

Abstract: To have a high quality experimental growth of crystals, understanding the equilibrium crystal shape (ECS) in different thermodynamic growth conditions is important. The factor governing the ECS is usually the absolute surface formation energies for surfaces (or edges in 2D) in different orientations. Therefore, it is necessary to obtain an accurate value of these energies in order to give a good explanation for the observation in growth experiment. Historically, there have been different approaches proposed to solve this problem. This paper is going to review these representative literatures and discuss the pitfalls and advantages of different methods.

Key words: surface; first principle; morphology

1. Introduction

The fundamental thermodynamic theory of surfaces, initialized by the American scientist Josiah Willard Gibbs, is one of the most practical tools for the study of surface-related phenomena[1-4]. For this approach, the key quantity is the surface free energy (or surface energy),which is equal to the excess free energy per unit area on account of the creation of surfaces, compared with the bulk structure, plus the energy change due to deformation in liquid or reconstruction in solid[5]. For a liquid surface, the phrase "surface tension" was widely used for the surface free energy. This phrase had also been extended to the surfaces of solids previously, which,however, produces a confusion between the surface free energy and the intrinsic surface stress of solids, revealed by Gibbs[6]. Therefore, currently, the term "surface tension" is rarely used for solids. Instead, we use the term surface energy, which depends on the temperature in the experiment or computational simulation bywhereF,U,TandSare the surface (free) energy, internal energy, temperature and entropy respectively.

For surfaces of solids, especially metals and semiconductors, the surface energy is important in many related fields, determining the equilibrium shape of monocrystals, brittle fracture, or the rate of sintering. Wulff construction[7]is a prevalent technique being used in the prediction of the morphologies because it is able to formulate the relation between the surface energy of surfaces (and edges of 2D lattices) and the relative feasibility of their formations in different directions.Hence, the principle of Wulff construction will be introduced as an appetizer before any estimation on the surface energy.

The anisotropy of surface free energybrings about the idea of equilibrium shape, which is known as the Wulff shape,attributed to the pioneering works by Curie and Wulff[6,8],which is the shape that minimizes its total surface energy. It is conventional to visualize the anisotropy of surface energies throughwhich are polar plots ofversus the polar anglesand the azimuthal angleφof the unit vector perpendicular to the corresponding surfaces. In the 2D case,is only related to the polar anglemeasured from a certain appropriate crystallographic direction. The process of Wulff construction transforms theinto equilibrium shape. According to Wulff construction, the equilibrium shape is the inner convex hull bounded by planes (perpendicular lines for 2D case) drawn perpendicular to each direction at a distanceis a normalization constant) from the origin, shown in Fig. 1. As a result, orientations with planes that lie outside of the inner hull are unimportant in the construction of the final equilibrium shape because their energies are too high.This is indicated by planes that draw the "ears" at the corner of the Wulff shape in Fig. 1. Such an exclusion of unimportant planes results in an equilibrium shape with edges and/or sharp corners. Practically, we do not have an explicit functionγ(θ,φ) in an arbitrary direction, we may however obtain the surface energy of surfaces in 3D or edges in 2D in certain directions. Therefore, direct calculation of equilibrium shape from the experimental or simulation data of surface energies is more favorable. According to Ref. [9], the distancerfrom the origin of the crystal shape in arbitrary direction is given by

wherehis the unit vector of the arbitrary direction, andis the surface energy of a plane or an edge with unit normal vectorm. The minimum value ofis chosen among different surfaces or edges in directionm. Even though the energy of surfaces or edges are obtained only in limited directions, we can calculate the radius of crystal shape at every angle and hence the equilibrium crystal shape(ECS). It is therefore crucial to obtain an accurate value of sur-face free energy so ECS can be directly estimated.

Fig. 1. (Color online) Workflow of Wulff construction: (I) draw a (black line) in which is used as a normalization constant; (II) draw planes(green line) at every point on the that are perpendicular to the line drawn from the origin to that point; (III) Wulff shape is obtained as the inner convex hull and "ears" appear as indications of missing angles of the Wulff shape.

Historically, among solids the surface energy of elemental crystals, mostly metals, was the earliest to be studied by researchers[10]. The surface energies of semiconductors, such as silicon and germanium were later measured[11]. Yet, to date, accurate experimental values for specific facets are still rarely available. Mostly, it was conducted by estimating the solid-liquid interfacial energy and liquid surface tension of the material to calculate the isotropic solid surface energy at finite temperatures, like the melting temperatures, which is then extrapolated to 0 K under isotropic approximations[12,13]. Besides,there were some techniques, such as zero-creep and cleavage techniques, being used to estimate quantitative values of surface energy, which were applied to a limited number of solids, usually metals[14]. Moreover, the 3D equilibrium crystal shapes were also used to estimate the absolute surface energy of a solid[14,15]. These experimental values have been summarized and reviewed elsewhere[12,16-18]. All of these known values are regarded as averages over a range of crystallographic orientations and these values yield large discrepancies by different experimental approaches. In order to precisely control the surface structure and composition, theoretical calculations by first-principles or semi-empirical methods have been conducted to determine such quantities, especially the anisotropy of that for specific facets[19-27]. Due to the recent success of density function theory (DFT) in the investigation of electronic structure of many-body systems in the field of condensed matter physics and quantum chemistry, a high throughput database of surface energies of elemental crystals based on DFT has been published[28].

In the framework of DFT[29,30], slab geometry (nanoribbon geometry for the case of 2D lattices) is widely used to model the surface of metal or semiconductor thin films[28,31,32].Surfaces in semiconducting compounds are different and more complicated. These surfaces can be divided into three types: nonpolar, polar symmetric, and polar or semi-polar and non-symmetric surfaces. For the first two types, they are composed of symmetric surfaces on the top and bottom of the thin film, hence their surface energies can be obtained easily by assuming that the formation energies are identical on both surfaces. The third type of surfaces are asymmetric with a significant dipole moment perpendicular to the surface planes. Among these types of surface, the formation energies are difficult to obtain. Strictly speaking, the formation energies of polar surfaces may diverge for highly ionic crystal,such as alkali halide crystals with a strong dipole field, and charge transfer among surfaces are significant[32-41]. However,for semiconducting compounds, the covalent nature may still make the formation energy of surfaces converge when the thickness of the thin film is large enough. For simplicity, it is therefore preferable to base our discussion on surfaces created from relatively covalent binary compound crystals. In this review article, we focus on the computational aspect of surface energies of semiconducting binary compounds,briefly on non-polar surfaces and edges and mainly on their polar surfaces, semi-polar surfaces, and polar edges for 2D materials.

Before any algorithm is discussed, let us briefly have an overview on the history of the algorithm developments so as to have a general understanding of key advancements in the algorithm design of semiconductor surfaces and edges.

For non-polar surfaces, Feibelman, in 1983, used metallic Al and Mg crystal as examples to demonstrate the calculation of non-polar surface energies[42], which has been also widely applied in semiconductors until now.

For polar surfaces, Chetty and Martin, in 1991, made the first attempt to calculate the absolute surface energy as an integral of local energy density[43]. Even though this preliminary method is not as accurate as the prevalent strategies nowadays, it provided the basic idea of creating an adapted surface unit cell bounded by high symmetry planes for the calculation. Then in 2004, Zhang and Wei proposed a method using different sizes of wedges and a slab for the calculation of polar surface energies, using GaAs as an example[44]. Making use of the subtraction between different wedge sizes, edge(or called ridge) contribution can be cancelled out. The polar surface energy can be estimated by algebraic operation on the energies obtained from both wedges and slab. Later, Rempelet al., in 2005, using CdSe as an example, proposed a similar method to Zhang and Wei but Rempel’s method made use of wedges terminated by both cat- and an-ion at each size in the cancellation procedure[45]. After that, Jenichen, in 2013,proposed a method of using a heterojunction supercell, constructed by both ZB- and WZ-GaAs, to approximate unknown surface energies from known ones[46]. Nevertheless,this approach is still not free of the coupling between conjugated surfaces, which is the major obstacle. In 2016, a tetrahedral cluster was invented by Zhanget al. to calculate the chemical potential of pseudo-H atoms at different surface locations[47], which satisfy the ECM to minimize the stress induced by the passivation. These pseudo atoms were then used to passivate the bottom surface of the slab model so as to mimic the semi-infinite structure. The energies of polar surfaces can easily be obtained. This method can be directly applied to both wurtzite and zinc-blende semiconductors.

For semi-polar surfaces, Liet al., in 2015, proposed a wedge scheme to find the energy difference between surfaces and a reference polar surface in the system of GaN[9].Wedges of different sizes and slabs of both polar and semi-polar surfaces were involved. However, individual surface energies were still not found. Besides, the discrepancy between the Wulff construction and the experimental observation indicated that errors existed in the calculation. Later, Zhanget al.,in 2018, suggested a slab model of GaN with a bottom surface being cut into a zigzag shape composed of non-polar and polar surfaces which can be easily passivated[48]. The semi-polar surface at the top side can then be calculated individually. The steric effect caused by passivation can also be included by simulating its contribution in a slab calculation with a “well” that has similar steric corners. This approach is able to remove the effect of unphysical charge transfer and steric effect in the steps, yielding a high accuracy. Currently,Setaet al., in 2019, proposed a scheme including both passivated wedges and slab of AlN to estimate the surface energies[49], which is a similar method to that proposed by Liet al.but has an improvement on the removal of unphysical charge transfer by passivation. In addition, the temperature effect had been taken into account by calculating the partition function of translational, vibrational, and rotational motions of the atoms. However, direct passivation of semi-polar surface was involved, which may induce the steric effect between pseudo-H atoms and finally affect the accuracy of surface energy estimation.

For polar edges of 2D materials, Mukherjeeet al., in 2011,presented a strategy with pure bare triangular clusters of h-BN in single size, each of which has only one type of exposed edge, to calculate the edge energies individually[50]. It is now known that the unphysical charge transfer substantially affects the accuracy of the algorithm. In 2015, Caoet al.offered a scheme as an extension to Mukherjeeet al.’s method[51]. Bare triangular MoS2clusters of different sizes were used so that the error from corner contribution is reduced by energy subtraction between nanoclusters of different sizes.This algorithm can only predict the morphology of MoS2at the S-rich limit while discrepancy with the experimental observation in the Mo-rich limit indicates there are pitfalls in the calculation. There could be unphysical charge transfer at the corner of triangular clusters, which was later shown to be removed by passivation, and the omission of temperature effect. Later, in the study of morphology of h-BN nanoclusters,half-passivated nanoribbons were created to estimate the energies of polar edges individually by Zhanget al. in 2018[3]. Ordinary H atoms were employed in the passivation, whose chemical potential was obtained from the quadratic fitting of the equation formulating the total energy of fully-passivated triangular nanoclusters. In addition, the temperature effect was included by extrapolating the extra Gibbs free energy of H atoms from the experimental data. The morphologies obtained from the theoretical calculation are consistent with the experimental works.

After the historical review, the algorithm designs for the calculation of different surface and edge energies are going to be discussed separately as follows.

Before diving into polar and semi-polar surfaces and polar edges, it is good to have a brief review on the non-polar surface or edges, and compare them with polar and semipolar ones. For slab (ribbon) exposing non-polar surfaces(edges), the constituent elements are in the stoichiometric ratio within the whole structure. Therefore, there is no need to consider the chemical potential contribution from individual elements. Also, since the top and bottom surfaces are identical, it is possible to obtain the formation by assuming both surfaces have the same contribution to the total formation energy. This idea was first proposed by Feibelman in metallic materials[42]. We can directly apply the following equation to obtain the surface energy ?

Fig. 2. (Color online) A slab created by cleaving a zinc-blende structure in (111) plane, grey and yellow atoms represent atom species A and B. Note the resultant upper and lower surface is of different termination.

Fig. 3. (Color online) Wedge structure of size n = 4 composed of two equivalent (111) and one (001) surface used in the calculation scheme of Zhang et al. [44].

For symmetric surface (edge for 2D lattice) pairs, the surface (edge) energyis given by

2. Polar surfaces

At the earliest stage of polar surface study of semiconductors, fractionally charged pseudo hydrogen (pseudo-H) has been suggested to passivate one of the surfaces of the slab model so as to remove the charge transfer between surfaces[53]. This approach allows us to study the polar surfaces individually. However, at the time of publication, the removal of the electrostatic effect between surfaces by this kind of passivation has been shown, but has not been implemented in the surface energy calculation. The first attempt to calculate absolute surface energy defined a local energy densityand integrated over the surface region[54]. Chetty and Martin[43]generalized the symmetry argument of Appelbaum,Baraff and Hamann[55]to low symmetry surface by creating a symmetry-adapted surface unit cell bound by high symmetry planes, and computed the surface energy as an integral of local energy density. This method is first applied to the calculation of ideal GaAs (111) andabsolute surface energy, bridging previous studies of relative surface energy on these surfaces to show that Ga-trimer and As-trimer reconstructedsurfaces are always more stable than reconstructions on (111) or (110) surfaces[56]. However, calculations by Mollet al. using the same procedure failed to agree on the splitting of slab energy between the (111) andsurface,leading to a discrepancy that Ga-trimer is significantly less stable[57], possibly due to a nontrivial approximation in the local energy density[44].

Zhang and Wei[44]first attributed the origin of surface energy to the local bonding environment of individual surface atoms, and proposed a more sophisticated geometry to be constructed to compute the absolute surface energy of polar surfaces from surfaces of known energies, applicable to surfaces that are inaccessible by the standard slab method from surfaces of known energies. In their work[44], an infinite wedge geometry with 2 equivalent surfaces and one distinct surface (Fig. 3) is used. Under this geometry, the total surface energy of the wedge can be attributed to energies from three surfaces and three edges. The unknown surface energy can be solved by taking the energy difference between two reasonably large wedges of different sizes so that the converged edge contribution cancels out. Rempelet al.[45]also considered the convergence with respect to wedge size, and proposed to base the cancellation on the difference of edge energies by considering both cat- and an-ion termination at each size. Using this wedge method, Zhang and Wei[44]shows a result similar to Moll[57]in GaAs, and the method was applied to the CdSe polar surfaces[45,58]and GaN with surface passivation to reduce charge transfer between surfaces[59].

Fig. 4. An illustration of a slab containing an interface and two passivated surfaces.

Fig. 5. (Color online) A ZB(111)/WZ(001) heterojunction supercell consists of 6 WZ and ZB layers used in the calculation scheme of Tang et al.[4]. Note the two interfaces indicated by dashed lines are inequivalent in that ion termination at the interface exchanged.

A heterojunction supercell (Fig. 4) was also used to approximate unknown surface energies from known ones[4,46]. The total energy of the slab is assumed to consist of additive components from two surfaces, interface and bulk[46],

In Eq. (7), chemical potential has to be obtained from strained bulk simulation using the slab lattice constant to account for stress due to lattice mismatch. The heterojunction scheme is employed to study the (001) andsurface of wurtzite GaAs[46], and subsequently ZnO[4]. The heterojunction geometry is similar to the slab approach, thus is not free of the averaged energy problem by construction, it can be seen in ZB(111)/WZ(001) interfaces as shown in Fig. 5 that the heterostructure would consist of two inequivalent interfaces. An averaged interface energy of two interfaces may substitute for exact interface energy, smaller interface energy compared to surface, interface being coherent, similarity in lattice parameters and local environment may justify such approximation in this work[4], however, a small lattice mismatch and coherent interface is not always possible.

Fig. 6. (Color online) Tetrahedral cluster of zinc-blende structure of size n = 4 composed of identically passivated (111) surfaces, passivated edges and corners. The figure is adapted from Ref. [60].

With the key assumption of localized energy contributions, the energy of an isolated crystal can be computed[44,58]as

The passivated surface energy can then be systematically obtained by cancellation of bulk, edge and corner contributions. The transferability of such an estimation scheme mainly depends on 1) the ability of the proposed structure to capture the local bonding environment of the surface to be estimated, and 2) good size convergence of the proposed structure so that systematic cancellation of other contributions is possible. Zhanget al. proposed a tetrahedral cluster reproducing the same symmetry as the zinc-blende[47]. The definition ofis taken to be the energy associated with the surface (edge, or corner) divided by the number of passivation hydrogen attached to A (or B) involved at the surface (edge, or corner). Pseudo-chemical potential (PCP)andrepresenting the atomic chemical potential plus the passivation contribution can then be normalized by simple surface atom counting. A crude estimation of pseudo-chemical potential can be obtained from the smallest cluster of 4 pseudo-hydrogen passivating a host atom A as

This estimation was shown to yield an acceptably accurate surface energy estimation[47]. More physical estimations of pseudo-chemical potential can be constructed by observing that the nature of passivating hydrogen is mostly influenced by its local bonding environment, which can be classified according to the hydrogen attachment to a host atom A at surface, edge or corner. For a zinc-blende structured cluster as shown in Fig. 6, the total energy of a cluster of sizecan be expressed

Calculation of 4 clusters of different sizes allows solutions for all unknownsIt should be noted that althoughis physically interpreted as the bulk energy which should be equal to that from bulk calculation, solving for the termavoids error due to bulk energy differences across bulk and cluster calculations.

With a fractional hydrogen[53]passivation scheme, in the case of multiple surface dangling bonds, more than one hydrogen is usually required per surface atom to satisfy the electron counting model (ECM) requirement and to preserve surface symmetry. Passivating hydrogen atoms can become too densely placed and induce a strong steric effect. The steric effect is particularly severe in systems with small lattice constants, and in the directions where periodicity is absent so uniform distribution of stress is not guaranteed. The significant non-uniform stress produces visible distortion and longranged variation in bond lengths and bond angles and results in a slow size convergence, which deteriorates cancellation of edge contribution in the wedge method between small sized wedges[60]. Zhanget al.[60]introduced a passivation by other atoms to compute the absolute surface energy of ZB/WZ ZnO and GaN, a pseudo-chemical potential is defined per surface passivating atom similar to the procedure in the tetrahedral cluster method. The choice of surface passivating atom should a) satisfy ECM same as the case of hydrogen; and b) have a correct size that minimizes the stress induced by the passivation, to achieve better accuracy and convergence with respect to wedge size.

Absolute interface energy is of equivalent physical interest including the determination of the wetting condition[61]and band offset[62]. It is also interesting to note that the connection between surface and interface energy established in a slab construction (Fig. 4) allows the determination of the absolute interface energy. To determine the absolute stability of interfaces, a common energy reference among all interfaces at different terminations and orientations must be known. In traditional superlattice approaches[4,63,64], an interface cannot be created independently of the other complementary interface, which prohibits the determination of the absolute interface energy. The superlattice approach is thus only capable of comparing relative stability of reconstructions at a single interface[63], or the calculation of the averaged interface energy of the two complementary interfaces[4]. Coupling between two interfaces would also induce unphysical charge transfer and dipole-dipole interaction. Akiyamaet al.[62]and Zhang[61]independently demonstrated the viability of the slab construction to determine the absolute interface energy, with the former utilizing the wedge method with passivation and the latter utilizing the pseudo-hydrogen method for surface. Similar methodology is used to compare the stability between different terminations and reconstructions of{112} grain-boundary insuggesting the engineering importance of this method.

3. Semi-polar surfaces

Fig. 7. (Color online) GaN crystal with 3 different types of surface cut.The semi-polar one is highlighted pink.

Over the past few decades, even though the technologies in industrial application, like the quantum dot light-emitting diodes (LEDs)[66,67], of WZ based semi-conductors have been becoming mature[33-38], it is still a challenge to have high quality crystal growths as well as control of their morphologies[39-41,60]. On account of the substantial achievements in InGaN based LEDs, Group-III nitrides have drawn a lot of attention[36]. To date, there is a major limitation to GaN-based optoelectronic devices that only blue emitters have been produced by polar GaN grown on ac-plane (0001) sapphire[33-36,68]. It it also unfavorable to manufacture greenand yellow-light LEDs with high efficiency, given that high quality InGaN with a high indium concentration is a prerequisite[41,68-71], because of the miscibility gap led by significant atom mismatch among indium and gallium and the piezoelectric effect upon polar surfaces. The growth along non-polar or semi-polar plane in GaN thin film could be a suitable strategy to overcome this critical barrier, especially in a semipolar direction[72-79], it is because the recruitment of relatively larger indium atoms can be assisted by any site that is under tensile stress[41,75]. Besides, owing to the fact that there is a weaker piezoelectric effect among semi-polar surfaces[80],leading to both intensified indium recruitment[68,75,76]and weakened quantum confined Stark effects[80-82]. Nevertheless, different from polar surfaces, there was a lack of elementary knowledge about semi-polar surfaces because no workable algorithm was present for the calculation of absolute surface energies of individual semi-polar surfaces.

The definition of a semi-polar surface, with an example shown in Fig. 7, was initially made by Bakeret al.[83]as surfaces being cut in the planes with one of theh,k, andiindex andlindex being nonzero in the (hkil) Miller-Bravais indexing convention, crossing the hexagonal unit cell diagonally and forming a non-orthogonal angle with thec-plane.

A comprehensive understanding of the absolute surface energies of all possible GaN surfaces is crucial to the estimation of equilibrium shape in the thermodynamic stability study, leading to important factors that can be used to modulate the crystallographic growth of GaN, which are regarded as the key issue in the realization of broadband and multi-color emission[84-91].

There is an early method proposed by Jindal and Shahedipour-Sandvik in 2009[92], that tries to find the energies of different GaN surfaces so as to create the ECS. However, only the average value of two conjugated surface energies can be acquired. In order to calculate the ECS in a more precise manner, the WZ wedge scheme was employed to find the energy differences between surfaces and a reference polar surface in 2014[9]. They proposed that since individual surface energies cannot be directly obtained, Eq. (6) is not applicable.Therefore, they first found the difference in surface energies between semi-polar and polar surfaces, then the surface energies of semi-polar surfaces can be obtained by further addition and subtraction of polar surface energies. The detailed procedures are illustrated below.

Fig. 8. (Color online) Wulff construction of one of the 2D cross-sections of GaN. The yellow shaded area is a quarter of ECS in the cross-section. This strategy is from Ref. [9].

Fig. 9. (Color online) Workflow of finding the difference in crystal plane radii. Blue and black notations correspond to unrelaxed and relaxed surface structures respectively. This strategy is from Ref. [9].

A 2D scheme of Wulff construction is applied to one of the cross-sections of GaN as indicated in Fig. 8. The length ofis fixed but the origin can be varied in position between A and B. The position of origin, which is unknown now, depends on the individual energies between surfaces (0001)andis the difference in crystal plane radii of planeand the plane (0001) in [0001] direction. We can do the same trick for planeto obtainHence, we obtainFor differentandvalues, the blue lines can be shifted vertically. After obtainingand, the origin can be located by the normal vectors of planesanda quarter of the Wulff construction is done.

However, based on the experimental observation from[93], the nanocrystals did not show the appearance of asurface while the simulated results contain this surface for a large range ofThis discrepancy indicates that their results are rather approximated and also it is difficult to judge by individual surface energies since only relative energies were found. Besides, passivation was not performed on the wedge structures so that unphysical charge transfer should be present, being especially severe at the corner. The relative stability of different surfaces could be affected so as to alter the ECS. Therefore, this algorithm may not be able to give accurate predictions on the equilibrium shape in the growth study.

The acquisition of accurate energies of semi-polar surfaces is particularly difficult for three reasons: firstly, the conventional slab method cannot be used to deal with individual semi-polar surfaces resulting from the structural asymmetry; secondly, large computational input in the form of wedges, usually with high index planes, are involved which leads to high computational cost; thirdly, as semi-polar surfaces are sometimes of a step nature, it is not always feasible to passivate the bottom surfaces of slabs with pseudo-H atoms in the absence of significant unphysical charge transfer and steric effect, which deteriorates the result accuracy[47,60].

To overcome these mentioned difficulties, Zhanget al. introduced a fundamentally different algorithm in 2018, using GaN as an example[48]. The practice of passivation on polar surfaces is generally not applicable to semi-polar surfaces because a substantial steric effect may appear among the passivating agents and ECM may not be easily satisfied[94]. A feasible approach is to cut the semi-polar surface of the slab with a step-like structure so that the surfaces exposed are instead non-polar and polar which can be passivated by suitable pseudo-H much more easily, as shown in Fig. 10(a). Nevertheless, solely applying the above modification may not fully solve the problem because the steric effect may still be present at the location of included angles (or corners)between non-polar and polar surfaces, indicated in Fig. 10(b).Thus, a unique treatment can be implemented to take the steric effect into account. This new strategy is essentially different from the early treatments in the calculation of surface energies and is generally applicable to the study of other high indexed surfaces which are difficult to handle.

According to Eq. (3),

Fig. 10. (Color online) (a) and (b) are slabs with upper semi-polar surfaces of m- and a-family, respectively, and with bottom side cut into step-structure in which the non-polar and polar surfaces are passivated by either H or pseudo-H. These figures are adapted from Ref. [48].

Fig. 11. (Color online) Slab with a well being cut with width and depth as w and d, respectively, that mimic the steric effects between pseudo hydrogen at the concave corner between the polar and non-polar plane. This figure is adapted from Ref. [48].

However, taking the semi-polar surfaceas an example, the passivating H atom may experience the steric effect at the concave corner of bottom surface as indicated in Fig. 10(b). There was another treatment proposed to include this steric effect[48]. As indicated in Fig. 11, a slab with a well is created in which all surfaces are passivated with pseudo hydrogen. The steric effect on the pseudo-H passivating the Ga atoms can be calculated by

Zhang has used a slab with both sides cut with a zigzag structure to implement the convergence test to give a residual error less than 1.5 meV/Å2, indicating the high accuracy of the method. This new algorithm to estimate the absolute energy of semi-polar surface is completely different from the traditional methods which are based on wedges and slabs. It is because this new method is applicable to an arbitrary surface as long as we can passivate the polar and non-polar planes at the bottom surface with zigzag structure.

Later on another Japanese group, Setaet al., published some literature in 2019 using both slabs and wedges to estimate the absolute energy of semi-polar surfaces[49]. His scheme is similar to that proposed by Li in 2014[9]except that Seta modified the surface of the wedges by passivation of hydrogen so as to remove unphysical charge transfer, indicated in Fig. 12.

After calculating the energies of wedges in Fig. 12 for different sizes, the same procedures were repeated by interchanging the position of Al and N atoms. Therefore, the energy difference (here we use the notationfor surface energy instead of) between passivated surfaceand) is given by

Fig. 12. (Color online) Cross section view of AlN triangular wedge with surface and (0001) which are passivated by hydrogen. Orange, silver and red spheres represent Al, N and H atoms respectively. The area bounded by a black line demonstrates that the removal of the area creates a smaller wedge. The strategy is from Ref. [49].

Seta has given one more improvement on the temperature dependence on the estimation of surface energies by including the translational, rotational and vibrational motion of atoms in gaseous phase into the chemical potential

Nearly at the same time, Akiyama, from the same Japanese group, proposed another algorithm to calculate the energy of polar and semi-polar surfaces simultaneously[96]. In the literature, it was claimed that Zhang’s method may not be general because he has used the ZB tetrahedral cluster to obtain the PCPs of passivating hydrogen, which were then applied to the WZ structure of GaN. The error coming from the deviation between ZB and WZ could be large if the materials have large ionicity. Also, he claimed that bond length between cation and pseudo-H could be too long to form the tetrahedrally coordinated atomic configuration so that the method is not applicable to compounds with small atomic size like BN and AlN. Therefore, he continued to use the scheme of creation of ideal slabs and wedges without any passivation of pseudo-H to formulate multiple equations describing the relation between different polar and semi-polar surfaces. Under the assumption of orientational independence of twofold or threefold coordinated Ga and N surface atoms,the energies of the twofold or threefold surface atoms can be obtained and used to estimate the energy of semi-polar surfaces by counting the number of twofold or threefold surface atoms on the corresponding plane. However, the avoidance of using pseudo hydrogen may induce a very severe error. Without the passivation on the surface atoms of the semi-polar plane, the degree of distortion will be different in different orientation so that the twofold or threefold surface atoms will experience a different local electronic and strain environment. To avoid such distortions and keep the orientation independence, the author, instead, calculated the unrelaxed structures to obtain the absolute energies. This eventually may lead to a large error due to the unphysical stress within the unrelaxed structure. It can be easily observed by the average of the absolute energy of polar surfaces (0001) andwhich is significantly larger than the average values obtained from the conventional slab. In addition, the fundamental assumption of orientation independence may not be correct because when the structures are allowed to relax, the number and angle of bondings depend on the local environment. Therefore, the idealization made in the literature may not be applicable in the description of the real situation. In general, the key success in the high accuracy algorithm in the cluster or modified slab method[48]is largely due to the localization of charge density by pseudo-H passivation. The passivation energy is largely transferable across different microstructures. Therefore, the stability of the algorithm will be sacrificed if the passivation is removed.

4. Polar edges of 2D materials

In the world of 2D materials, graphene is the system that was first discovered and the most intensively studied[97-101].Its growth is relatively simple compared compounds, because it consists of one element. Thermodynamically, the ECS is relatively easy to be investigated because no matter how it is cut, the edges are always non-polar and hence the calculation of the edge energy is straightforward[101].

Fig. 13. (Color online) (a)/(c) and (b)/(d) are the (top view/side view) of h-BN and nanoribbon with edges of opposite polarities, respectively.

In recent years, the family of 2D materials has grown larger, including compounds such as hexagonal boron nitride(h-BN) that is an insulator, and molybdenum dichalcogenideswhere X can be S, Se or Te, which are semiconductors[102]. The latter are often considered quasi-2D materials because they consist of several layers of atoms, unlike graphene or h-BN. Different morphologies of h-BN, MoS2and MoS2were observed in experiments. Triangle, truncated triangle and hexagon shapes were observed in both h-BN and MoS2while fractal, three-point stars and multi-apex triangles were observed in MoS2[103-105]. Usually convex edges suggest that the growth is near the equilibrium limit, while concave edges suggest a relatively non-equilibrium limit. Therefore,for triangle, truncated triangle, and hexagon shapes, the equilibrium shape can be investigated by the energy calculations under various growth conditions. However, other shapes with concave edges are usually due to the growth highly limited by kinetics and may often be off equilibrium[106], leading to a high complexity in the simulation of the reaction. In the following, only the shapes in thermodynamic equilibrium will be reviewed.

To understand various equilibrium shapes in experiments, obtaining the energies of different edges are crucial.Due to the structural asymmetry of the nanoribbon of these compounds, as shown in Fig. 13, edges with polarities emerge. The direct calculation of polar edge energy is no longer achievable because the early method used in graphene can only estimate the average energies of two opposite zigzag edges in the ribbons[101].

In the past few years, several groups have developed their own methods to estimate the absolute energy of polar edges. Most of them are based on the creation of triangular nanoclusters terminated at the edges with the same polarity[3,51,107]. Therefore, the excess Gibbs free energy[108]can be obtained by deducting the bulk chemical potential contributions, and the energy of the polar edge can be obtained after the excess energy is divided by three. Even for similar approaches, there are still some technical details that may significantly affect the accuracy. In the following, three examples based on first principle calculations using DFT will be examined and compared so that we can see the pitfalls of some early methods.

In all methods, polar edge energies can be obtained under different chemical potentials. The range of the chemical potentials can be obtained from a standard calculation of phase diagram of various secondary phases of the compound, which has been widely applied in the energy calculations of point defects, surfaces, and interfaces[48,65,109].

Besides chemical potential, passivation of edges is the key to the calculation of edge energy. For early methods, passivation was not taken into the consideration[51,107]. The omission of edge passivation may lead to severe errors in the estimation of edge stability and equilibrium shape. Experiments and theory[110,111]both suggest that a graphene sheet grown next to the transition metal step edge has lower energy for both zigzag and armchair edges, as a result of charge transfer from the step edge of the transition metal substrate to the edge of nano-structures. In addition, boron nitride grown in the direction with highest substrate atom packing has higher stability[112]. It was suggested that it is easy for the electron transfer from the highest packing plane of transition metal substrate to the edges of nano-structures[112]. Also, to model the asymmetric polar edges properly, it is essential to remove their interaction by passivating one side. During the growth, decomposition of precursors or carrier gases may lead to different passivation configurations that may strongly affect the final equilibrium shapes. Therefore, it is important to include the edge passivation effect in the estimation of edge and shape stability in 2D or quasi-2D materials.

The first theoretical attempt was made by estimating the average energies of h-BN edges terminated by B and N[50].However, such an estimation is too rough to derive accurate equilibrium shapes. Later, triangular clusters were created so that only one type of zigzag (ZZ) edge, which is polar, was investigated at a time[107]. In that paper, the energy of the polar edges was formulated by the following equations

Fig. 14. (Color online) (a) The computational setup for triangular clusters with green dots as boron atoms and silver dots as nitrogen atoms. (b) The result of equilibrium shapes at different chemical potential ranges[107] in which blue, red and black are ZB, ZN and AC edges, respectively.

The results suggest that it is possible for armchair-edged hexagon to exist at the mid-range of chemical potential. In addition, other literature gives the same computational prediction on the stability of the armchair edge[50,113]. However,there has been no experimental observation of the existence of armchair hexagons or truncated triangles with armchair corners, while the alternating B- and N-terminated hexagons was observed experimentally[103]. Therefore, this calculation scheme may contain pitfalls in simulating the physical system under real conditions.

Fig. 15. (Color online) (a) The triangular clusters of different sizes enclosed by two triangles with length (b) Structure of four main type of zigzag edges. Other reconstructions of edges were also studied but they are not important in the final equilibrium shape. For details please refer to Ref. [51].

Fig. 16. (Color online) Equilibrium shape of under different chemical potential of S atoms constructed by Wulff construction. The figure is adapted from Ref. [51].

The bare triangular cluster was reported to contain corner distortions and inter-edge couplings[114], which may produce extra systematic errors and lower the calculation accuracy. Usually, the edges of the clusters have a different electronic structure from the ribbons, with significant charge transfers and atomic reconstructions, which are especially severe near the corners. The most ideal situation is that all the edges in the triangles have the same electronic environment as those in the ribbon. To mimic the electronic environment of edges in the ribbons, it is important (1) to passivate one side of it and remove the inter-edge interaction; (2) to estimate the passivation energy of hydrogen by constructing triangles with passivated edges. However, this literature does not include detailed procedures in the passivation. Only in a later example[3], a more detailed scheme was provided. During the CVD growth of BN, a large amount of hydrogen may exist due to the decomposition of the precursors and passivate the edges, the passivation energy and related reconstructions were not systematically explored. In addition, the temperature effect of the passivation may play an important role in the stability[3]and was not considered.

Another method proposed by Caoet al.to find the equilibrium shape of MoS2[51]is better than the previous one in handling the corner problem in the polar-edged triangular cluster. The structural inputs and different types of polar edges being studied are shown in Fig. 15. The energies of the triangular clusterwith side lengthare

Fig. 17. (Color online) B, N and H atoms are denoted by pink, blue and white spheres respectively. (a) Passivated and unpassivated zigzag and armchair edges. (b) Reconstruction of seven- and five- rings on the ZZN and ZZB edges, respectively. (c) Ribbon of bottom zigzag edged passivated with hydrogen and arbitrary configuration on the upper zigzag edge. (d) N-terminated passivated triangular cluster of size m = 5. (d) Bare N-terminated triangular cluster with corner distortion as indicated by red circle. (a) Ribbon with fully passivated zigzag edges. The figure is adapted from Ref. [3].

From the Wulff construction, the S-terminated triangular shape can be observed in an S-rich condition in which the shape matches the experimental results[104,115,116]. Besides,this S-rich edge is the ZZ-S2 edge which is a Y-shape rather than a pure zigzag structure, confirming the experimental result in Ref. [116]. However, when the chemical potential of S was reduced, a truncated triangular or hexagonal structure with ZZ-Z and ZZ-Z2 were predicted. This does not match the result of[116]in which the hexagonal structure is terminated by an alternating S-monomer attached Mo-edge and hydrogen-passivated S-edge. This discrepancy could probably be attributed to the exclusion of calculation that the S-monomer attached Mo-edge and also to the omission of consideration of the temperature effect. Therefore, the calculations failed to yield the proper experimental equilibrium shapes in the whole chemical potential range.

Before entering the last example, it is good to mention that the method of Cao's example is suitable for obtaining a fast calculation. Yet, there is another method proposed by Zhanget al.[3]dealing with h-BN that can reduce the error fromby bare triangular toby passivated triangular cluster. In his calculation, the importance of restoring the original bulk electronic configuration was emphasized.Theoretically, if we want to calculate the energy of an arbitrary edge, we have to create a semi-infinite crystal with only one edge exposed. This condition is imitated by passivation on one of the edges on a nanoribbon as shown in Fig. 17(c).The absolute energy of, for example, a zigzag boron (ZZB)edge is given by

Therefore, instead of a direct calculation of edge energy from the bare triangular cluster, the chemical potentials of passivating hydrogens have to be first estimated from the passivated cluster. The reason for the more reliable calculation of the chemical potential of passivating hydrogen than the edge energy of the bare triangle is that the passivation helps to reduce corner distortion and the unphysical charge transfer[3]which can be compared in Figs. 17(d) and 17(e). In the cal-culation of chemical potential of passivating hydrogen in a triangular cluster, different size (m) of passivated triangular clusters are calculated so as to obtain their total energy

Fig. 18. (Color online) Total energy of H-passivated triangular clusters with different size (m) and the corresponding non-linear fitting. The figure is adapted from Ref. [3].

wheremis the cluster size andis the chemical potential of hydrogen atoms at corners, as shown in Fig. 17(d). Practically, we can obtainunder different values ofandby non-linear fitting. In the thermodynamic equilibrium between the edges and the bulk h-BN, from Eq. (3),

There are three red-colored parameters to be fittedTherefore, a quadratic fitting is needed afterof differentmare obtained. The fitting graph is shown in Fig. 18. A direct check for the accuracy of the fitting is to compare the fittedwith a separate calculation of the 2D monolayer.

After the estimation of hydrogen chemical potential, the half-passivated ribbon with arbitrary configuration on the upper edge (Fig. 17(c)) can be calculated to obtain the absolute energy of the particular polar edge by Eq. (22). Also, Zhang has proposed a self-consistency check to ensure the accuracy of the algorithm by calculating the residual errorEr.Ercan be calculated by Eq. (26) after the calculation of total energy of both sides passivated ribbonEp. Zhang has shown the error is reduced from 3.43% to 0.12% when compared with the bare triangular method.

After that, Zhang had shown in Fig. 3 of their literature[3]that bare armchair (ARM) is the most stable, which matches the computational results given by[50,113]but does not match the experimental results[103]. Therefore, Zhang included the calculation of a passivated edge (upper edge) in which the chemical potential of the hydrogenis half of the total energy of H2moleculesAlso, he had included the consideration of the temperature effect by

After the discussion of several methods of calculating the polar edge energy, the last one proposed by Zhang, is able to capture more physical pictures and also gives the highest accuracy because it includes the temperature effect and passivation which leads to stabilization effect to all type of edges. It is also capable of revealing the important role played by hydrogen atoms in the growth of 2D h-BN monolayer.

5. Conclusion

Fig. 19. (Color online) Equilibrium shapes of h-BN nanocluster under different chemical potentials at 1300 K, consisting of H-passivated edges. Yellow, green and black lines are of ZZBH, ZZNH and ARMH edges respectively. The figure is adapted from Ref. [3].

We have reviewed some important historical algorithms on the assessment of surface and edge stability of various semiconducting compounds. The key concept for a successful algorithm is to eliminate the long range charge transfer and interaction of different surfaces or edges by passivating dangling bonds and mimicking the electronic environment of the desired surfaces or edges. In addition, not all passivation can yield a reliable result because an electron counting model has to be satisfied and steric effects should be avoided. To estimate the localized steric effects, it is possible to perform further simulation that can mimic the stressed local configuration. Still, further investigations of quasi-2D structures are highly important, yet largely missing because they lack effective passivation schemes on the edges. With all the technological advancements, we can safely conclude that a highly accurate algorithm combining a reasonable analysis of passivation and temperature effects can have strong predictive power in the equilibrium shape under various growth conditions and the dawn of a highly effective collaboration between theoreticians and experimentalists may largely improve the field of crystal growth and device fabrication of semiconductors.

Acknowledgements

The research is supported by HKRGC, GRF with the Project Codes of 14307219, 14307018, 14301318, and 14319416; and by direct grant from CUHK.