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

Electric field, uniform within the slab, emerging due to Fermi level pinning at its both sides is analyzed using DFT simulations of the SiC surface slabs of different thickness. It is shown that for thicker slab 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 enough number of the Si-C layers. We show that using 12 diatomic layers leads to satisfactory results. It is also demonstrated that the change of the opposite side slab termination, both by different type of atoms or by their location, can be used to adjust electric field within the slab, creating a tool for simulation of surface properties, depending on the doping in the bulk of 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 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($000 \bar{1}$) are analyzed and compared using field dependent DFT simulations.


Introduction
Silicon carbide can exist in very large number of polytypes, of which hundreds were already identified [1]. They are created due to different stacking sequences, still preserving 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 a 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 building material was developed [3]. These structures are based on the material having its properties, extremely attractive in many applications in electronics [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 cm 2 · V −1 · s −1 and 720 cm 2 · V −1 · s −1 for 6H-SiC and 4H-SiC respectively), and finally, high blocking voltage observed in typical MOSFET devices [3]. It may be used in 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 applied 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 huge number of publications.
The most successful method of the growth of 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, though technically difficult due to very high temperatures required, was developed to the stage in which the large size SiC crystals are grown and are distributed commercially. Yet, despite considerable success, the method is still plagued by difficulties most notably by creation of micro-and nano-pipes [10], that deteriorate or even block functioning of SiC-based devices. The density of these defects was substantially reduced which allowed to use of Lely grown SiC substrates in technology of electronic devices. The pipe problem still remains essential, thus basic research and development of the technology are both needed to alleviate it completely.
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 the chemical instabilities, which affect the chemical state of SiC surface, leading to structural and chemical transformation which may contribute to the generation and preservation of the pipes during the growth. It is therefore understandable that the atomistic structures of various SiC surfaces, belonging to different polytypes, were intensively studied. Here we limit our discussion to the stoichiometric polar surfaces of hexagonal polytypes, i.e. 2H-, 4H-and 6H-SiC (0001) and SiC(0001) surfaces.
In fact, the 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. Despite the fact that the DFT studies of these surfaces have considerable value as the reference point for the investigations of the 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 surface 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 ideal 0.63Å distance. The second layer, i.e. C and Si layers move in opposite direction by 0.04Å and 0.11Å which is much smaller fraction of the ideal distance 1.88Å. Third layer relaxation is small, moving the atoms by 0.02 A and 0.03Å, respectively. Again the direction is reversed, the third layer atoms are moved towards SiC interior.
The above described relaxation type indicates that the ionic contribution in this effect is essential and the subsurface fields are important factors. Despite that these authors calculated surface band structure of both surfaces without reference to this field [11]. Using DFT LDA, they obtained indirect bandgap of 6H-SiC bulk equal to 1.97 eV, which reflects well known under-estimate 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 was attributed to difference in ionic character of both surfaces. A charge transfer towards C-surface was predicted and proved by the charge density distribution plots. It was also shown that both 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: C-state is quasisymmetric with some dominance of the internal part, while in Si-related state the charge is located exclusively in the region above the Si atom. The dispersion of both states is generally small which is due to 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 large variety of different reconstructions such as ( 3) and (9 × 9) [12]. The plethora of different symmetries was studied and several different models were proposed [13,14,15]. The subject remains under intense dispute. A new impetus to the studies of 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 favor 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 study the methodology of thick slab simulations, the role of electric field ( Fig. 1) and the influence of stacking sequence on electronic and structural properties of polar SiC(0001) and SiC(0001) surfaces. The paper is organized as follows. In Sec. 2 we give a brief overview of employed computational methods. Sec. 3 describes the construction of surface models. We analyze there the influence of finite number of layers on 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. Fig. 1 explains the subsurface field and its connection with surface acceptor and donor states together with its modeling in the slab simulations. As described in detail in Sec. 4, there is significant difference between Shockley and Tamm surface states. As showed in Fig. 1, Shockley states can be transformed into bulk states under sertain electric field conditions. On the other hand, Tamm states, irrespective of electirc field, are placed in the band gap. In Sec. 4 we present obtained structural and electronic properties of examined surfaces. Sec. 5 summarizes our work.

Method of calculation
We have carried out self-consistent total-energy calculations based on density functional theory. The code VASP developed at the Institut für Materialphysik of Universität Wien was employed [18,19,20]. The projector augmented wave (PAW) approach [21] was used in its variant available in the VASP package [22]. For the exchange-correlation functional generalized gradient approximation (GGA) according to Perdew, Burke and Ernzerhof (PBE) [23] was selected. A plane wave cutoff energy was set to 500 eV. In bulk calculations 11x11x11 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 Fig. 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 number of Si-C double layers. We used standard notation i.e. SiC (0001) surface is the one terminated by Si atoms. The slabs having 6-18 atomic layers were used, i.e. the slabs consisting of up to 36 atoms. Due to the fact that SiC(0001) and SiC(000-1) surfaces can be prepared without surface reconstruction [26], (1x1) supercells were employed. The bottommost C atoms had their broken bonds either terminated by hydrogen atoms, or left uncovered. A particular attention was paid to k-point sampling and convergence. In surface calculations Monkhorst-Pack scheme was applied, using 9x9x1 k-point mesh for bare and covered surfaces. Selected parameters were sufficient to obtain relative energies converged to few meV. To prevent artificial electrostatic field in vacuum region between repeated unit cells a dipole correction was used [27]. Three topmost Si-C layers were relaxed using a conjugate-gradient (CG) algorithm in Sec. 4.1 (Structural properties). 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 to simulate the influence of the electric fields and the doping of bulk semiconductor allows to study the interplay between the ionicity of the surface, the existence of electric field and the position of Fermi level [28,29]. Mayer and Marx [30] claimed that ZnO surface can be described using purely ionic model, where each Zn-O double layer exhibits a dipole moment perpendicular to the surface. In such model we obtain spontaneous polarization which is independent of the thickness of the slab. Top part of Fig. 3 shows that potential difference at both sides is independent of slab thickness, which means that field inside the slab is smaller for larger thickness and there is divergence for large thickness. Thus Mayer and Marx prediction is incorrect. Moreover, electron density at the ends of the slab decreases ( ∆ρ < 0) for thicker surface.
The energy of quantum states of H termination atom can be controlled by their distance to the nearest Si/C atoms. In accordance to the stability condition of fermion systems, the Fermi energy should be positioned halfway between the empty and occupied states. The energy of the surface states may be changed by electric field which is 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 slab is constant which leads to divergence for large slab thickness. This was repeated later by Noguera [32]. This prediction also turned out to be incorrect. In fact, the condition of the stability of fermion system leads to smaller field in thicker slab as it was shown in Fig. 3. From these calculations we can learn about influence of 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 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 Fig. 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 the equality of the Valence Band Maximum (VBM) on the negative charge side and Conduction Band Minimum (CBM) at the Fermi energy level. Using ionic arguments in investigation of ZnO surfaces, they claimed that the charge separation leads to uniform field within the slab, as is observed in the central part of potential diagrams in Fig. 3. This argument was based on earlier argument by Tasker [31] that charge separation creates 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 Ref. [29], this prediction is incorrect, the field is not kept constant but the difference in electric potential between both sides of the slab is preserved. This is also demonstrated in the potential plots in Fig. 3, where the difference of potentials is identical for all three slabs. Thus for thicker slab, 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 Fig. 4. These BD and DOS are calculated for the same slab models as electric field and charge diagrams presented in Fig 3. As it is shown for all these three cases, the Fermi level is kept at the energy of the quantum states associated with both    sides of the slab. Our BD and DOS prove that the Fermi level is attached to the states, independently on 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 Fig. 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 between. This charge was plotted in Fig. 3 for systematic variation of the slab thickness and it was also shown that average charge density decreases exponentially (Fig. 5). Here average charge is defined as average charge in the center of slab, i.e., interval with a width of 2.5Å from the center of the slab. In fact, the interaction between both sides of the slab have two principal contributions: electrostatic and quantum overlap. Electrostatic contribution is in fact useful for 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 how thick slab should be used for simulation of SiC surfaces. The linear part in potential profile and the negligible overlap (Fig. 5) are the necessary conditions for accurate modelling. 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 the electronic properties of the semiconductor surfaces. The extent of this modification needs direct modelling by numerical procedures. According to our best knowledge, no such awareness or the methodology can be found in the literature. However, suitable approach was developed lately [28,29]. It was demonstrated that the change of the opposite slab termination can be used for simulation of the subsurface field in polar GaN(0001) surfaces. It was also shown that the change of the distance between the hydrogen termination atoms and the layer can affect the electric field, allowing to change its magnitude and also its direction. The same approach was applied to SiC slabs. In Fig. 6 the results of such approach are shown. It is demonstrated that the approach is universal: in all three polytypes the procedure allows us to obtain different magnitude and direction of electric field within the slab. In Fig. 7 we show points in which we get potential cross sections. Potential profiles presented in Fig. 3 (2H-SiC) comes from point A in the middle of hexagon. 4H-and 6H-SiC polytypes have atoms in 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 Fig. 6.

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 cases the 12 DL slab models were used to ensure good accuracy of the calculations. We also studied the dependence of the properties on the different values of subsurface electric field, which was adjusted according to our findings described in Sec. 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 (c.f. Fig.  8). The relaxation data is collected in Table 2. From these data it follows that the significant relaxation is limited to the first two layers. The most significant changes are observed in the [0001] and [0001] direction. No significant difference between the polytypes was observed. Also there is no dependence on the electric field which indicates weakly ionic character of the polar SiC surfaces. The data is in good agreement with Ref. [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 Bond lengths related to relaxed atoms are showed in Table 2 .  in collection of diagrams plotted in Fig. 9. In Fig. 9 the band diagram for the Si-surface of 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 these 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, 4H two bands associated with bottom atoms (i.e. hydrogen and carbon) and Si-top atoms are  Fig. 6. However, for 6H-SiC (distance of 1.5Å) field is much smaller than for 2H-and 4H-SiC.
Similar analysis was made for C-terminated SiC(0001). Band structures and DOS are presented in Fig. 10. As in the previous case overlap of H-terminated and Sibottom wave function strongly affects electronic structure. Unlike SiC(0001) surface electric field changes from negative to positive values. Band associated with C-top atom states initially are located near the VBM, but under influence of electric field are moved across the entire band gap. It is dispersionless. On the other hand, band associated with bottom atoms states (H-and Si-bottom) resemble 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 atoms states) shifted outside of the conduction band, i.e., with a change of electric field within the slab it is moved to the Fermi level. Moreover, dispersion of this band disappears with increasing distance of H atoms. Fermi level is pinned to the termination atom states.
The surface states are divided into Shockley [33] and Tamm [34] states. It is often claimed that there is no real physical distinction between these two types. In fact there is 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 Fig. 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 role in the electrical conductivity. In contrast to Shockley states, Tamm states do not disappear with the potential well. This difference provides clear procedure for 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 recent literature describing way in which this states can be distinguished [35]. Zak showed 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 Fig. 3 we can observe such 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.  [36]. To denote type of states band structure and DOS have been analyzed for different methods of termination. The states which disappear in the reversed electric field are Shockley states. The states which move across the entire band gap preserving their localization at the particular surface are Tamm states.
Following changes in the band structure and DOS and comparing results in SiC(0001) and SiC(0001) side it is possible to identify surface states associated with Si and C edge atoms. This method can be used as a criterion for identification of Tamm and Shockley surface states.
Moreover, changing position of H-termination atoms we can move Fermi level which means that we can control doping type of SiC surface. In Fig. 9 and Fig. 10 doping type is changed from n-type and p-type respectively to almost undoped SiC surface.

Summary
Extensive density functional study of bare and H-terminated SiC surfaces was presented. Different SiC(0001) and SiC(0001) polytypes: 2H, 4H, 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) surface and no significant differences between them were found. It was shown that for all cases (H-terminated 2H-, 4H-, 6H-SiC) electric field within the slab can be adjusted. Field intensity dependence on the H positions for SiC(0001) and SiC(0001) have been presented (Table 3 and 4). Figures 9 and 10 prove that Fermi level is pinned to the surface states. The energy of surface related states of both Si-and C-terminated surfaces are moved by the subsurface field with respect to VBM. In case of Si-terminated surface the energy is shifted by about 2 eV while in 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 level.