I. INTRODUCTION
Some group IV semiconductors (carbon, silicon, and SiC) and group III–V compounds crystallize in a cubic and/or a hexagonal form and exhibit polytypism, in which the nearest neighbor atoms are tetrahedrally coordinated and have sp 3 hybrid characters. The polytypes are characterized by a stacking sequence with a given repeat unit along the hexagonal c-axis direction and are theoretically possible to have endless permutations of the sequences.Reference Estevez-Rams and Martinez-Mojicar1 Among the polytypes for these elements and compounds, the silicon carbide (SiC) systems have been the most attractive and motivated tremendous amounts of theoretical and experimental investigations of their fundamental and technological properties for many years.Reference Bernstein, Gotsis, Papaconstantopoulos and Mehl2–Reference Rutter and Heine9 The SiC polytypes differ only in the stacking order of double layers of carbon (C) and silicon (Si), and more than 200 polytypes have experimentally been determined to date.Reference Jepps and Page10
From the theoretical view, the problem of mathematically generating all nonequivalent polytypes for a given length (L) shares common problems with the counting problem. Some mathematicians and/or crystallographers have ardently been challenging this issue.Reference Estevez-Rams and Martinez-Mojicar1 Because of the exponential explosion with L, smart counting procedures that avoid the necessity of generating all the sequences have been necessary to effectively count the number of nonequivalent polytypes for a given length of L. Estevez-Rams and Martinez-MojicarReference Estevez-Rams and Martinez-Mojicar1 have recently presented an efficient polytype generation algorithm and have released a computational source code generating the hk notationsReference Jagodzinski11 for nonequivalent polytypes on the International Union of Crystallography electronic archives, so that we can use their method conveniently.
Using the first-principles calculations within the density functional theory (DFT), some papers have reported the correlations between energetics and stacking configurations in the polytypes for group IV semiconductors and group III–V compounds.Reference Bernstein, Gotsis, Papaconstantopoulos and Mehl2–Reference Rutter and Heine9 However, due to the computational cost, the data of energetics by the first-principles calculations, by all means, cannot help becoming discrete ones for possible configurations. The trends of material properties such as the correlations between energetics and hexagonality, therefore, have not fully been clarified. In addition, since there are few comparative studies on energetics and electronic structure for group IV semiconductors and group III–V compounds,Reference Raffy, Furthmüller and Bechstedt3, Reference Kobayashi and Komatsu4, Reference Yeh, Lu, Froyen and Zunger12, Reference Wright13 the chemical trends on the structural and energetic properties for these polytype compounds have not been comprehensively understood.
In this paper, also motivated by the work of Estevez-Rams and Martinez-Mojicar,Reference Estevez-Rams and Martinez-Mojicar1 we have systematically revaluated the energetics of polytypes for carbon, silicon, SiC, BN, AlN, GaN, and InN including some not previously calculated ones. For all these elemental systems, the atomistic geometry, energetics, and electronic structure for these compounds with up to the periodic stacking length of L = 8 have been carefully calculated based on the DFT within the generalized gradient approximation (GGA). Using the Axial Next-Nearest-Neighbor Ising (ANNNI) model extracted from the GGA calculations, we have also studied the energetics for more than 6 million kinds of nonequivalent stacking polytypes with up to L = 30, whose configurations have been deduced by the efficient polytype generation algorithm,Reference Estevez-Rams and Martinez-Mojicar1 and illustrated some trends of physical properties for these compounds.
II. THEORETICAL
A. Crystal structures and polytype notations
The polytypes considered in the present work can be basically characterized as a close-packed stacking sequence of double layers (bilayers) for IV–IV and III–V elements in the <111> direction of the cubic zincblende (or diamond) structure. In these close-packed structures, three different layers are often denoted by A, B, and C, corresponding to different atomistic positions in the planes perpendicular to the stacking direction (c-axis). These descriptions for nearest neighbor stacking result in the so-called hk notation,Reference Jagodzinski11 where a letter h is assigned to each layer with the same type of layer at both sides (e.g., AXA, BXB, CXC) and a letter k is assigned otherwise. In addition, using this hk notation, the “hexagonality” (σH) is essentially defined as σH = n h/(n h + n k), where n h and n k are numbers of layers in a hexagonal or cubic environment, respectively. The hexagonality is often used to classify the physical properties such as structural parameters, total energetics, phonon mode, and electronic strucure in the polytypes constructed of hexagonal close-packed stacking.Reference Bernstein, Gotsis, Papaconstantopoulos and Mehl2–Reference Rutter and Heine9, Reference Yeh, Lu, Froyen and Zunger12, Reference Wright13 Another notation commonly used in the representation of stacking sequences is the Hägg notation.Reference Sebastian and Krishna14 In this notation, the cyclic stacking orders of AB, BC, and CA are all denoted by up-spin (+), whereas the reverse stacking orders (BA, CB, and AC) are denoted by down-spin (−).
In this paper, the character of U(D) is used instead of the symbol of +(−). The Hägg and hk notations are bijective between them and also with respect to the ABC notation. Three nonequivalent stacking polytypes for the periodic length of L = 7 are depicted along with the ABC, Hägg, and hk notations in Fig. 1.
Using the computational code presented by Estevez-Rams and Martinez-Mojicar,Reference Estevez-Rams and Martinez-Mojicar1 we can easily get the nonequivalent polytypes with any periodic stacking length of L. There are only one possible polytype with up to L = 5 (2H, 3C, 4H, and 5H), while it is possible to consider two, three, and six kinds of nonequivalent crystal structures for 6H, 7H, and 8H, respectively. Fifteen kinds of crystal structures for each elemental system (i.e., carbon, silicon, SiC, BN, AlN, GaN, and InN), therefore, have been calculated by the first-principles method in the present work. We have tabulated the hk notations and space groups for the crystal structures considered by the first-principles method in Table I. The hexagonality, σH, for each structure is also listed in this table. In this paper, each structure is referred to using the symbol tabulated in the second column in Table I.
Schematic illustration of constructing nonequivalent stacking polytypes exemplified in the case of periodic length of L = 7. Three nonequivalent stacking polytypes for the periodic length of L = 7 are shown along with the ABC, Hägg, and hk notations.
Structural notations for nonequivalent polytypes calculated within the GGA. The hk notations and space groups are listed for the polytypes considered. In this paper, each structure is referred to using the symbol tabulated in the second column. Hexagonality is defined as σH = n h/(n h + n k), where n h and n k are numbers of layers in a hexagonal or cubic environment, respectively.
B. First-principles calculations
For carbon, silicon, SiC, BN, AlN, GaN, and InN systems, the atomistic geometry, energetics, and electronic structure for these compounds with up to the periodic stacking length of L = 8 have systematically been calculated based on the DFT. The first-principles calculations for the periodic boundary systems are carried out by the ab initio pseudopotential theory within the GGA. We use the CASTEP codeReference Clark, Segall, Pickard, Hasnip, Probert, Refson and Payne15 to solve the pseudopotential Schrödinger equation self-consistently. We adopt the ultrasoft pseudopotential introduced by Vanderbilt.Reference Vanderbilt16 The semicore d electrons have been treated as valence states for Ga and In. The wave functions are expanded on plane-wave basis sets with the kinetic energy cutoff of E cut = 1000 eV for all the systems studied. The Perdew–Burke–Ernzerhof functional of the GGAReference Perdew, Burke and Ernzerhof17 is used. As for the Brillouin zone sampling, the Monkhorst–Pack gridsReference Monkhorst and Pack18 of 13 × 13 × 8, 13 × 13 × 4, 13 × 13 × 4, 13 × 13 × 2, 13 × 13 × 2, 13 × 13 × 2, and 13 × 13 × 13 are adopted for 2H, 4H, 5H, 6H, 7H, 8H, and 3C configurations, respectively.
A high precision x-ray diffraction measurement for the 4H and 6H polytypes of SiC systems along with the first-principles calculations based on the local density approximation (LDA) by Bauer et al. have clearly shown that the internal degrees of freedom such as atomic relaxations are key ingredients for the structural and energetic properties in these systems.Reference Bauer, Kräußlich, Dressler, Kuschnerus, Wolf, Goetz and Käckell19 In the present work, we have, therefore, performed careful optimizations for the unit cells and internal coordinates in the systems considered. The criteria for optimizing the internal coordinates and unit cell structures are chosen so that the total energy change, maximum interatomic force, maximum stress acting on the unit cell, and maximum atomic displacement should be less than 5.0 × 10−6 eV/atom, 0.01 eV/Å, 0.02 GPa, and 5.0 × 10−4 Å, respectively. Since all polytypes considered are nonmetallic, the convergence of lattice parameters for the calculated systems is confirmed to be sufficient with the grids and criteria adopted. From the convergence tests of total energies and comparison with other DFT theoretical calculations, our choice of the k-space grids and E cut has been assumed to produce much smaller absolute error than the polytype energy differences for all the systems studied except SiC.
As the most severe case for estimating energetic differences among polytypes in the present work, Table II exemplifies the total energy difference for the 2H, 4H, and 6H(2) from 3C stacking configuration in the SiC system along with the recent DFT calculation results. All calculations including ours show that the 3C structure is significantly higher in energy than 4H by about 1 meV/atom, and the 2H is almost 3 meV/atom higher than 3C. This indicates that our calculations are sufficiently converged for the purposes of the present work. While claims for the ground state structure and ordering in energy for the SiC system have not been consistent among the DFT calculations,Reference Bernstein, Gotsis, Papaconstantopoulos and Mehl2, Reference Kobayashi and Komatsu4–Reference Limpijumnong and Lambrecht6, Reference Käckell, Wenzien and Bechstedt8, Reference Rutter and Heine9 as will be shown in this paper, we assume that the structures with moderate hexagonality such as 4H and/or 6H(2) are lowest in energy, 2H with perfect hexagonal stacking is probably highest, and 3C with zero hexagonality is in the middle, in the case of the SiC system.
Total energy difference for the 2H, 4H, and 6H(2) from 3C stacking configuration of SiC system in the unit of meV/atom.
C. ANNNI model
Since the polytypes considered here differ only in the stacking sequence along the <111> direction in the cubic 3C and/or the <0001> one in hexagonal (or rhombohedral) polytypes, the resulting energetic differences can be reasonably described by a one-dimensional Ising-type model, the ANNNI model.Reference Raffy, Furthmüller and Bechstedt3, Reference Limpijumnong and Lambrecht6, Reference Rutter and Heine9, Reference Lambrecht, Limpijumnong, Rashkeev and Segall20–Reference Heine, Cheng and Needs23 In this model, each bilayer, i, is characterized by two possible pseudospin states of σi = ±1, which have the same meaning as the Hägg notation corresponding to the orientation of triple bonds in a sp 3 tetrahedron.Reference Sebastian and Krishna14 The total energy per pair atoms for a system can be approximated as
$$E_{{\rm{total}}} = E_0 - {1 \over L}\mathop\sum\limits_{i = 1}^L {\mathop\sum\limits_{n = 1}^\infty {J_n } } {\rm{\sigma }}_i {\rm{\sigma }}_{i + n} ,$$where the layer index i is summed over the total number of layers L, E 0 is the self-energy of a layer, and Jn is the interaction strength for two layers with separation n. We restrict the n sum up to the third neighbors (n ≤ 3) assuming that the long-range interactions are small. In the present work, other effective Hamiltonian such as four spin interaction termsReference Limpijumnong and Lambrecht6, Reference Cheng, Needs and Heine21 have also been neglected. For all elemental systems considered, the interaction parameters of E 0, J 1, J 2, and J 3 are deduced by fitting the total energies of 2H, 3C, 4H, and 6H(2) calculated using the GGA. These interaction parameters can be obtained using the following equation:
$$\left( {\matrix{ {E_0 } \cr {J_1 } \cr {J_2 } \cr {J_3 } \cr } } \right) = {1 \over 4}\left( {\matrix{ 1 & 1 & 2 & 0 \cr 2 & { - 1} & 2 & { - 3} \cr { - 1} & { - 1} & 2 & 0 \cr 0 & { - 1} & { - 2} & 3 \cr } } \right)\left( {\matrix{ {E_{2H} } \cr {E_{3C} } \cr {E_{4H} } \cr {E_{6H(2)} } \cr } } \right),$$where E 2H, E 3C, E 4H, and E 6H(2) are the total energies for 2H, 3C, 4H, and 6H(2) configurations, respectively.
The ANNNI parameters extracted from the present GGA calculations for 2H, 3C, 4H, and 6H(2) polytypes are listed in Table III. As shown in Table III, the J 1 parameters are positive for carbon, silicon, SiC, and BN, while those for AlN, GaN, and InN are negative. The absolute values of |J 2| are smaller than |J 1| and those of |J 3| are much more for all elemental systems, so that the assumption that the long-range interactions in the polytypes considered are small is valid. In addition, the ANNNI parameters for carbon, silicon, SiC, and GaN extracted from the present GGA analyses are in good agreement with the values reported in Ref. Reference Raffy, Furthmüller and Bechstedt3 (for C and Si), Ref. Reference Limpijumnong and Lambrecht6 (for SiC), and Ref. Reference Wright13(for GaN), respectively. In addition, the ANNNI parameters for BN and AlN deduced from the present GGA calculations agree also well with those constructed from the LDA results by Kobayashi and Komatsu.Reference Kobayashi and Komatsu4
Parameters of the ANNNI model extracted from the present GGA calculations for 2H, 3C, 4H, and 6H(2) polytypes (in meV per pair). Other parameters reported in literature are also listed for comparison.
III. RESULTS AND DISCUSSION
A. Structural properties of polytypes
In the present GGA calculations, we have found some strong linear correlations between the structural properties (P) and the hexagonality (σH) in the polytypes considered. We have expressed these linear correlations found for the physical values of P against the hexagonality by the following linear equation:
$$P = k_P {\rm{\sigma }}_{\rm{H}} + P_0 .$$The coefficient parameters of kP and P 0 for some physical properties deduced from the present work are listed in Table IV.
Linear correlations for the structural properties against the hexagonality (σH) found in the present GGA calculations, where the associated properties (P) can be approximated by the equation of P = kPσH + P 0. The relations for the lattice parameters of a and c/n, the ratio of c/(na), and the atomic volume of V A are shown. We cannot find the linear correlation between V A and σH for Si and SiC systems (see Fig. 5).
Figure 2 shows the optimized lattice constants of a and c/n as a function of hexagonality for carbon, silicon, SiC, BN, AlN, GaN, and InN within the present GGA calculations. We have found that there is a strong linear correlation between the lattice constants (a and c/n) and the hexagonality (σH) in the polytypes for any elemental systems considered in the present work. The lattice constant parameter of a decreases with increasing hexagonality for carbon, silicon, SiC, and BN. The opposite trend is true for AlN, GaN, and InN, in which the 2H stacking structure (wurtzite) with perfect hexagonality is the ground state as described in the following section. The distortion in the close-packed plane due to the hexagonality change is much smaller for SiC and GaN than for other systems.
Optimized lattice constants of a (left panels) and c/n (right ones) as a function of hexagonality for carbon, silicon, SiC, and III–V nitrides within the present GGA calculations. For comparison with other structures, the structure of 3C is described by the notation with a rhombohedral supercell. The lattice parameter c is normalized by n, the number of bilayers per unit cell. The lines are deduced by the linear least squares fitting.
In contrast, the bilayer distance of c/n increases with rising hexagonality for carbon, silicon, SiC, and BN, while the opposite trend has again been confirmed in AlN, GaN, and InN systems. The values of lattice parameters of a and c/n for carbon and BN are nearly equal to each other and also for SiC and AlN. As shown in Table IV, the variation of lattice parameter c/n to hexagonality is again smaller in SiC and GaN systems than in other ones. These results indicate that the lattice parameters in SiC and GaN systems are more insensitive to the change of hexagonality compared with those in other systems. In addition, for each elemental system, the variation along the c-axis distortion against hexagonality is confirmed to be larger than that in close-packed plane by a factor from two to three.
In Fig. 3, we show the optimized ratio c/(na) as a function of hexagonality for carbon, silicon, SiC, and III–V nitrides within the present GGA calculations. We have also found that there is a strong linear correlation between the c/(na) value and the hexagonality in the polytypes for any elemental systems considered. The c/(na) increases with rising hexagonality for carbon, silicon, SiC, and BN, while the opposite trend has been detected in AlN, GaN, and InN systems. It is interesting to note that the slope constants of kc /(na) estimated for silicon and BN are nearly equal to each other.
The optimized ratio c/(na) as a function of hexagonality for carbon, silicon, SiC, and III–V nitrides within the present GGA calculations. The lines are deduced by the linear least squares fitting. The dotted horizontal line indicates the ideal value of 
$c{\rm{/}}(na) = \sqrt {2/3} $ in 3C structure.
In addition, we have also found an interesting relation between the polytype energetics and this slope constant kc /(na). Figure 4 shows the correlation between the GGA-calculated 2H − 3C energy difference, 
$\Delta E_{2H - 3C}^{{\rm{GGA}}} ( \equiv E_{2H} - E_{3C} )$, and the slope constants of kc /(na). From the location of AlN with large ionicity to that of carbon with large covalency, we can find a linear scaling between 
$\Delta E_{2H - 3C}^{{\rm{GGA}}} $ and kc /(na), which is approximated by a linear equation passing around the origin as shown in Fig. 4. From the systematic calculations for many compound systems using the LDA, Yeh et al. have also reported that the 
$\Delta E_{2H - 3C}^{{\rm{LDA}}} $ shows a linear scaling against an atomistic orbital-radii coordinate, 
$\overline R $, that depends only on the properties of the free atoms constructing the associated compounds, where the 
$\overline R $ is an orbital-dependent quantity and can be calculated from atomic pseudopotentials.Reference Yeh, Lu, Froyen and Zunger12 As mentioned in Sec. II. B., many first-principles calculations including ours show that the 2H structure is significantly higher in energy than 3C by about 3 meV/atom for SiC system. Since the slope constant kc /(na) for SiC is very small but has consistently positive value as shown in Figs. 3 and 4, this energetic ordering between 2H and 3C in SiC based on the first-principles calculations reported so far should be appropriate considering the correlation found in the present work.
The correlation between the GGA-calculated 2H − 3C energy difference and the slope constants of kc /(na) listed in Table IV. The lines are deduced by the linear least squares fitting.
The optimized atomic volumes of V A as a function of hexagonality for carbon, silicon, SiC, and III–V nitrides within the present GGA calculations are depicted in Fig. 5. It has been found that there is a relatively strong linear correlation between the atomic volume and the hexagonality for carbon, BN, AlN, GaN, and InN systems. The atomic volume of V A increases with rising hexagonality for carbon and BN and decreases for AlN, GaN, and InN. In silicon and SiC, however, the V A is assumed to have little correlation with the hexagonality. As will be described in the next section, the ground state is the 3C stacking for carbon and BN and the 2H for AlN, GaN, and InN. In addition, the V A increases with rising hexagonality for carbon and BN and decreases for AlN, GaN, and InN. The thermodynamic stability of the ground state under high pressure is, therefore, qualitatively predicted to be robust against the hexagonality change for these elemental systems.
Optimized atomic volume of V A as a function of hexagonality for (a) carbon, (b) silicon, (c) SiC, (d) BN, (e) AlN, (f) GaN, and (g) InN within the present GGA calculations. The lines are deduced by the linear least squares fitting.
B. Total energetics in polytypes
We have tabulated the total energies (meV/pair) of the polytypes calculated for carbon, silicon, SiC, and III–V nitrides using the GGA in Table V. In this table, the GGA total energies are measured from the equilibrium energy of 3C configuration for each elemental system. From the present GGA calculations, the stacking structures with the most stable in energy have been confirmed to be the 3C for carbon, silicon, and BN, the 4H for SiC, and the 2H for AlN, GaN, and InN, respectively. The lattice energy increments for the 2H with perfect hexagonality from the 3C structure have been calculated as 24.9, 11.0, 2.5, 17.4, −22.4, −5.2, and −8.7 meV/atom for carbon, silicon, SiC, BN, AlN, GaN, and InN, respectively. These values can directly be compared with other theoretical values of 25.3 (Ref. Reference Raffy, Furthmüller and Bechstedt3), 10.7 (Ref. Reference Raffy, Furthmüller and Bechstedt3), 2.3 (Ref. Reference Bernstein, Gotsis, Papaconstantopoulos and Mehl2), 18.8 (Ref. Reference Kobayashi and Komatsu4), −19.8 (Ref. Reference Kobayashi and Komatsu4), −5.5 (Ref. Reference Wright13), and −7.5 (Ref. Reference Zoroddu, Bernardini, Ruggerone and Fiorentini24) meV/atom in the literature.
Relative energies (meV/pair) of polytypes for carbon, silicon, SiC, and III–V nitrides. The symbols for nonequivalent polytypes are defined in Table II. The lattice energy increments from 3C lattice calculated with the present GGA are listed. The values in parentheses are defined as 
$\Delta E_{{\rm{ANNNI}}}^{{\rm{GGA}}} \equiv E_{{\rm{ANNNI}}} - E_{{\rm{GGA}}} $, where E ANNNI and E GGA are the total energies calculated with the present ANNNI model and the GGA, respectively.
In Table V, we have also tabulated the accuracy of the present ANNNI model for estimating total energies of nonequivalent polytypes in each elemental system, which is defined by 
$\Delta E_{{\rm{ANNNI}}}^{{\rm{GGA}}} \equiv E_{{\rm{ANNNI}}} - E_{{\rm{GGA}}} $. The maximum errors of total energies calculated using the ANNNI model for the polytypes with up to the L= 8 have been 0.17, 0.14, 0.09, 0.05, 0.65, 0.06, and 0.11 meV/atom for carbon, silicon, SiC, BN, AlN, GaN, and InN systems, respectively. The ANNNI models with the present parameters are, therefore, confirmed to be an excellent ones for estimating the total energetics in both homopolar and heteropolar polytypes with various covalent and ionic characteristics. Since the ANNNI model projected from the total energetics calculated with the appropriate optimization technique using the GGA is confirmed to be accurate one, we have performed the evaluations on the energetics for tremendous amounts of polytypes, whose calculations are impossible within the first-principles schemes.
Figure 6 depicts the ANNNI total energies as a function of hexagonality (left panels) and periodicity (right panels) for 6,209,704 kinds of nonequivalent polytypes in carbon, silicon and InN, whose configurations have been deduced using the efficient polytype generation algorithm reported by Estevez-Rams and Martinez-Mojicar.Reference Estevez-Rams and Martinez-Mojicar1 Energy is measured from the 3C configuration in these figures. We have also plotted the results of the present GGA calculations by open circles. The dotted horizontal lines are the energy levels for 2H, 3C, and 4H stacking sequences, respectively. In the left panels, as a guide for the eye, we show the triangular area where the energies for 2H, 3C, and 4H are connected. The solid lines in the right panels show the upper and lower limits of total energies of polytypes in each periodicity. The same figure as Fig. 6 for the SiC polytypes is shown in Fig. 7.
The ANNNI total energies of the nonequivalent polytypes for carbon [(a) and (b)], silicon [(c) and (d)], and InN [(e) and (f)] relative to that of the 3C structure are plotted as a function of hexagonality (left panels) and periodicity (right panels). The open circles denote the values using the GGA calculations. The dotted horizontal lines are the energy levels for 2H, 3C, and 4H stacking sequences, respectively. In the left panels, as a guide for the eye, we show the triangular area where the energies for 2H, 3C, and 4H are connected. The solid lines in the right panels show the upper and lower limits of total energies of polytypes in each periodicity.
Same figure as Fig. 6 for SiC polytypes.
Focusing on the left panels in Figs. 6 and 7, one of the most interesting universal trends for any elemental systems is that almost all the points of total energies as a function of hexagonality for the polytypes are located in the inner region of a triangle with vertices of the 3C (σH = 0), 4H (σH = 1/2), and 2H (σH = 1) points, while a few exceptional points can be seen especially in the silicon and SiC systems. This trend indicates that the energetic ordering among the 3C, 4H, and 2H configurations is a key ingredient for grasping the general view of energetics for the polytypes of semiconductor compounds considered.
Within the ANNNI model, the geometrical equations of the triangle edges for 3C − 2H, 3C − 4H, and 4H − 2H can be written by:
$$E_{3C - 2H} ({\rm{\sigma }}_{\rm{H}} ) = 2(J_1 + J_3 ){\rm{\sigma }}_{\rm{H}} ,$$
$$E_{3C - 4H} ({\rm{\sigma }}_{\rm{H}} ) = 2(J_1 + 2J_2 + J_3 ){\rm{\sigma }}_{\rm{H}} \,\,\rm and,$$
$$E_{4H - 2H} ({\rm{\sigma }}_{\rm{H}} ) = \,2(J_1 - 2J_2 + J_3 )\left( {{\rm{\sigma }}_{\rm{H}} - {1 \over 2}} \right) \cr+ (J_1 + 2J_2 + J_3 ),$$ respectively. In addition, we define the energy ratio of 
$R_{4H} \equiv \Delta E_{4H - 3C} /\Delta E_{2H - 3C} $, where the energy increments from the 3C are defined as 
$\Delta E_{2H - 3C} \equiv E_{{\rm{total}}}^{2H} - E_{{\rm{total}}}^{3C} $ and 
$\Delta E_{4H - 3C} \equiv E_{{\rm{total}}}^{4H} - E_{{\rm{total}}}^{3C} $, respectively. It should be noticed that the triangle of 3C–4H–2H becomes flat and can be seen as a straight line, when the R 4H is approaching the value of 1/2. Substituting the present ANNNI parameters for Eqs. (4), (5), and (6), we can get an energetics map of polytypes for the compounds as shown in Fig. 8.
Energetics map for polytypes deduced using the present ANNNI analyses.
The values of R 4H are estimated as 0.37, 0.31, and 0.41 for carbon, silicon, and BN, respectively, in which the ANNNI values of J 1 are positive and the ground state is the 3C stacking. Since the values of R 4H for these elemental systems are much less than 1/2, the 3C–4H–2H triangles in the energy maps are geometrically convex downward for these systems. One of the remarkable features in the systems with relatively strong covalent characteristic is that a clear energy folding point appears at σH = 1/2 in the energetics map as shown in Fig. 8. It should be also noticed here that the strong covalent characteristic for BN with sp3-bonding and its similarity to carbon diamond have theoretically been addressed using the first-principles analyses on their structural, lattice dynamical, and dielectric properties by Karch and Bechstedt.Reference Karch and Bechstedt25 In contrast, the values of R 4H are 0.48, 0.54, and 0.49 for AlN, GaN, and InN with strong ionicity, respectively, in which the J 1 values are negative and the ground state is the 2H stacking. Since the R 4H values for these systems are nearly equal to 1/2, the 3C–4H–2H triangles in the energy map are geometrically flat and can be seen as straight lines. In addition, it is interesting to note that the energy distribution region for the polytypes in Fig. 8 systematically rotate clockwise around the origin (σH = 0) with changing from the strong covalency systems to strong ionicity ones. The energies of the polytypes with hexagonality from σH= 0 to 1/2 are almost degenerate for the SiC system with intermediate characteristic between covalency and ionicity.
Within the full-potential linear muffin tin orbital method using the GGA, Limpijumnong and LambrechtReference Limpijumnong and Lambrecht6 have reported that the 15R polytype with the stacking of (UUDDD)3 (σH = 2/5) is the lowest in energy in the SiC system. In general, there are 203 kinds of nonequivalent polytypes with the stacking length of L = 15, in which the hexagonality is in the region from σH = 2/15 to 4/5. Using the present ANNNI parameters for SiC system, the most energetically favorable configuration among the polytypes with L = 15 is confirmed to be the stacking of (UUDDD)3, and the energy difference from the 3C is calculated by −1.48 meV/pair. The ordering in energy has, therefore, been predicted by 4H < 15R < 6H(2) with the associated energies of −1.56, −1.48, −1.43 meV/pair, respectively. While the energy ordering of these three polytypes in the SiC system seems to be controversial depending on the calculation schemes,Reference Bernstein, Gotsis, Papaconstantopoulos and Mehl2, Reference Limpijumnong and Lambrecht6 these energy differences are almost degenerate and within the estimated accuracy.
The right panels in Figs. 6 and 7 show the ANNNI total energies as a function of periodicity for 6,209,704 kinds of nonequivalent polytypes in each elemental system. In the case of periodic length of L = 30, there exist 2,990,968 kinds of nonequivalent polytypes with the hexagonality from σH = 1/15 to 14/15. Since there are few kinds of nonequivalent polytypes when the periodicity is short, the energies for the polytypes associated become discrete in the energetic maps of E total versus periodicity. It is, however, possible to consider stacking sequences with lower and higher hexagonality as the periodic length of L increases because the number of nonequivalent polytypes gives the exponential explosion with L. Therefore, the energies of polytypes become continuous with increasing periodicity in the energetic maps. The solid lines in the right panels of Figs. 6 and 7 show the upper and lower limits of total energies of polytypes in each periodicity. Depending on the ground state stacking and the highest energy one for each elemental system, these upper and lower limit lines approach asymptotically these associated energy levels with increasing periodicity.
In the ANNNI model, the total energy difference from 3C for each configuration can be generally expressed as 
$E_{{\rm{total}}} = \sum\nolimits_{i = 1}^3 {a_i J_i } $. Figure 9(a) shows the hexagonality as a function of periodicity for the SiC polytypes with the lowest energy in each periodicity. The hexagonality for the SiC polytypes with energetic stability approaches asymptotically, accompanying a step-like behavior, the value of 4H (σH = 1/2) with increasing periodicity. We have also depicted the behaviors of the ANNNI coefficients of the polytypes with the lowest energy in each periodicity in Fig. 9(b). All these coefficients are found to asymptotically approach those of 4H (
$\Delta E_{4H}^{3C} = J_1 + 2J_2 + J_3 $) with increasing hexagonality, also accompanying step-like behaviors.
(a) The hexagonality as a function of periodicity for the SiC polytypes with the lowest energy in each periodicity. (b) The coefficient parameters (ai) of ANNNI model for the SiC polytypes associated. The dotted horizontal lines other than zero indicate the associated values of 4H. All these values asymptotically approach those of 4H (
$\Delta E_{4H}^{3C} = J_1 + 2 J_2 + J_3 $) with increasing periodicity.
We have also found that most polytypes with the lowest energy in each periodicity have been experimentally observed and determined to date.Reference Jepps and Page10 This fact indicates that the thermodynamic stability is still one of the important factors in selecting the polytypic structures with a given stacking length, though the energy differences of polytypes are generally quite small in SiC system. Focusing on the polytype structures with long periodicity, it is interesting to note that the stacking sequences with energetic stability are basically constructed of the combinations of 4H (UUDD: σH =1/2), 6H(2) (UUUDDD: σH =1/3), and 15R [(UUDDD)3: σH=2/5], which are the polytypes most often observed in the experimental crystal growth and do not have energetic degeneracy among the polytypes with the same associated periodicity.
IV. CONCLUDING REMARKS
We have investigated correlations between energetics and hexagonal stacking configurations for carbon, silicon, SiC, BN, AlN, GaN, and InN polytypes with sp 3-bonded networks, using the comparative analyses based on the DFT and the ANNNI model. The atomistic geometry, energetics, and electronic structure for these compounds with up to the periodic stacking length of L = 8 have been carefully calculated within the GGA. Using the ANNNI model extracted from the GGA calculations, we have also studied the energetics for more than 6 million kinds of nonequivalent stacking polytypes with up to L= 30, whose configurations have been deduced by the efficient polytype generation algorithm, and illustrated some trends of structural and energetic properties for these compounds. Since we have performed the evaluations on the energetics for tremendous amounts of polytypes, for which calculations are impossible within the first-principles schemes, more comprehensive trends on the correlations have been stated in this work. Therefore, the energetic concepts presented here might also be helpful information for controlling crystal growth and/or band gap engineering based on the polytypism in wide-gap semiconductors.
ACKNOWLEDGMENTS
The authors sincerely thank the Director General of Corporate Research and Development Laboratories, T. Miyake for giving us an opportunity to carry out this research. The authors’ gratitude also goes to M. Igarashi, Y. Ohashi, A. Yauchi, and Y. Sumitomo (Sumitomo Metal Industries, Ltd.) for their encouragement in this work. The authors would like to also thank Prof. K. Suzuki for fruitful discussions on the SiC polytype systems.
