The epitaxy of 2D materials growth

Two dimensional (2D) materials consist of one to a few atomic layers, where the intra-layer atoms are chemically bonded and the atomic layers are weakly bonded. The high bonding anisotropicity in 2D materials make their growth on a substrate substantially different from the conventional thin film growth. Here, we proposed a general theoretical framework for the epitaxial growth of a 2D material on an arbitrary substrate. Our extensive density functional theory (DFT) calculations show that the propagating edge of a 2D material tends to align along a high symmetry direction of the substrate and, as a conclusion, the interplay between the symmetries of the 2D material and the substrate plays a critical role in the epitaxial growth of the 2D material. Based on our results, we have outlined that orientational uniformity of 2D material islands on a substrate can be realized only if the symmetry group of the substrate is a subgroup of that of the 2D material. Our predictions are in perfect agreement with most experimental observations on 2D materials’ growth on various substrates known up to now. We believe that this general guideline will lead to the large-scale synthesis of wafer-scale single crystals of various 2D materials in the near future.

In graphene CVD growth, the zigzag (ZZ) edge is generally the slowest propagating edge because of its highest barrier for edge propagation 23,24 and the alignment of graphene on a substrate has been broadly observed to be dependent on the symmetry of the substrate. For example, graphene islands grown on a Cu(111) or Cu(110) surface are well-aligned but those grown on a Cu(100) surface are observed to be along two perpendicular directions 5,11,12 . Among the three low-index Cu surfaces, highly robust alignment of graphene can be obtained on the Cu(111) surface and currently, the principal method of graphene single crystal production is by epitaxially growing graphene single crystals on Cu(111) surface [5][6][7] .
Because of the lower C 3V symmetry of hBN, the alignment of hBN on a substrate is different from that of graphene. Three of its six ZZ edges are nitrogen terminated and the other three are boron terminated (named as ZZN and ZZB edges hereafter). In most experiments, the ZZN edge has been proven to be the slowest propagating and kinetic Wulff construction leads to triangular hBN islands enclosed by three ZZN edges [25][26][27] . In contrast to epitaxial graphene growth, well-aligned hBN islands have rarely been observed. When grown on Cu(111) or Cu(110) surfaces, triangular hBN islands aligned along two opposite directions were found [13][14][15]28 , while those grown on Cu(100) surfaces had four different orientations 14 . Recent works have shown that well-aligned hBN islands can be successfully achieved by using a Cu substrate with tailored step edges, thus enabling epitaxial growth of wafer-scale hBN single crystals 8,9 . Similar to hBN, most transition metal dichalcogenides (TMDCs) possess three-fold symmetry and present very similar epitaxial behavior on substrates; for example, with two possible alignments on Au(111) 16 , Al 2 O 3 (0001) [17][18][19] , and GaN(0001) 20 surfaces. Well-aligned WS 2 islands have been grown on hBN surface 21 and nearly well-aligned WSe 2 islands have been grown on a vicinal Al 2 O 3 (0001) surface 22 . Most recently, centimeter scale single-crystalline MoS 2 was obtained by the coalescence of well-aligned MoS 2 grains on a vicinal Au(111) surface 10 .
All these experimental observations strongly indicate that the alignment of a 2D material on a substrate depends on both its symmetry and that of the substrate and a general theory for 2D materials epitaxy that helps to predict the alignment of various 2D materials on different substrates is highly desirable to serve as a guideline for experimental design.
Here, based on extensive density functional theory (DFT) calculations, we present a general theory to explain how the alignment of a 2D material on a substrate is intimately related to its interaction with the substrate and how the epitaxial growth of the 2D material on a substrate is critically dependent on the interplay of the symmetries of the 2D material and the substrate. Our theory explains most known experimental observations on 2D material epitaxial growth and hence can serve as a guideline for the experimental synthesis of various 2D single crystals, as well as 2D polycrystalline materials with designed grain boundaries.

Results
2D material-substrate interaction and the alignment of a 2D material on a substrate. There are hundreds of important 2D materials and the possible substrate types are also of the same order of magnitude. So, it is impossible to calculate the interactions of all possible combinations of 2D materials and substrates, which is greater than 100,000. Without losing the generality, we can classify the interactions between 2D materials and various substrates into two sceneries: (i) The edge of the 2D material is terminated by the substrate, such as graphene or hBN on an active metal substrate, where the strong interaction between the edge of the 2D material and the pristine substrate facet determines the alignment of the 2D material and its epitaxial growth behavior 29,30 ; (ii) The edge of the 2D material is self-passivated or terminated by active atoms from the environment of its growth, such as H or OH groups [31][32][33][34] , where the weak interaction between the bulk of the 2D material and the pristine substrate facet dominates the alignment of the 2D material.
To establish an epitaxial relationship between a 2D material and a substrate for scenery (i), we firstly explore the interaction between the edges of graphene or hBN with the three low-index Cu surfaces, namely Cu(111), Cu(100), and Cu(110). The calculated binding energies between a graphene (hBN) ZZ (ZZN) edge on the three low-index Cu surfaces as a function of the angle of edge alignment that is defined as the angle between the edge and a Cu〈110〉 direction of the substrate, are shown in Fig. 1 (please refer to method section and Supplementary Figs. 1-6 in Section 1 of supplementary information (SI) for more details on the calculation and modeling). We clearly see that on each of the three low-index surfaces, the strongest binding energy appears when the ZZ(ZZN) edge of graphene (hBN) is along a Cu〈110〉 direction of the surface; the difference between the binding energy minimum and maximum is significant, >0.2 eV per edge atom. Hence, on a Cu surface, a well-aligned small graphene or hBN island of~2 nm (which has only~200 atoms of which~40 are at the edge) has an energy advantage of >8 eV over misaligned ones. This binding energy difference is large enough to maintain a growing graphene or hBN island in a well-aligned configuration on a Cu surface.
To elucidate the reason behind the strongest binding of a ZZ (ZZN) edge along the 〈110〉 direction on a Cu surface, we have plotted the electron density distributions of Cu(111), Cu(100), and Cu(110) surfaces, respectively, in Fig. 2a-c. It can be seen that on all the three surfaces, the isosurface fluctuation in electron density is the lowest along the 〈110〉 direction, indicating that the close-packed 〈110〉 atomic rows form a pattern with alternative ridges and valleys of uniform height on the surface. In a straight ZZ(ZZN) edge of graphene(hBN), the less stable edge atoms form a straight line and this straight edge is preferentially passivated by either a ridge or a valley of the Cu surface instead of crossing over ridges and valleys on the surface, which results in distortion of the edge. To further illustrate preferential passivating of the graphene ZZ edge by a close-packed atomic row, we compare the atomic structures of the interfaces and the charge density differences of the graphene edge along both 〈110〉 and other directions on the three low-index Cu surfaces ( Fig. 2d-i). When the graphene ZZ edge is aligned along the Cu〈110〉 direction, all the edge atoms are well passivated by a Cu〈110〉 atomic row and the edge remains straight. In contrast, if the graphene edge is along another direction, some of the edge atoms are poorly passivated and the edge is no longer straight because of the fluctuating ridge-valley pattern of the surface. The above analysis clearly shows the superiority of the close packed direction of a substrate in passivating a high-symmetric edge in a 2D material.
Our analysis of a straight graphene/hBN edge preferring the direction of valley or ridge of the isosurface of Cu substrate can be generally applied to most combinations of 2D materials on various substrates. It is noted that the lattice constant of graphene/hBN ZZ edge matches that of the 〈110〉 direction of Cu substrate well. Thus, it is worth to consider a system without perfect lattice-match. To address the effect of lattice-mismatch, we consider the interaction between a ZZ edge of graphene and the Pt(111) surface, where the lattice constant of graphene ZZ edge is about 12.6% smaller than that of Pt〈110〉 direction. The binding energies of a graphene ZZ edge on the Pt(111) surface as a function of the alignment angle of the graphene edge are shown in Supplementary Fig. 7. As expected and similar to that of a graphene ZZ edge on Cu(111) surface, the graphene ZZ edge prefers the alignment along a 〈110〉 direction of the Pt (111) surface. As an example, above calculation suggests that the conclusion that the slowest propagating (also highsymmetric) edge of a 2D material prefers to align along the high-symmetric direction of an active metal substrate, regardless of the lattice-match between the 2D material and the substrate.
If the edge of a 2D material is self-terminated, such as the selenium-terminated edges of TMDC materials 33 , or terminated by the H or OH functional groups, such as the edges of graphene or hBN grown on a less active metal substrate 31,32,34 . In such a scenery, the interaction between the 2D material edge and the substrate is no longer very strong and the dominating interaction is the weak interaction between the bulk of the 2D material and the substrate. To elucidate the alignment of a 2D material on a substrate for scenery (ii), we calculate the interaction between hBN and Cu(111) (Fig. 3a, b, and Supplementary Fig. 8) and Au(111) (Supplementary Fig. 9) surfaces, respectively, by using the periodic boundary condition models and the interaction between a triangular WS 2 cluster and the hBN surface (Fig. 3c, d and Supplementary Fig. 10). All these calculations, together with previous studies of graphene on Cu(111) surface 35 and MoS 2 on Al 2 O 3 (0001) surface 36 , show that a high-symmetric direction of a 2D material (ZZ directions of graphene, hBN and TMDCs) prefers to align along the highsymmetric directions of a substrate, such as the 〈110〉 directions of Cu(111) and Au(111) surfaces, and 〈11 20〉 direction of hBN and Al 2 O 3 (0001) surfaces.
The above results allow us to draw a conclusion of the alignment of a 2D material on an arbitrary surface, i.e., a highsymmetric direction of the 2D island prefers to align along a high-symmetric direction of the substrate. Although we just explored a very limited systems of 2D materials grown on various substrates, as will be seen later, this rule is in accordance with most experimental observations on the epitaxial growth of various 2D materials, such as the CVD synthesis of graphene, hBN, and TMDCs on various transition metals or nonmetallic substrates (please refer to Section 2 in SI). Thus, we believe that this rule can be applied for the epitaxial growth of various 2D materials on different substrates.  material on a substrate, we proceed to discuss the interplay between the symmetry of a 2D material and that of the substrate in epitaxial growth. Let us consider a 2D material with a G 2D symmetry group on a substrate with a symmetry group of G sub . The symmetry group of the whole system, G 2D@Sub , must be a subgroup of either G 2D or G Sub because any symmetry operation of G 2D@Sub will not change the alignment of the 2D material or the substrate. As shown in SI, we have proved that the number of equivalent but different directions of a 2D material on a substrate can be calculated by where G sub j jand G 2D@sub j jare the orders, or numbers of different symmetry operations, of G sub and G 2D@sub , respectively.
According to the principle of 2D materials alignment discussed in the previous paragraph, the presence of a highsymmetry edge of a 2D material along a high-symmetry direction of the substrate ensures that the symmetry group of the whole system, G 2D@Sub , is the largest subgroup of both G sub and G 2D . We have considered various combinations of the symmetries of the 2D material and the substrate and the numbers of equivalent but different alignments of various 2D materials are shown in Table 1.
Without loss of generality, we have used fcc(111), fcc(100), fcc (110), and hBN(0001) surfaces as different examples of a substrate with 6-, 4-, 2-and 3-fold symmetries to illustrate the alignment of various 2D materials on them. Figure 4 presents the various ways in which 2D materials with 2-, 3-, 4-and 6-fold symmetries are aligned on these substrates. From the figure, we can deduce that in order to keep the whole system with the highest symmetry, there are: 1, 2, 1 and 1 equivalent but different alignments for a 6-fold symmetric 2D material on 6-, 4-, 2-and 3-fold symmetric substrates; 2, 4, 2, and 1 equivalent but different alignments for a 3-fold symmetric 2D material on 6-, 4-, 2-and 3-fold symmetric substrates; 3, 1, 1 and 3 equivalent but different alignments for a 4-fold symmetric 2D material on 6-, 4-, 2-and 3-fold symmetric substrates; 3, 2, 1 and 3 equivalent but different alignments for a 2-fold Table 1 The number of equivalent but different orientations of a 2D material on a substrate based on the interplay between their symmetries.  1  2  1  1  3  1  3  1  2  4  1  2  3  2  3  1 Here, G Sub , G 2D and G 2D@Sub are symmetry groups of the substrate, the 2D material and the 2D material-substrate system, |G sub | and |G 2D@sub | are the orders of G Sub and G 2D@Sub ; and N 1 = |G sub |/ |G 2D@sub | is the number of equivalent but different directions of the 2D material on the substrate.

BN-(0001)
[1120] 0°30°6 0°0°3  Table 1. Besides the number of equivalent but different alignments of a 2D material on a substrate, Fig. 4 also gives the misalignment angles of equivalent islands of 2D materials. On the substrates with 6-, 4-, 3-, 2-fold symmetries, the misalignment angles are i 3 π; i 2 π; 2i 3 π; iπ; i ¼ 1; 2; 3;… respectively. It is important to note that a high-symmetric edge of a 2D material along a high-symmetric direction of a substrate is critical for the above analysis. If the high-symmetric edge of the 2D material is along a low-symmetric direction of the substrate with mirror symmetry, the symmetry group G 2D@Sub will have no mirror symmetry and the least number of equivalent but different alignments of the 2D material will be 2, which makes orientational uniformity impossible.
Comparison with experimental observations. We have summarized most of known experimental obversions on 2D materials epitaxial growth, including graphene growth on low-index Cu surfaces (Supplementary Table 1), hBN growth on low-index Cu surfaces (Supplementary Table 2), TMDCs growth on various low-index substrates, including Al 2 O 3 , Au, GaN, hBN (Supplementary Table 3), and various 2D material grown on different high-index substrates (Supplementary Table 4). It is interesting to note that there is a perfect agreement between the experimentally observed numbers of alignments and misalignment angles of 2D materials on these substrates with those predicted by our theoretical analysis. Hence, we believe that the epitaxial relationship is valid for the epitaxial growth of most 2D materials grown on different substrates.
Strategies toward the epitaxial growth of various 2D single crystals. From Eq. (1), we can see that once the symmetry group of the substrate, G Sub , is a subgroup of the 2D material, G 2D , or in other words the symmetry group of the system G 2D@Sub could be same as G Sub . This means that there is only one most preferential alignment of the 2D material on the substrate and the orientational uniformity of a large number of islands of the 2D material is possible. As shown in Table 1 and Fig. 4, this can be applied to 2D materials with 6-fold symmetry, such as graphene on a 6-, 3-or 2-fold symmetric substrate, 2D materials with 4-fold symmetric on a 4-or 2-fold symmetric substrate, 3-fold symmetricy 2D materials on a substrate with 3-fold symmetry, and a 2-fold symmetric 2D material on a 2-fold symmetric substrate. Up to now, the seamless stitching of well-aligned graphene islands have been realized on Cu(111) and Ge(110) surfaces, both of which are in accordance with the above described analysis [5][6][7]37,38 . Among the thousands of known 2D materials, most of them have the 3-fold symmetry, such as the most explored hBN and TMDCs. As shown above, a substrate with 3-fold symmetry is expected to be suitable for the epitaxial growth of 2D materials with a three-fold symmetry, but it is difficult to find proper lowindex substrates with 3-fold symmetry. Although some low-index substrates has the 3-fold symmetry, such as the Cu(111) surface and the Al 2 O 3 (0001) surface, the atoms of the top layer of the substrate generally have a higher symmetry, such as the atoms of the top Cu(111) layer have the 6-fold symmetry. So, for a C 3V 2D material on a C 3V substrate, there are generally two deep local minima in the formation profile which corresponds two high symmetric configurations, such as the 0˚and 60˚configurations of a WS 2 on the hBN surface shown in Fig. 3d and Supplementary  Fig. 10, and the 0˚and 60˚configurations of MoS 2 on the Al 2 O 3 (0001) surface in Fig. 3h of ref. 34 (please refer to Supplementary  Fig. 11 for the configuration difference). Experimentally, antiparallel TMDC islands on 3-fold symmetric substrates are generally seen and parallel aligned TMDC islands could be realized by precise control of the experimental condition (Supplementary Table 3) [17][18][19][20][21] . From Eq. (2), we can see that if the substrate has the C 1 symmetry group, the condition for epitaxial growth can be satisfied for any 2D materials, implying that on a substrate with no symmetry, we may be able to achieve orientational uniformity for any 2D materials. In practice, substrates with very low symmetry could be a high-index surface or a vicinal surface of a low-index surface. As illustrated in Fig. 5, a high-index surface has a large number of low-index terraces connected by parallel step edges. These step edges of the substrate can serve as nucleation sites to initiate the growth of the 2D material. Furthermore, these step edges interact preferentially with an edge of the 2D material to promote the orientational uniformity of the 2D material. In this manner, the orientational uniformity of various 2D materials has been widely observed. As listed in Supplementary Table 4, epitaxial growth of well-aligned hBN islands have been observed on Cu(102), Cu(103), and vicinal Cu(110) surfaces 39,40 , where one of the three edges of the triangular hBN island is aligned along the step edge of the substrate. In addition, well-aligned WSe 2 islands were also observed on Al 2 O 3 (0001) surface with step edges 22 . Recently, such a strategy has been used to grow wafer scale single- crystalline hBN on vicinal Cu(110) surface and Cu(111) surface with step egdes 8,9 , and centimeter scale single-crystalline MoS 2 on an Au(111) surface with aligned step edges 10 . DFT calculations in these studies have shown that the weak interaction between the bulk of the 2D material and the substrate and/or the strong interaction between the edge of 2D material and step edge of the substrate lead to the favorable alignment of the 2D islands along the step edges of the substrate [8][9][10]39,40 . Since high-index surfaces can be easily obtained by miscutting a single crystal, we believe that this could be a general strategy for the synthesis of various 2D materials in the future.

Discussion
We would like to note that the fact that a high symmetric direction of a 2D material prefers to align along a high symmetric direction of a substrate, which is revealed by our extensive DFT calculations in this study, is the foundation of our main conclusion. Otherwise, orientational uniformity of a 2D material on a substrate, even the symmetry group of the substrate is a subgroup of that of the 2D material, cannot be realized. Currently, most previous studies on the synthesis of large-sized single crystalline hBN and MoS 2 employed vicinal (111) or (110) surfaces with parallel step edges [8][9][10] , where 2D materials tend to align along these step edges and the general study on the epitaxial 2D materials growth on various high index surface is very rare. Besides the vicinal surfaces which are close to one of the lowindex surfaces, our study also predicts that the high-index surfaces that are largely deviated from all the low-index surfaces, such as the (123) surface of an fcc material, are also ideal for the epitaxial growth of large-scale single-crystalline 2D materials. Furthermore, our theory provides a principle to determine the alignment of a 2D material on any given substrates and, thus, it offers a strategy of synthesizing 2D materials with well-defined grain boundaries, for instance, polycrystalline graphene with grain boundaries of a 30 o misorientation angle can be synthesized on an fcc(100) surface, similar to the case of graphene growth on a liquid Cu surface 41 .
In conclusion, our study clearly demonstrates that the interplay of the symmetries of the 2D material and the substrate is critical for the epitaxial growth of 2D materials. Both theoretical analysis and experimental observations show that a high-symmetric direction of a 2D material tends to be aligned along a high-symmetric direction of a substrate, so that the 2D material-substrate system has the highest possible symmetry. Based on the symmetry analysis and the rules for the preferential alignment of a 2D material on a substrate, we established a library of the different possible alignments of various 2D materials on different substrates to serve as a guideline for experimental design. Furthermore, we theoretically proved that the epitaxial growth of a 2D single crystal can be realized only if the symmetry group of the substrate is a subgroup of that of the 2D material. To meet the requirement for single-crystalline 2D material growth, we propose using substrates with high-index surfaces which have lower symmetry to template the epitaxial growth of various 2D materials; this strategy has been successfully demonstrated (Nature 570, 91 (2019); Nature 579, 219 (2020); ACS nano 14, 5036 (2020)) and is in agreement with experimental observations. After submission of this manuscript, we have noticed that same strategy has been employed for the epitaxial synthesis of singlecrystalline nanoribbons of TMDCs (Nat. Mater. DOI 10.1038/ s41563-020-0795-4 (2020)), which further validated the proposed approach of synthesizing large area single-crystal 2D materials. Our study thus provides a theoretical foundation for the synthesis of wafer-scale 2D single crystals.

Methods
DFT calculations of edge binding energies. All DFT calculations were carried out via the Vienna ab initio simulation Package (VASP) 42,43 . The exchange-correlation effect was treated by the Perdew-Burke-Ernzerhof generalized gradient approximation (GGA) 44 . The interaction between valence electrons and ion cores was described by the projected augmented wave (PAW) method 45 and the k-point mesh was sampled by a separation of 0.03 Å −1 .
To calculate the binding energy of a graphene ZZ edge to Cu (111), (100), and (110) surfaces, Cu slabs consisting of three (111), (100), or (110) atomic layers were constructed to mimic the Cu substrates. Graphene nanoribbons along the ZZ direction, which were three hexagons wide and one of the two ZZ edges passivated by hydrogen, were constructed. Because of the incommensurate lattice constants between the graphene ZZ nanoribbon and the Cu substrates, only a small number of periodic structures can be constructed with a graphene ZZ nanoribbon adsorbed on low-index Cu surfaces along different directions, of which the number of atoms is not too large and can be handled by DFT calculations, as shown in Supplementary Figs. 1-3. The initial distance between the graphene ZZ nanoribbon and the Cu substrate is set to 3.06 Å, which was estimated from DFT-D2 calculations 46 , and it is the typical equilibrium distance between a graphene layer and the underlying Cu substrate surface. Structure optimization was conducted with the atomic positions of the lowest Cu atomic layer fixed. In addition, the vertical positions of carbon atoms that are passivated by hydrogen were also fixed during structure optimization. To eliminate the imaginary interaction between periodic images along the vertical direction, a vacuum layer with a 15 Å thickness was used to separate the Cu slabs. All the structures were relaxed until the force on each unfixed atom was <0.01 eV/Å, with an energy convergence of 10 −4 eV.
In a similar manner, hBN nanoribbons along the ZZ direction were also constructed with their ZZB edges passivated by hydrogen. The structures with hBN ZZ ribbons adsorbed on low-index Cu substrate surfaces are shown in Supplementary Figs. 4 and 6. The distance between hBN nanoribbon and the Cu substrate is set to be 3.10 Å, which was obtained by optimizing a hBN sheet on Cu(111) surface by DFT-D2 calculations 46 . During structure optimization, the vertical positions of boron atoms passivated by hydrogen were fixed. In addition, the atomic positions of the lowest Cu atomic layer were also fixed.
The method similar to above calculations was also employed to calculate the binding energies between the graphene ZZ edge and Pt(111) surface. The calculated structures are shown in Supplementary Fig. 7.
The binding energy of a ZZ edge of graphene or hBN to the substrate is defined as where E ribbon , E sub , and E total are, respectively, the energies of the nanoribbon, the substrate and the total energy of the nanoribbon adsorbed on the substrate; L is the edge length of the nanoribbon.
To calculate the weak interaction between hBN wall and Cu(111) or Au(111) surface, a hBN layer was stacked to a Cu(111) or Au(111) slab consisting of three atomic layers under periodic boundary condition and with different alignment angles. The calculated structures are provided in Supplementary Figs. 8 and 9. During structure optimization, the bottom atomic layer of the metal slab was fixed. The binding energy between the hBN wall and the metal substrate is defined as where E total , E hBN , and E sub are the energies of the whole system, the hBN layer and the substrate, respectively. N BN is the number of BN pairs of the hBN layer in the unit cell of the whole system. To calculate the weak interaction between WS 2 and hBN layer, a triangular WS 2 cluster consisting of 60 S atoms and 27 W atoms was placed on a hBN layer with different orientations (Supplementary Fig. 10). Because the edges of the WS 2 cluster are passivated by S and the hBN layer is chemically inert, the interaction between the WS 2 cluster and the hBN layer should be dominated by the WS 2 wall and the hBN layer. The binding energy between the WS 2 cluster and the hBN layer is defined as where E total , E hBN , and E WS2 are the energies of the whole system, the hBN layer and the WS 2 cluster, respectively. N W is the number of W atoms in the WS 2 cluster.

Data availability
The data that support the findings of this study are available from the corresponding authors on reasonable request.