Chinese Journal of Chemical Physics  2019, Vol. 32 Issue (6): 693-700

#### The article information

Lei-lei Li, Shuo Feng

Influence of Neighboring Layers on Interfacial Energy of Adjacent Layers

Chinese Journal of Chemical Physics, 2019, 32(6): 693-700

http://dx.doi.org/10.1063/1674-0068/cjcp1812291

### Article history

Accepted on: February 13, 2019
Influence of Neighboring Layers on Interfacial Energy of Adjacent Layers
Lei-lei Lia , Shuo Fengb
Dated: Received on December 25, 2018; Accepted on February 13, 2019
a. CAS Key Laboratory of Mechanical Behavior and Design of Materials, Department of Modern Mechanics, University of Science and Technology of China, Hefei 230027, China;
b. Hefei National Laboratory for Physical Sciences at the Microscale, Collaborative Innovation Center of Chemistry for Energy Materials, CAS Center for Excellence in Nanoscience, School of Chemistry and Materials Science, University of Science and Technology of China, Hefei 230026, China
Abstract: The binding energy and generalized stacking-fault energy (GSFE) are two critical interface properties of two dimensional layered materials, and it is still unclear how neighboring layers affect the interface energy of adjacent layers. Here, we investigate the effect of neighboring layers by comparing the differences of binding energy and GSFE between trilayer heterostructures (graphene/graphene/graphene, graphene/graphene/boron nitride, boron nitride/graphene/boron nitride) and bilayer heterostructures (graphene/graphene, graphene/boron nitride) using density functional theory. The binding energy of the adjacent layers changes from -2.3% to 22.55% due to the effect of neighboring layer, with a very small change of the interlayer distance. Neighboring layers also make a change from -2% to 10% change the GSFE, depending on the property of the interface between adjacent layers. In addition, a new simple expression is proven to describe the GSFE landscape of graphene-like structure with high accuracy.
Key words: Generalized stacking-fault energy    Binding energy    Neighboring layers    Adjacent layers    Density functional theory
Ⅰ. INTRODUCTION

Van der Waals (vdW) heterostructures are ideal platforms for exploring exciting physical phenomena [1, 2]. The types of materials, stacking configurations, and twisting angles can be controlled, enabling the physical properties to be further manipulated [3-14]. These excellent properties are related to the structure characteristics of two-dimensional materials (2DMs); 2DMs consist of strong intralayer bonds and weak VdW interlayer interactions [10, 11, 15]. Due to these weak force, layers in 2DMs can easily slide relative to each other and form the domain wall structure [16-18]. Recent studies reported that domain walls can cause topological valley transport phenomena etc. [19-24]. To better understand the formation mechanism of these physical phenomena, deep research into the interface properties between layers is needed. There are two critical energies of interfaces, the binding energy and generalized stacking-fault energy (GSFE), which determine the energy dissipation path of a system. Binding energy is defined as the energy per area required to separate the interaction structure into two parts [25], and indicates the ability of a structure to resist delamination. GSFE is the difference in energy between the stable stacking configuration and uniformly disregistered configurations [17], and GSFE is an indication of anti-sliding ability.

To date, many studies about the interface energy of bilayer heterostructures are reported [17, 25, 26]. Studies on multilayer 2DMs reported that the properties of heterostructures depend on the number of layers [27-29], which implies a long-range interaction within 2DM layers. Since that VdW interactions are long-range [30] and neighboring layers influence the interface energy of adjacent layers [31], multi-layer heterostructures are expected to have different interfacial properties from bilayer heterostructures. However, the interface energy of multilayer heterostructures are rarely reported. It is still unclear how neighboring layers affect the interface energy of adjacent layer. The definitions of the neighboring layer and adjacent layer are shown in supplementary materials. The heterostructures consisting of alternating stacking configuration of graphene (G) and boron nitride (BN) exhibited tunable metal-insulator transitions [32]. Bilayer G is an excellent platform for studying physical phenomena, such as topological edge states [12, 19]; when the physical phenomena of bilayer G are measured or when a new generation of electronic devices is created, bilayer G inevitably contacts the substrate. Also, BN has a structure similar to G with a 1.8% lattice mismatch, and BN has a wide band gap; these excellent properties make BN an ideal substrate [33]. Moreover, trilayer G has more adjustable degrees of freedom and exhibits many meaningful physical phenomena, such as electronic properties depending on stacking structure [34, 35].

Therefore, in this work, to investigate the interface energy of multilayer structures, we select trilayer heterostructures of G-G-G, G-G-BN, and BN-G-BN. Firstly, we prove a new expression that can describe the GSFE landscape of G-like structures. Then, the influence of the neighboring layer on the interface energy of the adjacent layers is investigated.

Ⅱ. CALCULATION METHODS

Atomic structures and interface energy are investigated using density functional theory (DFT) calculations. The Vienna ab initio Simulation Package (VASP) with projector augmented wave (PAW) potentials is employed. The energy cutoff of 500 eV is chosen [36]. The force and energy convergence criterion are set to 0.01 eV/Å and 10$^{-5}$ eV, respectively. The periodic boundary condition is set to 20 Å with a vacuum region. For the hexagonal heterostructure, we adopt a 30$\times$30$\times$1 Gamma pack mesh.

Full optimization is performed to obtain the stable stacking configuration. The PBE-D3 method can reproduce stacking configuration correctly [26]. Thus, the PBE-D3 method is chosen to describe VdW interactions.

Ⅲ. RESULTS AND DISCUSSION A. Interface energy of G/G and G/BN

Although the binding energy and equilibrium interlayer spacing of G/G and G/BN have been widely studied, there are still no uniform values. For example, the binding energy of G/G is calculated to be -22.8 meV/atom using LDA method, and is -91.35 meV/atom using RPA method [17, 26]. To compare the differences between bilayer heterostructures and trilayer heterostructures, we firstly calculate the interface energy of two bilayer heterostructures. Their optimized structures are shown in FIG. 1. The equilibrium interlayer spacing $d_{\rm{eq}}$ is defined as the distance between two layers in the optimized structure, as shown in FIG. 1 (c) and (d). The value of binding energy $\Delta E_{\rm{b}}$ can be obtained using Eq.(1),

 $\begin{eqnarray} \Delta {E_{\rm{b}}} = \frac{{ E({\rm{A}}) + E({\rm{B}}) - E({\rm{A}} + {\rm{B}})}}{S} \end{eqnarray}$ (1)
 FIG. 1 Stacking configurations of (a) top view for G/G, (b) top view for G/BN, (c) side view for G/G, and (d) side view for G/BN. The red dotted line indicates the corresponding positions of an atom in the top view and in the side view. The interlayer distance is 3.51 and 3.43 Å for G/G and G/BN, respectively.

where $E$(A+B) indicates the energy of the interaction structure. $E$(A) and $E$(B) are the energies of the isolated structures of A and B, respectively. $S$ refers to the area of the interface between two layers. The values of $d_{\rm{eq}}$ and $\Delta E_{\rm{b}}$ are listed in Table Ⅰ. $\Delta E_{\rm{b}}$ and $d_{\rm{eq}}$ of G/G are 277.89 mJ/m$^2$ and 3.51 Å, respectively, which are consistent with previously reported calculations [17, 26, 37]. $\Delta E_{\rm{b}}$ and $d_{\rm{eq}}$ of G/BN are 221.26 mJ/m$^2$ and 3.43 Å, respectively, which also agree with calculational and experimental values [17, 38]. The values of in-plane lattice parameter $a_0$ for G/G and G/BN are 2.47 Å and 2.49 Å, respectively. These results are in agreement with existing reports [17, 26, 37].

Table Ⅰ The stacking fault energy $\gamma_{\rm{sf}}$ (in mJ/m$^2$), maximum stacking fault energy $\gamma_{\rm{peak}}$ (in mJ/m$^2$), binding energy $\Delta E_{\rm{b}}$ (in mJ/m$^2$), equilibrium interlayer spacing $d_{\rm{eq}}$ (in Å), in-plane lattice parameter $a_0$ (in Å), amplitude $W_h$ (in mJ/m$^2$) in Eq.(3), and ratio $\delta$ in Eq.(3) for bilayer G/G and G/BN.

GSFE is introduced to describe atomic interactions across a crystal face [39], and the expression was improved by Xiang et al. to describe the $\gamma_{\rm{sf}}$ surface of metal [40]. Subsequently, Zhou et al. found that GSFE landscapes of G/G and G/BN can also be described well [17], the expression is as follows,

 $\begin{eqnarray} F_1(\phi , \varphi ) \hspace{-0.12cm}&=&\hspace{-0.12cm} c_0\hspace{-0.06cm} +\hspace{-0.06cm} c_1\bigg[\cos\frac{2\pi }{a_0}\bigg(\phi \hspace{-0.06cm} +\hspace{-0.06cm} \frac{\varphi }{\sqrt 3 }\bigg) \hspace{-0.06cm} +\hspace{-0.06cm} \cos\frac{2\pi }{a_0}\bigg(\phi \hspace{-0.06cm} -\hspace{-0.06cm} \frac{\varphi }{\sqrt 3 }\bigg) \nonumber\\ &&{}\hspace{-0.12cm} + \cos\frac{4\pi \varphi }{\sqrt 3a_0}\bigg] + c_2\bigg[\cos\frac{2\pi }{a_0}(\phi + \sqrt 3 \varphi )\nonumber\\ &&{}\hspace{-0.12cm} + \cos\frac{2\pi }{a_0}(\phi - \sqrt 3 \varphi ) + \cos\frac{4\pi \phi }{a_0}\bigg]+\nonumber\\ &&{}\hspace{-0.12cm}c_3\bigg[\cos\frac{2\pi }{a_0}\bigg(2\phi + \frac{2\varphi }{\sqrt 3 }\bigg) + \cos\frac{2\pi }{a_0}\bigg(2\phi -\frac{2\varphi }{\sqrt 3 }\bigg)\nonumber\\ &&{}\hspace{-0.12cm} + \cos\frac{8\pi \varphi }{\sqrt 3 a_0}\bigg]+ c_4\bigg[\sin \frac{2\pi }{a_0}\bigg(\phi -\frac{\varphi }{\sqrt 3 }\bigg)\nonumber\\ &&{}\hspace{-0.12cm} - \sin \frac{2\pi }{a_0}\bigg(\phi + \frac{\varphi }{\sqrt 3 }\bigg) +\sin \frac{4\pi \varphi }{\sqrt 3 a_0}\bigg]+\nonumber\\ &&{}\hspace{-0.12cm}c_5\bigg[\sin \frac{2\pi }{a_0}\bigg(2\phi - \frac{2\varphi }{\sqrt 3 }\bigg) -\cos\frac{2\pi }{a_0}\bigg(2\phi +\frac{2\varphi }{\sqrt 3 }\bigg) \nonumber\\ &&{}\hspace{-0.12cm} +\cos\frac{8\pi \varphi }{\sqrt 3a_0}\bigg] \end{eqnarray}$ (2)

where $c_0$, $c_1$, $c_2$, $c_3$, $c_4$, $c_5$ are constants. $\phi$ and $\varphi$ are the disregistries along the zigzag and armchair directions, respectively. The periodic length of the armchair and zigzag directions are $\sqrt{3}a_0$ and $a_0$ respectively. $a_0$ refers to the in-plane lattice parameter for G/G and G/BN, the values of $a_0$ are 2.47 and 2.49, respectively. Note that $F_1(\phi, \varphi)$ is the excess energy per area (in reference to the stable stacking configuration) at $(\phi, \varphi)$=(0, 0).

The stacking fault energy ($\gamma_{\rm{sf}}$), the unstable stacking fault (saddle point) energy in the armchair direction ($\gamma_{\rm{sp}}$), and the maximum energy in the GSFE landscape ($\gamma_{\rm{peak}}$) are marked in FIG. 2. The unstable stacking fault energy in the zigzag direction ($\gamma_{\rm{sp}}$) is also indicated in FIG. 3. The definition of these energies is proposed by Wang et al. [25].

 FIG. 2 GSFE vs. disregistry along the armchair direction for (a) G/BN and (b) G/G. Reference configurations are both set to be stable stacking configurations. The positions and corresponding disregistries of $\gamma_{\rm{sf}}$, $\gamma_{\rm{sp}}$, and $\gamma_{\rm{peak}}$ are marked.
 FIG. 3 GSFE landscapes derived using Eq.(2) for (a) G/BN and (b) G/G. GSFE landscapes derived using Eq.(3) for (c) G/BN and (d) G/G. The unit of energy is mJ/m$^2$. The arrows represent the preferred sliding directions (SDs). The locations of $\gamma_{\rm{sf}}$, $\gamma_{\rm{sp}}$, $\gamma_{\rm{peak}}$, and $\gamma_{\rm{sp}'}$ are indicated.

For G/G and G/BN, GSFE versus disregistry along the armchair direction is shown in FIG. 2 (a) and (b), respectively. The values of $\gamma_{\rm{peak}}$ and $\gamma_{\rm{sf}}$ for G/G are 32.86 and 0.00 mJ/m$^2$, respectively. The values of $\gamma_{\rm{peak}}$ and $\gamma_{\rm{sf}}$ for G/BN are 40.54 and 34.19 mJ/m$^2$, respectively. These results are close to the values obtained from calculations using ACFDT-PRA [17], which indicates GSFE can be computed accurately by DFT-D3. For both FIG. 2 (a) and (b), six typical points are selected to fit Eq.(2). A GSFE landscape can be determined from the fitted parameters in the same way as proposed by Zhou et al. [17]. Selected disregistries and their corresponding energies are given in supplementary materials (Table S1), and the fitted parameters of $c_0$, $c_1$, $c_2$, $c_3$, $c_4$, $c_5$ are listed in Table Ⅱ. GSFE landscapes of G/BN and G/G are shown in FIG. 3 (a) and (b), respectively.

Table Ⅱ The fitted $c_0$, $c_1$, $c_2$, $c_3$, $c_4$, $c_5$ (in mJ/m$^2$) in Eq.(2) for G/G and G/BN.

Eq.(2) expresses well the GSFE landscape of G-like structure, but the expression is relatively complex, and the fitting method needs to be selected according to the symmetry of structure [17]. Marom et al. also proposed the registry index model [41], but the parameters do not directly reflect the characteristic of GSFE. On the other hand, the GSFE of a triangular lattice structure has been well described by Eq.(3) [42]. The parameter in Eq.(3) directly reflects the ratio of $\gamma_{\rm{peak}}$ to $\gamma_{\rm{sf}}$. Because G and BN are hexagonal lattices, which are similar to a triangular lattice, it is believed that Eq.(3) can also describe the GSFE landscape of a G-like structure. The expression of Eq.(3) is as follows,

 $\begin{eqnarray} F(\phi , \varphi ) \hspace{-0.12cm}&=&\hspace{-0.12cm} W_h\bigg[\frac{3}{2} + \cos 2\pi \bigg(\frac{\phi }{a_0}+ \frac{\varphi }{\sqrt 3a_0} + \frac{2}{3}\bigg) +\nonumber\\ &&{}\cos 2\pi \bigg(\frac{2\varphi }{\sqrt3a_0} + \frac{1}{3}\bigg) +\cos 2\pi \bigg(\frac{\phi }{a_0}- \nonumber\\ &&{} \frac{\varphi }{\sqrt3a_0}+\frac{1}{3}\bigg)\bigg]+W_h\delta\bigg[\frac{3}{2} +\cos 2\pi \bigg(\frac{\phi }{a_0}+\nonumber\\ &&{} \frac{\varphi }{\sqrt3a_0}+\frac{1}{3}\bigg) +\cos 2\pi \bigg(\frac{2\varphi }{\sqrt3a_0} + \frac{2}{3}\bigg) + \nonumber\\ &&{} \cos 2\pi \bigg(\frac{\phi }{a_0} - \frac{2\varphi }{\sqrt3a_0}- \frac{1}{3}\bigg)\bigg] \end{eqnarray}$ (3)

where $W_h$ is the amplitude, and its value can be determined by $\gamma_{\rm{peak}}$:

 $\begin{eqnarray} W_h=0.22\gamma_{\rm{peak}} \end{eqnarray}$ (4)

$\delta$ is the ratio of $\gamma_{\rm{sf}}$ to $\gamma_{\rm{peak}}$:

 $\begin{eqnarray} \delta =\frac{\gamma _{\rm{sf}}}{\gamma _{\rm{peak}}} \end{eqnarray}$ (5)

It should be noted that we firstly directly relate $W_h$ to $\gamma_{\rm{peak}}$, and $W_h$ is only a reference value in the original equation. The values of $W_h$ for G/G and G/BN are 7.23 and 8.92 mJ/m$^2$, and $\delta$ for G/G and G/BN are 0.00 and 0.84, as listed in Table Ⅰ. The GSFE landscapes of G/BN and G/G described by Eq.(3) are shown in FIG. 3 (c) and (d), respectively.

FIG. 3(c) gives the same values and positions of $\gamma_{\rm{sf}}$ and $\gamma_{\rm{peak}}$ for G/BN shown in FIG. 3(a). The sliding direction and energy barrier can be determined by comparing the values of $\gamma_{\rm{sp}}$ and $\gamma_{\rm{sp}'}$ [17]. The preferred sliding direction for G/BN is zigzag, and the associated energy barrier is 32.81 mJ/m$^2$, which are consistent with FIG. 3(a). FIG. 3(d) also accurately describes the values and positions of $\gamma_{\rm{sf}}$ and $\gamma_{\rm{peak}}$ for G/G shown in FIG. 3(b). The preferred sliding direction for G/G is armchair, and the associated energy barrier is 3.62 mJ/m$^2$; the information agrees well with FIG. 3(b). Therefore, Eq.(3) can also accurately describe the GSFE landscapes of G-like structures. Compared with Eq.(2), Eq.(3) is easier to fit; the form is more compact, and it directly reflects the characteristic of the GSFE from $\delta$. The relatively large value of $\delta$ indicates that the system has high-level stacking-fault energy, and the preferred sliding direction is zigzag. Otherwise, the system has low-level stacking-fault energy, and the armchair direction is preferred, which is similar to the previous conclusion [42].

B. Interface energy of G-G-G, G-G-BN, and BN-G-BN 1. Stable stacking configurations

The stable stacking configurations of G-G-G, G-G-BN, and BN-G-BN are shown in FIG. 4(a)-(f). The values of $d_{\rm{eq}}$ and $\Delta E_{\rm{b}}$ are listed in Table Ⅲ. For G-G-G, the energetically preferred stacking configuration is ABA, which is consistent with the previously reported results [43]. The upper and lower Gs are symmetric about the middle layer, and thus, the interface has only one type. G-G-G can be seen as the structure of a monolayer G adsorbed on bilayer AB-stacked G, and the interface is G/G. This structure is denoted as G-G/G, whose $\Delta E_{\rm{b}}$ is 297.05 mJ/m$^2$. Compared with G/G, this value increases by 6.89%. The value of $d_{\rm{eq}}$ is still 3.51 Å.

 FIG. 4 Stable stacking configurations for the (a) top view of G-G-G, (b) top view of G-G-BN, (c) top view of BN-G-BN, (d) side view of G-G-G, (e) side view of G-G-BN, and (f) side view of BN-G-BN. The red dotted lines indicate the corresponding positions of an atom in the top view and in the side view.
Table Ⅲ The stacking fault energy $\gamma_{\rm{sf}}$ (in mJ/m$^2$), maximum stacking fault energy $\gamma_{\rm{peak}}$ (in mJ/m$^2$), binding energy $\Delta E_{\rm{b}}$ (in mJ/m$^2$), equilibrium interlayer spacing $d_{\rm{eq}}$ (in Å), in-plane lattice parameter $a_0$ (in Å), amplitude $W_h$ (mJ/m$^2$) in Eq.(3), and ratio $\delta$ in Eq.(3) for trilayer G-G/G, BN-G/G, G-G/BN, and BN-G/BN.

The energetically preferred stacking configuration of G-G-BN is ABC. There are two types of interfaces. The first type consists of monolayer G adsorbed onto bilayer G-BN, and the interface is G/G. This structure is denoted as BN-G/G, and $\Delta E_{\rm{b}}$ is 273.21 mJ/m$^2$, which is 1.68% lower than $\Delta E_{\rm{b}}$ for G/G. The value of $d_{\rm{eq}}$ is still 3.51 Å. The second type consists of monolayer BN adsorbed on bilayer AB-stacked G-G, and the interface is G/BN. We denote this as G-G/BN. The values of $\Delta E_{\rm{b}}$ and $d_{\rm{eq}}$ are 220.75 mJ/m$^2$ and 3.39 Å, respectively. These values have changed by -2.30% and -1.17%, respectively, of the corresponding values for G/BN. This weakening trend of $\Delta E_{\rm{b}}$ was also reported by Yelgel [31].

Similar to G-G-G, the stable stacking configuration of BN-G-BN is ABA. The upper and lower BN layers are symmetric about the middle layer G, and thus, there is only one type of interface. BN-G-BN can be viewed as monolayer BN adsorbed on bilayer G-BN with an interface of G/BN; we denote this as BN-G/BN. The values of $\Delta E_{\rm{b}}$ and $d_{\rm{eq}}$ for BN/G-BN are 271.15 mJ/m$^2$ and 3.39 Å, respectively. The changes in these data are 22.55% and -1.17%, respectively, of the corresponding values for G/BN.

The neighboring layer has little influence on the value of $d_{\rm{eq}}$ for the adjacent layers, and this is consistent with previously calculated trends [26, 31]. However, the neighboring layer has a relatively large influence on $\Delta E_{\rm{b}}$, and this may further explain the difference between theoretical and experimental values [25]. In addition to the limitations of the calculations and experimental method themselves, the measured value is unavoidably affected by the neighboring layer or substrate. The extent of the impact is related to the level of polarity of the neighboring layer. Polar materials, such as BN, have a complex influence. However, nonpolar materials, such as G, have simpler influence. In most cases, the binding energy between layers can be well estimated using the Lennard-Jones (LJ) pair potential [44, 45]. The functional derived from LJ potential relationship between binding energy and layer spacing is

 $\begin{eqnarray} \Delta E_{\rm{b}}(d) = \frac{U_0}{d^4} \end{eqnarray}$ (6)

where $U_0$ is the constant associated with the interaction structure. $d$ is the layer spacing. (See supplementary materials for the detailed derivation of this equation.) Then we studied the relationship between $\Delta E_{\rm{b}}$ and $d$, as shown in FIG. 5.

 FIG. 5 Relationship between binding energy $\Delta E_{\rm{b}}$ and layer spacing $d$, and $d_{\rm{eq}}$ is the equilibrium interlayer spacing.

As seen in FIG. 5, when $d$ is increased to $2d_{\rm{eq}}$, $\Delta E_{\rm{b}}$ decreases to 6.25% of the original value $U_0$, therefore, if only the LJ potential is considered, the binding energy of adjacent layers will be about 6% larger due to the effect of neighboring layer. Also, the effect is negligible when $d$ continues to expand. The relative deviation of $\Delta E_{\rm{b}}$ between G-G/G and G/G confirms this trend. To further confirm our conclusion, we compare the theoretical and experimental values. The experimental value of graphite is about 390.00 mJ/m$^2$ [25], and the value we calculate for G/G is 277.89 mJ/m$^2$. Because graphite is a multilayer material, in which the upper and lower layers of G/G both have adjacent layers, the calculated value should be increased by 12% if the experimental value and the theoretical value are to be compared. We can find that $\Delta E_{\rm{b}}$ of four-layer graphene is 311.24 mJ/m$^2$, which is closer to the experimental value. But it does not explain why the interaction structure with BN is different from this trend. In fact, the influence of VdW forces is very complicated. For example, VdW forces can cause redistribution of charge [30, 46]. The analysis of Bader charge shows that, the electron distributions of G-G-BN and BN-G-BN are obviously more uneven than that of G-G-G, and redistribution of charge also causes variation in energy. This may be a simple explanation for why the heterojunction with BN is different, and more detailed mechanisms require further research. The Bader charge analyses of G-G-G, G-G-BN, and BN-G-BN are reported in Table S2 in supplementary materials.

As a result of our calculation, $\Delta E_{\rm{b}}$ is dominated by the adjacent layer, and the neighboring layer can modify its value.

2. GSFE of G-G/G, BN-G/G, G-G/BN, and BN-G/BN

In this section, we investigate the effect of the neighboring layer on the GSFE of adjacent layers. FIG. 6 (a)-(d) shows GSFE versus disregistry along the armchair direction for G-G/G, BN-G/G, G-G/BN, and BN-G/BN. The values of $\gamma_{\rm{peak}}$ and $\gamma_{\rm{sf}}$ are given in Table Ⅲ. For G-G/G, $\gamma_{\rm{peak}}$ is 36.35 mJ/m$^2$, and $\gamma_{\rm{sf}}$ is 0.00 mJ/m$^2$. For BN-G/G, $\gamma_{\rm{peak}}$ is 32.18 mJ/m$^2$, and $\gamma_{\rm{sf}}$ is still 0.00 mJ/m$^2$. For G-G/BN, $\gamma_{\rm{peak}}$ and $\gamma_{\rm{sf}}$ are 45.20 and 36.87 mJ/m$^2$, respectively. For BN-G/BN, $\gamma_{\rm{peak}}$ and $\gamma_{\rm{sf}}$ are 46.15 and 36.79 mJ/m$^2$, respectively. The values of $a_0$ for G-G/G, BN-G/G, G-G/BN, and BN-G/BN are 2.47, 2.48, 2.48, and 2.49 Å, respectively. These data, together with $W_h$ and $\delta$ calculated from Eq.(4) and Eq.(5), are reported in Table Ⅲ.

 FIG. 6 GSFE versus disregistry along the armchair direction for (a) G-G/G, (b) BN-G/G, (c) G-G/BN, and (d) BN-G/BN. Reference configurations are set as the respective stable stacking configurations. The positions and corresponding disregistries of $\gamma_{\rm{sf}}$, $\gamma_{\rm{sp}}$, and $\gamma_{\rm{peak}}$ are marked.

Based on Eq.(3), we also obtained the GSFE landscapes of G-G/G, BN-G/G, G-G/BN, and BN-G/BN, respectively, as shown in FIG. 7(a)-(d). The layers in these vdW heterostructures can slide relative to each other, and the sliding direction is indicated in the GSFE landscape. The sliding directions are the armchair direction for G-G/G and BN-G/G and the zigzag direction for G-G/BN and BN-G/BN. By comparing the sliding directions of G-G/G and BN-G/G with that of G/G, we find that the sliding direction of G/G is not changed by the neighboring layer. Similarly, comparation of the sliding directions of G-G/BN, BN-G/BN, and G/BN shows that, the sliding direction G/BN is not altered by the neighboring layer. In short, the GSFE landscape of adjacent layers is hardly changed by neighboring layers.

 FIG. 7 GSFE landscapes derived using Eq.(3) for (a) G-G/G, (b) BN-G/G, (c) G-G/BN, and (d) BN-G/BN. Units are mJ/m$^2$. The sliding direction (SD) and the locations of $\gamma_{\rm{sf}}$, $\gamma_{\rm{sp}}$, $\gamma_{\rm{peak}}$, and $\gamma_{\rm{sp}'}$ are indicated.

To quantitatively analyze the influence of the neighboring layer on the GSFE of adjacent layers, we compared the value of GSFE of the trilayer heterostructure with that of the bilayer heterostructure. The calculation method is:

 $\begin{eqnarray} \frac{F_{\rm{tri}}(\phi , \varphi ) -F_{\rm{bi}}(\phi , \varphi )}{F_{\rm{bi}}(\phi , \varphi )}\nonumber \end{eqnarray}$

$F_{\rm{tri}}(\phi, \varphi)$ refers to the GSFE of a trilayer heterostructure. $F_{\rm{bi}}(\phi, \varphi)$ indicates the GSFE of a bilayer heterostructure. The interfaces for comparison are the same type. The calculated results are shown in FIG. 8. Compared with the GSFE of G/G, the GSFE of G-G/G increases by 10.65%, and the GSFE of BN-G/G decreases by 2.07%. Compared with the GSFE of G/BN, the GSFE of G-G/BN increases by 9.00%-11.00%, and the GSFE of BN-G/BN increases by 8.50%-10.50%. The neighboring layer has a relatively large influence on the value of the GSFE of adjacent layers. On the other hand, as seen in FIG. 8 (a) and (b), the influence of the neighboring layer G or BN on the GSFE of the G/G interface is uniform. However, the influence of the neighboring layer G or BN on the GSFE of the G/BN interface is not uniform. Thus, the changed value of GSFE depends on the property of the interface between adjacent layers.

 FIG. 8 Variation values of the GSFE for (a) G-G/G versus G/G, (b) BN-G/G versus G/G, (c) G-G/BN versus G/BN, and (d) BN-G/BN versus G/BN.

To summarize, the influence of neighboring layers on the GSFE of adjacent layers is that it does not alter the features of a GSFE landscape, but the neighboring layers can modify the values. Also, the distribution of changed values depends on the interface property.

Ⅳ. CONCLUSION

Using first-principles calculations, we have compared the differences in the interface energies between trilayer heterostructures (graphene/graphene/graphene, graphene/graphene/boron nitride, boron nitride/graphene/boron nitride) and bilayer heterostructures (graphene/graphene, graphene/boron nitride). It is found that the neighboring layer has little effect on the equilibrium layer spacing of adjacent layers, but it has a relatively large influence on the binding energy. The neighboring layer changes the value of GSFE of the adjacent layer, and doesn't change the GSFE landscape. We also found that the relative deviation of the GSFE depends on the interface properties of the adjacent layer. In addition, we proved a new expression that can be used to describe the GSFE landscape of graphene-like structures.

Although different two-dimensional layered materials have their own characteristics, they share some common features. In particular, almost all of the layers are coupled via weak by VdW forces. Deep understanding of any member of this family of materials will have universal value for the larger family.

Supplementary materials: The definitions of the neighboring layer and adjacent layer, the detailed derivation of the relationship between the binding energy and interlayer spacing, the values used to fit GSFE landscapes, and the results of Bader charge analysis are given in the supplementary materials.

Supporting Information
 FIG. S1 The schematic diagram of the neighboring layer and adjacent layer.

2. Derivation of the relationship between the binding energy and interlayer spacing

 FIG. S2 The schematic diagram of bilayer vdW heterostructure

The atom interaction between two adjacent layers is complex, and herein we take the standard form of the Lennard-Jones potential:

 ${U_{LJ}}(r) = \frac{A}{{{r^{12}}}} - \frac{B}{{{r^6}}}$ (1)

where r is the distance between the two atoms, A is the constant for the repulsive interaction, and B is the constant for the attractive interaction. Then, we integrate the energy between one atom in the upper layer and all the atoms in the lower layer to obtain the density of the interaction energy Uvdw :

 $\begin{array}{l} {U_{vdw}}({d_{eq}}) = 2\pi \rho \int\limits_0^\infty {[\frac{A}{{{{(d_{eq}^2 + {x^2})}^6}}} - \frac{B}{{{{(d_{eq}^2 + {x^2})}^3}}}]} xdx\\ = \pi \rho (\frac{A}{{5d_{eq}^{10}}} - \frac{B}{{2d_{eq}^4}}) \end{array}$ (2)

where $\rho$ is the number of atoms per unit area of the lower layer, deq refers to the equilibrium interlayer spacing.

And the derivative of Uvdw to deq is 0, namely,

 $\frac{{\partial {U_{vdw}}}}{{\partial {d_{eq}}}} = \pi \rho (\frac{{ - 10A}}{{5d_{eq}^{11}}} - \frac{{ - 4B}}{{2d_{eq}^5}}) = 0 \to A = d_{eq}^6*B$ (3)

Formula (3) is brought into Equation (2):

 $U_{v d w}\left(d_{e q}\right)=\frac{-3 \pi \rho B}{10} \frac{1}{d_{e q}^{4}}=\frac{-U_{0}}{d_{e q}^{4}}\left(U_{0}=\frac{3 \pi \rho B}{10}\right)$ (4)

The binding energy is that:

 $\Delta E_{b}=-U_{vdw}$ (5)
Table S1 The dataset of selected disregistries along the armchair direction and corresponding energies (mJ/m2) that can be employed to fit the GSFE landscapes.
Table S2 Bader charge analysis for G-G-G, G-G-BN, and BN-G-BN. The average number of transferred electrons per layer (10-4 e). '+' indicates atoms obtain electron, and '-' indicates atoms lose electron.
Reference
 [1] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43(2018). DOI:10.1038/nature26160 [2] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, and E. Kaxiras, Nature 556, 80(2018). DOI:10.1038/nature26154 [3] D. K. Bediako, M. Rezaee, H. Yoo, D. T. Larson, S. F. Zhao, T. Taniguchi, K. Watanabe, T. L. Brower-Thomas, E. Kaxiras, and P. Kim, Nature 558, 425(2018). DOI:10.1038/s41586-018-0205-0 [4] J. S. Alden, A. W. Tsen, P. Y. Huang, R. Hovden, L. Brown, J. Park, D. A. Muller, and P. L. McEuen, Proc. Natl. Acad. Sci. USA 110, 11256(2013). DOI:10.1073/pnas.1309394110 [5] R. Bistritzer, and A. H. MacDonald, Proc. Natl. Acad. Sci. USA 108, 12233(2011). DOI:10.1073/pnas.1108174108 [6] T. Ohta, J. T. Robinson, P. J. Feibelman, A. Bostwick, E. Rotenberg, and T. E. Beechem, Phys. Rev. Lett. 109, 186807(2012). DOI:10.1103/PhysRevLett.109.186807 [7] J. He, K. Hummer, and C. Franchini, Phys. Rev. B 89, 075409(2014). DOI:10.1103/PhysRevB.89.075409 [8] Y. Shimazaki, M. Yamamoto, I. V. Borzenets, K. Watanabe, T. Taniguchi, and S. Tarucha, Nat. Phys. 11, 1032(2015). DOI:10.1038/nphys3551 [9] G. Li, A. Luican, J. L. Dos Santos, A. C. Neto, A. Reina, J. Kong, and E. Andrei, Nat. Phys. 6, 109(2010). DOI:10.1038/nphys1463 [10] A. K. Geim, and I. V. Grigorieva, Nature 499, 419(2013). DOI:10.1038/nature12385 [11] Y. Liu, N. O. Weiss, X. Duan, H. C. Cheng, Y. Huang, and X. Duan, Nat. Rev. Mater. 1, 16042(2016). DOI:10.1038/natrevmats.2016.42 [12] G. M. Rutter, S. Jung, N. N. Klimov, D. B. Newell, N. B. Zhitenev, and J. A. Stroscio, Nat. Phys. 7, 649(2011). DOI:10.1038/nphys1988 [13] X. Chen, C. Tan, Q. Yang, R. Meng, Q. Liang, J. Jiang, X. Sun, D. Yang, and T. Ren, Phys. Chem. Chem. Phys. 18, 16229(2016). DOI:10.1039/C6CP01083K [14] J. B. Wu, Z. X. Hu, X. Zhang, W. P. Han, Y. Lu, W. Shi, X. F. Qiao, M. IjiaÌŝ, S. Milana, and W. Ji, ACS Nano 9, 7440(2015). DOI:10.1021/acsnano.5b02502 [15] J. Wang, G. Liang, and K. Wu, J. Chem. Phys. 30, 649(2017). [16] I. V. Lebedeva, A. A. Knizhnik, A. M. Popov, Y. E. Lozovik, and B. V. Potapkin, Phys. Chem. Chem. Phys. 13, 5687(2011). DOI:10.1039/c0cp02614j [17] S. Zhou, J. Han, S. Dai, J. Sun, and D. J. Srolovitz, Phys. Rev. B 92, 155438(2015). DOI:10.1103/PhysRevB.92.155438 [18] S. Dai, Y. Xiang, and D. J. Srolovitz, Phys. Rev. B 93, 085410(2016). DOI:10.1103/PhysRevB.93.085410 [19] L. J. Yin, H. Jiang, J. B. Qiao, and L. He, Nat. Commun. 7, 11760(2016). DOI:10.1038/ncomms11760 [20] Y. Zhang, T. T. Tang, C. Girit, Z. Hao, M. C. Martin, A. Zettl, M. F. Crommie, Y. R. Shen, and F. Wang, Nature 459, 820(2009). DOI:10.1038/nature08105 [21] G. W. Semenoff, V. Semenoff, and F. Zhou, Phys. Rev. Lett. 101, 087204(2008). DOI:10.1103/PhysRevLett.101.087204 [22] X. Li, F. Zhang, Q. Niu, and A. MacDonald, Phys. Rev. Lett. 113, 116803(2014). DOI:10.1103/PhysRevLett.113.116803 [23] L. Jiang, Z. Shi, B. Zeng, S. Wang, J. H. Kang, T. Joshi, C. Jin, L. Ju, J. Kim, and T. Lyu, Nat. Mater. 15, 840(2016). DOI:10.1038/nmat4653 [24] L. Ju, Z. Shi, N. Nair, Y. Lv, C. Jin, J. Velasco Jr., C. Ojeda-Aristizabal, H. A. Bechtel, M. C. Martin, and A. Zettl, Nature 520, 650(2015). DOI:10.1038/nature14364 [25] W. Wang, S. Dai, X. Li, J. Yang, D. J. Srolovitz, and Q. Zheng, Nat. Commun. 6, 7853(2015). DOI:10.1038/ncomms8853 [26] I. V. Lebedeva, A. V. Lebedev, A. M. Popov, and A. A. Knizhnik, Comput. Mater. Sci. 128, 45(2017). DOI:10.1016/j.commatsci.2016.11.011 [27] W. Gao, and A. Tkatchenko, Phys. Rev. Lett. 114, 096101(2015). DOI:10.1103/PhysRevLett.114.096101 [28] M. Tosun, D. Fu, S. B. Desai, C. Ko, J. S. Kang, D. H. Lien, M. Najmzadeh, S. Tongay, J. Wu, and A. Javey, Sci. Rep. 5, 10990(2015). DOI:10.1038/srep10990 [29] T. Korhonen, and P. Koskinen, Phys. Rev. B 92, 115427(2015). DOI:10.1103/PhysRevB.92.115427 [30] A. Ambrosetti, N. Ferri, R. A. DiStasio, and A. Tkatchenko, Science 351, 1171(2016). DOI:10.1126/science.aae0509 [31] C. Yelgel, J. Appl. Phys. 119, 065307(2016). DOI:10.1063/1.4941552 [32] L. Ponomarenko, A. Geim, A. Zhukov, R. Jalil, S. Morozov, K. Novoselov, I. Grigorieva, E. Hill, V. Cheianov, and V. Fal'Ko, Nat. Phys. 7, 958(2011). DOI:10.1038/nphys2114 [33] A. Pakdel, Y. Bando, and D. Golberg, Chem. Soc. Rev. 43, 934(2014). DOI:10.1039/C3CS60260E [34] S. Hattendorf, A. Georgi, M. Liebmann, and M. Morgenstern, Surf. Sci. 610, 53(2013). DOI:10.1016/j.susc.2013.01.005 [35] M. Yankowitz, I. Joel, J. Wang, A. G. Birdwell, Y. A. Chen, K. Watanabe, T. Taniguchi, P. Jacquod, P. San-Jose, and P. Jarillo-Herrero, Nat. Mater. 13, 786(2014). DOI:10.1038/nmat3965 [36] S. Chen, J. Wang, J. Huang, and Q. Li, Chin. J. Chem. Phys. 30, 36(2017). DOI:10.1063/1674-0068/30/cjcp1605113 [37] Y. Baskin, and L. Meyer, Phys. Rev. 100, 544(1955). DOI:10.1103/PhysRev.100.544 [38] S. J. Haigh, A. Gholinia, R. Jalil, S. Romani, L. Britnell, D. C. Elias, K. S. Novoselov, L. A. Ponomarenko, A. K. Geim, and R. Gorbachev, Nat. Mater. 11, 764(2012). DOI:10.1038/nmat3386 [39] V. Vitek, Philos. Mag. 18, 773(1968). DOI:10.1080/14786436808227500 [40] S. Dai, Y. Xiang, and D. J. Srolovitz, Acta Mater. 61, 1327(2013). DOI:10.1016/j.actamat.2012.11.010 [41] N. Marom, J. Bernstein, J. Garel, A. Tkatchenko, E. Joselevich, L. Kronik, and O. Hod, Phys. Rev. Lett. 105, 046801(2010). DOI:10.1103/PhysRevLett.105.046801 [42] J. Snyman, and H. Snyman, Surf. Sci. 105, 357(1981). DOI:10.1016/0039-6028(81)90168-0 [43] M. Aoki, and H. Amawashi, Solid State Commun. 142, 123(2007). DOI:10.1016/j.ssc.2007.02.013 [44] Z. Ye, A. Otero-De-La-Roza, E. R. Johnson, and A. Martini, Nanotechnology 26, 165701(2015). DOI:10.1088/0957-4484/26/16/165701 [45] M. Neek-Amal, and F. Peeters, Appl. Phys. Lett. 104, 041909(2014). DOI:10.1063/1.4863661 [46] M. Sadhukhan, and A. Tkatchenko, Phys. Rev. Lett. 118, 210402(2017). DOI:10.1103/PhysRevLett.118.210402

a. 中国科学技术大学近代力学系，中国科学院材料力学行为和设计重点实验室，合肥 230026;
b. 中国科学技术大学化学与材料科学学院，合肥微尺度物质科学国家研究中心，能源材料化学协同创新中心，中国科学院纳米科学卓越创新中心，合肥 230026