Large-area epitaxial growth of curvature-stabilized ABC trilayer graphene

The properties of van der Waals (vdW) materials often vary dramatically with the atomic stacking order between layers, but this order can be difficult to control. Trilayer graphene (TLG) stacks in either a semimetallic ABA or a semiconducting ABC configuration with a gate-tunable band gap, but the latter has only been produced by exfoliation. Here we present a chemical vapor deposition approach to TLG growth that yields greatly enhanced fraction and size of ABC domains. The key insight is that substrate curvature can stabilize ABC domains. Controllable ABC yields ~59% were achieved by tailoring substrate curvature levels. ABC fractions remained high after transfer to device substrates, as confirmed by transport measurements revealing the expected tunable ABC band gap. Substrate topography engineering provides a path to large-scale synthesis of epitaxial ABC-TLG and other vdW materials.

V an der Waals (vdW) materials composed of individual atomic layers held together by vdW interactions have attracted considerable attention due to their unique physical properties that can be tuned by manipulating interlayer twisting angles [1][2][3] , creating heterostructures 4,5 , and controlling stacking configurations 6,7 . One example is trilayer graphene (TLG) 8 , which has two stacking configurations in its natural form. The ABA or Bernal configuration is a semi-metal with tunable band overlap 9,10 , while the ABC or rhombohedral configuration is a semiconductor with a gate-tunable band gap 7,11,12 . In fact, few-layer graphene with ABC stacking is an interesting system for all thicknesses, as it is predicted to have topologicallyprotected surface states [13][14][15] and high-temperature surface superconductivity 16,17 . Although ABA appears to be the ground state and is observed more often in exfoliated trilayer graphene, the thermodynamic preference over ABC is small, and a mixture is usually found. There are indications that the choice of substrate or the number of layers can affect this delicate ordering balance 18,19 .
To date, electronic transport studies of ABC-TLG have relied heavily on exfoliated samples, a process that is not scalable and limits the range of experimental investigations. A critical experimental challenge is to develop large-area synthesis of highquality ABC-TLG [20][21][22] , which would catalyze a wide range of studies of the unexplored physics and applications of this unique material system. To answer this need, we develop a high-yield chemical vapor deposition (CVD) process to grow epitaxial TLG on Ni-Cu gradient alloy substrates. Angle resolved photoemission spectroscopy (ARPES) is used to confirm the epitaxial growth process and the characteristic band structures of ABC-TLG and ABA-TLG. Transmission electron microscopy is used to assess the atomic structure of the material, and it reveals the existence of quasi-lamellar patterns of ABA-TLG and ABC-TLG. Based on these observations, we develop a model for how substrate topography and, more specifically, the curvature of surface corrugations, plays a principal role in stabilizing ABC-TLG during growth. This curvature-based stacking selection (CBSS) mechanism is confirmed by an experiment correlating CVDgrown TLG material with the original growth substrate. We then demonstrate that the substrate curvature can be engineered through CVD growth time for controllable ABC-TLG yields of 59%. We further find that transfer of TLG from the corrugated growth substrate to flat SiO 2 leads to annihilation of ABA-ABC domain walls and increased ABC and ABA domain size, as revealed by infrared scanning near-field optical microscopy (IR-SNOM). Electron transport measurements demonstrate that the band gap of the CVD-grown ABC-TLG could be tuned by a perpendicular electric field. Our finding that substrate topography can be used to stabilize the ABC stacking structure may enable development of related CBSS synthesis strategies for other vdW epitaxies with specific interlayer stacking configurations and study concerning their physical properties and applications.

Results
Growth and characterization of ABC-TLG. As shown schematically in Fig. 1a, TLG was synthesized on a Ni-Cu gradient alloy substrate in back-diffusion growth mode 23 . Catalytic substrates were prepared by depositing a 100-nm thick Ni film on one side of a flat Cu foil. The Ni-Cu foils were then annealed at 1050°C for 5 min to create a Ni-Cu gradient alloy substrate with a Ni-rich side and a Ni-poor side. During CVD growth, the top TLG layer forms first over the Ni-rich side of the substrate, and then the middle and bottom TLG layers grow beneath the top layer via carbon back-diffusion from the Ni-poor side 23 . Additional information about the growth process is provided in the Methods section. As shown in Fig. 1b, the process resulted in high-yield TLG growth, with an areal coverage of~30% over the entire substrate and flake sizes of 15-50 µm. All as-grown TLG domains displayed a compact hexagonal shape, attesting to their highly oriented, crystalline nature.
We performed systematic micro-ARPES and nano-ARPES measurements of the band dispersions of the material on the growth substrate, with 140 µm and 120 nm spot size respectively, to evaluate the electronic properties of the as-grown material. We directly probed the binding energy vs. momentum (E vs. k == ) curves close to the Fermi level; these are strikingly distinctive for TLG with ABA, ABC, and AAA stacking configurations 24 . The orientation, size, and distribution of graphene grains, as well as their registry with the substrate were investigated using micro-ARPES. Figure 1c shows a 3D ARPES plot, including the Fermi surface plane of the multi-grain graphene film, which discloses the complete electronic band structure of the graphene flakes and the metallic substrates. Strong π Dirac graphene bands are visible at the K high symmetry points of the reciprocal lattice. Moreover, well-defined d bands from Cu and/or Cu-Ni alloy are evident with binding energies of 2-4 eV. Distinctive d and sp substrate bands in the 3D ARPES plots indicate a high degree of crystallinity in the metallic substrate foils after CVD growth.
The graphene-substrate registry and the orientation of the graphene flakes are revealed by the multi-grain Fermi surface shown in Fig. 1d. A set of copper sp bands are observed dispersing in reciprocal space as expected for a copper single crystal (111) plane. Three well-defined Dirac cones are observed at k == ¼ 1.703 Å −1 , corresponding to the π graphene bands that cross the Fermi level. In contrast to earlier reports on few-layer graphene grown by CVD 24 , no arc-like shape states around the cones are observed, indicating that our graphene flakes are oriented along the ΓΚ graphene [ 110] direction of copper, in registry with the substrate lattice 25 . These results demonstrate that our process leads to epitaxial growth of TLG on the Ni-Cu gradient alloy substrate.
The ARPES plot of Fig. 2a was used to select those states exclusively belonging to the graphene flakes and then record the nano-ARPES image of Fig. 2b. This image displays the intensity variations of the Dirac states at the Fermi level (yellow rectangle in Fig. 2a) throughout a large portion of the sample and reveals the distribution, size, and morphology of the graphene flakes. All the flakes exhibit identically-oriented hexagonal shapes, some more intense than the others in Fig. 2b (indicated by arrows, the intensity line scan profile of the green line is shown in Fig. 2c). To precisely and unequivocally resolve the electronic structure of each flake, the electronic band dispersion close to the Fermi level (E vs. momentum) was recorded for individual graphene grains using high energy and momentum resolution nano-ARPES 26 . The most intense flakes in Fig. 2b show TLG electronic dispersion characteristics; the flake indicated by the black arrow has two bands with quadratic dispersion and one band with linear dispersion (Fig. 2f), as expected for ABA-TLG. The flake indicated by the blue arrow has a cubic band that is particularly flat at the Fermi level and two quadratic bands at similar energies that are shifted along k == . (Fig. 2g), as expected for ABC-TLG 27,28 . TLG flakes containing a mixture of both stacking types were also identified (Fig. 2h). All TLG flakes were n-doped, with the Dirac point located below the Fermi level due to the combined effects of the electrostatic surface potential and the "pillow effect" caused by surface adsorbates [29][30][31] . These nano-ARPES results are consistent with our density-functional theory (DFT) calculations ( Fig. 2i-j, Supplementary Figs. 8,9), which qualitatively confirm the observed band dispersions for the two TLG configurations, as well as the flat substrate bands around −1.5 eV (see Supplementary Fig. 4). Due to differences in the density of states around the Dirac point, the induced charge transfer from the substrate to graphene leads to a greater band shift in ABA-TLG as compared with ABC-TLG, in agreement with the data in Fig. 2f-g.
Dark-field transmission electron microscopy (DFTEM) was used to probe the stacking morphology and crystallinity of the asgrown TLG. This method can provide layer number contrast for oriented graphene 32 (see Methods section for details). Figure 2d shows a DFTEM image generated by selecting one of the six 1 210 f g. diffracted beams with the objective aperture. For oriented TLG, these diffracted beams from each layer constructively interfere, and the contrast between monolayer, bilayer, and trilayer graphene is evident in the image. The uniform appearance of the trilayer hexagon in the middle suggests that the TLG is crystalline and highly oriented. To obtain contrast between different stacking configurations, we instead selected one of the six 0 110 f g diffracted beams 33 , which with ABC stacking completely destructively interferes and with ABA stacking only partially interferes. The 0 110 f g DFTEM image (Fig. 2e) of the same region shown in Fig. 2d reveals a quasi-lamellar pattern of TLG regions with ABA and ABC stackings. Approximately 30% of the TLG flake is in the ABC configuration, twice that reported for exfoliated samples 34 and larger than usual levels in graphite. We observed flakes with up to 100% ABC content (see Supplementary  Fig. 2), as also observed in exfoliated samples 35 . AAA stacking was not observed, consistent with the ARPES results, while some twisted TLG was obtained (see Supplementary Fig. 3).
CBSS mechanism of ABC-TLG formation. The large-area fraction of ABC phase and the quasi-lamellar ABA/ABC domain morphology (Fig. 2e) suggest that the growth substrate properties and/or growth kinetics have an ABC-stabilizing effect. We performed ab initio calculations on flat, idealized substrates to first determine whether substrate chemistry contributes to this effect. We found that the relative energetics of multilayer graphene adsorbate/substrate systems are dominated by the interaction strength of the bottom-most graphene layer with the substrate, while differences due to variations in stacking of the upper layers are small. Energy differences computed between the most stable ABA and ABC arrangements on both Cu(111) and Ni-doped Cu(111) are therefore also small (~0.1 meV/atom) and barely within the accuracy of DFT (see Supplementary Note 1). We concluded that ABA and ABC are nearly isoenergetic on flat Cu(111) and Ni-doped Cu(111), in line with the coexistence of both structures in graphite and exfoliated TLG 34 . Importantly, this near equivalence between ABA and ABC energies suggests that the abundance of the two states can be tuned in finite-sized flakes by modulating other factors that slightly bias their relative energetics.
One such factor is substrate topography, which we discovered can be exploited to energetically stabilize the (typically) less prevalent ABC phase. Our growth substrates had a quasi-periodic pattern of surface corrugations (discussed further below) with typical curvatures of jκ 0 j % 10 7 m À1 at maxima and minima (Fig. 3a). When adjacent graphene layers retain atomic registry, substrate curvature induces an in-plane interlayer strain ε κ ðxÞ % κðxÞd % 0:5%, where d is the distance between graphene layers (Fig. 3b). Graphene must stretch (compress) in regions of negative (positive) curvature to maintain registry with adjacent layers. Above a critical curvature κ c , interlayer slip (i.e., interlayer dislocation formation) becomes energetically preferable to maintaining atomic registry between layers at all points. For weakly bound atomic layers of the type studied here 36 , κ c is $10 6 À 10 7 m À1 . Growth on a corrugated substrate with κ 0 > κ c will result in arrays of interlayer dislocation lines, several types of which induce stacking shifts between ABA and ABC (i.e., domain walls = curvature-induced interlayer dislocations, see Fig. 3c and Supplementary Note 2).  To quantify these effects and identify the factors that control interlayer dislocation 37 formation and domain size, we formulated a first principles-informed continuum mechanics description of the structure of vdW materials on non-flat surfaces (see Methods section for details). The dislocation character of a domain wall is denoted ξ À X u l Y, where ξ is the domain wall line direction (z = zigzag or a = armchair), X and Y specify the stacking configuration (B Bernal/ABA, R rhombohedral/ ABC) on the left and right side of the interface, respectively, and u and l specify the magnitude of b e , the edge component of the Burgers vector b (in units of a/2, where a = 0.14 nm is the atomic nearest neighbor distance) between the upper two and lower two graphene layers, respectively (see Fig. 3c and Supplementary Note 2 for further details of the notation).
The results demonstrated that curvature stabilizes interlayer dislocations and that certain types of stable dislocations in turn imply the presence of ABC domains. Specifically, our model indicates that the equilibrium state changes with increasing curvature from perfect ABA to states containing interlayer dislocations (Fig. 3e). A substrate with appropriate surface curvature (red-orange region in Fig. 3e) could fully stabilize arbitrarily large ABC domains containing only z À R 2 2 R interfaces. However, the experimental values of κ 0 span from ABA equilibrium (darker turquoise region) to z À B 4 4 B equilibrium (lighter turquoise region), such that some flakes may be perfect ABA, some ABC/ACB (z À R 2 2 R) or ABA/ACA (z À B 4 4 B), and some a mixture of ABA and ABC. The computed energies of the z À B 0 2 R (mixed ABA/ ABC), perfect ABC, and z À B 1 2 B (ABA/ACA) states are only slightly higher than those of the equilibrium states within the region of interest, indicating that these may also be observed. This is in agreement with the identification of pure ABA, pure ABC, and flakes with both stackings via nano-ARPES ( Fig. 2f-h) and DFTEM ( Fig. 2e and Supplementary Fig. 2). This stabilization of the ABC phase results from the difference in geometrically allowed partialdislocation content between ABA and ABC stacking variants, as detailed in Supplementary Note 2 and Supplementary Figs. 10-14. Principally, ABC/ACB configurations (z À R 2 2 R) accommodate 1D curvature strain more efficiently than ABA/ACA (z À B 1 2 B) due to their larger edge character (larger b e , see Fig. 3d). The ABC phase is therefore preferred over an intermediate range of curvatures despite its larger bulk energy.
Our description of CBSS indicates that rich internal stacking structures will be commonplace in CVD multilayer graphene systems ( Supplementary Fig. 25) and provides a framework for understanding how the interplay between corrugations and domain walls can stabilize different structures. The more prevalent views that assume internally structureless flakes (e.g.,  ABAB…-stacked multilayer graphene islands 38 ) and/or structural domains with only simple periodic stacking sequences (e.g., a mixture of ABAB… and ABCABC… 33 ) are not likely representative of actual states.
To directly test the CBSS mechanism, we marked a CVD-grown sample such that DFTEM images of domain wall morphologies could be mapped back to AFM images of the substrate topography at the precise location where growth occurred (see Supplementary Note 3 and Supplementary Fig. 19 for experimental details). Figure 3f shows that the dominant domain wall direction is precisely aligned with that of the substrate corrugations, and further that individual domain walls are strongly correlated with lines of maximum curvature (indicated by white arrows in Fig. 3f). Imperfect correlation is expected due to substrate corrugation irregularities and changes in domain wall configuration induced during cooling and transfer. Importantly, we found that substrate regions with lower than average corrugation curvatures (e.g., at grains with atypical surface crystal orientations) yielded significantly less ABC coverage than those with curvatures in or near the ABC/ACB-stabilizing range shown in Fig. 3e. This is seen clearly in the IR-SNOM stacking map (Supplementary Fig. 6) for a single TLG flake that spans three different Cu grains. The map shows 40% ABC coverage over one grain that has near-optimal curvature (~1-2 × 10 7 m −1 ) and only 7% ABC coverage over the other two grains which have significantly lower curvatures (~1-2 × 10 6 m −1 ). Further details of the IR-SNOM technique are provided below.
Previous studies have shown that ABC-TLG domain size can be manipulated using an electric field 39 or mechanical force 40 Fig. 3 Curvature-stabilized ABC-TLG. a AFM topography image showing the corrugated substrate surface, 1D topography profile with a Fourier series fit, and Gaussian fits to the corresponding curvature profiles at corrugation extrema. The AFM image of the substrate was taken after growth and transfer of the graphene film. The measured curvatures were within the resolution of our AFM (see Supplementary Fig. 5). b Schematic of TLG on curved Cu demonstrating how geometric curvature and interlayer-interactions lead to in-plane strain, ε κ . C i denotes graphene where i is the layer index. Substrate and graphene curve lengths are L sub ¼ R 0 θ and L i ¼ ½R 0 þ 4 À i ð Þd i θ, respectively, such that ε κ % ðL i À L iþ1 Þ=L i % d i =R 0 ¼ κd i . c Schematics of interlayer dislocation configurations in TLG. ? and ? denote full and partial dislocations, respectively. The notation is described in the text. d Top views of the z À B 1 2 B. and z À R during synthesis. Substrate corrugations of the type seen here, long observed on graphene growth substrates but poorly understood, have recently been shown to form and coarsen beneath graphene during growth 41 . We propose (and provide evidence) that the primary driving force for this phenomenon is the reduction of the total step/curvature-induced interlayer disregistry/dislocation energy (E disreg þ E disloc , see Methods section) as substrate step edges merge into increasingly coarse corrugations (see Supplementary Note 5 and Supplementary Fig. 22). The resulting corrugation evolution during growth feeds back into graphene stacking domain morphology in a co-evolving process. Incorporating this overlayer-driven corrugation kinetics into our continuum model, we obtain a time-dependent description of domain wall energetics during simultaneous corrugation coarsening and graphene growth (see Supplementary Note 5 and Fig. 3e).
These insights indicate that ABC yield can be controlled and optimized by engineering curvature through tailoring synthesis parameters (growth time, temperature, substrate chemistry, crystallography). Accordingly, increasing growth time from 1 to 6 h led to a monotonic increase in the final average peak corrugation curvatures (measured on our Cu (111) substrates) from~7.0 × 10 6 to~1.4 × 10 7 m −1 (Fig. 3g and Supplementary Fig. 23). ABC yield also increased monotonically from~22% to~59% between 1 and 3 h, before decreasing to~36% at 6 h ( Fig. 3g and Supplementary Fig. 24). The increase in ABC yield to a maximum at substrate curvatures 1.0-1.2 × 10 7 m −1 is in good agreement with the predictions of Fig. 3e, and the decrease at 6 h is consistent with our expectation of a finite range of optimal curvatures. These results substantiate our physical description of CBSS and demonstrate that the growth process can be designed to control and optimize ABC yield. More precise control and higher yields may be achievable with further exploration of the effects of growth temperature, substrate chemistry, and substrate crystallography (see Supplementary Note 5).
We performed additional experiments to confirm that curvaturestabilized ABC domains were preserved and exhibited appropriate electronic properties after transfer and conformation to flat dielectric substrates commonly used for device fabrication. TLG flakes were transferred to a SiO 2 substrate and imaged with infrared scanning near-field optical microscopy (IR-SNOM; Fig. 4a). In this measurement, infrared light (red beam in the figure) is focused onto a metal-coated AFM tip, and the backscattered light is collected by a detector in the far field. Local infrared conductivities are thus obtained across the TLG flake, generating contrast between ABA (green) and ABC (purple) regions 40,42 . IR-SNOM is uniquely suited for the task of analyzing the domain structure of the TLG synthesized here because its tip-enhanced signal offers~20 nm spatial resolution (diameter of the metallic AFM tip), greatly superior to the~1 µm typically achieved by Raman spectroscopy 34 . Figures 4b, d show TLG flake topographies as measured by AFM. The flake height is uniform across the hexagonal domains, except at the 3-7 nm high wrinkles visible as white lines. The corresponding IR-SNOM image, however, exhibits distinct contrast between ABA and ABC regions (Fig. 4c, e). Both ABA-TLG and ABC-TLG were observed via IR-SNOM after transfer to SiO 2 , and three striking phenomena were identified.
These phenomena can be understood in light of the physical picture described by our continuum model. First, the difference in ABC-TLG yield between samples with zigzag and armchair domain walls is consistent with the prediction of our continuum model that maximum ABC stabilization occurs at the zigzag direction (see Fig. 3e and Supplementary Fig. 10). Second, the growth substrate possesses a quasi-periodic array of peaks and valleys with nearly equal but opposite curvature (Fig. 3a), which favors the creation of TLG with a quasi-periodic array of domain walls with alternating sign. After removal of the domain wallstabilizing Ni-Cu substrate and transfer to flat SiO 2 , many adjacent interlayer dislocation pairs will be attracted toward each other, and those pairs whose Burgers vectors sum to zero will fully annihilate upon contact, as described in Supplementary  Figs. 15-18. The average ABA and ABC domain sizes are therefore expected to increase upon transfer to flat SiO 2 . Strikingly, we found that large-area CVD ABC-TLG can be reliably retained after transfer onto an hBN substrate (Supplementary Fig. 26). This stands in sharp contrast to the case for exfoliated graphene, where it is reported 43 that ABC-TLG is almost always converted to ABA-TLG after transfer onto hBN. Our CVD approach is therefore expected to have a significant impact on research directed at the Mott-insulator state 43 and tunable superconductivity 44

in TLG.
Transport measurement of ABC-TLG. The large ABC-TLG and ABA-TLG domain sizes on SiO 2 as identified by IR-SNOM allowed us to fabricate FET devices and perform transport measurements on TLG regions with different stacking configurations. Confirmation of their expected electronic behavior is critical for potential applications. As illustrated in the inset of Fig. 4f, FET devices were created in dual-gate configurations using 40 nm thick Al 2 O 3 dielectric layers deposited by atomic layer deposition as the top-gate dielectric and 250 nm thick thermally grown SiO 2 as the bottom-gate dielectric. We measured the resistance of TLG channels as a function of the top and bottom-gate voltages at a temperature of 1.8 K. Figure 4f shows resistance vs. top-gate voltage (R-V TG ) curves of an ABC-TLG channel. The existence of a tunable energy band gap in the material is evidenced by the increase in the on/off ratio as the strength of the out-of-plane electric field is increased. The same measurement showed that a sizeable ABC band gap was also achieved at room temperature (Supplementary Note 4 and Supplementary Fig. 21), further evidence of the high quality of our epitaxially-grown TLG. At 1.8 K, an on/off ratio of more than 3000 was observed, comparable to exfoliated samples 35,45,46 . No such increase was observed for the ABA FET (Fig. 4g), in good agreement with its predicted band structure and measurements by others on exfoliated samples 7,45 .

Discussion
We have demonstrated a scalable and controllable approach to CVD growth of TLG with an enhanced yield of large ABCstacked domains. The growth strategy was based on a Ni-Cu gradient alloy substrate that facilitates the back-diffusion growth mode. The key mechanism introduced and underlying the approach is curvature-based stacking selection (CBSS); the use of substrate curvature on the nanoscale to generate interlayer dislocations and locally stabilize ABC domains. The quasi-periodic corrugated topography of our growth substrate (commonly observed) typically generates alternating lamellar regions of ABC and ABA, as revealed by DFTEM and IR-SNOM measurements. An optimal range of curvatures, tunable through growth time, leads to maximized ABC yield, consistent with our continuum CBSS model predictions. The coarsening of substrate topography beneath growing graphene overlayers is attributed to the reduction of interlayer disregistry/dislocation energy achieved by increasing corrugation scale. ARPES results on as-grown TLG confirmed the epitaxial nature and electronic hallmark of ABC-TLG and ABA-TLG, while the existence of a sizable, gate-tunable band gap for ABC-TLG regions was confirmed by electron transport measurements. Significant increases in ABC and ABA domain size were observed following transfer from growth substrates to flat SiO 2 substrates, consistent with the expectation that interlayer dislocation pairs annihilate upon transfer to flat SiO 2 . CBSS, optimized through substrate topography engineering, is broadly applicable to epitaxial synthesis of vdW materials beyond ABC-TLG, and should enable further exploration of their promising physical properties and provide opportunities for applications in next-generation nanoelectronic devices.

Methods
TLG synthesis and transfer. The TLG growth was carried out in an atmospheric pressure CVD system with a 1-inch furnace (Lindberg Blue M, Thermo Scientific Co.). Cu foils of 25 µm thickness (Item #46365, Alfa Aesar) were cleaned with 5.4% HNO 3 for 40 s and two DI water baths for 2 min each. A nickel film of 100 nm thickness was sputtered onto one side of the Cu foil (Lesker PVD75 DC/RF Sputter System). Ni-Cu foils were cut into 1 cm × 4 cm pieces and then loaded into the CVD furnace. The furnace was ramped to 1050°C at a rate of 60°C/min in a flow of 500 sccm Ar without H 2 . The growth substrate was annealed for 5 min at 1050°C in a flow of 500 sccm Ar + 30 sccm H 2 . TLG was grown using 1.8 sccm CH 4 (1% in Ar) + 30 sccm H 2 for 3 h. The reactor was then rapidly cooled to room temperature in a flow of 10 sccm H 2 and 1000 sccm Ar.
PMMA (MicroChem Corp., PMMA950, A4) was spun over the surface of the TLG as grown on Ni-Cu foils and then baked at 100°C for 2 min The substrate was removed by an etcher solution (Transene Company, Inc. CE-100). The PMMA supported TLG was cleaned in 10% HCl solution and two water baths, and was transferred onto a 250 nm SiO 2 /Si substrate with prefabricated gold makers. After air drying and baking at 150°C for 3 min, PMMA was removed by soaking the sample in acetone overnight and cleaning with IPA. TLG on SiO 2 /Si substrates were annealed in H 2 /Ar forming gas at 400°C for 1 h before IR-SNOM characterization.
IR-SNOM. Scattering-type scanning near-field optical microscopy is based on tapping mode atomic force microscopy (AFM). A metallic AFM tip is illuminated with a focused CO 2 laser beam with wavelength 10.6 µm and the enhanced local field under the tip interacts strongly with the sample underneath it. The elastically ARTICLE scattered light containing the local optical properties of the sample substrate is collected by a MCT (HgCdTe) detector in the far field. Near-field optical images with spatial resolution of sub 20 nm can be achieved simultaneously with topography by recording the scattered light while scanning the sample. ABA/ABC graphene have different optical conductivities at 10.6 µm due to their different band structures and thus can be differentiated in the near-field images. In bilayer graphene, the contrast of domain walls in near-field images stems from surface plasmon reflection at the domain walls and feature a single or double bright line depending on the domain wall type 47 . Operation in tapping mode rather than contact mode 40 prevents imaging-induced domain wall motion, and no such motion was observed in these experiments.
AFM. An atomic force microscope (AFM, Icon Bruker) equipped with a probe with a tip radius of 10 nm (HQ-300-Au, Oxford Instruments) was used to characterize the topography of the growth substrate, in a different facility from that in which IR-SNOM measurements were performed. DFTEM. TLG was transferred onto a Cu grid with an amorphous carbon support. An aperture was placed in the diffraction plane of an electron microscope to create a real space image from electrons scattered only at a specified angle. With oriented stacked graphene, the 1 210 f gplanes of A, B, and C layers are aligned, and the diffracted beams from each layer interfere constructively. Selecting one of the six f1 210g diffraction spots with the objective aperture give layer number contrast for oriented graphene. In ABC-TLG, the first order 0 110 f gplanes of B and C layers are displaced by 1/3 and 2/3 of the plane spacing, respectively, from the same planes in the A layer. The 0 110 f gdiffracted beams from each layer therefore show complete destructive interference in ABC-TLG. In ABA-TLG there is only partial destructive interference of the 0 110 f gbeams 33 . The 0 110 f gdiffraction spots therefore provide the desired stacking contrast in oriented TLG. In BLG, when selecting a 0 110 f gdiffracted beam, a non-zero sample tilt can break the symmetry between the AB and BA layers 32 . This can be seen in Fig. 2e, where although the layer normal of the sample was aligned to the optic axis, the not perfectly flat TEM support leads to contrast between AB and BA stackings in some areas of the bilayer region.
ARPES. The ARPES and nano-ARPES experiments were carried out at the ANTARES beamline of synchrotron SOLEIL. It is equipped with a Fresnel zone plate (FZP) to focalize the beam and order selection aperture (OSA) to eliminate higher diffraction orders. The sample was mounted on a nano-positioning stage which was placed at the coincident focus point of the electron analyzer and the FZP. Nano-ARPES can be operated in two ways: (1) the point mode where a spectrum is collected from a specific point of the sample and (2) the imaging mode, where photoelectron intensity from a specific electron energy range is mapped over the sample to create a two dimensional image of the electronic states of interest. All photoemission measurements were performed at a temperature of ∼60 K. The photoelectron spectra were obtained using a hemispherical MBS analyzer. The instrumental energy resolution is ∼12 meV and practical momentum resolution is ∼0.02 1/Å for an emission angle ∼20°D FT calculations. All-electron density-functional theory (DFT) calculations were carried out with the Fritz-Haber-Institute-Ab initio-Molecular-Simulations 48 package using the Perdew-Burke-Ernzerhof (PBE) generalized-gradient approximation 49 supplemented by the TS method 50 to account for dispersive interactions. We use the suggested tight setup and accompanied basis functions throughout. A Cu (111) surface was generated from an optimized copper slab using six layers, where the bottom three layers were kept fixed to their bulk positions. A vacuum layer of 30 Å was introduced to separate periodic images along z, and a finite-size correction for the slab was employed. The Brillouin zone was sampled by a 11 × 11 × 1 Monkhorst-Pack k-point grid 51 , and k-point convergence was improved by filling the electronic states with a Gaussian broadening of 0.1 eV. Models including nickel doping of 3.7% were calculated in 3 × 3 surface supercells of the above-described copper surface, reducing the number of layers to three and replacing one copper atom with nickel. Binding energies were defined throughout as , and E 0 n À graphene ð Þ are the total energies of the combined Cu(111) plus n-graphene sheet system, the defined Cu(111) surface, and n-graphene sheets, respectively. Transition states were searched with the climbing-image nudged-elastic-band method using five images. 52 Continuum model for CBSS in vdW materials. The total energy of an elastic multilayer system (weakly bonded sheets) was divided into contributions from inplane disregistry strain, interlayer dislocations, and stacking phase bulk energies; For trilayer graphene on corrugated Cu (111), the substrate topography was assumed to vary only along x such that the disregistry energy over a region of size L Cu (per unit length y) is where ε i ¼ b i N i À ε 0;i is the elastic strain and ε κ;i x ð Þ ¼ d i =½κ À1 x ð Þ þ P 3 j¼iþ1 d j is the geometric strain. E 2D ¼ 340 N=m 53 is the 2D Young's modulus of graphene, i is the interlayer index (i = 1, 2, and 3 correspond to C 1 -C 2 , C 2 -C 3 , and C 3 -Cu, respectively; see Fig. 3b), b i is the dislocation Burgers vector magnitude along x (edge component), N i is the integer dislocation density between x = 0 and x = L Cu , ε 0,i is the in-plane misfit strain between flat layers in i, and d i = 0.335 nm 54 is the interlayer separation. We separate geometric strains from elastic strains in E disreg to emphasize the physical competition between these two effects. Geometric strain ε κ.i is that introduced via curvature under the condition of epitaxial registry between layers (as illustrated in Fig. 3b), and elastic strain ε i is that introduced by interlayer dislocations. Without dislocations, the entire geometric strain is stored as elastic energy within the system. With a particular amount of dislocations, the introduced elastic strain exactly offsets the geometric strain, and the strain energy goes to zero. With this form of E disreg , the strain reference state is that in which a given graphene layer rests strain-free on either the substrate (i = 3) or the graphene layer below it (i = 1 or 2). This form also invokes the assumption that dislocations relieve strain uniformly over each local peak or valley. The only non-material-parameter input into the model, the substrate curvature profile κ(x), was determined by fitting AFM-measured substrate topographies to a Fourier sine series and parameterizing the resulting curvatures at each corrugation extremum as κ x ð Þ ¼ κ 0 e Àx 2 =2σ 2 (Fig. 3a). We further employed ε 0;1 ¼ ε 0;2 ¼ 0, and neglected i = 3 (the substrate). Bending energy differences between states can also be neglected (the bending modulus of graphene is extremely small). The total dislocation line energy (per unit length y, adapted from ref. 37 ) was expressed as where a = 0.14 nm is the nearest neighbor distance in graphene, E edge = 0.318 × 10 −10 J/m and E screw = 1.091 × 10 −10 J/m are the core energy coefficients, θ i is the angle between the dislocation line and Burgers vector directions, and B = 22.08 × 10 −20 J is the bending modulus of graphene. The first term in the square brackets is the core energy and the second term is the contribution from long-range interactions between interlayer dislocations. The bulk energy of the trilayer is E phase ¼ E i L Cu P 3 j¼1 d j for single phase states, where E i is the bulk energy of ABA or ABC per unit volume, and E phase ¼ 1 2 ðE ABA þ E ABC ÞL Cu P 3 j¼1 d j for two phase states.
We considered ten possible TLG states; perfect ABA, perfect ABC, and eight configurations with periodic dislocation arrays (Fig. 3c). There are five configuration categories: ABA-ABA, ABA-ACA, ABC-ABC, ABC-ACB, and ABA-ABC. The existence of the latter three was experimentally confirmed (see Supplementary Fig. 1; ABA and ACA cannot be distinguished by DFTEM 32 ). The energy of each state was determined as a function of maximum local curvature κ 0 (at σ = 16 nm, the average measured feature width), as shown in Fig. 3e for interlayer dislocations along the zigzag direction (see Supplementary Fig. 10 for the armchair direction).