A comparative DFT study of electronic properties of 2H-, 4H- and 6H-SiC(0001) and SiC() clean surfaces: significance of the surface Stark effect

The electric field, uniform within a slab, emerging due to Fermi level pinning at both sides of the slab, is analyzed using DFT simulations of SiC surface slabs of different thicknesses. It is shown that for thicker slabs the field is nonuniform and this fact is related to the surface state charge. Using the electron density and potential profiles, it is proved that for high-precision simulations it is necessary to take into account a sufficient number of SiC layers. We show that the use of 12 diatomic layers leads to satisfactory results. It is also demonstrated that the change of the opposite side slab termination, both by different types of atoms or by their location, can be used to adjust the electric field within the slab, creating a tool for simulation of surface properties, depending on the doping in the bulk of the semiconductor. Using these simulations, it was found that, depending on the electric field, the energy of the surface states changes in a different way than the energy of the bulk states. This criterion can be used to distinguish Shockley and Tamm surface states. The electronic properties, i.e. energy and type of surface states of the three clean surfaces: 2H-, 4H-, 6H-SiC(0001) and SiC(), are analyzed and compared using field-dependent DFT simulations.


Introduction
Silicon carbide can exist in a very large number of polytypes, of which hundreds have already been identified [1]. They are created due to different stacking sequences, still preserving the tetrahedral configuration of silicon and carbon atoms. These polytypes belong to wide bandgap semiconductors, having its gaps, ranging from 2.403 eV (3C), 3.101 eV (6H), 3.285 eV (4H) to 3.300 eV (2H) [2]. Also, SiC stacking faults, being planar large area defects, could be obtained in a controlled fashion. These structures, chemically identical with the material of the locally adjustable bandgap, are new, potentially important, quantum structures [3]. In addition, a number of quantum structure-based devices, using SiC as the building material, were developed [3]. These structures are based on the material having its properties, extremely attractive for many applications in electronics (see e.g. [4]), that is, a wide-bandgap semiconductor having superior electronic properties, displaying high values of several important parameters like breakdown electric field (2.2 × 10 6 V cm −1 ), forward current density (up to 1 kA cm −2 ), saturated electron drift velocity (2 × 10 7 cm s −1 for 4H-SiC and 2.5 × 10 7 cm s −1 for 3C-SiC) and room temperature mobility (370 and 720 cm 2 V −1 s −1 for 6H-SiC and 4H-SiC, respectively) and, finally, a high blocking voltage observed in typical MOSFET devices [3]. It may be used in the high-temperature environment because it is thermally resilient and has high thermal conductivity [5]. Being mechanically extremely hard and chemically inert [6], SiC can be used in a number of important applications, mostly in high-power, high-temperature and high-frequency electronics [7]. Therefore SiC and the structure based on this compound attracted tremendous interest, resulting in a huge number of publications.
The most successful method of growing large-size high-quality crystals is based on physical vapor transport (PVT), known as a modified Lely method first developed by Tairov and Tsvetkov [8,9]. The method, although technically difficult due to the very high temperatures required, was developed to the stage in which large size SiC crystals were grown and distributed commercially. Yet, despite considerable success, the method is still plagued by difficulties, most notably by the creation of micro-and nano-pipes [10], which deteriorate or even block the functioning of SiC-based devices. The density of these defects was substantially reduced, which allowed the use of Lely-grown SiC substrates in the technology of electronic devices. The pipe problem still remains essential; thus basic research and development of the technology are needed to alleviate it completely. 3 It is well known that the modified Lely method relies on the chemical decomposition of the source SiC material into Si, Si 2 C and SiC 2 volatile species. The transport process is, therefore, prone to chemical instabilities, which affect the chemical state of the SiC surface, leading to structural and chemical transformations that may contribute towards the generation and preservation of pipes during growth. It is therefore understandable that the atomistic structures of various SiC surfaces, belonging to different polytypes, were intensively studied previously.
In fact, clean hexagonal SiC(0001) surfaces are difficult to achieve, as they are typically covered by Si or C atoms, or other adsorbed species like oxygen, hydrogen, etc. However, density functional theory (DFT) studies of these surfaces have considerable value as a reference point for the investigation of adsorption phenomena and their influence on the electronic and structural properties of SiC(0001) and SiC(0001) surfaces. In an early study of these systems, Sabisch et al [11] showed that both Si-and C-terminated surfaces undergo strong relaxation of the topmost layer by surface reconstruction. The top Si and C layers move 0.15 and 0.25 Å towards the SiC interior, a considerable fraction (about 23 and 40%) of the ideal distance 0.63 Å. The second layer, i.e. C and Si layers, moves in the opposite direction by 0.04 and 0.11 Å, which is a much smaller fraction of the ideal distance 1.88 Å. Third-layer relaxation is small, moving the atoms by 0.02 and 0.03 Å, respectively. Again the direction is reversed; the third layer atoms are moved towards the SiC interior.
The above-described relaxation type indicates that the ionic contribution to this effect is essential and the subsurface field is an important factor. However, Sabisch et al [11] calculated the surface band structure of both surfaces without reference to this field. Using DFT LDA, they obtained an indirect bandgap of bulk 6H-SiC equal to 1.97 eV, which reflects the well-known underestimate of bandgap energy in this approach. The band structures of both surfaces were compared with the bulk data. They identified the flat band in the gap, originating from dangling bonds of surface Si and C atoms. The Si-related surface state is about 1.5 eV above the C-related states, which is attributed to the difference in ionic character of both surfaces. A charge transfer towards the C-surface was predicted and proved by using charge density distribution plots. It was also shown that both the states are half-filled, i.e. the surface is metallic. It is interesting to note that the charge density distribution of the dangling bonds is totally different: the C-state is quasi-symmetric with some dominance of the internal part, whereas in the Si-related state, the charge is located exclusively in the region above the Si atom. The dispersion of both the states is generally small, which is due to a small overlap of the neighboring atom wavefunctions. As expected, larger dispersion is observed for Si-related states [11].
The adsorbate-covered SiC(0001) and SiC(0001) surfaces display a wide large variety of different reconstructions such as ( 3) and (9 × 9) [12]. A plethora of different symmetries was studied and several different models were proposed [13]- [15]. The subject remains under intense dispute. New impetus to the study of the nonstoichiometric SiC surface was given by the discovery of graphene on both Si-and C-terminated SiC surfaces [16]. Since the subject is vast and not in the mainstream of the present work, it will not be discussed here.
It is interesting to note that the surface structure, defined by the positions of the adsorbed atoms, is very sensitive to external factors. This effect was investigated by Righi et al [17]. They have found that, at the surface, the two-dimensional effect favors the formation of 3C-SiC. In contrast to adsorbate covered surfaces, the clean polar surfaces of other hexagonal polytypes were expected to have essentially similar properties, as the differences in atom arrangement are observed deep in the bulk. Still this effect could be investigated along with the role of the subsurface electric field.
The present paper is devoted to the study of the methodology of thick slab simulations, the role of electric field (figure 1) and the influence of stacking sequence on the electronic and structural properties of polar SiC(0001) and SiC(0001) surfaces. The paper is organized as follows. In section 2, we give a brief overview of employed computational methods. Section 3 describes the construction of surface models. We analyze there the influence of finite number of layers on the reliability of calculations. We also study the effect of surface termination by hydrogen atoms and show that they can be used to adjust electric field inside a slab. Figure 1 explains the subsurface field and its connection with the surface acceptor and donor states, together with its modeling in the slab simulations. As described in detail in section 4, there is a significant difference between Shockley and Tamm surface states. As shown in figure 1, Shockley states can be transformed into bulk states under certain electric field conditions. On the other hand, Tamm states, irrespective of the electric field, are placed in the bandgap. In section 4, we present the obtained structural and electronic properties of examined surfaces. Section 5 summarizes our work.

Method of calculation
We have carried out self-consistent total-energy calculations based on DFT. The VASP code developed at the Institut für Materialphysik of Universität Wien was employed [18]- [20]. The projector augmented wave (PAW) approach [21] was used in its variant available in VASP package [22]. For the exchange-correlation functional, generalized gradient approximation (GGA) according to Perdew, Burke and Ernzerhof (PBE) [23] was selected. The plane wave cutoff energy was set as 500 eV. In bulk calculations a 11 × 11 × 11 Monkhorst-Pack k-point mesh was used [24]. Lattice parameters obtained for bulk SiC polytypes (2H, 4H, 6H), summarized in table 1, are in good agreement with the experiment [4, 6,25]. The band structures of the hexagonal polytypes are presented in figure 2. Generally, these band structures are in good agreement with the previously obtained DFT results [11]. As it follows from these data, the bandgaps are much lower than the experimental findings. This is a well-known systematic error of DFT LDA/GGA approximations.
The surface slabs representing the 2H-, 4H-and 6H-SiC(0001) surfaces are built using different numbers of SiC double layers (DL). We used standard notation, i.e. the SiC(0001) surface is the one terminated by Si atoms. Slabs having 6-18 atomic layers were used, i.e. slabs consisting of up to 36 atoms. Due to the fact that SiC(0001) and SiC(0001) surfaces can be prepared without surface reconstruction [26], (1 × 1) supercells were employed. The bottommost C atoms had their broken bonds either terminated by hydrogen atoms or left uncovered. Particular attention was paid to k-point sampling and convergence. In surface calculations, the Monkhorst-Pack scheme was applied, using a 9 × 9 × 1 k-point mesh for the bare and covered surfaces. The selected parameters were sufficient to obtain relative energies converged to a few meV. To prevent artificial electrostatic field in the vacuum region between repeated unit cells, a dipole correction was used [27]. The three topmost SiC layers were relaxed using a conjugate-gradient (CG) algorithm in section 4.1. Due to the fact that significant relaxation is limited to the first two layers, in other cases only two layers were relaxed.

Preparation of surface models
Recently developed formalism allowing us to simulate the influence of the electric field and the doping of a bulk semiconductor allows us to study the interplay between the ionicity of the surface, the existence of electric field and the position of the Fermi level [28,29]. Meyer and Marx [30] claimed that the ZnO surface can be described using the purely ionic model, where each ZnO DL exhibits a dipole moment perpendicular to the surface. In such a model we obtain spontaneous polarization, which is independent of the thickness of the slab. The top part of figure 3 shows that the potential difference at both the sides is independent of slab thickness, which means that the field inside the slab is smaller for larger thickness and there is divergence for large thicknesses. Thus, the purely ionic model is unstable and has to be corrected by the rearrangement of the charges. Moreover, the electron density at the ends of the slab decreases ( ρ < 0) for thicker surface.
The energy of the quantum states of the H termination atom can be controlled by their distance to the nearest Si/C atoms. In accordance with the stability condition of fermion systems, the Fermi energy should be positioned halfway between the empty and the occupied states. The energy of the surface states may be changed by the electric field, which is the well-known Stark effect. Meyer and Marx [30] recognized that fact, but they erroneously assumed that these are always conduction and valence band states. In fact, more frequently these are surface states. Also following Tasker [31], they deduced that the field in the slab is constant, which leads to divergence for large slab thicknesses. This was repeated later by Noguera [32]. This prediction also turned out to be unstable. In fact, the condition of the stability of the fermion system leads to smaller field in the thicker slab as it was shown in figure 3. From these calculations, we can learn about the influence of the subsurface electric field on the equilibrium and dynamic properties of surfaces. This in turn depends on the energy and concentration of principal point defects in the bulk. At present, this is a completely unexplored area.
Electric potential energy profiles are useful quantities to examine when analyzing the influence of the finite size of the slab. As it follows from figure 3, these profiles inside the slab consist of three parts: the linear part in the middle and nonlinear parts at both ends. The linear part is in accordance with the scenario predicted first by Meyer and Marx [30], where they claimed that the electric field is created due to charge transfer within the slab, leading to equality of the valence band maximum (VBM) on the negative charge side and the conduction band minimum (CBM) at the Fermi energy level. Using ionic arguments in the investigation of ZnO surfaces, they claimed that the charge separation leads to a uniform field within the slab, as is observed in the central part of potential diagrams in figure 3. This argument was based on an earlier argument by Tasker [31] that charge separation creates a uniform field within the slab, which leads to divergence for thicker slabs. Later on, this argument was reformulated by Noguera for the case of oxide crystals [32]. As it was discussed in [29], this prediction is unstable; the field is not kept constant but the difference in electric potential between both the sides of the slab is preserved. This is also demonstrated by the potential plots in figure 3, where the difference of potentials is identical for all three slabs. Thus for thicker slabs, the electric field is smaller.
The identification of the physical factor creating the field within the slab by Meyer and Marx [30] is confirmed by the band diagram (BD) and density of states (DOS), plotted in figure 4. The BD and DOS are calculated for the same slab models as the electric field and charge diagrams presented in figure 3. As it is shown for all the three cases, the Fermi level is kept at the energy of the quantum states associated with both the sides of the slab. Our BD and DOS prove that the Fermi level is attached to the states, independently of the slab thickness. Thus, the diagram proves that the Fermi level is pinned at the bare Si face and C face. The data presented in figure 3 show that the central part of the potential diagram is flat, i.e. this part of the slab is electrically neutral. On the contrary, the edge parts are not, which proves that there is an uncompensated charge at both sides of the slab, creating the uniform field in between. This charge was plotted in figure 3 for the systematic variation of the slab thickness and it was also shown that the average charge density decreases exponentially ( figure 5). Here, average charge is defined as the average charge at the center of the slab, i.e. an interval with a width of 2.5 Å from the center of the slab. In fact, the interaction between both the sides of the slab has two principal contributions: electrostatic and quantum overlap. The electrostatic contribution is in fact useful for the simulation of subsurface fields. The quantum overlap should be avoided as it could cause spurious effects. Hence, both diagrams can be used for the determination of how the thick slab should be used for the simulation of SiC surfaces. The linear part in the potential profile and negligible overlap (figure 5) are the necessary conditions for accurate modeling. From these diagrams, it follows that the slab having 12 DL is sufficient to simulate SiC polar surfaces.
An electric field is an additional important component affecting the structure and electronic properties of semiconductor surfaces. The extent of this modification needs direct modeling  by numerical procedures. To the best of our knowledge, no such awareness or methodology is found in the literature. However, a suitable approach was developed lately [28,29]. It was demonstrated that a change of the opposite slab termination can be used for the simulation of the subsurface field in polar GaN(0001) surfaces. It was also shown that a change of the distance between the hydrogen termination atoms and the layer can affect the electric field, allowing us to change its magnitude and also its direction. The same approach was applied to SiC slabs. In figure 6, the results of such an approach are shown. It is demonstrated that the approach is universal: in all three polytypes the procedure allows us to obtain different magnitudes and directions of the electric field within the slab. In figure 7, we show points where we get potential cross sections. Potential profiles presented in figure 3 (2H-SiC) come from point A in the middle of the hexagon. 4H-and 6H-SiC polytypes have atoms at point A. In this case, point B is used. In order to achieve consistency with other polytypes, point B was used for each potential profile in figure 6.
In the classical theory of semiconductor surfaces, the band and surface state energies are identically affected by the subsurface field, i.e. the difference of the (bulk) band state and surface state energies is independent of the band bending, equal to e 0 V s . Since the Fermi level is, as a rule, located in the bandgap, Schottky barrier approximation, assuming complete ionization of the defects, can be used. The density of surface charge, ρ sur = q/S, where q is the surface charge for a single lattice site covering area S, is related to electric field E as follows ρ sur = εε 0 E.
(1) The Schottky barrier is characterized by the parabolic potential distribution that can be expressed as where E, the magnitude of the electric field at the surface, is given by    [40]. Points where we get potential profiles are marked. and z 0 denotes the charged layer width. The width z 0 can be obtained from the neutrality condition, z 0 = ρ sur /(N i e 0 ), giving In the classical approach, the surface band bending is obtained assuming that the energy of the surface states measured with respect to the band states is independent of the field intensity. Then, from DFT, the Fermi energy position with respect to the surface state energy is determined, from which the band bending is obtained, i.e. e 0 V s . The magnitude of band bending e 0 V s , derived from DFT calculations, is equal to the difference between the surface state energy (for the case of pinning) and the dominant point defect energy. The ionization energies of the principal impurities in SiC are: N (shallow donor): -E ∼ = 0.06 eV; and Al (shallow acceptor): E ∼ = 0.06 eV [37,38]. These defect ionization energies are small and slightly vary from type to type and therefore they can be treated as corrections to the VBM and CBM energies. From the relations (equations (1)-(4)), a complete description of the electric state of the surface is therefore obtained, i.e. E, ρ sur and z 0 are determined, provided that the density of the donor/acceptor defects and their energy levels are given. Thus, the point defect density can be set arbitrarily from which the remaining parameters could be obtained. The results obtained here paint a different picture: the energies of the band states and surface states change differently, leading to a relative shift of their energies, caused by the subsurface electric field. Generally, the change of the quantum state energies is denoted as the Stark effect. As the surface states are differently localized, the electric field changes their energies with respect to band states, a phenomenon called the surface states Stark effect (SSSE). Using the electric field from figure 3 for 18 DL, which is about 2 MV cm −1 (0.02 V Å −1 ), and SiC dielectric permittivity ε = 9.7 [39], the planar charge density ρ sur = 1.71 × 10 −2 C m −2 is obtained that translates into the charge for a single lattice site: q sur = 8.71 × 10 −3 e 0 . The charge transfer is therefore relatively small; still it affects the electric properties of the surface considerably. Assuming the benchmark density of the fully ionized point defects, N i = 10 19 cm −3 , for the simulated field intensity, the width of the charged layer is z 0 = 10.7 nm. Using equation (2), this translates to band bending e 0 V s ∼ = 1 eV. This result proves that a relatively small charge is able to induce band bending of about 30% of the bandgap. In the more practical approach, for a given electric field E, derived from the DFT simulation results, the surface charge density is obtained. From the difference of the energies of the surface slab and in the bulk, the magnitude of band bending e 0 V s could be obtained, neglecting or accounting for small acceptor/donor ionization energies. In more precise calculations, however, this has to be accounted for by subtraction of the energy difference between the Fermi energy levels at the surface, obtained from the above DFT calculations, and the energy of dominant point defects. Then the depletion width could be obtained from equation (3), accounting for the electric field intensity E. This provides a complete description of the electric conditions at the surface in the new model. The difference is that the energy bending has to be obtained for each field intensity independently as the Fermi energies change, from which the charge density and the width of the charged layer is determined. Thus, the point defect density for which the simulation is carried out is obtained from equation (3) and cannot be set independently. For a given doping level, the DFT calculations have to be adjusted. In practice, it is more convenient to simulate the surface properties and then derive the corresponding point defect density. That completes the determination of the electric state of the surface. In summary, the difference between the classical approach and the proposed new one is that in the first, the density of the point defect can be adjusted arbitrarily, and the remaining values could be obtained from equations (1)-(4). In the new approach, the density of charged defects has to be adjusted to the derived surface potential barrier V s , and the field intensity E, from which the corresponding defects density is obtained from equations (1), (3) and (4).

Surface properties
The simulation results include the structural and electronic properties of clean SiC(0001) and SiC(0001) surfaces of all three polytypes. Following the results from the previous section, in all the cases the 12 DL slab models were used to ensure good accuracy of the calculations. We also studied the dependence of the properties on different values of the subsurface electric field, which was adjusted according to our findings described in section 3.

Structural properties
For bare SiC (0001) and SiC (0001) surfaces, no surface reconstruction was obtained. The surface modification was limited to relaxation of several atomic layers (cf figure 8). The relaxation data is collected in table 2. From these data are follows that significant relaxation is limited to the first two layers. The most significant changes are observed in the [0001] and [0001] directions. No significant difference between the polytypes was observed. Also there  is no dependence on electric field, which indicates weakly ionic character of the polar SiC surfaces. The data are in good agreement with the paper by Sabisch et al [11].

Electronic properties
In contrast to structural properties, the electronic properties of bare SiC(0001) and SiC(0001) surfaces are severely affected by the subsurface electric field. This effect is shown in the collection of diagrams plotted in figure 9. In figure 9, the BD for the Si surface of the 2H polytype is presented. These results indicate that the surface states associated with the Si atoms are moved down with respect to VBM. The two surface states are those associated with Si-top, C-bottom and hydrogen termination atoms. The band associated with the H atoms has considerable dispersion arising from the overlap with C dangling bonds. These states are essentially empty. The band related to Si-top atoms is almost dispersionless. The energy of the states changes strongly with the electric field so that for these fields the states move across the entire bandgap. For 2H and 4H, two bands associated with the bottom atoms (i.e. hydrogen and carbon) and Si-top atoms are close to each other. Yet there are some differences between SiC polytypes. For 6H SiC(0001), the states are initially more separated than for 2H and 4H. When the distance between hydrogen and carbon is increased, the band related to the H-and C-bottom atoms moves to the Fermi level. The Fermi level moves in a valence band direction in all cases. Moreover, two additional bands associated with C-bottom atoms are shifted into the valence band when the H-termination atoms are sufficiently close to the surface. With the diversion of hydrogen atoms, the case is more and more similar to the bare surface, cf figure 4. In all these cases, it is observed that the Fermi level is pinned to the surface states related to dangling Si bonds of the SiC(0001) surface. Field intensities have been shown in table 3. For three polytypes (2H, 4H and 6H), the electric field changes from positive to negative values, cf figure 6. However, for 6H-SiC (distance of 1.5 Å), the field is much smaller than for 2H-and 4H-SiC.   A similar analysis was made for C-terminated SiC(0001). Band structures and DOS are presented in figure 10. As in the previous case, overlap of H-terminated and Si-bottom wavefunctions strongly affects electronic structure. Unlike the SiC(0001) surface, the electric field changes from negative to positive values. Bands associated with C-top atom states are initially located near the VBM but, under the influence of the electric field are moved across the entire bandgap. It is dispersionless. On the other hand, band associated with bottom atom states (H-and Si-bottom) resembles the shape of the valence bands. However, when the hydrogen atoms are moved, it bends to the Fermi level. Additionally, there is one more band (also associated with bottom atom states) shifted outside of the conduction band, i.e. with a change of  [33] and Tamm [34] states. It is often claimed that there is no real physical distinction between these two types. In fact, there is a deep difference between them. Shockley surface states are those associated with band states, which are localized at the surface by the potential well created by band bending (see figure 1). Thus, they disappear when there is no potential well. They are localized only in the direction perpendicular to the surface. Tamm states of surface atoms arise from dangling bonds at the surface. They are fully localized and they do not play a role in the electrical conductivity. In contrast to Shockley states, Tamm states do not disappear with the potential well. This difference provides a clear procedure for the determination of the type of surface states. Generally, it is difficult to precisely indicate which states are Shockley and Tamm states. There are several works in the recent literature describing ways in which these states can be distinguished [35]. Zak [35] showed a graphical description of two kinds of crystal terminations. For different termination types, we expect Shockley or Tamm states. Following this classification, surface states induced by surface perturbation should be called Tamm states. In figure 3, we can observe such a perturbation at the edge of each bare surface, i.e. shallower wells at the ends of the slabs. Because of this perturbation Tamm surface states can be induced.
Moreover, it was found that Tamm and Shockley surface states can coexist in a semiconducting superlattice [36]. To denote the type of states, the band structure and DOS have been analyzed for different methods of termination. The states that disappear in the reversed electric field are Shockley states. The states that move across the entire bandgap preserving their localization at the particular surface are Tamm states.
Following the changes in the band structure and DOS and comparing the results for SiC(0001) and SiC(0001) sides it is possible to identify the surface states associated with Si and C edge atoms. This method can be used as a criterion for the identification of Tamm and Shockley surface states.
Moreover, changing the position of H-termination atoms, we can move the Fermi level, which means that we can control the doping type of the SiC surface. In figures 9 and 10, doping type is changed from n-type and p-type, respectively, to almost undoped SiC surface.

Summary
An extensive density functional study of bare and H-terminated SiC surfaces is presented in this paper. Different SiC(0001) and SiC(0001) polytypes: 2H, 4H and 6H were considered and compared. It was determined that twelve Si-C layers is a sufficient number for high-precision calculations. Structural optimization was performed for 2H-, 4H-and 6H-SiC(0001) as well as SiC(0001) surfaces and no significant differences between them were found. It was shown that for all cases (H-terminated 2H-, 4H-and 6H-SiC) the electric field within the slab can be adjusted. Field intensity dependence on the H positions for SiC(0001) and SiC(0001) has been presented (tables 3 and 4). Figures 9 and 10 prove that the Fermi level is pinned to the surface states. The energy of surface-related states of both Si-and C-terminated surfaces is moved by the subsurface field with respect to VBM. In the case of the Si-terminated surface, the energy is shifted by about 2 eV, whereas for the C-terminated surface, this shift amounts to about 1 eV. This effect can be used as a simulation tool, e.g. for the identification and subsequent analysis of surface states as well as changing n-type and p-type doping levels.