Stacking in incommensurate graphene/hexagonal-boron-nitride heterostructures based on ab initio study of interlayer interaction

The interlayer interaction in graphene/boron-nitride heterostructures is studied using density functional theory calculations with the correction for van der Waals interactions. It is shown that the use of the experimental interlayer distance allows to describe the potential energy surface at the level of more accurate but expensive computational methods. On the other hand, it is also demonstrated that the dependence of the interlayer interaction energy on the relative in-plane position of the layers can be fitted with high accuracy by a simple expression determined by the system symmetry. The use of only two independent parameters in such an approximation suggests that various physical properties of flat graphene/boron-nitride systems are interrelated and can be expressed through these two parameters. Here we estimate some of the corresponding physical properties that can be accessed experimentally, including the correction to the period of the Moir\'{e} superstructure for the highly incommensurate ground state of graphene/boron-nitride bilayer coming from the interlayer interaction, width of stacking dislocations in slightly incommensurate systems of boron nitride on stretched graphene and shear mode frequencies for commensurate graphene/boron-nitride systems, such as a flake on a layer. We propose that the commensurate-incommensurate phase transition can be observed in boron nitride on stretched graphene and experimental measurements of the corresponding critical strain can be also used to get an insight into graphene/boron-nitride interactions.


I. INTRODUCTION
Following the discovery of 2D analogues of graphene a new research field has emerged on layer-by-layer design of van der Waals heterostructures 1,2 opening fresh possibilties to observe unusual properties and physical phenomena. Absolute leaders among such heterostructures are those based on graphene and hexagonal boron nitride (h-BN). Chemically inert, flat and wide gap boron nitride substrate allows to avoid graphene rippling, supress the charge inhomogeniety and improve the carrier mobility providing a quality comparable to suspended graphene. [3][4][5][6] The advanced performance of graphene/boron-nitride interfaces has been already utilized in a number of nanoelectronic devices where fewlayer and monolayer boron nitride serves as a gate dielectic 7-10 or tunnel barrier, 11,12 including novel fieldeffect tunneling transistors. 13 Nanocapacitors based on graphene/hexagonal-boron-nitride heterostructures have been proposed theoretically 14 and fabricated 15 recently. Multi-layer graphene intercalated by boron nitride has been also suggested for the use in ultra-scaled interconnects of integrated circuits. 16 Double-layer graphene with an ultrathin boron nitride spacer has been considered to study strong Coulomb drag 17 and tunable metal-insulator transition. 8 It has been also predicted that high-temperature superfluidity can be observed in such a heterostructure. 18 Though hexagonal boron nitride and graphene have very close lattice constants, the small lattice mismatch of 1.8% turns out sufficient for creation of the Moiré superstructure even when the layers are perfectly aligned. [19][20][21] A weak periodic potential of this superstructure perturbs the electronic spectrum of graphene and leads to emergence of superlattice minibands manifested through the Hofstadter's butterfly and Dirac points near the edges of the superlattice Brillouin zone. [22][23][24][25][26] The electronic and optical properties of graphene 27-32 on boron nitride modified due to the formation of the Moiré superstructure can serve as a basis for development of new technologies.
While the interlayer interaction causes significant adjustment of aligned graphene and boron nitride layers and commensurate domains arise, [19][20][21]33 relative rotation of the layers brings the system to the fully incommensurate state in which the interlayer interaction landscape is extremely smooth. Drastic changes in the structure, electronic and optical properties of the system are observed upon such a rotation. 19,30 A new phenomenon of macroscopic self-reorientation of graphene towards crystallographic directions on the underlying boron nitride crystal has been recently demonstrated. 34 On the other hand, robust superlubricity that can be used to reduce the friction in nanoelectromechanical systems has been predicted for interfaces between graphite and boron nitride bulk. 35 As long as energetic characteristics of interlayer interaction in graphene/boron-nitride heterostructures have not been measured yet experimentally, theoretical studies hold the key to understanding the properties related to the interlayer interaction and elaboration of nanodevices based on these properties. Here we present an approach that allows accurate and detailed calculations of the dependence of interlayer interaction energy on the relative position of the layers in graphene/boron-nitride heterostructures and consider a wide set of experimentally measurable physical properties related to this interaction.
To analyze the interaction of graphene and boron nitride layers responsible for the phenomena mentioned above first-principles studies have been performed. 16,32,33,[36][37][38][39][40] Though density functional theory (DFT) calculations 33,36,39,40 predict the qualitatively correct dependence of the interlayer interaction energy on stacking of the layers, they give wrong magnitudes of relative energies of states with different stacking as can be deduced from comparison with the more accurate but computationally expensive random phase approximation (RPA) approach. 32,40 This failure of the DFT method is associated with its inability to describe the equilibium interlayer distance since, as known from publications for bilayer graphene, 41,42 the relative energies of states with different stacking decrease exponentially in magnitude upon increasing the interlayer distance. In the present paper we perform DFT calculations of the potential surface of interlayer interaction energy in graphene/boron-nitride heterostructures, i.e. the dependence of the interlayer interaction energy on the relative in-plane displacement of the layers, with the help of the vdW-DF2 functional. 43 We show that the use of the experimental interlayer distance instead of the optimized one allows to get the results close to the data obtained in the more accurate RPA approach 32,40 at a lower computational cost.
Since first-principles methods are restricted to small simulation times and length scales, semiemirical models able to reproduce the potential surface of interlayer interaction energy are invoked for efficient modeling of phenomena taking place in heterostructures. Several such models have been proposed for the interaction between graphene and boron nitride layers including the registrydependent 38 and Morse-type 33 interatomic potentials, registry index model 35 and approximations based on the first spatial Fourier harmonics. 31,40 Atomistic models using semiempirical potentials were employed to analyze structure 33,38,44 and van der Waals interactions 45 in Moiré patterns of graphene/boron-nitride systems. Superlubricity of heterostructures was studied qualitatively on the basis of the simple registry index model. 35,39 The approximations of the dependence of the interlayer interaction energy on the relative in-plane position of graphene and boron nitride layers by the first spatial Fourier harmonics were applied to analyze structural relaxation, tribological behavior and band gap landscapes. 28,31,40 However, the adequacy of such approximations has not been addressed. Here we examine the deviation of such an expression determined by the system symmetry from the potential energy surface obtained by the DFT calculations at the experimental interlayer distance.
The approximation of the potential energy surface at the fixed interlayer distance by the first spatial Fourier harmonics includes only two independent parameters and thus a number of properties of flat graphene/boron-nitride systems related to the interlayer interaction are determined to a large extent by these two parameters. A similar observation has been made previously for pure graphene and boron nitride systems 41,46,47 and the barrier to relative rotation of the layers, shear mode frequency, width of stacking dislocations and critical strain for the commensurate-incommensurate phase transition have been estimated for such materials. Here, in addition to these properties, we use the approximation of the potential energy surface to analyze the Moiré superstructure of the incommensurate ground-state graphene/boron-nitride heterostructure and to study how the interlayer interaction affects the superstructure period. The possibility of experimental measurements of the calculated properties is discussed.
The paper is organized in the following way. In section II we present the results of our DFT calculations of the potential energy surface of graphene/boron-nitride heterostructures. In section III this surface is approximated using the first spatial Fourier harmonics and the accuracy of this approximation is studied. In section IV the approximation is applied to evaluate physical properties related to the interlayer interaction for flat heterostructures with a different degree of incommensurability including the barrier to relative rotation of the layers, shear mode frequency, period of the Moiré pattern, width of stacking dislocations and critical strain for the commensurate-incommensurate phase transition. Finally conclusions are summarized.

II. DFT CALCULATIONS
The DFT calculations are performed using the nonlocal vdW-DF2 functional 43 taking into account van der Waals interactions as implemeted in the VASP code. 48 The projector augmented-wave method (PAW) 49 is applied. The rectangular unit cell including 4 atoms of each layer and having height of 20Å is considered under periodic boundary conditions. Integration over the Brillouin zone is performed using the Monkhorst-Pack method 50 with the 28 × 36 × 1 k-point grid (in the armchair and zigzag directions, respectively). The maximum kinetic energy of plane waves is 600 eV. The convergence threshold of the self-consistent field is 10 −7 eV. These parameters provide well-converged values of the barriers to relative sliding of graphene layers. 41 The calculations of the optimal size of the commensurate unit cell of the graphene/boron-nitride heterostructure and potential energy surface are performed at a fixed interlayer distance d = 3.33Å, which is close to the experimentally measured interlayer distances in bulk boron nitride 51-59 and graphite. 53,[60][61][62][63][64][65] It has been shown previously 47,66 that our appoach based on the combination of the vdW-DF2 functional 43 and the use of the experimental interlayer distance allows to describe adequately such properties of purely graphene or boron ni-tride systems related to in-plane relative motion of the layers as the shear mode frequency, shear modulus, barrier to relative sliding of the layers and width of stacking dislocations. The DFT calculations using the experimental interlayer distance were found to be much more accurate than the ones using the optimized interlayer distance, while at the experimental interlayer distance the vdW-DF2 functional was demonstrated to perform better than the PBE-D2, PBE-D3, PBE-D3(BJ), PBE-TS, PBE-TS/HI, PBE-TS+SCS, optPBE-vdW functionals. 66 The calculations show that the lowest energy stacking for the commensurate unit cell of the graphene/boronnitride heterostructure is AB1 in which boron atoms are located on top of carbon atoms and nitrogen atoms are on top of centers of hexagons, in agreement with previous findings 16,32,33,[36][37][38][39][40] (Fig. 1 The potential energy surface of the graphene/boronnitride heterostructure calculated for the optimized commensurate unit cell (Fig. 1a) is in qualitative agreement with the previous studies. 32,33,36,[38][39][40] Two types of maxima on the potential energy surface correspond to the AA and AB2 stackings (Fig. 1b). In the first of these stackings, all boron and nitrogen atoms are on top of carbon atoms and the energy of this stacking relative the ground-state AB1 stacking is 12.35 meV/atom, which is within the range of the values obtained at the optimized interlayer distance using different methods of 5.63 meV/atom (DFT-TS), 39 ∼5.4 meV/atom (vdW-DF2), 33 10.86 meV/atom (HSE+MBD), 38 15.68 meV/atom (DFT-D2), 5.35 meV/atom (vdW-DF2), 9.92 meV/atom (RPA) 40 and 10.5 meV/atom (RPA) 32 (note that all energies in the present paper are given in meV per atom in the upper (adsorbed) layer). The second smaller maximum with the relative energy 9.73 meV/atom corresponds to the AB2 stacking in which nitrogen atoms are located on top of carbon atoms and boron atoms are on top of centers of hexagons. For comparison, the liter-ature values for the relative energy of the AB2 stacking are 4.5 meV/atom (DFT-TS), 39 ∼4.6 meV/atom (vdW-DF2), 33 12.96 meV/atom (DFT-D2), 4.68 meV/atom (vdW-DF2), 8.67 meV/atom (RPA) 40 and 9 meV/atom (RPA). 32 The saddle point (SP) for the transition between adjacent energy minima (Fig. 1a) lies on the straight line in the armchair direction connecting the AA and AB2 stackings on the potential energy surface, 0.427Å away from the AB2 stacking and 1.015Å away from the AA stacking (Fig. 1b). The relative energy of this SP stacking corresponding to the barrier to relative in-plane motion of the graphene and boron nitride layers is 9.46 meV/atom, i.e. only 0.27 meV/atom smaller than the relative energy of the AB2 stacking. The previously reported values of the barrier obtained at the optimized interlayer distance include 3.9 meV/atom (DFT-TS), 39 ∼4.4 meV/atom (vdW-DF2), 33 ∼13 meV/atom (DFT-D2), ∼4.5 meV/atom (vdW-DF2) and ∼8.5 meV/atom (RPA). 40 It should be noted that the calculated barrier for the graphene/boron-nitride heterostructure is several times greater than the barriers to relative sliding of two graphene layers or two boron nitride layers. The reported barriers to relative sliding of graphene layers range from 0.5 to 2.4 meV/atom (Refs. 40-42, 66, 67, 70-72). There are also estimates of this barrier from the experimental measurements of the shear mode frequency and width of dislocations in few-layer graphene of 1.7 meV/atom (Ref. , respectively, i.e. six and two and a half times smaller than the value for the graphene/boron-nitride heterostructure.
The comparison of the relative energies of symmetric stackings for the graphene/boron-nitride heterostructure obtained here and in previous papers 32,39,40 shows that the approach applied in the present paper, i.e. the combination of the vdW-DF2 functional with the use of the experimental interlayer distance in DFT calculations, overestimates the characteristics of the potential surface of interlayer interaction energy by 10% -20% with respect to more accurate RPA calculations. 32,40 This error can be partially attributed to the fixed interlayer distance in our study since according to the RPA calculations, 32,40 the optimal interlayer distance varies by 0.2Å for different stackings. Nevertheless, this error is small compared, for example, to the scatter of the experimental data on the total binding energy of graphene layers [75][76][77][78] or more than 50% deviation of the relative energies of symmetric stackings following from DFT calculations at the optimized interlayer distance 39,40 from the RPA data. These observations are in agreement with the previous comparison of performance of different methods for pure graphene and boron nitride. 66 Therefore, the combination of the vdW-DF2 functional with the use of the experimental interlayer distance in standard DFT calculations is a computationally cheap but sufficiently accurate alternative to expensive methods, such as RPA.

III. APPROXIMATION OF POTENTIAL ENERGY SURFACE
The possibility to accurately approximate potential energy surfaces by expressions containing only the first spatial Fourier harmonics determined by the system symmetry has been previously demonstrated for various types of interacting layers, including purely graphene [40][41][42]46,72,79 and boron nitride 40,47 systems as well as carbon nanotube walls in the case when the walls are infinite and commensurate [80][81][82][83][84] and when the corrugations of the potential energy surface are determined by edges 85 or defects. 80 An expression for the potential surface of interaction energy of graphene and boron nitride layers including 5 fitting parameters has been also proposed. 40 In the present paper we show that this potential energy surface can be approximated using only two independent fitting parameters. A similar approximation of the potential energy surface with the first spatial Fourier components has been used for studies of strain distribution and band gap opening in a graphene layer on the boron nitride crystal. 28,31 However, the adequacy of this approximation has not been addressed.
We consider the graphene layer as adsorbed on the boron nitride one and use that the potential energy surface of an atom adsorbed on a 2D trigonal lattice can be approximated by the first Fourier harmonics as 86 where x and y axes are chosen in the armchair and zigzag directions, respectively, k x = 2π/( √ 3a), k y = 2π/a, a is the lattice constant, u describes the relative position of the atom with respect to the 2D lattice and point u = 0 corresponds to the case when the atom is located on top of one of the lattice atoms. In the case of graphene/boron-nitride heterostructures, we should sum up interactions of carbon atoms with the boron (CB) and nitrogen (CN) lattices. Thus, we get where u describes the relative position of the graphene and boron nitride layers, point u = 0 corresponds to the AA stacking and U 0,tot = 2(U 0CB + U 0CN ). It is seen that corrugations of the potential energy surface in flat graphene/boron-nitride heterostructures are described by two parameters U 1CB and U 1CN , which determine all physical properties related to the interlayer interaction, such as the barrier to relative sliding of the layers, shear modulus and shear mode frequency. 46,47,66,67,79 The parameters of the approximation can be found from the relative energies of the symmetric stackings AA, AB1 and AB2 as The values obtained are U 1CB = 0.58056 meV/atom and U 1CN = 2.7436 meV/atom. The greater value of U 1CN compared to U 1CB can be explained by the stronger repulsion between carbon atoms and negatively charged nitrogen ions compared to positively charged boron ions. The standard deviation of the approximated potential energy surface from the one obtained by the DFT calculations is only 0.0305 meV/atom (Fig. 1b,c and 2), which is within 0.25% of the energy difference between the AA and AB1 stackings. The largest deviations of the approximation from the DFT results of up to 0.07 meV/atom are observed for the regions around the saddle-point stacking SP (Fig. 2).

IV. PROPERTIES OF GRAPHENE/BORON-NITRIDE HETEROSTRUCTURES
In this section we consider graphene/boron-nitride systems characterized by different stacking patterns and degree of incommensurability. Using the approach introduced by Porovskiȋ and Talapov in their pioneer paper 87 on the commensurate-incommensurate phase transition and used for description of such a transition in multilayer films on surfaces (see, for example, Ref. 88) we distinguish the following states: (1) commensurate system in which the layers have the same lattice constant and are in the AB1 stacking corresponding to the minimal interlayer interaction energy, (2) slightly incommensurate system in which there is a small lattice mismatch of the layers leading to formation of large commensurate domains separated by narrow incommensurate boundaries corresponding to stacking dislocations, (3) highly incommensurate system in which there is a considerable mismatch in the lattice constants of the layers but the positions of atoms in the layers are perturbed by the interlayer interaction and (4) fully incommensurate system in which there is a considerable mismatch in the lattice constants of the layers and the interlayer interaction virturally does not disturb the positions of atoms in the layers. The ground state of the graphene/boron-nitride bilayer with the aligned layers is the highly incommensurate system (3) with the period of the Moiré pattern comparable to the width of boundaries between commensurate domains. [19][20][21]33 The relative rotation of the graphene and boron nitride layers leads to the transition to state (4), in which the potential energy surface is extremely smooth and no deformations of the layers due to the interlayer interactions are observed. 19 A commensurate system (1) can be realized for a flake on a large substrate layer when the flake size is smaller than the period of the Moiré pattern or in the case when the graphene layer is stretched so that its lattice constant reaches the one for boron nitride. Biaxial stretching is required for bilayer geometry, while uniaxial stretching is sufficient for a ribbon of the width smaller than the period of the Moiré pattern on a wide substrate layer or on another ribbon. A slightly incommensurate system (2) can be achieved in this case when the lattice constant of graphene is slightly different from the one for boron nitride.
Our model of the interlayer interaction energy described by Eq. 2 does not take into account changes in the interlayer distance and thus we limit our consideration to the case of flat systems. The experimental measurements gave the height variation for supported graphene/boron-nitride systems of about 0.4Å (Ref. 21). This is comparable to the 0.2Å variation of the optimal interlayer distance for different stackings according to the RPA calculations. 32,40 Therefore, the deviation of 10% -20% for the relative energies of symmetric stackings obtained here at the fixed interlayer distance from the RPA data 32,40 for the variable interlayer distance (see section II) provides an estimate of the accuracy of our model for supported graphene/boron-nitride bilayers. The varia-tion in the interlayer distance should be also insignificant for commensurate or fully incommensurate systems. On the other hand, very strong surface corrugation with the characteristic amplitude of 8Å was predicted for freestanding graphene/boron-nitride bilayers on the basis of atomistic simulations. 33,38 To be able to describe quantitatively the properties of such corrugated structures the model should be supplemented by the terms describing the dependence of the interlayer interaction energy on the interlayer distance similar to approximations of the potential energy surface proposed for graphene/boronnitride heterostructures 31,40 and graphene bilayer 42 and interlayer atomic potentials 41,70 for graphene. Nevertheless, even in the present form the model can still provide a qualitative insight into the phenomena coming from the interlayer interaction and serve as an important step towards development of more complicated models with account of variations in the interlayer distance.
The possibility to describe the potential energy surface of graphene/boron-nitride heterostructures with only two independent parameters within Eq. 2 means that a large number of properties of such systems are interrelated and expressed through these two parameters, similar to the cases of pure graphene and boron nitride systems. 41,46,47 Below we estimate the following physical properties of systems of different incommensurability. First the barrier to relative rotation of graphene and boron nitride layers, which corresponds to the transition from state (1) to state (4), is obtained for commensurate flakes on a large substrate layer. The shear mode frequencies are evaluated for various commensurate systems (1). The period of the Moiré pattern is analyzed for the highly incommensurate state (3) corresponding to the ground state of the graphene/boron-nitride bilayer with the aligned layers. The phase transition from the commensurate state (1) to the slightly incommensurate state (2) with a low density of stacking dislocations and characteristics of these dislocations are studied in the last part of this section.

A. Barrier to rotation
We start consideration of physical properties of graphene/boron-nitride systems from the barrier to relative rotation of the layers from the commensurate state. As mentioned above, though the ground state of the graphene/boron-nitride bilayer is incommensurate and it is difficult to realize the rotation in macroscopic systems, such a rotation is relevant for commensurate flakes on a large substrate layer. From studies for graphene it is known that superlubric behavior [89][90][91][92] and diffusion 93,94 of flakes on a periodic substrate occurs via rotation, which brings the system to the fully incommensurate state 19 with an extremely smooth potential energy surface. Correspondingly, the diffusion coefficient of flakes is determined by the barrier to rotation from the commensurate state to the fully incommensurate one. 93,94 The potential energy of the latter can be found as the aver-age interlayer interaction energy given by Eq. 2 Therefore, the barrier to relative rotation of the graphene and boron nitride layers from the commensurate state with the AB1 stacking can be estimated as U in =Ū − U (AB1) = −3U 1CB /2 + 3U 1CN = 7.36 meV/atom. This value is in agreement with the estimate of the relative energy of the fully incommensurate state from the RPA calculations 32 of ∼7 meV/atom and is somewhat greater than the reported barriers to relative rotation of two graphene layers of ∼5 meV/atom (Ref. 46

B. Shear mode frequency
One of the experimentally measurable quantities that can be used to probe the potential surface of interlayer interaction energy is the frequency of the shear mode in which commensurate layers rigidly slide with respect to each other parallel to the plane. 41,46,47,66,79 To estimate shear mode frequencies for commensurate systems of the graphene/boron-nitride bilayer, a graphene flake on a boron nitride layer and a boron nitride flake on a graphene layer the calculations of the curvature ∂ 2 U/∂x 2 of the potential energy surface around the AB1 minimum are performed for the commensurate unit cells with different lattice constants corresponding to the optimized lattice constants of the bilayer and single graphene and boron nitride layers. The graphene and boron nitride layers are rigidly displaced within 0.05Å in the armchair direction and the obtained energy curves are approximated by parabolas.
The calculations show that the curvatures of the potential energy surface ∂ 2 U/∂x 2 for the lattice constants a BN and a C corresponding to the isolated boron nitride and graphene layers differ only by 1.9% and 2.6%, respectively, from the value for the lattice constant a of the bilayer. The curvature of the potential energy surface for the bilayer can be also estimated from Eq. 2 as ∂ 2 U/∂x 2 = 4π 2 (2U 1CN − U 1CB )/a 2 = 0.031 eV/Å. This estimate is only 1.3% greater than the value obtained directly from the DFT calculations.
The shear mode frequencies f E are then found as 46,47,66 Here for the graphene flake on the boron nitride layer  Table I. As seen from this table, the results for the graphene flake on the boron nitride layer and the boron nitride flake on the graphene layer are virtually indistinguishable due to the close reduced masses and curvatures of the potential energy surface. It should be also noted that the estimated shear mode frequency for the graphene/boron-nitride bilayer of 37 cm −1 is somewhat higher than the frequencies for bilayer graphene 66 of 29 cm −1 and boron nitride 47,66 of 34 cm −1 obtained within the same computational approach. This value is also greater than the results of experimental measurements for bilayer graphene of 28 ± 3 cm −1 (Ref. 95) and 32 cm −1 (Ref. 96).

C. Moiré pattern
Let us now obtain the correction to the period of the Moiré pattern in the incommensurate ground state of graphene/boron-nitride heterostructures induced by the interlayer interaction. First we consider the relationship between the elastic and interlayer interaction energies in the fully incommensurate and commensurate states of the aligned graphene and boron nitride layers. Maintaining the structure of the layers the same as in the absence of the interlayer interaction is related to excess in the average density of the interlayer interaction energy compared to the commensurate state of w int = U in /σ = 44 mJ/m 2 , where σ = √ 3a 2 /4 = 2.7Å 2 is the area per one carbon atom. Transformation of the incommensurate structure to the commensurate one with the uniformly stretched graphene layer and compressed boron nitride layer, on the other hand, is associated with the penalty in the elastic energy of w el = kδ 2 /(1 − ν), where δ = (a BN − a C )/a = 1.78% is the relative lattice mismatch, k is the reduced elastic constant of the layers and ν is the Poisson ratio. In the present paper we use the values of the elastic constants of graphene and boron nitride layers k C = 331 J/m 2 and k BN = 273 J/m 2 calculated previously 67 using the same computational approach. Although slightly different values of the Poisson ratio were obtained for graphene and boron nitride, ν C = 0.174 and ν BN = 0.201, in this section we assume that the Poisson ratios of the layers are close and approximately equal to that for graphene ν = ν C (in section IV D the difference in the Poisson ratios of the graphene and boron nitride layers is taken into account). In the case when both of the graphene and boron nitride layers are free to relax in the plane (this can be achieved by using an incommensurate substrate such as a rotated graphene or boron nitride layer), the reduced elastic constant is k = k C k BN /(k C + k BN ) = 150 J/m 2 and the density of the elastic energy in the commensurate state w el = 56 mJ/m 2 . Comparison of energies w el and w int suggests that in the ground state the heterostructure should be incommensurate. For a free graphene layer on the boron nitride crystal (similar to the experimental studies [19][20][21], the interface boron nitride layer can be considered as nearly rigid. 32,44 In this case the reduced elastic constant of the layers k ≈ k C is approximately twice greater, the density of the elastic energy in the commensurate state is w el = 124 mJ/m 2 and the incommensurate state is even more preferred. The same behavior with w el = 102 mJ/m 2 should also be observed for a free boron nitride layer on the fixed graphite substrate.
Incommensurability of the graphene/boron-nitride systems is manifested through the formation of Moiré superstructures. [19][20][21]33 In the case when the change in the bond lengths of the layers due to the interlayer interaction is neglected and the layers are completely aligned we can estimate that the period of the Moiré pattern is L 0 = a/δ = 14.2 nm, in agreement with the experimental measurements of ∼14 nm 19,20 69 Nevertheless, since the characteristic interlayer interlayer interaction energy w int is comparable to the characteristic elastic energy w el , significant modulation of positions of atoms in the interacting layers is possible. Such a modulation has been previously observed in atomistic simulations 33,38,44 and calculations using continuum models. 28,31 However, the correction to the period of the Moiré pattern induced by the interlayer interaction has not been considered explicitly. Below we quantify this correction for aligned graphene and boron nitride layers.
In heterostructures formed by uniformly deformed graphene and boron nitride layers (with the constant in-plane strains) the relative displacement of atoms in the layers varies uniformly with the absolute position of atoms within the layers u 0,x = xa/L, u 0,y = ya/L, where L is the period of the Moiré pattern, axes x and y are chosen in the armchair and zigzag directions and the point x = 0, y = 0 corresponds to the AA stacking following the notations from sections II and III. The Moiré pattern has the same hexagonal symmetry as the potential energy surface, which is imposed by the geometry of the layers, with the only difference that the unit cell of the Moiré pattern is scaled by a factor of L/a as compared to the unit cell of the potential energy surface. Obviously this symmetry is not changed by the interlayer interaction. However, the interlayer interaction modifies the period of the Moiré pattern L = L 0 + ∆L and the relative displacement of atoms in the layers u = u 0 + ∆ u. Here we assume that the perturbations ∆L and ∆ u are small so that ∆L L 0 and ∆u a. Since the symmetry of the pattern is not affected by the interlayer interaction, the correction ∆ u to the relative displacement of atoms should vanish at the symmetry planes. This is automatically satisfied for the correction of the form i.e. in the case when the correction is proportional to the force of interlayer interaction. The interlayer interaction energy can then be expanded as It is sufficient to stop here at the first-order terms in α as it can be easily checked that the second-order terms vanish upon the integration over the unit cell of the Moiré superstructure and do not contribute to the total energy. Integrating Eq. 7 over the unit cell of the Moiré superstructure, we get the average interlayer interaction energy per unit area where 1/2 (9) and this parameter corresponds to the characteristic corrugation of the potential energy surface. Tensile ( ) and shear (τ ) strains associated with the corrections to the period of the Moiré pattern and relative displacement of the layers are given by Using notation for the total strain the elastic energy of the heterostructure can be written as The average elastic energy per unit area is therefore given by The total energy of the heterostructure ∆w = W el + U −Ū relative to the incommensurate state with no in-plane deformation of the layers due to the interlayer interaction is minimal for Let us now check that the resulting correction to the relative displacement of atoms in the layers is actually much smaller than the lattice constant ∆u a and thus the perturbation approach applied here is valid. From Eqs. 6 and 14 and relation L ∼ a/δ it can be estimated that the characteristic magnitude of the correction to the relative displacement of atoms of the graphene and boron nitride layers is This means that the approach used here makes sense as long as the characteristic magnitude of corrugation of the potential energy surface U 1 is less or comparable to the characteristic elastic energy kδ 2 . In our case U 1 = 15 mJ/m 2 , kδ 2 = 46 mJ/m 2 when both of the graphene and boron nitride layers are free to relax in the plane, 103 mJ/m 2 for the graphene layer on the fixed boron nitride substrate and 85 mJ/m 2 for the boron nitride layer on the fixed graphite substrate, which ensures ∆u a and adequacy of the approach used. One of the consequences of the perturbation approach used is that in the limit of L ≈ L 0 the strain in each free layer does not depend on whether the second layer is free or fixed. This is clear from Eq. 14 as the coefficient α, which determines the relative displacement ∆u of atoms in the layers, is inversely proportional to the reduced elastic constant k of the system. In the case when both of the layers are free to relax in the plane, absolute displacements of atoms in the graphene and boron nitride layers correspond to the fractions k BN /(k C +k BN ) and k C /(k C + k BN ) of the relative displacement ∆u, respectively. As a result, the same displacement as on the fixed substrate with k = k C or k = k BN , respectively, is obtained.
The distribution patterns of strain in free layers for L ≈ L 0 are shown in Fig. 3 in units of 0 = −4π 2 αU 1 /L 2 0 = 1.94 · 10 −3 for the graphene layer and 2.35·10 −3 for the boron nitride layer. Such a strain distribution is qualitatively similar to the distribution of bond lengths observed in atomistic simulations 33,38,44 and elastic energy maps obtained within continuum models. 28,31 As can be expected, the maximal strains of about 3 0 are observed in the regions with the AB1 stacking, where the layers tend to be commensurate in order to minimize the interlayer interaction energy (Fig. 3). These maximal strains of 0.6 -0.7% are in agreement with variations of the bond length of about 1% measured experimentally for the graphene layer on the boron nitride crystal. 19 The minimal strains of about 0.6 0 are achieved for the stacking that is intermediate between the AA and AB1 stackings. Adjusting the period L of the Moiré superstructure allows to reduce the strain though the¯ term (see Eq. 10).
With account of Eq. 14, the energy of the heterostructure relative to the incommensurate state with no inplane deformation of the layers due to the interlayer interaction takes the form The optimal change of the period of the Moiré superstructure is then roughly When one of the graphene or boron nitride layers is fixed, the correction to the period of the Moiré superstructure is rather small. For the fixed graphene layer, ∆L = 0.32 nm, i.e. 2.3% of the period L 0 of the Moiré superstructure without account of in-plane deformation of the layers due to the interlayer interaction. For the fixed boron nitride layer, ∆L = 0.21 nm, i.e. 1.5% of L 0 .
In the latter case, the estimated period of the Moiré superstructure is L = 14.4 nm. The same as L 0 , this is close to the experimentally measured period of the Moiré superstructure for the aligned graphene layer on the boron nitride crystal of about 14 nm 19,20 and 15±1 nm. 21 However, these experimental data are insufficient to confirm the effect of the interlayer interaction. A more significant correction to the period of the Moiré superstructure of ∆L = 1.30 nm, i.e. 9.2% of L 0 , with the resulting period of L = 15.5 nm, is expected when both of the graphene and boron nitride layers are free to relax in the plane, which can be the case for the graphene/boron-nitride bilayer on an incommensurate substrate, for example, on a rotated layer of graphene or boron nitride. Experimental measurements for such a system could help to distinguish the effect of the interlayer interaction.

D. Stacking dislocations in slightly incommensurate systems
Let us now consider stacking dislocations in slightly incommensurate graphene/boron-nitride systems. Stretching graphene so that its lattice constant increases up to the one for boron nitride a BN results in formation of the commensurate structure with the minimal interlayer in-teraction energy. Futher stretching or compressing the graphene layer should lead to a competition in the elastic energy of the free boron nitride layer and interlayer interaction energy resolved through the transition to the slightly incommesurate state with stacking dislocations. In other words, a commensurate-incommensurate phase transition 87 should take place at some critical strain associated with the difference in the average lattice constants of the graphene and boron nitride layers.
As discussed in the beginning of section IV, a slightly incommensurate system can be realized by biaxial stretching of the graphene layer for bilayer geometry or by uniaxial stretching if the boron nitride layer is a ribbon of the width smaller than the period of the Moiré superstructure. In the latter case stretching should take place along the ribbon axis and we disregard the effect of the ribbon edges on the potential surface of interlayer interaction energy, which is reasonable if the ribbon width is much greater than the lattice constant and the edges are properly terminated, e.g. by hydrogen. We also assume that the density of stacking dislocations is low so that the interaction between the dislocations can be neglected. These assumptions allow us to use the formalism of the two-chain Frenkel-Kontorova model, 97 in which two layers are represented as chains of particles connected by harmonic springs and coupled by van der Waals interactions. The model has been already applied to study the commensurate-incommensurate phase transition in double-walled carbon nanotubes, 83,97 bilayer graphene 67,98 and boron nitride 47,67 as well as edge stacking dislocations in bilayer graphene. 99 Following the approach from papers, 47,67,98 it is first necessary to choose the adequate approximation of the dislocation path, i.e. the curve on the potential energy surface described by the dependence of relative displacement u of the layers on the coordinate in the direction perpendicular to the boundary between commensurate domains that minimizes the formation energy of dislocations. In the case of graphene/boron-nitride heterostructures, the minimum energy path between the AB1 minima passing through the SP stacking only slightly deviates from the straight line between adjacent energy minima in the zigzag direction. The barrier along this straight path (corresponding to the SP' stacking) is 9.81 meV/atom (Fig. 1c). This is only 4% higher than the barrier along the minimum energy path through the SP stacking and within the accuracy of the DFT calculations with corrections for van der Waals interactions. The SP' stacking is obtained from the SP stacking by shifting the layers by only 0.29Å. Therefore, it can be safely assumed that the dislocation path in graphene/boron-nitride heterostructures lies along the straight line between adjacent energy minima AB1.
On the basis of Eq. 2, the relative potential energy along the straight dislocation path is then approximated as where k 0 = 2π/b, u is the relative displacement of the layers along the dislocation path, which changes from 0 to the magnitude of the Burgers vector b = a BN , and V max is the barrier per unit area. The barrier along this path calculated from Eq. 2 is U max = 2(2U 1CN −U 1CB ) = 9.81 meV/atom, which is exactly equal to our DFT value. The variation of the barrier due to changes in the lattice constant from a to a BN lies within 2% and is neglected. Therefore, we estimate the barrier per unit area as V max = 4U max /( √ 3a 2 BN ) = 0.0571 J/m 2 . The dislocation width depends on the angle β between the Burgers vector and the normal to the boundary between commensurate domains 47,67 and is given by where K(β) = E cos 2 β + G sin 2 β describes the dependence of the reduced elastic constant of the layers on fractions of tensile and shear character in the dislocation, E and G are the reduced tensile and shear elastic constants. When the boron nitride and graphene layers are of similar width, both of the layers participate in the formation of the dislocation and the reduced constants are given by , respectively. For the boron nitride ribbon on the wide graphene layer, the dislocation appears only in the boron nitride layer and, correspondingly, E ≈ E BN and G ≈ G BN .
As seen from Fig. 4a, the dislocation width l D decreases with increasing the angle β, i.e. changing the dislocation character from tensile to shear. The estimated dislocation widths l D in the case of similar widths of the graphene and boron nitride layers are 6 -9 nm, which is much smaller then the width of full dislocations in bilayer boron nitride of 12 -16 nm (Ref. 47 and 67) evaluated within the same approach and only a little smaller than the corresponding results for partial dislocations in bilayer graphene of 8 -14 nm (Ref. 67 and 98). These values also lie within the range of experimental data on dislocation widths in bilayer graphene of 6 -11 nm (Ref. 73, 100, and 101). While the barrier V max along the dislocation path in graphene/boron-nitride heterostructures is several times greater than that in bilayer graphene, the length of this path is also greater (the dislocation path length is equal to the lattice constant a BN in the heterostructure and to the bond length l C ≈ a BN / √ 3 in bilayer graphene) and the resulting dislocation widths are close.
The dislocation energy W 0 , i.e. the energy of the state with one stacking dislocation relative to the commensurate state with equal lattice constants of the layers per unit length of the boundary between commensurate domains, follows the same dependence on the angle β as the dislocation width W 0 (β) = 2b 2K(β)V max /π (Refs. 47 and 67). For the system of two layers of similar width this energy ranges from 0.42 meV/Å for tensile dilsocations to 0.27 meV/Å for shear dislocations, which is about four times greater than for bilayer graphene 67,98 and comparable to the estimates for bilayer boron nitride. 47,67 For the boron nitride ribbon on the wide graphene layer, the corresponding values are 0.57 meV/Å and 0.36 meV/Å, respectively.
Let us now consider the commensurateincommensurate phase transition in the system in which the average lattice constant of graphene in one of the directions (along the ribbon axis in the case of the boron nitride layer of the ribbon shape) is slightly smaller or larger than the average lattice constant of boron nitride. The commensurate phase should be observed for strains C in the graphene layer within the interval C c,− = δ − c < C < C c,+ = δ + c , where c is the critical strain associated with the average lattice mismatch of the graphene and boron nitride layers at which the formation energy of stacking dislocations goes to zero. The commensurate-incommensurate phase transition in similar systems of pure graphene and boron nitride with the uniaxially stretched substrate layer was considered in our previous paper. 67 In that paper we derived the expressions for c and the optimal angle β c between the boundary separating commensurate domains and the Burgers vector which minimizes the formation energy of stacking dislocations at the critical strain (or at any strain in the substrate layer if the free layer is a ribbon). According to these derivations, the optimal angle β c between the boundary separating commensurate domains and the Burgers vector depends on the angle α between the Burgers vector and the direction in which the average lattice constants of the layers are slightly different (Fig. 4b) as In the case of the graphene/boron-nitride systems the critical strain c associated with the average lattice mismatch between the graphene and boron nitride layers 67 is given by while the absolute critical strain in the graphene layer takes the values of C c,± (α) = δ ± c (α). As seen from Fig. 4c, the critical strain c for α = 0 • is about 0.98% when the layers are of similar width and 1.33% for the boron nitride ribbon on the wide graphene layer. The dependence of the critical strain on the angle α is rather weak for angles within 30 • . A shallow minimum in the critical strain c of about 0.96% is reached at α = 23.5 • for the layers of similar width. For the boron nitride ribbon on the wide graphene layer, the minimal c is about 1.31% and is reached at α = 21.3 • . With increasing α beyond 30 • the critical strain c grows fast and tends to infinity at α → 90 • . Such a behavior is explained by the fact that the formation energy of stacking disloca-tions cannot be reduced by stretching the substrate layer in the direction perpendicular to the Burgers vector.
As the Burgers vector is aligned along one of the zigzag directions, for a given geometry of the system, only six orientations of the Burgers vector are possible. Among these six types of possible dislocations, the dislocations with the angle α = α 0 , where 0 ≤ α 0 ≤ 30 • is the smallest angle between one of the zigzag directions and the direction in which the average lattice constants of the graphene and boron nitride layers are slightly different, are characterized by the lowest critical strain c and should appear first upon changing the lattice constant of the graphene layer relative to the boron nitride layer. Therefore, c (α 0 ) corresponds to the critical strain associated with the average lattice mismatch of the layers at which the commensurate-incommensurate phase transition takes place in graphene/boron-nitride heterostructures (Fig. 4c). It should be mentioned that this critical strain for the layers of similar width is only slightly greater than the corresponding value for bilayer boron nitride of 0.75% and almost three times greater than the critical strain for formation of the first dislocation in bilayer graphene of about 0.36% (Ref. 67).

V. CONCLUSIONS AND DISCUSSION
The potential energy surface of commensurate graphene and boron nitride layers has been calculated within the DFT approach based on the vdW-DF2 functional. 43 It is shown that such an approach combined with the use of the experimental interlayer distance provides reliable characteristics of the potential energy surface that are in agreement with more accurate but much more expensive computational methods, such as RPA. The calculated barrier to relative motion of graphene and boron nitride layers of 9.5 meV/atom is several times greater than the values for two graphene or boron nitride layers.
It has been checked that the approximation of the calculated potential energy surface by the first spatial Fourier components, an approach which is widely used in literature 28,31,40 to analyze electronic properties, is highly accurate. Such an approximation allows to describe the potential energy surface of flat graphene/boron-nitride heterostructures using only two independent parameters. This means that a large number of physical properties of graphene/boron-nitride systems are interrelated and are expressed through these two parameters as well. In the present paper we have used the approximation of the potential energy surface to estimate the shear mode frequency in the commensurate state, the barrier to relative rotation of graphene and boron nitride layers from the commensurate state to the fully incommensurate one, the width of stacking dislocations in slightly incommensurate systems and the period of the Moiré superstructure of the highly incommensurate ground state of the graphene/boron-nitride bilayer. It is shown that the interlayer interaction in the highly incommensurate ground state of the graphene/boronnitride systems with the layers aligned leads to modulation of atomic positions in the layers, which results in the increase of the period of Moiré superstructure by 9% in the case of free graphene and boron nitride layers on an incommensurate substrate and about 2% in the case when one of the layers is fixed. The maximal strains in the graphene and boron nitride layers are on the order of 0.6 -0.7%, in agreement with the experimental observations. 19 Commensurate or slightly incommensurate graphene/boron-nitride systems can be realized by stretching the graphene layer or restricting the system geometry to flakes of the size smaller than the period of the Moiré superstructure. The shear mode frequencies for the commensurate graphene/boron-nitride heterostructure, a graphene flake on a boron nitride layer and a boron nitride flake on a graphene layer have been calculated. The estimated barrier to relative rotation of the layers from the commesurate to fully incommensurate state of 7.4 meV/atom is comparable to the previously reported values for pure graphene or boron nitride systems.
The characteristrics of stacking dislocations have been studied for slightly incommensurate systems of boron nitride and graphene layers of similar width with the biaxially stretched graphene layer or a boron nitride ribbon of the width smaller than the period of the Moiré superstructure on a wide graphene layer stretched along the ribbon axis. The estimated dislocation widths for these two types of systems lie in the ranges of 6 -9 nm and 8 -13 nm, respectively, where the lower bound corresponds to shear dislocations and the higher one to tensile dislocations. It is suggested that changing the lattice con-stant of graphene with respect to boron nitride in one of the directions (in the system with the boron nitride ribbon along the ribbon axis) should result in observation of the commensurate-incommensurate phase transition. The estimated strain intervals for the graphene layer corresponding to the commensurate phase are about 0.8 -2.7% in the case when the graphene and boron nitride layers are of similar width and 0.4 -3.1% for the boron nitride ribbon on the wide graphene layer.
Let us now discuss the possibility of experimental measurements of the calculated properties related to the interaction of graphene and boron nitride layers. The estimated barrier to relative rotation of the layers is relevant for dynamics of flakes on periodic substrates, i.e. for such processes as superlubricity [89][90][91][92] and diffusion. 93,94 In the latter case the diffusion coefficient of flakes is exponentially dependent on this barrier. 93,94 The shear mode frequency in bilayers can be measured by Raman scattering 96 and coherent phonon spectroscopy, 95 as demonstrated for graphene. 95,96 The transmission electron microscopy 73,100 and scanning tunneling microscopy 101 allow direct measurements of the dislocation width by analogy with papers 73,100,101 on graphene. The experiments for heterostructures with the stretched graphene layer should make it possible to observe the commesurate-incommensurate phase transition. The period of the Moiré pattern for graphene on boron nitride crystal has been already determined by scanning tunneling microscopy and atomic force measurements. [19][20][21] The consideration of flat heterostructures with both free layers, which can be realized in the case of an incommensurate substrate, such as a rotated graphene or boron nitride layer, would allow to distinguish the change in the period of the Moiré pattern coming from the interlayer interaction. Such studies could provide a valuable experimental insight into the graphene/boron-nitride interlayer interaction due to the link between these properties and the potential energy surface established in the present paper.