of graphite by active elements First-Principles Study of Graphene-6H SiC Surface Interactions

— Interactions of graphene with 6H-SiC {0001} surfaces are numerically investigated from first principles. In order to describe the bulk structure and its 6 bilayer thick surfaces correctly, bare and dipole-corrected atomic relaxations are considered. The obtained lattice parameters and bulk modulus are in good agreement with experimental values. The calculated indirect band gap width of 2.10 eV is smaller than the experimental value due to the nature of the computational method. Geometrical optimization of the surfaces, where dipole correction is applied, reveals that the first two bilayers displace significantly, where the relaxations of the very top bilayer is more pronounced. Band structures of the {0001} surfaces possess two flat bands around the Fermi level due to unsaturated bonds on opposite faces. When one layer of C atoms are introduced on the Si-terminated surface, it behaves as a tightly-bound buffer layer. This is also the case for the C-terminated surface when van der Waals interactions are taken into account. In contrast, disregarding these interactions yields free-standing graphene like behavior for the first C overlayer. On both surfaces, the second C overlayer is free-standing where the corresponding band structures incorporate Dirac-cone like features. high-temperature superconductors, photonic and phononic crystals, quantum well structures and green synthesis of metallic nanoparticles.

Graphene exhibits different physical properties on opposite surfaces of hexagonal SiC terminated by Si or C atoms [20,21]. Thus, understanding of graphene-substrate interactions is vital for design of graphene based nanodevices and investigation of their physical properties. For instance, the first graphene layer grown on Si-terminated face of SiC is tightly bound to surface, which has a complex reconstruction, while free-standing graphene layers are obtained starting from the second layer [20]. In contrast, the first layer grown on the Cterminated face of SiC is free standing, where there exist rotational defects between graphene layers [22].
Considering the complex surface reconstruction of hexagonal SiC, first numerical studies employing densityfunctional theory (DFT) calculations on graphene-SiC interactions involved relatively small surface reconstructions [8,23,24]. Although these studies suggested that the very first layer on both Si-terminated and C-terminated surfaces were tightly-bound buffer layers, more recent studies indicate that the buffer layer does not appear on the C-terminated face [25,26]. Thus, it is important to understand the correct behavior of graphene on opposite faces of hexagonal SiC.
In this work, graphene-substrate interactions for the commonly employed 6H-SiC{0001} surfaces are numerically investigated from first principles via DFT calculations. Structural and electronic properties of graphene on opposite surfaces are compared to understand the formation mechanism and physical properties of graphene on these surfaces.

II. MATERIAL AND METHODS
Crystal structure of the 6H silicon carbide substrate in the S3 configuration [27] is seen in Fig. 1(a). Here, each carbon (silicon) atom is tetravalently bonded to four silicon (carbon) atoms, where the bond length is 1.89 Å and the C-Si-C bond angle is 109.57 degrees. It is designated 6H since its unit cell

First-Principles Study of Graphene-6H SiC Surface Interactions
Ahmet Cicek and Bulent Ulug G in the hexagonal crystal class comprises 6 bilayers (BLs) in the ABCACB… stacking order [28]. The Si-C bonds follow a zig-zag pattern over the (1120) plane shown in Fig. 1(b).
Positions of the atoms on this plane are denoted as zi (i=1…12) (Fig. 1(c)), whereas the vertical distances in a BL and between two BLs are denoted as di and hi, respectively. Thus, 6(di+hi)=c, lattice constant in the [0001] (z) direction. The first Brillouin zone for 6H SiC and the high-symmetry points over the {0001} surfaces are given in Fig. 1(d). Graphene is grown on 6H SiC {0001} surfaces. Here, the (0001) surface is terminated by a sub-layer composed of only Si atoms, whereas the (000 1) surface is terminated by a C sublayer. Growth on 6H SiC takes place through a thermal process where Si atoms sublime and the remaining C atoms on the very top of the surface form surface reconstructions [29]. The most common surface reconstruction for epitaxial graphene growth is denoted as 6 3 6 3 30 R  [30]. However, in this study a more computationally afordable surface reconstruction of 3 3 30 R  is considered [8,23]. The top view of the (0001) surface of 6H SiC is seen in Fig.  2(a), where the enumeration of top BL atoms is also shown. Also in Fig. 2(a) is the graphene ML considered on 6H SiC surface and the enumeration of C atoms in this layer.
Geometry of a graphene ML on 3 3 30 R  reconstruction of 6H SiC is on the right of Fig. 2(a). The side view of the geometry under study for the case of (0001) surface is given in Fig. 2(b). Here, a BL graphene is demonstrated, where dGA and dGG are graphene-substrate and graphene-graphene distances, respectively. In the case of graphene BL, the ABAB… stacking order is considered.
The numerical problems in Fig. 1 and Fig. 2 are three dimensional (3D) in nature. However, the study of solid surfaces requires termination of an otherwise infinite system along a given direction. While the unit cell in Fig. 1(a) defined with Bloch periodic boundary condition is sufficient for the study of bulk crystals, study of surfaces requires supercell structures, where the finite-size crystal is surrounded by a vacuum layer of thickness dvac (h0 and h6 in Fig. 1(c)). In this study, an asymmetric supercell, as seen in Fig. 3, due to electrical polarity of the 6H SiC {0001} surfaces is considered. Since the Si and C atoms on opposite faces have unequal electron affinities, a net surface charge density develops [31]. For large number of cells, this polarization can be ignored. However, for small supercells as in Fig. 3, ignoring it may result in significant error. Graphene-6H SiC interactions are numerically investigated from first principles (ab-initio) via DFT calculations implemented in Quantum Espresso software [32]. This software employs a plane-wave basis set for the calculation of density functionals from Kohn-Sham equations. The local density approximation (LDA) and generalized gradient approximation (GGA) for the exchange-correlation functionals are utilized. Perdew-Zunger (1981) parametrization [33] in case of LDA and Perdew-Wang 91 (PW91) parametrization [34] for GGA are employed.
Brillouin zone integrals are discretized through Monkhorst-Pack (MP) [35] k-point meshes. In case of bulk 6H SiC and its {0001} surfaces, a 12x12x2 MP grid is used. In free-standing graphene calculations, 6x6x6 MP grids are employed. In addition, the cut-off energy in all calculations is 50 Ry.
The calculations are conducted in four stages. First, the equilibrium geometry for bulk 6H SiC in Fig. 1(a) is obtained through unconstrained geometric relaxation. Here, the convergence criterion is reduction of Helmann-Feynman forces below 2.0×10 -4 Ry/Bohr (5.1×10 -3 eV/Å). In addition, the band structure of bulk 6H SiC is obtained. In the second stage, 6H SiC {0001} surfaces are geometrically optimized.
Here, the 2D model in Fig. 1(c) is considered. In each case, 3 of 6 BLs are relaxed, while the rest are fixed. The vacuum layer thickness in these calculations is adopted as dvac=20 Å.

III. RESULTS AND DISCUSSION
Numerically-obtained lattice constants (a and c) from the Birch-Murnaghan equation of state fitting [36] are given in Table 1. They are in good agreement with both experimental [37][38][39] and previous numerical [28,[40][41][42] data. However, LDA approach resulted in a slightly closer value to the experimental value of a=3.081 Å [37]. In addition, the bulk modulus (B0) values are also close to the experimentallyobtained values [38,39], while GGA approach yields a closer value. Besides, the derivative of bulk modulus ( 0 B  ) is considerably larger than the experimental value [39]. Results of geometrical relaxations of the bulk 6H SiC indicates that the relaxations exhibit repeating behavior with 3 BL periods. Therefore, di+3=di and hi+3=hi (i=1,2,3). The values of di and hi in Table 2 are in good agreement with the results of Käckell et al [40]. Overall, the atoms experience relaxations around 1 pm. The band structure of geometrically optimized 6H SiC is given in Fig. 4. Here, an indirect band gap of EG=2.10 eV between the  and L points is obtained through LDA calculations, which is considerably smaller than the experimental value of 3.08 eV [43]. This is a well-known underestimation in LDA band structure calculations [44]. A similar result is obtained for GGA calculations. To calculate the band gap correctly, hybrid functionals are utilized, which are not considered in this work.  By employing the supercell in Fig. 3 and dipole correction, geometrical optimization of 6H SiC (0001) (Si terminated) and (000 1) (C terminated) surfaces are carried out. In both cases, only the first 3/6 BLs are relaxed and the variations of the di and hi values with respect to the corresponding values in bulk 6H SiC are calculated. The results are presented in Table 3. It is clearly seen that the relaxation is largest in the top BL for both cases. Percent variations of di and hi for the second BLs are considerably smaller, Table 3. Moreover, relaxations in the third BLs are negligible. Besides, for relaxation from the Siterminated surface, LDA calculations yield values almost twice of the values obtained through GGA calculations. Closer inspection of Table 3 reveals that relaxation of the first C sublayer on the (000 1) surface is relatively higher. This is related to the charge transfer towards the (000 1) side, whose polarity is higher.
Band structures of the 6H SiC {0001} surfaces overlaid on the projected bulk bands are given in Fig. 6. The presented band structures are obtained through LDA calculations by considering the relaxed surfaces. The Fermi level (EF) is above the valence band maximum of the bulk (0 eV) by 0.5 eV and 0.3 eV for the Si-and C-terminated surfaces, respectively. In both cases, two surface bands localized around EF are observed, Figure 6. These bands originate from the dangling p states of the Si1 and C6 atoms on the opposite faces. The band due to C6 is flat since this side is more planar [45]. In addition, a single band in the stop band around -10 eV is also seen for both surfaces, Fig. 6. The band structures of ML and BL graphene are given in Fig. 7. Since results of LDA and GGA calculations agree well, only the results of LDA calculations will be given in the following band structures. Fig. 7(a) shows that the EF and Dirac point (ED) coincide, as expected. Moreover, the Dirac cone is observed around these energies. In case of BL graphene in ABAB.. stacking, a slight band gap opening is observed, Fig. 7(b).  Fig. 2(a). It contains 18 Si and 18 C atoms belonging to the substrate and 8 C atoms belonging to the graphene superstrate. Fig. 8 shows that the first C layer on the (0001) surface does not exhibit free-standing graphene behavior, where it forms a buffer layer. When the first C layer initial distance to the (0001) surface is 2.6 Å or 3.0 Å, it approaches the surface upon relaxation, where the CG12 and CG17 atoms form mostly covalent bonds with the SiA12 and SiA13 atoms, Fig. 8(a). Besides, the SiA11 atom, whose dangling bond is not saturated, experiences less upwards displacement thant the SiA12 and SiA13 atoms. In case of the first C (buffer) layer, the value of the London s6 parameter (0 or 1) does not have significant effect on final relaxations on the (0001) surface.
The buffer layer in Fig. 8 is, on the average, separated from the top Si layer of the (0001) surface by 2.23 Å and exhibits a buckling of 0.27 Å (11% of graphene lattice constant). These results are close to the value (2.58 Å) reported by Mattausch and Pankratov [8]. Thus, the first C layer on the 6H SiC (0001) surface is a buffer layer.
The second C layer on the (0001) surface behaves like freestanding graphene, Fig. 8(b). The distance of this layer to the buffer layer is calculated as 3.39 Å and 3.02 Å, respectively, for s6=0 and 1. Furthermore, the buckling in the second C layer is much less significant (around 2 pm).
Band structures for one and two C overlayers on 6H SiC (0001) surface are given in Fig. 9. Inclusion of dispersion forces by varying the London s6 parameter does not have a significant effect on the obtained bands. The flat bands around 0.5 eV in each case are related to the unsaturated dangling states of the substrate. Fig. 9(a) reveals that the graphene Dirac cone does not appear for the buffer layer. On the other hand, although overlapping with the substrate bands, a Diraccone like behavior is observed in Fig. 9(b) for two C overlayers. The final atomic positions for the geometrically optimized system of one C overlayer on the (000 1) surface of 6H SiC are given in Fig. 10(a). In contrast to Fig. 8(a), the first C layer is not bonded to the surface when London s6=0. Thus, inclusion of the dispersion forces has a significant effect on the resulting interactions. When s6=1, the first C overlayer is covalentlybonded to the (000 1) surface. The second case is similar to the results of Mattausch and Pankratov where local spin density approximation (LSDA) without explicitly considering dispersion forces [8,23]. In contrast, calculations by Magaud et al [25] on C-rich (2x2)C surface reconstruction of 4H SiC suggest that the first C overlayer behaves like free-standing graphene. The latter is in agreement with experimental results and the DFT results for the 13 13 46.1 R  surface reconstruction of 6H SiC [46].
When s6=1, The relaxations for the (000 1) surface are significantly larger than those on the (0001) surface for both one C overlayer and two overlayers. The vertical displacement in the buffer layer is as high as 0.37 Å. Due to buckling in the buffer layer, the minimum bond angle is found to be 115.03°. The average distance between the buffer layer and the (000 1) decreases to 2.02 Å. In contrast, the second C layer behaves like free-standing graphene, Fig. 10(b), as in the case of (0001) surface in Fig. 8(b). The fluctuations in the first C overlayer when s6=0 are as low as 2.0 pm, where the C layers behave like free-standing graphene. The bond length in this layer does not significantly change from the initial value of 1.52 Å. Moreover, the graphene-substrate is measured as 2.94 Å.
The covalent bonding to surface on the right of Fig. 10(a), where the van der Waals forces are included in LDA calculations, stems from the nature of the numerical method and the van der Waals radii (r0=1.45 Å) of the C atoms. While, the LDA calculations usually tend to over-bonding, the initial distance between the released C atoms is close to r0. Thus, a tightly-bound atoms are obtained in these calculations.
The sharp contrast in Fig. 10(a) is also clearly visible in Fig.  11(a), where the band structure of one C overlayer on the (000 1) surface is presented. While, a Dirac-cone like behavior, hallmark of graphene band structure, is observed for s6=0, only four flat bands due the the buffer layer C atoms and unsaturated substrate C atoms appear for s6=1. The discrepancy is also visible for two C overlayers, Fig. 11(b).

IV. CONCLUSION
In conclusion, density-functional theory calculations of graphene on 6H silicon carbide {0001} surfaces reveal that the free-standing graphene behavior is observed starting from the second carbon overlayer on either surface. In case of the Siterminated (0001) surface, the first C layer is a buffer layer, covalently bonded to the surface. This is the case irrespective of dispersion forces are accounted for. In contrast, the first C overlayer acts like a free-standing graphene layer on the Cterminated 6H SiC (000 1) surface when the van der Waals interactions between the overlayer and surface are neglected. However, intrinsic nature of LDA calculations and the van der Waals radii of C atoms lead to covalently-bonded first C overlayer (i.e. a buffer layer) when van der Waals interaction is taken into account. In all cases, both first C overlayers and the top substrate bilayer experience significant displacement from their equilibrium positions. However, the fluctuations on the second C overlayer is negligible. The results of this work can be used in estimating behavior of epitaxially-grown graphene few layers on SiC surfaces and developing their applications.