Stacking domains in graphene on silicon carbide: a pathway for intercalation

Graphene on silicon carbide (SiC) bears great potential for future graphene electronic applications because it is available on the wafer-scale and its properties can be custom-tailored by inserting various atoms into the graphene/SiC interface. It remains unclear, however, how atoms can cross the impermeable graphene layer during this widely used intercalation process. Here we demonstrate that, in contrast to the current consensus, graphene layers on SiC are not homogeneous, but instead composed of domains of different crystallographic stacking. We show that these domains are intrinsically formed during growth and that dislocations between domains dominate the (de)intercalation dynamics. Tailoring these dislocation networks, e.g. through substrate engineering, will increase the control over the intercalation process and could open a playground for topological and correlated electron phenomena in two-dimensional superstructures.

Graphene on silicon carbide (SiC) bears great potential for future graphene electronic applications because it is available on the wafer scale and its properties can be custom tailored by inserting various atoms into the graphene/SiC interface. It remains unclear, however, how atoms can cross the impermeable graphene layer during this widely used intercalation process. Here we demonstrate that in contrast to the current consensus, graphene layers grown in argon atmosphere on SiC are not homogeneous, but instead are composed of domains of different crystallographic stacking as they have been observed in other systems. We show that these domains are intrinsically formed during growth and that dislocations between domains dominate the (de)intercalation dynamics. Tailoring these dislocation networks, e.g., through substrate engineering, will increase the control over the intercalation process and could open a playground for topological and correlated electron phenomena in two-dimensional superstructures. DOI: 10.1103/PhysRevMaterials. 2.104005 Graphene can routinely be produced on the wafer scale by thermal decomposition of silicon carbide (SiC) [1][2][3][4]. Due to the direct growth on SiC(0001) wafers, epitaxial graphene (EG) naturally forms on a wide-band-gap semiconductor, providing a doped or insulating substrate compatible with standard complementary metal-oxide semiconductor (CMOS) fabrication methods. Hence, EG is a contender for future graphene electronic applications such as power electronics [5,6], high-speed transistors [7], quantum resistance standards [8], and terahertz detection [9]. In EG, the first hexagonal graphene layer resides on an electrically insulating monolayer of carbon atoms that are sp 3 bonded to silicon atoms of the SiC(0001) surface [1][2][3]10]. The presence of this so-called buffer layer strongly affects the graphene on top, e.g., by pinning the Fermi level. Consequently, the graphene properties can be tuned via intercalation of atoms into the bufferlayer/SiC interface. The intercalation of hydrogen is most widely used and results in the conversion of the buffer layer to a quasifreestanding graphene (QFG) layer by cutting the silicon-carbon bonds and saturating silicon dangling bonds with hydrogen. This treatment reverses the graphene doping from n type to p type and improves the mobility [11,12]. Intercalation of heavier atoms is used to further tailor the graphene properties, e.g., to form pn junctions [13,14], magnetic moments [15], or potentially superconducting [16] and topologically nontrivial [17] states.
It has been realized that the quality of the grown graphene can be greatly improved by reducing the desorption rate of silicon atoms (which allows for a raised growth temperature), for example, by encapsulation of the SiC chip [3] or by growth in a low-pressure silane environment [18,19] or in an ambient-pressure argon atmosphere [1,2]. Graphene (EG and QFG) grown on SiC using these methods appears homogeneous with low-defect concentration in most techniques [1][2][3][4]11]. Together with the fact that layers span virtually unperturbed over SiC substrate steps [20][21][22], this has led to the consensus of perfectly crystalline graphene. On the other hand, two observations point to a less perfect sheet. First, the charge carrier mobility is generally low, even at cryogenic temperatures [2,12]. Second, an ideal graphene sheet is impermeable even to hydrogen [23,24], whereas a wide variety of atomic and molecular species has been intercalated into EG [11,[14][15][16][17][18]. Stacking domains as they have been observed in vacuum-grown graphene on SiC [25] and freestanding bilayer graphene [26,27] could explain these contradictions. In this paper, we demonstrate that graphene grown on SiC in argon atmosphere is, in fact, less homogeneous than widely believed but is fractured into domains of different crystallographic stacking order. We use advanced low-energy electron microscopy (LEEM) methods and ab initio calculations to show that those domains are naturally formed during growth due to nucleation dynamics and built-in strain. Their presence is thus intrinsic to all graphene-on-SiC materials, including high-quality graphene grown in an argon atmosphere.  show bright-field LEEM images of two QFG samples (see the Appendix for details on sample growth and hydrogen intercalation) with areas of different graphene thickness. Bright-field images are recorded using specularly reflected electrons that leave the sample perpendicular to the surface [see Fig. 1(c)]. The main contrast mechanism in this mode is the interaction of the imaging electrons with the thickness-dependent, unoccupied band structure of the material, which is used to unambiguously determine the number of graphene layers [28][29][30]. Large, homogeneous areas of bilayer, trilayer, and four-layer graphene can thus be distinguished in Figs. 1(a) and 1(b), supporting the notion of perfect crystallinity.
In stark opposition to this generally accepted view, the dark-field images in Figs. 1(d) and 1(e) clearly reveal that all areas are actually fractured into domains of alternating contrast. The symmetry breaking introduced in dark-field imaging, where the image is formed from one diffracted beam only [cf. Fig. 1(f) and the Appendix], leads to strong contrast between different stacking types of the graphene layers [25,31]. In fact, the contrast between different domains inverts [ Figs At first glance, the observation of different stacking orders might seem surprising, as it is known that graphene layers grown on SiC(0001) are arranged in Bernal stacking [3]. However, two energetically equivalent versions of Bernal stacking (AB and AC) exist and have been observed in other graphene systems [25][26][27]32]. The AC stacking order can be thought of either as AB bilayer, where the top layer is translated by one bond length, or, alternatively, as a full AB bilayer rotated by 60 degrees [Figs. 2(a) and 2(b)]. Consequently, AB and AC stacking are indistinguishable in bright-field imaging. Subsequent layers can be added in either orientation, generating more complicated stacking orders for trilayer and beyond.
In order to identify the exact stacking in each area, we simulate bilayer and trilayer graphene slabs in different stacking orders and compare their reflectivity with measured lowenergy electron reflectivity spectra. The latter are extracted from the intensity of an area in a series of spectroscopic LEEM images recorded at different electron landing energy (see Movies 1 and 2 in the Supplemental Material [33] for such measurements of the area in Fig. 1(b) in bright-field and dark-field geometry, respectively). While different domains show identical bright-field reflectivity (cf. Supplemental Material, Fig. 1 [33]), dark-field spectra extracted from different bilayer domains [marked blue and orange in Figs. 2(c) and 1(e)] are clearly distinguishable. Moreover, four distinct reflectivity curves are observed for trilayer graphene [ Fig. 2(d)]. Figures 2(e) and 2(f) show theoretical dark-field spectra, obtained by ab initio calculations (see the Appendix for computational details), of different bilayer and trilayer stacking orders, respectively. The excellent agreement of theoretical and experimental data in Figs. 2(c) and 2(e) is clear evidence that the assignment of Bernal AB and AC stacking orders for different bilayer domains is correct. Moreover, the comparison of Figs. 2(d) and 2(f) shows that by using these dark-field LEEM methods, we can distinguish the more complicated trilayer stacking orders: Bernal, ABA (cyan) and ACA (pink) versus rhombohedral, ABC (purple) and ACB (brown). Due to the small electron penetration depth in LEEM, however, the spectra fall into two families (ABA and ABC versus ACA and ACB) dominated by the stacking order of the top two layers.
In addition to their stacking orders, bilayer graphene and thicker areas differ in the morphology of the stacking domains [cf. Figs. 1(d) and 1(e)], which indicates two distinct formation mechanisms. Most notably, bilayer domains are smaller, triangular, and relatively regular. Similar morphologies have been observed in freestanding bilayer graphene, both etched from graphene-on-SiC [27] and transferred from copper [26,34], where they were linked to strain between the layers introduced during sample growth or fabrication. While uniform strain causes a Moiré reconstruction [ Fig. 3(a)], it is often energetically favorable to form domains of commensurate, optimal Bernal stacking. In this case, all strain is concentrated into the domain walls, thus forming dislocation lines [26,27], as sketched in Fig. 3(b). Upon close examination of Fig. 1(b), the network of these dislocations is visible as dark lines in our bright-field measurements. These lines correspond to the patterns observed by Speck et al. [18]. Such domain walls were predicted to host topological edge states [35,36], which has also been experimentally confirmed recently [32,37,38].
The size of the triangular domains shrinks for increasing uniform strain, while anisotropic strain causes domains elongated perpendicular to the strain axis. The observed average domain diameter of ∼100−200 nm coincides well with relaxation of the 0.2% lattice mismatch between buffer layer and first graphene layer [39] (see calculation in the Supplemental Material [33]). We thus conclude that the triangular domains in bilayer graphene result from strain thermally induced during growth and from the lattice mismatch with the SiC substrate and are intrinsic to the growth process. The presence of elongated triangular domains indicates nonuniform strain due to pinning to defects and substrate steps.
The larger, irregularly shaped domains that dominate trilayer and four-layer areas [Figs. 1(d) and 1(g)] can be explained by nucleation kinetics. To test this hypothesis, we study EG samples where the growth was stopped shortly after the nucleation of bilayer areas to prevent their coalescence (see the Appendix). The resulting small bilayer islands on monolayer terraces are shown in bright-field and dark-field conditions in Figs. 3(c) and 3(d), respectively. We observe that bilayer areas with a diameter below ∼300 nm form single domains of constant stacking order [either bright or dark in Fig. 3(d)] and that AB and AC stacked bilayer islands occur in roughly equal number. This indicates that new layers nucleate below existing ones in one of the two Bernal stacking orders randomly [2,3,10]. At the elevated growth temperature, dislocations in the existing layers can easily move to the edge of the new island where they annihilate. As islands of different stacking grow and coalesce, new dislocation lines are formed where they meet [cf. Fig. 2(a)]. This opens the interesting possibility to engineer the dislocation network by patterning the SiC substrate before graphene growth.
Notably, we observe strain-induced domains also in monolayer EG [ Fig. 3(d)] and between the bottom two layers in trilayer QFG (visible only for some energies, e.g., 33 eV in the Supplemental Material, Fig. 2 [33]). The prevalence of these triangular domains in all EG and QFG samples between the two bottommost layers demonstrates that stacking domains are a direct consequence of the epitaxial graphene growth and, consequently, are a general feature of this material system. The resulting dislocation network explains the linear magnetoresistance observed in bilayer QFG [40] and might be an important culprit for the generally low mobility in EG and QFG [12].
The presence of these strain-induced domains in EG as well as QFG raises the question of their role during (hydrogen) intercalation. Since the high hydrogen pressures necessary for intercalation are not compatible with in situ imaging, we investigate the inverse process. Figure 4(a) shows a time series of bright-field LEEM images of the area shown in Fig. 1(b) recorded at ∼1000 • C (cf. Movie 3 in the Supplemental Material [33]). At this temperature, hydrogen slowly leaves the SiC-graphene interface [11,12] and the n-layer QFG is transformed back to (n − 1)-layer (+ buffer-layer) EG. The change in the reflectivity spectrum accompanied with this conversion (cf. Supplemental Material, Fig. 1 [33]) yields a strong contrast [e.g., dark in the bilayer in Fig. 4(a)] and enables capture of the full deintercalation dynamics. Deintercalation starts at distinct pointlike defect sites where hydrogen can escape and proceeds in a highly anisotropic fashion. An overlay of the half-deintercalated state (15 min) with an image of the dislocations in the initial surface [ Fig. 4(b)] shows that deintercalation happens preferentially along dislocation lines. Although the dislocation lines are slightly mobile at higher temperatures [cf. Figs. 4(c) and 4(d) before and after deintercalation, respectively], their overall direction and density is preserved during the process. The local deintercalation dynamics reveal details of the underlying microscopic mechanism. Figures 4(e) and 4(f) show that deintercalation fronts move roughly linearly in time both perpendicular and parallel to dislocation lines. The velocity of the deintercalation fronts, however, is much larger parallel to dislocation lines (up to v = 95 nm s −1 ) than perpendicular to them (v ⊥ ≈ 0.1 nm s −1 ). The linear movement rules out that deintercalation is limited by hydrogen diffusion, but indicates that hydrogen desorption at the deintercalation front is the limiting factor. The nonlinear growth of the fraction of deintercalated area A EG [ Fig. 4(g)] demonstrates that deintercalation is also not capped by the venting of hydrogen from the defects where deintercalation starts [7 min in Fig. 4(a)]. We conclude from these observations that the hydrogen-desorption barrier is smaller within the dislocation lines than within the Bernalstacked domains, possibly triggered by the higher lattice strain in the former. While v ⊥ is the same for all areas, v varies from 0.2 to 95 nm s −1 [marked yellow and white in Fig. 4(a), respectively], suggesting that the deintercalation process is strongly affected by the precise atomic details of the dislocations. These findings indicate that not only the deintercalation, but also the intercalation of hydrogen and other species, all of which cannot penetrate graphene, is dominated by the presence of stacking domains. Consequently, their manipulation, e.g., by patterning the substrate, will open a route towards improved intercalation and tailored QFG on the wafer scale.
We conclude that graphene on SiC is a much richer material system than has been realized to this date. Specifically, we show that domains of AB and AC Bernal-stacking orders are always present in this material, even for high-quality argongrown samples, and even though its layers appear perfectly crystalline to most methods. We deduce that these domains are formed between the two bottommost carbon layers (either graphene and buffer layer for EG or bilayer QFG) by strain relaxation. In addition, the nucleation of grains of different stacking order during growth causes larger domains in thicker layers. We show that dislocation lines between domains dominate hydrogen-deintercalation dynamics, highlighting their importance for intercalation as well. By engineering these dislocation networks, we foresee wide implications for customized QFG for electronic applications. Moreover, the dislocation networks observed here could yield a wafer-scale platform for topological [32,[35][36][37][38] and possibly strongly correlated electron [41][42][43] phenomena when tailored into periodic structures.
The raw data of the findings presented here is available online [44].

Sample fabrication.
Graphene growth is carried out on commercial 4H-SiC wafers (semi-insulating, nominally on axis, RCA cleaned) at ∼1700 • C and 900 mbar Ar pressure for ∼30 min as described in Ref. [2]. To convert EG to bilayer QFG via hydrogen intercalation, the sample is placed in a carbon container and heated to 970 • C for 90 min at ambient hydrogen pressure, as described in Refs. [11,12]. Samples with small bilayer patches on large substrate terraces are achieved in a three-step process. First, SiC substrates are annealed at ∼1700 • C and 900 mbar Ar pressure for 30 min in a SiC container to enable step bunching. Second, unwanted graphitic layers formed during this process are removed by annealing the sample at 800 • C in an oxygen flow for 30 min. Third, graphene growth is carried out as described above.
Low-energy electron microscopy. The LEEM measurements are performed using the aberration-correcting ESCHER LEEM facility [45], which is based on a commercial SPECS P90 instrument and provides high-resolution imaging. Limitations on the angles of the incident and imaging beams make dark-field imaging in the canonical geometry, where the diffracted beam used for imaging leaves the sample along the optical axis, impossible. Instead, we use a tilted geometry where the incident angle is chosen such that the specular beam and the refracted beam used for imaging leave the sample under equal, but opposite, angles [illustrated in Figs. 1(f) and 1(i)]. The tilted incidence yields an in-plane k vector, which influences the reflectivity spectrum [30,46]. This is taken into account in our calculations, but needs to be considered when comparing to other LEEM and low-energy electron diffraction (LEED) data. Microscopy is performed below 2 × 10 −9 mbar and at 600 • C to prevent the formation of hydrocarbon-based contaminants under the electron beam. Images are corrected for detector-induced artifacts by subtracting a dark-count image and dividing by a gain image before further analysis. Figure 3 is corrected for uneven illumination by dividing by the beam profile. Additionally, the minimum intensity in images shown is set to black and maximum intensity is set to white to ensure visibility of all details. All dark-field images and images showing dislocation lines are integrated for 4 s; all other images for 250 ms.
Computations. All calculations were performed with a full-potential linear augmented plane-wave method based on a self-consistent crystal potential obtained within the local density approximation, as explained in Ref. [47]. The ab initio reflectivity spectra are obtained with the all-electron Bloch-wave-based scattering method described in Ref. [48]. The extension of this method to stand-alone two-dimensional films of finite thickness was introduced in Ref. [49]. Here, it is straightforwardly applied to the case of finite incidence angle to represent the experimental tilted geometry. An absorbing optical potential V i = 0.5 eV was introduced to account for inelastic scattering: the imaginary potential −iV i is taken to be spatially constant over a finite slab (where the electron density is non-negligible) and to be zero in the two semi-infinite vacuum half spaces. In addition, a Gaussian broadening of 1 eV is applied to account for experimental losses.