The article information
 Leilei Li, Shuo Feng
 李磊磊, 冯硕
 Influence of Neighboring Layers on Interfacial Energy of Adjacent Layers
 近邻层对相邻层界面能的影响
 Chinese Journal of Chemical Physics, 2019, 32(6): 693700
 化学物理学报, 2019, 32(6): 693700
 http://dx.doi.org/10.1063/16740068/cjcp1812291

Article history
 Received on: December 25, 2018
 Accepted on: February 13, 2019
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
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 [314]. These excellent properties are related to the structure characteristics of twodimensional 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 [1618]. Recent studies reported that domain walls can cause topological valley transport phenomena etc. [1924]. 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 stackingfault 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 antisliding 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 [2729], which implies a longrange interaction within 2DM layers. Since that VdW interactions are longrange [30] and neighboring layers influence the interface energy of adjacent layers [31], multilayer 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 metalinsulator 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 GGG, GGBN, and BNGBN. Firstly, we prove a new expression that can describe the GSFE landscape of Glike structures. Then, the influence of the neighboring layer on the interface energy of the adjacent layers is investigated.
Ⅱ. CALCULATION METHODSAtomic 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
Full optimization is performed to obtain the stable stacking configuration. The PBED3 method can reproduce stacking configuration correctly [26]. Thus, the PBED3 method is chosen to describe VdW interactions.
Ⅲ. RESULTS AND DISCUSSION A. Interface energy of G/G and G/BNAlthough 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
$ \begin{eqnarray} \Delta {E_{\rm{b}}} = \frac{{ E({\rm{A}}) + E({\rm{B}})  E({\rm{A}} + {\rm{B}})}}{S} \end{eqnarray} $  (1) 
where
GSFE is introduced to describe atomic interactions across a crystal face [39], and the expression was improved by Xiang et al. to describe the
$\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
The stacking fault energy (
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
Eq.(2) expresses well the GSFE landscape of Glike 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
$ \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
$ \begin{eqnarray} W_h=0.22\gamma_{\rm{peak}} \end{eqnarray} $  (4) 
$ \begin{eqnarray} \delta =\frac{\gamma _{\rm{sf}}}{\gamma _{\rm{peak}}} \end{eqnarray} $  (5) 
It should be noted that we firstly directly relate
FIG. 3(c) gives the same values and positions of
The stable stacking configurations of GGG, GGBN, and BNGBN are shown in FIG. 4(a)(f). The values of
The energetically preferred stacking configuration of GGBN is ABC. There are two types of interfaces. The first type consists of monolayer G adsorbed onto bilayer GBN, and the interface is G/G. This structure is denoted as BNG/G, and
Similar to GGG, the stable stacking configuration of BNGBN is ABA. The upper and lower BN layers are symmetric about the middle layer G, and thus, there is only one type of interface. BNGBN can be viewed as monolayer BN adsorbed on bilayer GBN with an interface of G/BN; we denote this as BNG/BN. The values of
The neighboring layer has little influence on the value of
$\begin{eqnarray} \Delta E_{\rm{b}}(d) = \frac{U_0}{d^4} \end{eqnarray} $  (6) 
where
As seen in FIG. 5, when
As a result of our calculation,
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 GG/G, BNG/G, GG/BN, and BNG/BN. The values of
Based on Eq.(3), we also obtained the GSFE landscapes of GG/G, BNG/G, GG/BN, and BNG/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 GG/G and BNG/G and the zigzag direction for GG/BN and BNG/BN. By comparing the sliding directions of GG/G and BNG/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 GG/BN, BNG/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.
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} $ 
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.
Ⅳ. CONCLUSIONUsing firstprinciples 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 graphenelike structures.
Although different twodimensional 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 Information2. Derivation of the relationship between the binding energy and interlayer spacing
The atom interaction between two adjacent layers is complex, and herein we take the standard form of the LennardJones 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 U_{vdw} :
$ \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
And the derivative of U_{vdw} to d_{eq} 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) 
[1] 
Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. JarilloHerrero, 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. SanchezYamagishi, 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. BrowerThomas, E. Kaxiras, and P. Kim, Nature
558, 425(2018).
DOI:10.1038/s4158601802050 
[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. OjedaAristizabal, 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. SanJose, and P. JarilloHerrero, 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/16740068/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/00396028(81)901680 
[43] 
M. Aoki, and H. Amawashi, Solid State Commun.
142, 123(2007).
DOI:10.1016/j.ssc.2007.02.013 
[44] 
Z. Ye, A. OteroDeLaRoza, E. R. Johnson, and A. Martini, Nanotechnology
26, 165701(2015).
DOI:10.1088/09574484/26/16/165701 
[45] 
M. NeekAmal, 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 
b. 中国科学技术大学化学与材料科学学院，合肥微尺度物质科学国家研究中心，能源材料化学协同创新中心，中国科学院纳米科学卓越创新中心，合肥 230026