First-Principles Calculation Study of Epitaxial Graphene Layer on 4H-SiC (0001) Surface∗

We performed first-principles energy-band calculations based on the generalized gradient approximation (GGA) in the Perdew-Burke-Ernzerhof (PBE) functional to determine the electronic structure of the graphene/SiC (0001) interface with a √ 3× √ 3 SiC periodicity. In the calculations, a structural model of graphene on a 4H-SiC substrate was constructed using a slab that contained four bilayers with dangling bonds terminated by hydrogen on the another surface of the slab. It is well known that the GGA method tends to overestimate the graphene layer spacing (GLS) when compared with experimental values reported previously. We included van der Waals interactions in the present calculations with the PBE-GGA exchange-correlation functional using the Grimme and TkatchenkoScheffler methods. The former scheme was observed to well reproduce the experimental value of the GLS compared with the latter one. Finally, the influence of Si and C vacancies on the electronic structure of the graphene/4H-SiC (0001) interface was investigated using these methods. [DOI: 10.1380/ejssnt.2016.107]


I. INTRODUCTION
Graphene is a unique material that can maintain a hexagonal honeycomb shape on any scale, even with a thickness of one atomic layer. Recently, the structural distinctiveness of graphene has attracted considerable attention in various fields such as electrochemical capacitors [1], field-effect transistor and inductors [2], and so on [3,4]. The energy band of single-layer graphene is characterized by the presence of linear dispersion around the K point in the two-dimensional (2D) Brillouin zone (BZ). This dispersion structure is referred to as the Dirac cone, and the mobility of electrons around the Dirac cone is very high along the line of symmetry of the 2D BZ [5,6]. Experimental values of this mobility have been estimated to be as high as 7 × 10 4 cm 2 /Vs at room temperature [7].
Thermal decomposition of SiC is a useful method for growing epitaxial graphene layers on a SiC substrate. For Si-rich (Si-terminated) surfaces such as 4H-SiC (0001) or 6H-SiC (0001), the first graphene layer deposited in contact with the SiC surface inhibits most of the effect of * This paper was presented at the 10th International Symposium on Atomic Level Characterizations for New Materials and Devices ' 3 SiC periodicity appears on the SiC surface [8], which is then reconstructed to obtain a structure with 6 √ 3 × 6 √ 3 periodicity, as confirmed by scanning tunneling microscopy and low-energy electron diffraction observations [9][10][11][12][13][14][15][16]. In contrast, C-rich (C-terminated) surfaces such as 4H-SiC (0001) and 6H-SiC (0001) display the formation of structures with 2 × 2 and 3 × 3 SiC periodicities on the SiC surface in the early stages of thermolysis [17].
To determine the electronic structure of the interface between epitaxial graphene and SiC substrate, an understanding of the sublimation of Si atoms from the SiC substrate and the epitaxial growth process of the graphene layer is required at a microscopic level. However, there is no report in the literature describing the electronic structure of a graphene/SiC interface that includes lattice defects. In this study, we report for the first time the electronic structure of a graphene/SiC interface with defects, including Si and C vacancies, using first-principles energyband calculations. In the present calculation, 4H-SiC (0001) √ 3 × √ 3 surface was used as a theoretical model, since a 6 √ 3 × 6 √ 3 structure model contains many atoms and copious computational resources are required. Van der Waals (vdW) forces were also considered using the Grimme method [18] and the Tkatchenko-Scheffler (TS) method [19] for determining interaction between graphene layers.

II. CALCULATION METHOD
SiC crystallizes into a hexagonal structure with space group P63mc and contains four formula units for a total of eight atoms [20]. This structure was firstly optimized by relaxing all the lattice constants and atomic positions using the CASTEP code [21,22] within the framework of a generalized gradient approximation (GGA) proposed by Perdew, Burke, and Ernzerhof (PBE) [23]. A slab model for investigating the growth of an epitaxial graphene layer on a SiC substrate was constructed with reference to one reported by Mattausch and co-workers [24]. Figure 1 shows top and side views of a 4H-SiC (0001) √ 3 × √ 3 structure on which a graphene bilayer with AB stacking was deposited, and the opposite surface was terminated by hydrogen. Here, we refer to this structure as model A. Model A was also optimized by relaxing all lattice constants and atomic positions using the CASTEP code. The optimization calculation was treated within the framework of the PBE-GGA considering vdW forces using the Grimme and TS methods to reproduce the interactions between graphene layers. The cutoff energy of the plane-wave expansions was set as 550 eV. Valence electrons were treated using Vanderbilt-type non-local ultrasoft pseudopotentials [25], which describe effective po-tentials of ions and tightly bound core electrons. Reciprocal space integration was performed with 25 irreducible kpoints, which correspond to a 7 × 7 × 1 grid of Monkhorst-Pack k-points for the integration of the BZ [26]. The slab in the present calculations had a vacuum region of 1 nm. Si and C vacancies were introduced by removing a Si atom and a C atom from the SiC (0001) surface, respectively. Slab models including Si and C vacancies were optimized by relaxing all the atomic positions using the CASTEP code and their lattice constants-a and c-were assumed to have the same values as ones in the model A. Table I lists the optimized lattice parameters of 4H-SiC structure, together with lattice constants of model A with and without Si and C vacancies. The optimizations of structure were performed for the ground state of the system at 0 K at 0 Pa within the framework of the PBE-GGA considering vdW forces using the Grimme and TS methods. The optimized lattice constants of 4H-SiC agree with the experimental values within 0.3% as shown in Table I. In the optimization calculation using the GGA scheme, the lattice constants of model A changed within 3.7% for the a axis and 2.9% for the c axis before the optimization. When vdW forces were added to the GGA potential, the lattice constants of model A was expanded slightly for the a axis and shortened for the c axis in comparison with the GGA calculation. In addition, lattice constant of graphene bilayer (2 × 2) deposited on 4H-SiC (0001) √ 3× √ 3 structure changed from 5.3407Å to 5.1332 A before and after the structure optimization, As a result, the mismatch of lattice constant between 4H-SiC and the graphene bilayer decreased from 7.29% to 4.20%. Table II shows the structural properties of 2-layer graphene on SiC (0001) obtained from various calculation schemes. Here, D 0 (or D ′ 0 ) and D 1 denote the distance between the top SiC (Si surface) and graphene 1-st layer and between the graphene 1-st and 2-nd layer, respectively. Varchon and co-workers performed a structure optimization calculation for the GGA framework and estimated the interlayer spacing between the first and second graphene layers to be 3.8Å [27]. This value is large in comparison with the experimental value reported for graphite crystal (3.5Å). Using the local spin density approximation (LSDA) method, Mattausch and co-workers calculated the interlayer spacing of graphene layers on a SiC (0001) surface that had undergone structure optimization and estimated it to be 3.3Å [24]. This value is closer to the experimental value for graphite when compared with the value obtained by the GGA approach, as mentioned above. We examined two methods considering vdW forces, GGA + Grimme and GGA + TS, to determine if they can reproduce the experimental interlayer distance of graphite. Table II summarizes the results of our calculations in comparison with reported values in the literature. It is clear that the GGA + TS method reproduces the interlayer distance of graphite very well. Therefore, we mainly used the GGA + TS scheme in the following calculations. Figure 2 shows the energy-band structure of model A obtained by the GGA + TS method. This struc- ture is very similar to those obtained by the GGA and LSDA methods, which suggests that the vdW force barely changes the electronic structure of model A. Further, using the GGA + TS method, we determined the electronic structure of a graphene/SiC interface that included Si or C vacancies. Three types of Si vacancies (Si1v, Si2v, and Si3v) and three types of C vacancies (C1v, C2v, and C3v) on the SiC (0001) surface were used for determining the electronic structure, as shown in Table III. Models Si1v, Si2v, and Si3v were constructed by removing Si1, Si2, and Si3 atoms from the SiC surface, respectively. Models C1v, C2v, and C3v were also constructed in the same manner. The lattice constants were assumed to have the same values as those in the original model (model A) without vacancies.

III. RESULTS AND DISCUSSION
We determined the relative stability of models with vacancy doping compared with the original model from the vacancy formation energy (E v ) for each type of supercell with a vacancy. The difference in total energy between undoped and vacancy-doped supercells is defined as the vacancy formation energy: where E host and E 1 are the total energies of undoped and vacancy-doped models with the same supercell size, respectively, and µ 1 is the energy per atom for Si atom (µ Si = −107.3376 eV) and C atom (µ C = −154.9610 eV). The chemical potentials are measured from particular reference systems and the E Si and E C refer the energies of reference systems per atom. We use Si or C crystal with diamond-type structure as reference systems. The vacancy formation energies were estimated to be 4.058 eV for Si1v, 3.324 eV for Si2v, and 4.050 eV for Si3v, whereas these were 1.599 eV for C1v, 1.590 eV for C2v, and 1.533 eV for C3v. These values are lower than those obtained for the SiC bulk reported by Ping and co-workers [31], who reported that the formation energies of Si and C vacancies are 7.859 eV and 2.042 eV, respectively. The present results suggest that Si and C vacancies are easily formed on a SiC surface on which a graphene layer is stacked in comparison with the SiC bulk. Figure 3 shows model Si1v and its energy-band struc- ture after structure optimization. When Si1 vacancy was doped into model A, the energy gap increased around the Dirac cone at the K point. Figure 4 indicates the displacement of each atom in model Si1v after structure optimization. When a Si vacancy was introduced into model A cell, C atoms occupying G1 and G2 sites (in the first layer of graphene as shown in Fig. 1) were shifted into the vacuum layer, and the topmost Si atoms on the SiC (0001) surface were displaced toward the bulk, although other atoms on the surface were barely displaced. As a result, the spacing between the first layer of graphene and topmost Si atom was 2.571Å, which is greater than that in model A without the Si vacancy, as listed in Table  III. Moreover, the distance between the first and second layers of graphene was 3.297Å, which is lower than that in model A. These results suggest that covalent bonding between the first layer of graphene and SiC (0001) surface weakened with the formation of the Si vacancy and the graphene layer exfoliated from the SiC surface. The weakening of the covalent interaction between the first layer of graphene and SiC (0001) surface is clear from a comparison of the charge density map before and after the introduction of the Si vacancy, as shown in Fig. 5. The energy-band structure and geometric optimization of model Si2v are similar to those of model Si1v because the Si2 atom has a local coordination environment that resembles that of the Si1 atom. On the other hand, model Si3v leads to results that are different from those of mod-els Si1v and Si2v. When a Si vacancy was doped at the Si3 site, the displacements of surface atoms on SiC (0001) were small in comparison with the cases of the Si1v and Si2v models. Figure 6 shows (a) model C1v, (b) its energy-band structure, and (c) the displacement of each atom in model C1v after structure optimization. This model has two notable characteristics. First, the Si atom at the Si3 site significantly changed its position upon doping with C1 vacancy in comparison with Si1 and Si2 sites. In other words, the Si3 atom approached closer to the C1 vacancy site than before vacancy doping. Similar behavior was observed for models C2v and C3v after structure optimization. Second, the energy-band structure of model A did not change significantly even when a C atom vacancy was doped into the model, especially around the K point of the 2D BZ. This behavior in reconfiguration that accompanies the formation of a C vacancy implies that the displacement of the Si3 atom absorbed the distortion energy with the formation of the C vacancy and was useful for stabilizing the energy of the system.

IV. CONCLUSION
In the present study, we performed first-principles energy-band calculations based on the GGA-PBE method including van der Waals interactions using the Grimme or TS methods to determine the electronic structure of the graphene/SiC (0001) interface. The former scheme was observed to well reproduce the experimental value of the GLS compared with the latter one. Then, the influence of Si and C vacancies on the electronic structure of the graphene/4H-SiC (0001) interface was calculated using the GGA-PBE + TS method. When Si1-or Si2-vacancies were introduced into model A, it was suggested that the covalent bonding between the first layer of graphene and SiC (0001) surface weakens and the graphene layer exfoliates from the SiC surface. On the other hand, the distortion energy with the formation of C1-, C2-or C3vacancies was relaxed by the displacement of Si3 atom.