Unusual role of epilayer–substrate interactions in determining orientational relations in van der Waals epitaxy Lei Liua,1, David A. Siegelb,1, Wei Chenc,d,1, Peizhi Liua,e, Junjie Guoe,f, Gerd Duschere,f, Chong Zhaog,h, Hao Wangg,h, Wenlong Wangh, Xuedong Baih, Kevin F. McCartyb, Zhenyu Zhangc,2, and Gong Gua,2 a Department of Electrical Engineering and Computer Science, University of Tennessee, Knoxville, TN 37996; bMaterials Physics Department, Sandia National Laboratories, Livermore, CA 94550; cInternational Center for Quantum Design of Functional Materials, Hefei National Laboratory for Physical Sciences at the Microscale, and Synergetic Innovation Center of Quantum Information and Quantum Physics, University of Science and Technology of China, Hefei 230026, China; dDepartment of Physics, The University of Tennessee, Knoxville, TN 37996; eDepartment of Materials Science and Engineering, University of Tennessee, Knoxville, TN 37996; fMaterials Science and Technology Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831; gInternational Center for Quantum Materials, School of Physics, Peking University, Beijing 100871, China; and hBeijing National Laboratory for Condensed Matter Physics, Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China

Using selected-area low-energy electron diffraction analysis, we showed strict orientational alignment of monolayer hexagonal boron nitride (h-BN) crystallites with Cu(100) surface lattices of Cu foil substrates during atmospheric pressure chemical vapor deposition. In sharp contrast, the graphene–Cu(100) system is well-known to assume a wide range of rotations despite graphene’s crystallographic similarity to h-BN. Our density functional theory calculations uncovered the origin of this surprising difference: The crystallite orientation is determined during nucleation by interactions between the cluster’s edges and the substrate. Unlike the weaker B– and N–Cu interactions, strong C–Cu interactions rearrange surface Cu atoms, resulting in the aligned geometry not being a distinct minimum in total energy. The discovery made in this specific case runs counter to the conventional wisdom that strong epilayer–substrate interactions enhance orientational alignment in epitaxy and sheds light on the factors that determine orientational relation in van der Waals epitaxy of 2D materials.

|

|

two-dimensional materials van der Waals epitaxy hexagonal boron nitride graphene orientational relation

|

|

R

esearch on van der Waals (vdW) heterostructures formed by stacking up various 2D crystals is an emerging field (1) because the numerous possible interfaces may lead to new physics not necessarily associated with the constituent 2D materials and, therefore, to novel applications. Most pioneering works on vdW heterostructures are based on mechanical placement of exfoliated 2D crystal flakes (2–4), where uncertainty in the orientational relation between layers is inevitable. For some purposes, e.g., to enhance graphene’s transport by placing it on hexagonal boron nitride (3) (h-BN, or simply BN hereafter), the orientational relation is not important. In other cases, however, a definitive orientational relation is necessary for certain physical properties to arise. Examples include the tunneling behavior of a graphene–insulator–graphene junction (5) and the band structure of graphene placed on top of BN (2). The only practical way to achieve certainty in orientational relation, or azimuthal order, in vdW heterostructures is epitaxial growth. Actually, vdW epitaxy has been studied for decades as a method to overcome lattice mismatch (6, 7) enabled by the relatively weak and flexible vdW interactions between the epilayer and the substrate. Although vdW epitaxy [or quasiepitaxy (8)] is typically incommensurate, there can be a well-defined orientational relation (8), which has been observed in many examples of vdW epitaxy (6, 7, 9–11). In some cases, however, a definitive epilayer–substrate orientational relation is lacking (12, 13). At the moment when vdW epitaxy is shaping into a new thrust area in 2D crystal research, it is imperative to uncover the mechanisms that determine the orientational relation. www.pnas.org/cgi/doi/10.1073/pnas.1405613111

Epilayer–substrate interactions lead to registry in commensurate epitaxy, which naturally guarantees orientational alignment between the epilayer and the substrate. In vdW epitaxy, one would intuitively expect that strong epilayer–substrate interactions enhance azimuthal order, and this trend has been observed in some experiments (14–16). Here we show by a case study that the conventional wisdom does not always apply and that the relationship between the azimuthal order and the epilayer–substrate interactions in vdW epitaxy is highly case-specific. This surprising finding sheds light on the determination of orientational relation in vdW epitaxy. One would also intuitively expect crystallography and symmetry similarities between the epilayer and substrate to be important, if not determining, factors to orientational alignment. Graphene has been grown with orientational alignment on crystallographically similar BN (10, 11). Lattice constant mismatched transition metal chalcogenides are orientationally aligned to each other in vdW epitaxial heterostructures (7), presumably because they share a threefold symmetry. In this work, we show that threefold symmetric BN exhibits definitive orientational alignments when grown on the fourfold symmetric Cu(100) surface. In stark contrast, the graphene–Cu(100) epitaxy exhibits a distribution of rotations (17–19). This difference occurs despite the close Significance van der Waals (vdW) heterostructures of dissimilar 2D material sheets held together by vdW interactions promise new physics and applications not necessarily associated with the constituent 2D materials. Only epitaxial growth can achieve the orientational alignment in vdW heterostructures needed to obtain certain novel phenomena. As a case study of vdW epitaxy, we experimentally find that hexagonal boron nitride strictly aligns to Cu(100), whereas crystallographically similar graphene is known to exhibit a wide spread of in-plane rotations. Theoretical investigation reveals that this stark difference occurs because the C–Cu interactions are stronger than the B–Cu and N–Cu interactions. This counterintuitive discovery sheds light on orientational relationships in vdW epitaxy and their case specificity. Author contributions: L.L. and G.G. designed research; L.L., D.A.S., and K.F.M. performed experiments; W.C. performed computations; L.L., D.A.S., W.C., P.L., J.G., G.D., C.Z., H.W., W.W., X.B., and K.F.M. contributed new reagents/analytic tools; L.L., D.A.S., W.C., K.F.M., Z.Z., and G.G. analyzed data; and L.L., D.A.S., W.C., K.F.M., Z.Z., and G.G. wrote the paper. The authors declare no conflict of interest. This article is a PNAS Direct Submission. 1

L.L., D.A.S., and W.C. contributed equally to this work.

2

To whom correspondence may be addressed. Email: [email protected] or [email protected].

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10. 1073/pnas.1405613111/-/DCSupplemental.

PNAS Early Edition | 1 of 6

APPLIED PHYSICAL SCIENCES

Edited by Paul L. McEuen, Cornell University, Ithaca, NY, and approved October 22, 2014 (received for review April 6, 2014)

crystallographic similarity between graphene and BN. First-principles calculations reveal that the difference between the two systems arises from epilayer–substrate interactions. Results Monolayer BN Single-Crystal Growth. Methods describes the growth by atmospheric-pressure chemical vapor deposition (APCVD) (20) on cold-rolled copper foils. Except for the precursor, the conditions were the same as we previously used to grow graphene crystallites (19). The wide thermal decomposition window and the violent decomposition behavior of the precursor ammonia borane (H3B−NH3) hinder the synthesis of monolayer BN crystallites with well-defined, energetically favored edges (21, 22). We leverage the diffusionlimited kinetics of APCVD to achieve a low feedstock arrival rate, which was further controlled by the precursor charge. Analogous to graphene APCVD using highly diluted methane (23, 24), the method results in monolayer single crystals with energetically favored edges, commonly believed to be nitrogen-terminated zigzag edges, and therefore exhibiting the distinctive equilateral triangle shape (22, 25). The scanning electron microscopy (SEM) image in Fig. 1A shows isolated, equilateral triangle-shaped BN domains grown with a low precursor charge. With increased precursor charge, crystallites coalesce into islands of complex shapes, such as butterflies and stars (Fig. 1B and SI Appendix, Fig. S1). Arrows overlaid on the BN domains in Fig. 1A suggest that these triangles are oriented only in several directions, a result that warrants an in-depth investigation. In addition, samples were characterized by X-ray photoemission spectroscopy and UV-visible spectroscopy (SI Appendix, Figs. S2 and S3) for chemical analysis and optical band-gap measurement, respectively. On the atomic scale, the element-contrast (Z-contrast) annular dark-field scanning transmission electron microscopy (ADF STEM) image and the corresponding intensity line profile (Fig. 1 C and D) show clear distinction between individual B and N atoms, which, considering the stacking of bulk hexagonal boron nitride (26), unambiguously confirm that the crystallites are monolayer.

Monolayer BN–Cu(100) Superstructure. To understand the limited number of allowed orientations of the BN crystallites suggested by Fig. 1A, low-energy electron microscopy (LEEM) and selectedarea low-energy electron diffraction (μ-LEED) were performed to determine the relative crystallography of the BN domains on the Cu foil substrate. Every diffraction pattern acquired across millimeter length scales of the surface indicated a Cu(100) surface termination, consistent with the fact that the cold-rolled Cu foil surface after thermal annealing consists nearly exclusively of large, (100)-oriented grains (17, 27, 28). Fig. 2A displays a representative μ-LEED pattern of a BN island, showing diffraction spots corresponding to Cu(100), BN, and a moiré formed between the BN and the underlying Cu(100). Because this pattern is rather complicated, we illustrate it schematically in Fig. 2B. The four first-order diffraction spots of the Cu(100) surface lattice are marked by brown arrows in Fig. 2A and are colored brown in Fig. 2B. Similarly, the six first-order diffraction spots corresponding to the BN overlayer are marked in blue. The moiré spots are circled and colored in green, and we note that many of the circles in Fig. 2A contain two or three diffraction spots; these are represented by single green dots in Fig. 2B. The BN island exhibits only one set of diffraction spots, which confirms its monocrystallinity. At the selected electron energy, the pattern in Fig. 2A shows a clear threefold symmetry: higher and lower intensity BN spots are marked by thicker and thinner blue arrows, respectively. To account for the moiré pattern, we first note that the overall diffraction pattern has two high-symmetry directions: a horizontal direction with many closely spaced diffraction peaks from the superstructure (moiré) and a vertical direction where the separation between rows of closely spaced superstructure peaks is the Cu (100) reciprocal lattice constant. In the vertical direction, the vertical components of the BN and Cu(100) reciprocal lattices roughly align, i.e., aCu* ∼ aBN*sin60°, where aCu* and aBN* are the reciprocal lattice constants of the BN and the Cu(100) 2D lattices, respectively, leading to each Cu spot being collinear with a pair of first-order BN spots (blue). The close match between the BN and nm) and Cu(100) real-space lattice constants (aCu ∼ aBN ∼ pffiffiffi0.25 ffi the relationships aCu* = 2π/aCu and aBN* = 4π/( 3 aBN) explain this coincidence. In the horizontal direction, the ratio of the Cu(100) to BN reciprocal lattice constants is very close to 5:6, which leads to a supercell that in real space is 5× the Cu unit cell or 6× the BN unit cell. In reciprocal space this match creates four moiré spots (marked by green circles) equally spaced at 1/5 the separation between the (00) spot and the first-order Cu spot, or five diffraction spots (including the first-order Cu spot) spaced at 1/6 the separation between the (00) beam and the first-order BN spot. An atomic model (Fig. 2C) is constructed based on these coincidences observed in the diffraction pattern: aCu ∼ aBN (within +2.1%) along the vertical direction, and 5aCu ∼ 6aBNsin60° (within –1.8%) in the horizontal direction. We point out that this model is an illustration of the orientational relation and periodicity; any lateral translation is allowable. As a consequence of small lattice mismatches (+2.1% and –1.8% in the two directions), each moiré spot marked by a green circle in Fig. 2A actually comprises three closely spaced spots. Orientational Alignment Between BN Crystallites and the Cu(100) Lattice. Unlike the graphene–Cu(100) system (17–19), μ-LEED

Fig. 1. Monolayer BN crystals obtained by APCVD. (A and B) Representative SEM images of BN crystals on Cu foils with different film coverages. Spatially isolated equilateral triangles (A) and complex structures (B) are achieved by controlling the amount of precursor. Red arrows in A indicate orientations of crystallites. (C) Z-contrast ADF-STEM image of a BN crystallite, showing singlelayer character. Pink and blue circles denote B and N atoms, respectively. (D) Intensity profile along the green line in C, identifying B and N atoms.

2 of 6 | www.pnas.org/cgi/doi/10.1073/pnas.1405613111

analysis on multiple BN crystallites surprisingly shows that only four equivalent orientations occur for BN on Cu(100), i.e., BN crystallites are well aligned to the underlying Cu(100) lattice. Fig. 3A shows a bright-field LEEM image, where BN islands are imaged as bright triangular regions. The entire field of view is a single Cu(100) grain. Fig. 3 C–F displays μ-LEED patterns acquired on four BN islands that are circled and labeled in Fig. 3A. These four diffraction patterns correspond to four orientations of the BN on the Cu(100) surface. In each pattern, one BN reciprocal primitive vector aligns with one reciprocal primitive vector of the Cu(100) surface. Thus, by symmetry there are four equivalent orientations of BN crystals on Cu (100). Furthermore, dark-field imaging (SI Appendix, Fig. S5) of the same region imaged in Fig. 3 shows that four and only four orientations account for all of the BN islands in Liu et al.

the entire field of view. Fig. 3B schematically illustrates the four equivalent BN orientations in real space. In addition, another set of LEEM–LEED data (SI Appendix, Fig. S4), with only the first-order Cu and BN diffraction spots visible, more simply documents the epitaxial relationships revealed in Fig. 3. Moreover, a better-resolved bright-field image (SI Appendix, Fig. S4A), which reveals the equilateral triangle shape of the BN crystallites, together with the diffraction patterns allows us to confirm that the edge orientations of the BN triangles are indeed zigzag. (In Fig. 3A the crystallites are incompletely imaged, possibly due to the nonplanarity of the Cu foil.) With the four allowed orientations accurately identified by LEEM and μ-LEED, SEM provides statistics over large areas. Fig. 4A shows a representative SEM image of BN islands grown on an as-received Cu foil. The orientations of triangular crystallites are carefully examined and analyzed, exploiting the distinctive equilateral triangle shape and the fact that their sides are zigzag edges. Four and only four possible orientations are observed, in good agreement with the LEEM–LEED measurement, as well as an earlier phenomenological observation (28) that the triangular domains tended to have one orientation. The uneven distribution between the orientations shown in the histogram (Fig. 4C) may be due to the vicinal surface of the Cu grain. Interestingly, on electrochemically polished foils (Fig. 4 B and D) one orientation dominates, perhaps because the polishing resulted in a surface whose steps and terraces were preferentially aligned along one direction. Oxygen Impurity Induced Misalignment. To test the sensitivity of the BN–Cu(100) orientational relation to the BN–Cu(100) interactions, we grew BN on oxygen-contaminated Cu foils without changing the recipe. Whereas surface oxygen has been associated with intriguing effects on graphene island morphology (29, 30), here we focus on its alteration to the cluster–substrate interactions. The LEED pattern (Fig. 5A) acquired from a region free of BN has four diffraction spots that are not present for a simple Cu(100) surface. These additional spots are attributed to the c(2 × 2) reconstruction of Cu(100) caused by adsorbed oxygen (31). We then placed the electron beam spot within individual BN islands to perform μ-LEED. Fig. 5B displays a diffraction pattern obtained in one island, showing the Cu, BN, and c(2 × 2) reconstruction spots; some of the BN–Cu(100) moiré spots present in Figs. 2 and 3 are faintly visible. More important, new moiré spots (circled in cyan) appear, corresponding to double diffraction from the c(2 × 2) and BN lattices. This moiré formed by the BN and the c(2 × 2) provides strong evidence that the oxygen impurity is underneath the BN island. Due to the oxygen, the BN crystallites no longer lock into the four aligned orientations on Cu(100). The large-area LEED pattern in Fig. 5C, along with the orientation statistics (SI Appendix, Fig. S6), shows multiple rotation angles on individual Liu et al.

Cu(100) grains. This wide angular distribution, along with the weak BN–Cu moiré spots, indicates that the intercalation of oxygen atoms between the BN crystallite and the Cu substrate alters the interaction and therefore the alignment between the two lattices, providing strong evidence that the BN–Cu interaction is essential to achieve the observed alignment of BN on pristine Cu(100). Density Functional Theory Calculations. To understand why BN rigorously aligns on Cu(100) whereas graphene does not, we performed a comparative density functional theory (DFT) study to calculate the total energy vs. rotation angle of BN and graphene clusters on Cu(100) (Fig. 6), which simulate the initial nucleation stage (32, 33). We used the Vienna ab initio simulation package (34) with the projector augmented wave (PAW) method (35, 36) and the Perdew-Burke-Ernzerhof parameterization of the generalized gradient approximation (PBE-GGA) (37), as well as DFTD2, a semiempirical approach that includes vdW interactions (38, 39). The investigated structures were first relaxed using PBEGGA functionals without vdW interactions to search for stable or metastable geometries around initial translational positions and

Fig. 3. Alignment between BN crystallites and Cu(100) surface lattice. (A) LEEM image (25-V, 15-μm-diameter field of view) of BN islands (bright) on a single Cu(100) grain. White circles denote regions where μ-LEED patterns (C–F) were obtained. (B) Real-space atomic model of four equivalent orientations of triangular BN crystallites on Cu(100). (C–F) μ-LEED patterns acquired at locations marked in A. Blue arrows mark the reciprocal primitive vectors of BN. Brown arrows in C mark the first-order diffraction spots of Cu(100). The BN crystals exhibit threefold symmetry in the diffraction pattern at the chosen electron energy.

PNAS Early Edition | 3 of 6

APPLIED PHYSICAL SCIENCES

Fig. 2. BN–Cu(100) superstructure. (A) μ-LEED pattern (63 V) obtained on a triangular BN crystallite on a Cu(100) grain. First-order diffraction spots of Cu are marked by brown arrows, and the brighter and dimmer BN spots are marked by thicker and thinner blue arrows, respectively. Reciprocal primitive vectors of BN and Cu are drawn as blue and brown dashed arrows, respectively. Selected moiré spots are marked by green circles. (B) Schematic of the diffraction spots arising from BN and Cu. (C) Atomic model of BN on Cu(100) derived from the diffraction pattern. A superstructure cell, marked by green arrows, is formed between BN and Cu(100) lattices. The two lattices match well at spacings of 5 Cu cells or 6 BN cells. Lateral registry shown is arbitrary. Primitive vectors of the BN and Cu lattices are also shown.

Fig. 4. Statistics of BN–Cu(100) rotational orientations. (A and B) Representative SEM images of BN islands on (A) an unpolished Cu foil and (B) an electrochemically polished Cu foil. Four possible orientations of BN crystals are denoted A1, A2, B1, and B2 shown by legends. (C and D) Orientation histograms of triangular BN crystallites within individual Cu grains. Cu grain boundaries are clearly identifiable in SEM, as delineated by dashed lines in A and B.

rotational orientations, followed by DFT-D2 calculations for a second-step relaxation. The relaxed BN cluster exhibits a domed structure, similar to graphene clusters on Cu(111) and Ir(111) (32, 40), indicating that the cluster interacts with the substrate predominantly at the periphery. Due to the three- and fourfold symmetries of BN and Cu(100), respectively, a rotation angle θ = Θ is equivalent to θ = Θ + 30°. Meanwhile, due to reflection symmetry, the total energy Etotal (Θ) = Etotal (−Θ). Therefore, we only need to consider θ ∈ [0°, 15°], and thus set the initial rotational angles to be θ = 0°, 5°, 10°, and 15°. Fig. 6A shows two of the stable or metastable geometries after relaxation by DFT-D2, θ = 0° and 6.2°. Fig. 6B displays the total energy vs. rotation angle plot for all of the stable–metastable configurations discovered by relaxation using PBE and DFT-D2. The total energy at θ = 0° was set to zero as the reference for both methods, which resulted in the same trend, suggesting that the effect of the net corrections due to vdW interactions between the cluster and the Cu substrate is quite close for all orientations and therefore plays a minimal role in determining the relative stability. One exception is that the θ = 9.1° geometry obtained by PBE relaxed into θ = 6.2° for DFTD2, probably due to vdW interactions between edge atoms and

Cu atoms underneath. The calculations support the experimentally observed alignment between the BN and Cu(100) lattices. The alignment of one edge of the cluster with the substrate lattice (Fig. 6A) substantially enhances the cluster–substrate interaction, and any rotation away from this high-symmetry configuration will make the cluster less stable. For comparison, we calculated the total energy vs. rotation angle for a graphene cluster with the same geometry. (Whereas graphene clusters should be equiangular hexagons, our clusters can be considered as six-sided with three edges of the shortest possible zigzag geometry.) The results are summarized in Fig. 6 C and D. Within the PBE, a few orientations near θ = 6.5° have nearly the same total energies as θ = 0°. By adding vdW interactions, some of these configurations even become more stable than θ = 0°. Moreover, other metastable orientations between θ = 0° and 6.5° are very likely to exist, although the present DFT study did not find any due to the stringent criteria used for force convergence in ionic relaxation. Therefore, in contrast with the strict alignment of BN on the Cu(100), the graphene cluster can rotate away from the highsymmetry orientation θ = 0°, consistent with the previous reported experiments (17), where the diffraction intensity is high for θ ∈ [0°, ∼6.5°] although nonzero elsewhere. Discussion The calculated narrow rotation angle window of low total energy for BN–Cu(100) and a much wider one for graphene–Cu(100) are consistent with experiments. Nevertheless, the DFT-obtained data points in Fig. 6 per se do not explain why BN–Cu(100) strictly aligns whereas graphene–Cu(100) exhibits nonzero diffraction intensity at all angles, because these data points are total energy “local minima” in structural relaxation. The question how stable these local minima (i.e., metastable configurations) are in the two systems might be answered by exploring the rotation barriers between these local minima at a significantly higher computational cost than in the present work. Instead, we show below that the knowledge of the exact energy landscape between these local minima is not necessary to gain insights into the alignment mechanism. By inspecting the Cu atoms in the top-view models of Fig. 6 A and C, we see that the graphene cluster induces more Cu atom rearrangement than does the BN cluster, making the Cu(100) substrate more “reconstructed” at the cluster periphery. (SI Appendix, Fig. S9 provides a more graphical illustration of the difference in substrate distortion.) At a growth temperature close to the melting point of bulk Cu, these Cu atoms are pulled by a rotated cluster from their original positions in the Cu lattice to maximize cluster-edge– substrate binding and therefore minimize the cluster–substrate system total energy. A graphene zigzag edge, if away from a Cu direction, induces appreciable rearrangement of surface Cu atoms, due to strong C–Cu interactions. This action blurs the

Fig. 5. Oxygen impurity induced misalignment. (A) Selected-area LEED pattern of a BN-free region of an oxygen-contaminated Cu foil, showing Cu(100) (1 × 1) and c(2 × 2) spots, distinguished by brown and red circles, respectively. (B) μ-LEED pattern from a single BN island. The blue circles indicate the first-order diffraction spots from BN. Moiré diffraction spots (circled in cyan) due to multiple scattering by a c(2 × 2) and BN reciprocal lattice vectors (as exemplified by yellow arrows) prove the presence of oxygen impurity underneath the BN. BN–Cu(100) moiré spots are circled in green. (C ) LEED pattern with multiple sets of BN diffraction spots, highlighted by blue circles, which reveal the misalignment of BN grains on the O–Cu(100) surface.

4 of 6 | www.pnas.org/cgi/doi/10.1073/pnas.1405613111

Liu et al.

distinctiveness of a high-symmetry, i.e., aligned, configuration as a total energy minimum. For BN–Cu(100), on the other hand, the surface Cu atom rearrangement is negligible, thus the high-symmetry orientation θ = 0° remains a distinctive total energy minimum. The origin of this difference is that the graphene zigzag edge–Cu interaction is much stronger than the BN zigzag edge–Cu interaction, as revealed by the present DFT calculations using the same method as in previous work (41). To compare the strengths of graphene–Cu and BN–Cu interactions, we define a cluster’s adsorption energy as the total energy of the isolated cluster and the bare Cu substrate minus that of the combined system. The adsorption energy of the graphene cluster is higher than that of the BN cluster by 4.19 eV, calculated by DFT-D2. This difference can be traced to the higher adsorption energy of C compared with B and N monomers. The carbon monomer adsorption energy on Cu(100) is larger than that of B and N by 0.71 eV and 1.65 eV, respectively. The vdW epilayer orientation is determined soon after a cluster nucleates (32, 33), when the edge atoms, with dangling bonds that interact strongly with the substrate, steer the cluster into the lowest energy orientation. As the crystallite grows, the vdW interaction between the cluster’s interior atoms and the substrate dominates over the edge–substrate interactions due to the increased area–perimeter ratio, and the orientation is nailed down. Indeed, we have shown that the energetic preference of BN toward alignment with the underlying Cu(100) is overridden when BN growth is templated by a fresh edge of a preexisting graphene seed crystallite that is not aligned to the Cu lattice (19). The choice of the BN– and graphene–Cu(100) systems for this study emphasizes the effect of interactions between the edge atoms and the metal substrate, which are stronger than, say, between cluster edges and a covalent 2D material substrate. van der Waals epitaxy is usually incommensurate, lacking registry, even for epilayer–substrate pairs with shared symmetry. Orientational alignment minimizes the total energy by maximizing epilayer– substrate interaction. Stronger cluster–substrate interactions in the case of graphene–Cu(100), however, cause substrate surface atom displacements. This obscures the aligned orientation as a distinctive total energy minimum and, therefore, reduces orientational order. The fact that BN and graphene have different Liu et al.

symmetries than Cu(100) signifies that the factors that control orientation go beyond merely shared symmetry between the epilayer and substrate. Interestingly, for the present case study comparing graphene and BN on Cu(100), the notion that stronger epilayer– substrate interactions lead to reduced orientational alignment runs counter to the conventional wisdom discussed earlier. Due to the complicated interplay between the cluster-edge– substrate interactions and other interactions such as those within the substrate, the orientational relation in vdW epitaxy is case-specific and does not always follow a rule of thumb. Contrary to the present case, graphene strictly aligns to Cu(111) (18) but BN does not (42). On the face-centered cubic Co(111) or the hexagonal close-packed Co(0001) surface, each crystallographically similar to Cu(111), however, BN exhibits strict alignment (43). In summary, we find that nuclei can be oriented by the nonvdW bonding between their edges and the substrate even in epilayer–substrate systems where vdW forces dominate the interplanar bonding. Thus, edge–substrate interactions must be generally considered when evaluating the vdW epitaxy of 2D crystals. We nonetheless point out that energetics is not the only factor. In early work by molecular beam epitaxy, alignment was observed in all experimented systems (6, 7, 9). In recent experiments where alignment is lacking (12, 13), kinetics might have led to misalignment for a relatively flat energy-rotation angle landscape that only slightly prefers alignment. Methods Synthesis and Transfer. BN crystallites were grown by APCVD on 25-μm-thick Cu foils (Alfa Aesar, 46986), which were sequentially washed by dilute nitric acid and deionized water, dried by a N2 gun, and immediately loaded into the furnace. For growth on oxygen-contaminated foils, the foils were kept in air for more than 1 d before loading. Before the growth, H3N−BH3 was loaded in a boat located in the upstream region of the tube. The furnace was gradually heated to 1,050 °C over 60 min in Ar:H2 (400 and 100 standard cubic centimeters, respectively), followed by annealing at the same temperature and gas flow rates for 10 min. To start the growth procedure, BN precursor was heated to sublime at ∼120 °C by a heating tape wrapped around the tube. The typical growth time was 10–30 min, followed by rapid cooling.

PNAS Early Edition | 5 of 6

APPLIED PHYSICAL SCIENCES

Fig. 6. DFT calculated stable and metastable BN–Cu (100) and graphene–Cu(100) geometries. (A and C) Top and side views of BN (A) and graphene (C) clusters on Cu(100) with rotation angles 0.0° and 6.2°. Graphene is shown to induce more distortion of Cu lattice than BN (SI Appendix, Fig. S9). (B and D) Total energy vs. rotation angle plots for BN (B) and graphene (D) clusters on Cu(100) calculated by PAWPBE (vdW not included; black) and DFT-D2 (vdW included; blue). Dashed line segments connecting data points for visual guidance signify that data points represent discrete stable or metastable configurations; kinetic pathways between these data points remain to be explored.

The electrochemical polishing treatment was accomplished in a homemade electrochemical cell using orthophosphoric acid electrolyte and a Cu cathode. A 2.0–2.5-V dc voltage was applied for 10–20 min. After polishing, the foils were washed by deionized water and dried by N2. BN crystallites were transferred onto transmission electron microscope (TEM) grids or other substrates for characterization. To facilitate the transfer, a thin layer of poly(methyl methacrylate) (PMMA, 2% in anisole) was spin coated on the Cu foil sample (4,000 rpm for 1 min), and the copper was etched away by floating the foil in a copper etchant (CE-100, Transene). The PMMA–BN membrane was transferred onto the surface of a 10% HCl solution to remove residual metal particles and then washed by deionized water several times. Finally, the film was scooped out by the desired substrate (Quantifoil TEM grid with 1.2-μm-diameter holes from Ted Pella, SiO2–Si, or optical quartz plate). The PMMA layer was removed by acetone vapor, followed by thermal annealing (350 °C for 2 h in an Ar:H 2 forming gas).

vacuum region of more than 11 Å was used to ensure decoupling between neighboring slabs. During structural relaxation, the atoms in the bottom two layers were fixed in their respective bulk positions, with all of the other atoms fully relaxed until the force on any given atom was smaller than 0.03 eV/Å. A 2 × 2 × 1 k-point mesh was used for a 7 × 7 Cu surface unit cell. The effect of spin polarization was examined. As a sanity check, we first placed a two-ring BN cluster on Cu(100). The PBE approach led to an unphysical picture, where the cluster stands up on the surface. After the vdW interaction was included, the physically reasonable picture was arrived at (SI Appendix, Fig. S7). Due to the dominance of mostly non-vdW edge–substrate interactions for clusters as depicted in Fig. 6, the accuracy of vdW modeling does not affect our conclusions despite its general importance. Clusters smaller than the ones in Fig. 6 are also examined, and the physical picture of the aligned BN–Cu (100) and rotated graphene on Cu(100) remains (SI Appendix, Fig. S8).

DFT Calculations. The lattice constant of Cu was obtained via structural optimization. The generic Cu(100) surface was modeled by a slab of four atomic layers, and the BN or graphene clusters were placed on top of the metal surfaces. A

ACKNOWLEDGMENTS. The experimental work was partially supported by National Science Foundation (ECCS-1231808) and Defense Advanced Research Projects Agency. Work at Sandia was supported by the Office of Basic Energy Sciences, Division of Materials and Engineering Sciences, US Department of Energy (DOE) under Contract DE-AC04-94AL85000. This research used resources of Oak Ridge National Laboratory’s Center for Nanophase Materials Sciences, which is sponsored by the Scientific User Facilities Division, Office of Basic Energy Sciences, US DOE. The theory work was partially supported by Natural Science Foundation of China (Grants 11034006 and 11204286) and National Key Basic Research Program of China (Grant 2014CB921103). The calculations were performed at National Energy Research Scientific Computing Center (NERSC) of the US DOE.

1. Geim AK, Grigorieva IV (2013) Van der Waals heterostructures. Nature 499(7459): 419–425. 2. Yankowitz M, et al. (2012) Emergence of superlattice Dirac points in graphene on hexagonal boron nitride. Nat Phys 8(5):382–386. 3. Dean C, et al. (2012) Graphene based heterostructures. Solid State Commun 152(15): 1275–1282. 4. Britnell L, et al. (2013) Resonant tunnelling and negative differential conductance in graphene transistors. Nat Commun 4:1794. 5. Feenstra RM, Jena D, Gu G (2012) Single-particle tunneling in doped grapheneinsulator-graphene junctions. J Appl Phys 111(4):043711. 6. Koma A, Saiki K, Sato Y (1989) Heteroepitaxy of a two-dimensional material on a three-dimensional material. Appl Surf Sci 41-42(2):451–456. 7. Koma A (1992) New epitaxial growth method for modulated structures using van der Waals interactions. Surf Sci 267(1-3):29–33. 8. Forrest SR (1997) Ultrathin organic films grown by organic molecular beam deposition and related techiques. Chem Rev 97(6):1793–1896. 9. Koma A (1999) Van der Waals epitaxy for highly lattice-mismatched systems. J Crystal Growth 201:236–241. 10. Yang W, et al. (2013) Epitaxial growth of single-domain graphene on hexagonal boron nitride. Nat Mater 12(9):792–797. 11. Tang S, et al. (2013) Precisely aligned graphene grown on hexagonal boron nitride by catalyst free chemical vapor deposition. Sci Rep 3:2666. 12. Shi Y, et al. (2012) van der Waals epitaxy of MoS₂ layers using graphene as growth templates. Nano Lett 12(6):2784–2791. 13. Hwang J, et al. (2013) van der Waals epitaxial growth of graphene on sapphire by chemical vapor deposition without a metal catalyst. ACS Nano 7(1):385–395. 14. Preobrajenski AB, Ng ML, Vinogradov AS, Mårtensson N (2008) Controlling graphene corrugation on lattice-mismatched substrates. Phys Rev B 78(7):073401. 15. Sutter P, Sadowski JT, Sutter E (2009) Graphene on Pt(111): Growth and substrate interaction. Phys Rev B 80(24):245411. 16. Emtsev KV, Speck F, Seyller T, Ley L, Riley D (2008) Interaction, growth, and ordering of epitaxial graphene on SiC{0001} surfaces: A comparative photoelectron spectroscopy study. Phys Rev B 77(15):155303. 17. Wofford JM, Nie S, McCarty KF, Bartelt NC, Dubon OD (2010) Graphene Islands on Cu foils: The interplay between shape, orientation, and defects. Nano Lett 10(12): 4890–4896. 18. Ogawa Y, et al. (2012) Domain structure and boundary in single-layer graphene grown on Cu(111) and Cu(100) films. J Phys Chem Lett 3(2):219–226. 19. Liu L, et al. (2014) Heteroepitaxial growth of two-dimensional hexagonal boron nitride templated by graphene edges. Science 343(6167):163–167. 20. Tay RY, et al. (2014) A systematic study of the atmospheric pressure growth of largearea hexagonal crystalline boron nitride film. J Mater Chem C 2(9):1650–1657. 21. Song L, et al. (2010) Large scale growth and characterization of atomic hexagonal boron nitride layers. Nano Lett 10(8):3209–3215. 22. Kim KK, et al. (2012) Synthesis of monolayer hexagonal boron nitride on Cu foil using chemical vapor deposition. Nano Lett 12(1):161–166. 23. Yu Q, et al. (2011) Control and characterization of individual grains and grain boundaries in graphene grown by chemical vapour deposition. Nat Mater 10(6):443–449.

24. Vlassiouk I, et al. (2011) Role of hydrogen in chemical vapor deposition growth of large single-crystal graphene. ACS Nano 5(7):6069–6076. 25. Liu Y, Bhowmick S, Yakobson BI (2011) BN white graphene with “colorful” edges: The energies and morphology. Nano Lett 11(8):3113–3116. 26. Geick R, Perry CH (1966) Normal modes in hexagonal boron nitride. Phys Rev 146(2): 543–547. 27. Cho J, et al. (2011) Atomic-scale investigation of graphene grown on Cu foil and the effects of thermal annealing. ACS Nano 5(5):3607–3613. 28. Guo N, et al. (2012) Controllable growth of triangular hexagonal boron nitride domains on copper foils by an improved low-pressure chemical vapor deposition method. Nanotechnology 23(41):415605. 29. Hao Y, et al. (2013) The role of surface oxygen in the growth of large single-crystal graphene on copper. Science 342(6159):720–723. 30. Zhou H, et al. (2013) Chemical vapour deposition growth of large single crystals of monolayer and bilayer graphene. Nat Commun 4:2096. 31. Harrison MJ, et al. (2006) Adsorbate-induced surface reconstruction and surface-stress changes in Cu(100)/O: Experiment and theory. Phys Rev B 74(16):165402. 32. Chen W, et al. (2012) Suppression of grain boundaries in graphene growth on superstructured Mn-Cu(111) surface. Phys Rev Lett 109(26):265507. 33. Zhang X, Xu Z, Hui L, Xin J, Ding F (2012) How the orientation of graphene is determined during chemical vapor deposition growth. J Phys Chem Lett 3(19):2822–2827. 34. Kresse G, Furthmüller J (1996) Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys Rev B Condens Matter 54(16): 11169–11186. 35. Blöchl PE (1994) Projector augmented-wave method. Phys Rev B Condens Matter 50(24):17953–17979. 36. Kresse G, Joubert D (1999) From ultrasoft pseudopotentials to the projector augmented-wave method. Phys Rev B 59(3):1758–1775. 37. Perdew JP, Burke K, Ernzerhof M (1996) Generalized gradient approximation made simple. Phys Rev Lett 77(18):3865–3868. 38. Grimme S (2006) Semiempirical GGA-type density functional constructed with a longrange dispersion correction. J Comput Chem 27(15):1787–1799. 39. Buˇcko T, Hafner J, Lebègue S, Ángyán JG (2010) Improved description of the structure of molecular and layered crystals: Ab initio DFT calculations with van der Waals corrections. J Phys Chem A 114(43):11814–11824. 40. Lacovig P, et al. (2009) Growth of dome-shaped carbon nanoislands on Ir(111): The intermediate between carbidic clusters and quasi-free-standing graphene. Phys Rev Lett 103(16):166101. 41. Chen H, Zhu W, Zhang Z (2010) Contrasting behavior of carbon nucleation in the initial stages of graphene epitaxial growth on stepped metal surfaces. Phys Rev Lett 104(18):186101. 42. Joshi S, et al. (2012) Boron nitride on Cu(111): An electronically corrugated monolayer. Nano Lett 12(11):5821–5828. 43. Orofeo CM, Suzuki S, Kageshima H, Hibino H (2013) Growth and low-energy electron microscopy characterization of monolayer hexagonal boron nitride on epitaxial cobalt. Nano Res 6(5):335.

LEEM and LEED Characterizations. Copper foils with APCVD-grown BN domains were transferred through air into an Elmitec LEEM III instrument, where they were outgassed in ultrahigh vacuum at about 300 °C. Bright-field LEEM images were formed from the specularly reflected (00) beam. MicroLEED patterns were obtained from areas either 0.5 or 2 μm in diameter.

6 of 6 | www.pnas.org/cgi/doi/10.1073/pnas.1405613111

Liu et al.