Interlayer interaction and related properties of bilayer hexagonal boron nitride: ab initio study

Alexander V. Lebedev*a, Irina V. Lebedevab, Andrey A. Knizhnikac and Andrey M. Popovd
aKintech Lab Ltd., Moscow 123182, Russia. E-mail: allexandrleb@gmail.com
bNano-Bio Spectroscopy Group and ETSF Scientific Development Centre, Departamento de Física de Materiales, Universidad del País Vasco UPV/EHU, San Sebastian E-20018, Spain. E-mail: liv_ira@hotmail.com
cNational Research Centre “Kurchatov Institute”, Moscow 123182, Russia
dInstitute for Spectroscopy of Russian Academy of Sciences, Troitsk, Moscow 142190, Russia. E-mail: popov-isan@mail.ru

Received 8th October 2015 , Accepted 7th December 2015

First published on 8th December 2015


Abstract

The principal characteristics of interlayer interaction and relative motion of hexagonal boron nitride (h-BN) layers are investigated by the first-principles method taking into account van der Waals interactions. The dependence of the interlayer interaction energy on the relative translational displacement of h-BN layers (potential energy surfaces) are calculated for two relative orientations of the layers, namely, for the layers aligned in the same direction and in the opposite directions upon the relative rotation of the layers by 180 degrees. It is shown that the potential energy surfaces of bilayer h-BN can be approximated by the first Fourier components determined by symmetry. As a result, a wide set of physical quantities describing the relative motion of h-BN layers aligned in the same direction including barriers to their relative sliding and rotation, shear mode frequency and shear modulus are determined by a single parameter corresponding to the roughness of the potential energy surface, similar to bilayer graphene. The properties of h-BN layers aligned in the opposite directions are described by two such parameters. The possibility of partial and full dislocations in stacking of the layers is predicted for h-BN layers aligned in the same and opposite directions, respectively. The extended two-chain Frenkel–Kontorova model is used to estimate the width and formation energy of these dislocations on the basis of the calculated potential energy surfaces.


1 Introduction

Although hexagonal boron nitride (h-BN) is one of the main components of nanodevices along with graphene, less is known about its physical properties. Nevertheless, a large number of phenomena1–27 that have been discovered for bilayer and few-layer graphene regarding the interaction of 2D layers can also be expected to occur in hexagonal boron nitride and can provide useful information for the development of new applications. These phenomena include, among others, the observation of Moiré patterns upon the relative rotation of graphene layers,1–3 atomic-scale slip-stick motion of a graphene flake attached to an STM tip on a graphene surface,4–6 self-retracting motion of the layers at their telescopic extension,8 diffusion and drift of a graphene flake on a graphene surface via rotation to incommensurate states9,10 and observation of dislocations in the stacking of graphene layers.7,18–26 The latter takes place when different parts of the system experience different relative shifts between the layers at the atomic scale and the system gets divided into a number of commensurate domains separated by incommensurate boundaries. Numerous experimental and theoretical studies have revealed that such dislocations can be used to tune electronic22–25,27 and optical26 properties of graphene and thus hold great promise for application in nanoelectronic devices. Variation of the electronic properties of graphene layers upon their relative in-plane displacement11–13 has already been proposed as a basis for various nanosensors.13–17 Understanding the interlayer interaction in h-BN could result in equally interesting observations and developments. It is not surprising that considerable efforts are being made to investigate the structural and energetic characteristics of h-BN.28–56

The structural and elastic properties of the ground state of bulk h-BN have been studied experimentally in detail. The interlayer distance,28–36 bulk modulus,31,32,35,36 shear modulus,35 shear mode frequency37,38 (E2g mode with adjacent layers sliding rigidly in the opposite in-plane directions) and frequency of relative out-of-plane vibrations (B1g ZO mode with adjacent layers sliding rigidly towards and away from each other)39 have been measured and stacking of the layers has been established.28,40 The interlayer distance has also been determined for few-layer h-BN,41–43 while for bilayer h-BN, only stacking is known from the experimental studies.42 Rotational stacking faults have been observed for few-layer h-BN.42 However, no experimental data on the interlayer binding energy, corrugation of the potential surface of interlayer interaction energy and barriers to in-plane relative displacement of the layers are available neither for bilayer nor for few-layer or bulk h-BN.

State-of-the-art first-principles methods have been applied to calculate the interlayer distance,39,43–55 interlayer binding energy,43,48,54 bulk modulus,44–48,50,51,53 shear modulus,53 shear mode frequency46,52,53 and frequency of out-of-plane relative vibrations39,45,46,52,53 for bulk h-BN. Much less theoretical data are available for bilayer h-BN, including the interlayer distance,43,48,54–56 interlayer binding energy,43,48,54,56 bulk modulus and frequency of out-of-plane vibrations.43 The dependence of the interlayer binding energy on the layer stacking has been analyzed both for bulk49,51,55 and bilayer h-BN54,55,57 and the barriers to relative sliding of the layers have been calculated.54,55,57 However, the results of these calculations on relative in-plane motion of h-BN layers are rather contradictory and obtained with the aim of studying the balance between van der Waals and electrostatic forces without relation to measurable physical quantities.

Considerably more detailed theoretical and experimental data are available for the interaction and relative motion of graphene layers. The approximation58 of the interaction energy of a single carbon atom in a graphene flake with a graphite layer by the first Fourier components determined by graphene symmetry has been used to derive a simple expression for the potential surface of the interaction energy of graphene layers.9,10,59 The latter has been confirmed both by density functional theory calculations with the dispersion correction (DFT-D) and calculations using empirical potentials with very good accuracy.10,60,61 According to this expression, the shape of the potential energy surface at a given interlayer distance is determined by a single parameter and thus links different physical quantities related to relative motion of the layers such as barriers to relative sliding and rotation, shear mode frequencies and shear modulus. Moreover, the approximation enabled the Frenkel–Kontorova model to be extended to the case of two interacting structures with a slightly different lattice constant62 to calculate the width and formation energy of dislocations in the stacking of graphene layers and to consider the commensurate–incommensurate phase transition in bilayer graphene upon stretching of one of the layers.63 Using the approximation, the barrier to relative motion of the layers in bilayer graphene was estimated, on the basis of the experimentally measured values of the shear mode frequency and width of dislocations in the layer stacking, to be 1.7 meV per atom (ref. 61) and 2.4 meV per atom (ref. 7), respectively (note that all energies for bilayers in the present paper are represented in meV per atom in the top (adsorbed) layer so that they are more or less equal to the exfoliation energy of the bulk material per atom and adsorption energy of a flake on the 2D layer per atom of the flake). These values are in good agreement with those predicted by theory: ∼1 meV per atom (ref. 64, DFT calculations for graphite within local density approximation (LDA) and generalized gradient approximation (GGA)), 2.6 meV per atom (ref. 59), 1.82 meV per atom (ref. 65, LDA calculations for bilayer graphene), 1–1.5 meV per atom (ref. 66, DFT-D calculations for polycyclic aromatic hydrocarbons adsorbed on graphene), 2.6 meV per atom (ref. 59), 2.07 meV per atom (ref. 60, DFT-D calculations for bilayer graphene), 1 meV per atom (ref. 59), 1.92 meV per atom (ref. 60, calculations for bilayer graphene using the vdW-DF functional67). It should also be noted that although there is a considerable discrepancy in the experimental estimates of the barrier to relative motion of graphene layers, its relative value is even less than that for the exfoliation energy of graphite ranging from 31 ± 2 to 52 ± 5 meV per atom (ref. 68–71).

The aim of the present paper is to predict experimentally measurable physical quantities associated with the potential surface of the interaction energy of h-BN layers. We use the vdW-DF2 approach72 and the adequacy of our calculations is confirmed by the agreement in the interlayer interaction energy for symmetric stackings of h-BN layers with the values obtained recently by the local second-order Møller–Plesset perturbation theory (LMP2),55 which is a high-level ab initio method that, different from DFT, fully describes vdW interactions. Analogously to the expression for the potential surface of interaction energy between graphene layers, simple expressions describing potential energy surfaces of h-BN layers aligned in the same and opposite directions are introduced with the parameters fitted to the results of the calculations. Based on these expressions a set of physical quantities related to the interlayer interaction and relative motion of h-BN layers is estimated, including the barriers to relative motion and rotation of the layers, shear mode frequencies and shear modulus for bilayer and bulk h-BN, length and formation energy of a dislocation in stacking of the layers in bilayer h-BN. The implications of the approximation introduced for the h-BN potential energy surface for further calculations of h-BN mechanical properties and simulations of dynamical phenomena related to oscillation of h-BN layers are discussed. The results obtained are also important for understanding the interactions in recently discovered73–76 graphene-h-BN heterostructures (see also ref. 77 for a review) showing the commensurate–incommensurate phase transition.76 The calculations performed here represent a necessary intermediate stage (along with the previous calculations for bilayer graphene60,61) before theoretical consideration of the properties of such heterostructures related to the interlayer interaction.

The paper is organized in the following way. Section 2 is devoted to first-principles calculations of elastic properties of a single h-BN layer and the potential surface of interlayer interaction energy in bilayer and bulk h-BN. The estimates of physical quantities associated with the interaction and relative motion of h-BN layers are presented in Section 3. Our conclusions are summarized in Section 4.

2 Methods and results

2.1 DFT calculations of potential energy surfaces

The DFT calculations have been performed using VASP code78 with the non-local vdW-DF2 functional72 to take into account van der Waals (vdW) interactions. The basis set consists of plane waves with a maximum kinetic energy of 600 eV. The interaction of valence electrons is described using the projector augmented-wave method (PAW).79 The periodic boundary conditions are applied to the rectangular unit cell with 2 boron and 2 nitrogen atoms in each h-BN layer. The height of the simulation cell is 20 Å for single-layer and bilayer h-BN and 6.66 Å for the bulk material. Integration over the Brillouin zone is performed using the Monkhorst–Pack method80 with a k-point grid of 24 × 20 × 1 for single-layer and bilayer h-BN and 24 × 20 × 4 for the bulk material (in Section 3, axes x and y are chosen in the armchair and zigzag directions, respectively). Convergence with respect to the number of k-points in the Brillouin zone and the maximum kinetic energy of the plane waves was tested previously for the energetic characteristics of bilayer graphene related to the interlayer interaction.60 The parameters chosen in the present study allow these properties to be evaluated within an accuracy of 2%. The convergence threshold of the self-consistent field is 10−6 eV.

The optimized bond length is l = 1.455 Å, which is in agreement with the experimental data for few-layer40–43 and bulk h-BN28–36 and previous calculations.39,43–53,55 The effect of the interlayer interaction on the structure of h-BN layers is neglected.55

The potential surface of the interlayer interaction energy U(x,y), i.e., the dependence of the interlayer interaction energy U on the relative displacement of two h-BN layers in the armchair (x) and zigzag directions (y) at a fixed interlayer distance d, has been calculated by shifting the layers with respect to each other as rigid. Account of deformations within the layers was shown to have a negligible effect on the potential surface of the interaction energy of the carbon nanotube walls81 and the shells of the carbon nanoparticles.82,83 Taking into account the mirror planes and three-fold rotational symmetry of a single h-BN layer, it is sufficient to explicitly calculate U(x,y) only for 1/6 part of the considered unit cell. The dependence U(x,y) has been obtained for 78 points uniformly distributed on this part of the unit cell. The interlayer distance has been set at the experimental value of d = 3.33 Å for the bulk material.28–36 The measurements for few-layer h-BN are less precise and include this value within the error bars 3.25 ± 0.10 Å (ref. 42). It should also be noted that the LMP2 calculations55 give the same interlayer distance of 3.34 Å for bilayer and bulk h-BN. Though significant efforts have been made recently to improve the description of the long-range vdW interactions in DFT,67,72,84,85 available DFT methods still fail to predict the equilibrium distance between 2D layers with sufficient accuracy.55,59,60,86 However, the magnitude of corrugations of the potential energy surface59,60,64 and the barriers to relative motion of the layers60 depend exponentially on the interlayer distance. Therefore, small deviations in the equilibrium interlayer distance predicted by different DFT methods result in significant errors in the properties related to interaction and relative motion of the layers.55,59,60,86 On the other hand, inclusion of vdW interactions almost does not affect the potential energy surface at a given interlayer distance.54,55,59,60,66 Thus, it is reasonable to study the potential energy surface at the equilibrium interlayer distance known from the experiments rather than to rely on the values provided by the DFT methods.

Let us now discuss the results of calculations of the potential energy surfaces for two orientations of h-BN layers in the simulation cell (Fig. 1). For h-BN layers aligned in the same direction (Fig. 1a and b), the potential surface of the interlayer interaction energy is very similar to that of graphene layers, which is in agreement with previous calculations.55 The energy minima correspond to the AB stacking in which nitrogen (boron) atoms of the top layer are located on top of boron (nitrogen) atoms of the bottom layer and boron (nitrogen) atoms of the top layer are on top of hexagon centers. It is clear that at this relative orientation of the layers all the configurations corresponding to the AB stacking with boron (AB2) or nitrogen (AB1) atoms of the top layer on top of the centers of the hexagons are equivalent. Shifting the layers by half of the bond length in the armchair direction (towards the nearest minimum) brings the systems to the saddle-point (SP) stacking which corresponds to the barrier to relative sliding of the layers. We find that this barrier is about 2 meV per atom (Tables 1 and 2), which is within the range reported for graphene7,59–61,64–66 and close to the LMP2 value ∼2.5 meV per atom for h-BN.55 Shifting the layers by one and a half of the bond length more in the armchair direction results in the AA stacking in which all the atoms of the top layer are located on top of the equivalent atoms of the bottom layer and the interlayer interaction energy reaches its maximum for the given interlayer distance. The energy of the AA stacking relative to the AB one is about 18 meV per atom (Tables 1 and 2), which is close to the DFT data for graphene7,59–61,64–66 and for bilayer55 and bulk49,51,55 h-BN as well as the LMP2 values of about 20 meV per atom55 for bilayer and bulk h-BN.


image file: c5ra20882c-f1.tif
Fig. 1 Calculated interlayer interaction energy of h-BN bilayer U (in meV per atom of the top layer) as a function of relative displacement of the layers in the armchair (x, in Å) and zigzag (y, in Å) directions at the interlayer distance of d = 3.33 Å: (a and b) h-BN layers aligned in the same direction and (c and d) h-BN layers aligned in the opposite directions. The energy is given relative to the AA′ stacking. (b and d) Black solid lines correspond to the calculated dependences of interlayer interaction energy U on displacement x in the armchair direction (y = 0) along the thick lines indicated in figures (a and c). Curves approximated according to eqn (2) (b) and (3) (d) are shown by red dashed lines. Structures of the symmetric stackings are indicated. Boron and nitrogen atoms are coloured in blue and magenta, respectively.
Table 1 Calculated relative energies of symmetric stackings of h-BN bilayer with respect to the AA′ stacking (in meV per atom of the top layer) for the interlayer distance of d = 3.33 Å. The LMP2 data are taken from ref. 55
Stacking vdW-DF2 LMP2
AB1′ 3.10 4.42
AB2′ 15.77 16.50
SP′ 3.57  
AB 0.40 0.24
SP 2.32  
AA 18.26 19.72


Table 2 Calculated relative energies of symmetric stackings of h-BN bulk with respect to the AA′ stacking (in meV per atom) for the interlayer distance of d = 3.33 Å. The LMP2 data are taken from ref. 55
Stacking vdW-DF2 LMP2
AB1′ 3.02 3.76
AB2′ 15.91 16.14
SP′ 3.57  
AB 0.38 0.44
SP 2.39  
AA 18.47 19.89


The potential energy surface for h-BN layers aligned in the opposite directions (Fig. 1c and d) is rather different from those for graphene and h-BN layers aligned in the same direction, which is in agreement with previous calculations.54,55 In this case there are two types of inequivalent AB′ stackings. In the AB1(2)′ stacking all of the boron (nitrogen) atoms of the top layer are on top of the boron (nitrogen) atoms of the bottom layer and the rest of the atoms of the top layer are on top of the hexagon centers. Also, the interlayer interaction energy reaches its maximum for the given orientation at the AB2′ stacking, while the local minima correspond to the AB1′ stacking along with the AA′ stacking in which all of the nitrogen and boron atoms of the top layer are on top of the boron and nitrogen atoms of the bottom layer, respectively. The saddle-point stacking corresponding to the barrier to relative sliding of the layers lies on the straight path between the nearest configurations corresponding to the AB1′ and AA′ but not in the middle (0.53 Å from AB1′ and 0.92 Å from AA′) since the minima are of different depth.

All the DFT methods predict that the energies of the AB, AB1′ and AA′ stackings are rather close and differ only by several meV per atom. However, the order of stability of these structures changes depending on the method and chosen interlayer distance.43,49,51,54,55,57 Our results indicate that the most stable stacking for both h-BN bulk and bilayer is AA′ (Tables 1 and 2) and these results are consistent with the experimental observations for bulk h-BN28 as well as with the LMP2 calculations55 and DFT calculations with account of nonlocal many-body dispersion (MBD)57 for bilayer and bulk h-BN. The AB stacking is higher in energy only by 0.4 meV per atom and the relative energy of the AB1′ stacking is 3 meV per atom, which is in agreement with the LMP2 results55 (Tables 1 and 2) and the experimental data42 that both the AA′ and AB stackings are observed for bilayer h-BN. The AB1′ minima in our case are very shallow and are close in energy to the transition state SP′ for relative sliding of the layers. The corresponding barrier for transition from the AA′ to AB1′ stacking is 3.6 meV per atom, somewhat higher than in the case of h-BN layers aligned in the same direction.

It should also be mentioned that switching the vdW contribution off (using only LDA correlation and rPW86 exchange72) does not change the order of the symmetric stackings (Table 3). Their relative energies with respect to the AA′ stacking are changed by 0.4–0.8 meV per atom, which is within 5% of the magnitudes of corrugation of the potential energy surfaces both for h-BN layers aligned in the same and opposite directions. An analogous conclusion has previously been made for graphene bilayer59,60,66 and h-BN bilayer54,55,57 for a set of functionals with the dispersion correction.

Table 3 Calculated relative energies of symmetric stackings of h-BN bilayer with respect to the AA′ stacking (in meV per atom of the top layer) for the interlayer distance of d = 3.33 Å with (vdW-DF2) and without (rPW86(x) + LDA(c)) account of van der Waals interactions
Stacking vdW-DF2 rPW86(x) + LDA(c) Difference
AB1′ 3.10 3.87 −0.76
AB2′ 15.77 15.30 0.48
SP′ 3.57 4.08 0.51
AB 0.40 0.85 −0.45
SP 2.32 2.68 0.36
AA 18.26 17.62 0.64


Another observation is that the calculated relative energies of different stackings for bilayer and bulk h-BN are nearly identical (Tables 1 and 2), which is in agreement with previous calculations for graphene.60,61 The relative deviation is within 0.2 meV per atom, i.e. about 1% of the magnitudes of corrugation of the potential energy surfaces for h-BN layers aligned in the same and opposite directions, and can be explained by interaction of non-adjacent layers in the bulk material.

2.2 Approximation of potential energy surfaces

To approximate the potential energy surface of bilayer h-BN we consider the top layer as being adsorbed on the bottom one and recall that the potential energy surface of an atom adsorbed on a 2D trigonal lattice can be approximated by the first Fourier harmonics as58:
 
U = 0.5(U0 + U1(2[thin space (1/6-em)]cos(k1x)cos(k2y) + cos(2k1x) + 1.5)), (1)
where k1 = 2π/(3l), image file: c5ra20882c-t1.tif l is the bond length and the point x = 0 and y = 0 corresponds to the case when the atom is located on top of one of the lattice atoms. In the case of h-BN, we should sum up interactions for boron atoms on the boron lattice (BB), nitrogen atoms on the nitrogen lattice (NN), boron atoms on the nitrogen lattice (BN) and nitrogen atoms on the boron lattice (NB). In the latter two contributions U1NB = U1BN and U0NB = U0BN. Then for h-BN layers aligned in the same direction we arrive at:
 
UI = U0,tot + 4.5U1NB + (U1NN + U1BBU1NB)×(2[thin space (1/6-em)]cos[thin space (1/6-em)](k1x)[thin space (1/6-em)]cos[thin space (1/6-em)](k2y) + cos[thin space (1/6-em)](2k1x) + 1.5), (2)
where the point x = 0 and y = 0 corresponds to the AA stacking and U0,tot = 2U0NB + U0NN + U0BB. It is seen that corrugations of the potential energy surface of h-BN layers aligned in the same direction are described by a single parameter (U1NN + U1BBU1NB), which determines all the physical properties related to the interlayer interaction for h-BN materials with adjacent layers aligned in the same direction.

For h-BN layers aligned in the opposite directions,

 
image file: c5ra20882c-t2.tif(3)
where the point x = 0 and y = 0 corresponds to the AA′ stacking. In this case corrugations of the potential energy surface are described by two parameters (2U1NB − 0.5U1NN − 0.5U1BB) and (U1NNU1BB) and these parameters are responsible for the properties of h-BN materials with adjacent layers aligned in the opposite directions.

The parameters of the approximation are found from the relative energies of the symmetric stackings AA′, AB1′, AB2′ and AB as

 
image file: c5ra20882c-t3.tif(4)

The values obtained are U1NB = −0.0888 meV per atom, U1BB = 3.328 meV per atom and U1NN = 0.512 meV per atom and allow straightforward physical interpretation. The negative value of U1NB corresponds to the attraction of ions of different sign, while positive values of U1BB and U1NN are related to repulsion of ions of the same sign. The greater value of U1BB compared to U1NN is associated with the larger size of the boron ion as compared to the nitrogen one.87

In addition to corrugations of the potential energy surface for h-BN layers aligned in the opposite directions and relative energy of the co-aligned configuration, which are reproduced exactly, the approximation is very accurate in relative energies of the SP and AA stackings for h-BN layers aligned in the same direction (Fig. 1b and d). The magnitude of corrugation for this orientation of the layers according to the approximation is EAAEAB = 4.5(U1NN + U1BBU1NB) = 17.68 meV per atom, which is only 0.18 meV per atom or 1% smaller than the DFT value obtained. The barrier to relative sliding of the layers aligned in the same directions according to the approximation is ESPEAB = 0.5(U1NN + U1BBU1NB) = 1.96 meV per atom. This is only 0.05 meV per atom or 2.6% greater then the DFT value. The standard deviation of the approximated potential energy surfaces from the ones obtained by the DFT calculations for h-BN layers aligned in the same and opposite directions are 0.056 meV per atom and 0.014 meV per atom, respectively, which is within 0.3% of their magnitudes of corrugation.

We should note that the possibility to accurately approximate the potential energy surfaces by expressions containing only the first Fourier harmonics has been previously demonstrated for interaction between graphene layers60,61,66,88 and between carbon nanotube walls, both for infinite commensurate walls81,89–92 and in the case where corrugations of the potential surface are determined by the contribution of edges93 or defects.81 Thus, we can expect that analogous expressions can describe potential energy surfaces for other layered materials with the van der Waals interaction between layers or for translational motion of large molecules physically adsorbed on crystal surfaces.

It should also be pointed out that the expression given by eqn (2) for approximation of the potential energy surface for h-BN layers aligned in the same direction is exactly the same as the one for graphene.9,10,59,61 The calculated barrier to relative sliding of h-BN layers aligned in the same direction and the magnitude of corrugation of the potential energy surface (Tables 1 and 2) are within the ranges reported for graphene.7,59–61,64–66 Therefore, it can be expected that these materials have very similar physical properties. This is confirmed in Section 3 by comparison of physical quantities estimated for h-BN on the basis of eqn (2) with the experimental data for graphene.

2.3 Mechanical properties of single-layer h-BN

To get the Young modulus and Poisson ratio of a single h-BN layer, strains of up to 0.8% in the armchair direction and up to 0.5% in the perpendicular zigzag direction have been considered. For each size of the unit cell, the positions of the atoms within the cell have been optimized using the quasi-Newton method94 up to a residual force of 6 × 10−3 eV Å−1. The energy dependence on the strain in the zigzag direction at a given strain in the armchair direction has been approximated by a parabola. The dependence of the energy corresponding to the parabola minimum on the elongation in the armchair direction gives the elastic constant k = 272.8 ± 0.5 J m−2 or Young modulus Y = k/d = 819.2 ± 1.5 GPa, where d = 3.33 Å is the interlayer distance in bulk h-BN, while the dependence of the position of the parabola minimum gives the Poisson ratio ν = 0.2011 ± 0.0010. These values are in excellent agreement with the experimental data for bulk h-BN of Y = 811 GPa and ν = 0.21 (ref. 35) and are close to the DFT results for bulk h-BN of Y = 860–900 GPa, ν = 0.21–0.22 (ref. 53), Y = 952 GPa, ν = 0.18 (ref. 46) and for single-layer h-BN of k = 279 J m−2, ν = 0.22 (ref. 95) and 291 J m−2, ν = 0.21 (ref. 96).

3 Properties related to interlayer interaction

As shown in the previous section, the potential surfaces of interlayer interaction energy for h-BN layers aligned in the same and opposite directions can be approximated by relatively simple functions involving only one or two independent parameters (eqn (2) and (3)), respectively. Therefore, all the properties of h-BN related to interaction and relative motion of the layers are basically determined by these one or two parameters for the layers aligned in the same and opposite directions, respectively. Our DFT calculations demonstrate that interaction of non-adjacent layers in bulk h-BN has a weak influence on the potential energy surface so that the same parameters are responsible for the physical behavior of bilayer, few-layer and bulk h-BN.

3.1 Shear mode frequency and shear modulus

The frequency of the shear mode E2g in n-layer and bulk (n → ∞) h-BN, in which adjacent layers slide rigidly in the opposite in-plane directions, is determined by the curvature of the potential energy surface in a given metastable state61:
 
image file: c5ra20882c-t4.tif(5)
where μ = p(mN + mB)/4 for n = 2p and μ = p(p + 1)(mN + mB)/2(2p + 1) for n = 2p + 1, mB and mN are masses of boron and nitrogen atoms, respectively, p is an integer, l is the bond length and we denote Ueff = (l/2π)22U/∂x2. Curvatures of the potential energy surfaces in the states corresponding to the AA′, AB1′ and AB stackings are described by the parameters Ueff(AA′) = (U1NN + U1BB − 4U1NB)/3 = 1.40 meV per atom, Ueff(AB1′) = (U1BB + 2U1NB − 2U1NN)/3 = 0.71 meV per atom and Ueff(AB) = (U1NN + U1BBU1NB)/3 = 1.31 meV per atom. The corresponding shear mode frequencies as functions of the number of layers n are shown in Fig. 2. The estimated shear mode frequency for the ground-state AA′ stacking of bulk h-BN is in good agreement with the experimental data37,38 and the results of previous DFT calculations46,52,53 (Fig. 2). The calculated dependences of the shear mode frequencies for the AA′ and AB stackings on the number of layers are similar to the experimental ones for graphene.97,98

image file: c5ra20882c-f2.tif
Fig. 2 Calculated shear mode frequencies f (in cm−1) as functions of the number of layers n for h-BN crystals with different stacking of the layers: AA′ — blue solid line, AB1′ — green dashed line, AB — red dash-dotted line. The experimental and theoretical data for bulk h-BN are included: □ — ref. 37 (experiment), • — ref. 38 (experiment), × — ref. 46 (LDA), ◊ — ref. 52 (LDA), ▲ — ref. 53 (LDA), ▼ — ref. 53 (revPBE). The experimental data for few-layer graphene and graphite are also shown for comparison: ■ — ref. 97 and ♦ — ref. 98.

The curvature of the potential energy surface also determines the shear modulus, which is expected to be weakly dependent on the number of h-BN layers and can be estimated as:

 
image file: c5ra20882c-t5.tif(6)
where image file: c5ra20882c-t6.tif is the area per atom in an h-BN layer and d = 3.33 Å is the interlayer distance. The shear moduli obtained for different stackings of the layers are given in Fig. 3. It is seen that the result for the AA′ stacking is within the range of the experimental35,99 and theoretical values53 for bulk h-BN.


image file: c5ra20882c-f3.tif
Fig. 3 Calculated shear moduli C44 (in GPa) for h-BN crystals (few-layer or bulk) with different stacking of the layers. The experimental data for bulk h-BN are taken from ref. 35 (Exp1) and ref. 99 (Exp2). The range of values from the previous DFT calculations53 for bulk h-BN is indicated.

3.2 Barrier to relative rotation of h-BN layers

As is the case for graphene,9,10 the potential energy surface of h-BN layers rotated to an incommensurate state should be smooth. Therefore, the interaction energy in such states can be estimated as an average of the interaction energy (eqn (2) and (3)) over in-plane relative displacements of the layers in the commensurate states:
 
Urot = 〈UIx,y = 〈UIIx,y = U0,tot + 1.5(U1NN + U1BB + 2U1NB) (7)

Thus, the barriers to relative rotation of h-BN layers from metastable states corresponding to the AA′, AB1′ and AB stackings are Erot(AA′) = UrotUII(AA′) = 1.5(U1NN + U1BB − 4U1NB) = 6.29 meV per atom, Erot(AB1′) = −1.5U1NN + U1BB + 3U1NB = 2.29 meV per atom and Erot(AB) = 1.5(U1NN + U1BB) − U1NB = 5.85 meV per atom, respectively. The values for the most stable AA′ and AB stackings are of the order of the data for graphene Erot ∼ 5 meV per atom (ref. 61) and 4 meV per atom (ref. 9 and 10). The obtained expressions for the potential energy surface (eqn (2) and (3)) can be also used to estimate the barriers to relative rotation of h-BN layers in configurations corresponding to Moiré patterns analogously to studies performed for graphene.100

3.3 Dislocations

Another manifestation of incommensurate states in two-dimensional layered materials is related to dislocations in stacking of the layers, which are observed when different parts of the system experience different relative shifts between the layers and correspond to incommensurate boundaries between commensurate domains.7,18–26 It has been demonstrated that such dislocations have a strong effect on the electronic22–25,27 and optical26 properties of graphene. The analysis of the structure of dislocations in few-layer graphene by transmission electron microscopy has already helped to get an experimental estimate of the barrier to relative displacement of the layers.7 Along with the properties discussed above the formation energy and width of dislocations in bilayer h-BN can be obtained from the potential surface of interlayer interaction energy.63

To consider dislocations in bilayer h-BN we follow the formalism of the two-chain Frenkel–Kontorova model62 which was applied previously to study dislocations in double-walled carbon nanotubes62,91 and graphene.63 Elementary dislocations in graphene correspond to partial dislocations with the Burgers vector [b with combining right harpoon above (vector)] equal in magnitude to the bond length.7,63 Similar dislocations (b = l, Fig. 4c and d) are also possible in h-BN layers aligned in the same direction as the potential surface of interlayer interaction in this case has degenerate energy minima AB1 and AB2 separated by a distance of one bond length l (Fig. 1a). However, for h-BN layers aligned in the opposite directions, the AB1′ energy minima are rather shallow and are much higher in energy than the AA′ minima (Fig. 1c). Therefore, in this case the minimal possible Burgers vector is equal in magnitude to the lattice constant (image file: c5ra20882c-t7.tif Fig. 4a and b), i.e. the elementary dislocations are full. Both full dislocations in h-BN layers aligned in the opposite directions and partial dislocations in h-BN layers aligned in the same direction are considered below. It should also be noted that even for graphene theoretical estimates of the formation energy and width have been limited to the tensile dislocation with the Burger vector [b with combining right harpoon above (vector)] normal to the boundary between commensurate domains,63 although experimentally dislocations with different angles between the Burgers vector and boundary between commensurate domains have been observed.7,19–22,24 Here we extend the Frenkel–Kontorova model to describe dislocations with the Burgers vector [b with combining right harpoon above (vector)] at an arbitrary angle β to the normal [n with combining right harpoon above (vector)] to the boundary between commensurate domains. In particular, the results for shear dislocations (Fig. 4b and d), where the Burgers vector is parallel to the boundary between commensurate domains β = ±π/2, are presented along with the results for tensile dislocations (Fig. 4a and c) with β = 0.


image file: c5ra20882c-f4.tif
Fig. 4 (a–d) Atomistic structures corresponding to full tensile (a) and shear (b) dislocations in h-BN bilayer with the layers aligned in the opposite directions and partial tensile (c) and shear dislocations (d) in h-BN bilayer with the layers aligned in the same direction. Boron and nitrogen atoms are coloured in blue and magenta, respectively. The directions of the unit vector [n with combining right harpoon above (vector)] normal to the boundary between commensurate domains (black dashed line) and Burgers vector [b with combining right harpoon above (vector)] are indicated. The magnitudes of the vectors are scaled up for clarity. (a and b) The directions of the Burgers vectors [b with combining right harpoon above (vector)]1 and [b with combining right harpoon above (vector)]2 corresponding to two straight pieces of the path of full dislocations ([b with combining right harpoon above (vector)] = [b with combining right harpoon above (vector)]1 + [b with combining right harpoon above (vector)]2) are shown. (e–g) Displacements u/l (l = 1.455 Å is the bond length) of atoms as functions of their position x (in nm) in the direction of the normal [n with combining right harpoon above (vector)] to the boundary between commensurate domains for full tensile (e) and shear (f) dislocations in h-BN bilayer with the layers aligned in the opposite directions and for partial (g) tensile and shear dislocations in h-BN bilayer with the layers aligned in the same direction. (e and f) The displacements along the Burgers vector (black solid line) and in the perpendicular direction (red dashed line) are shown. (g) The displacement along the Burgers vector is shown. The displacement in the perpendicular direction is zero. The results for tensile (black solid line) and shear dislocations (blue dashed line) are given. The atomistic structures of symmetric stackings across the dislocations are included with their position shown by vertical gray lines. The characteristic widths lD of the dislocations are indicated.

In the limit of large systems and low density of dislocations it can be assumed that dislocations are isolated. This means that the layers are commensurate at large distances from the dislocation center. Correspondingly, the boundary conditions for the relative displacement [u with combining right harpoon above (vector)] of atoms of the layers can be formulated as [u with combining right harpoon above (vector)] = 0 at x → −∞ and [u with combining right harpoon above (vector)] = −[b with combining right harpoon above (vector)] at x → +∞, where x is the coordinate along the normal [n with combining right harpoon above (vector)] to the boundary between commensurate domains and [b with combining right harpoon above (vector)] is the Burgers vector. According to the Frenkel–Kontorova model,62,63 the formation energy of isolated dislocations per unit length is determined by the sum of excessive elastic and van der Waals interaction energies originating from the deformation and incommensurability of the layers, respectively,

 
image file: c5ra20882c-t8.tif(8)
where the axis y is chosen along the boundary between commensurate domains, [u with combining right harpoon above (vector)]′ = d[u with combining right harpoon above (vector)]/dx describes the strain across the boundary, V([u with combining right harpoon above (vector)]) is the van der Waals interaction energy per unit area relative to the deepest energy minimum for the given alignment of h-BN layers (the AB and AA′ stackings for the layers aligned in the same and opposite directions, respectively), E = k/(1 − ν2) and G = k/2(1 + ν) are the tensile and shear elastic constants of the layers, respectively, k = Yd is the elastic constant under uniaxial stress and ν is the Poisson ratio. Introducing the coefficient (1 − ν2)−1 in the tensile elastic constant E we take into account that the tensile strain is kept zero in the direction along the boundary between commensurate domains (see also ref. 7). In the absence of an external load, one of the layers in the tensile dislocation is stretched in the direction across the boundary, while the other one is compressed.63 Therefore, in the perpendicular direction they tend to compress and stretch, respectively. The interlayer interaction, however, keeps the layers commensurate in the latter direction, providing the vanishing tensile strain along the boundary in the regions far from boundary crossings. It should also be noted that the coefficient (1 − ν2)−1 = 1.04 is only slightly different from unity and gives corrections to the formation energy and characteristic width of dislocations within 2%.

The dislocation path, i.e. the dependence of the relative displacement [u with combining right harpoon above (vector)] of h-BN layers on the coordinate x in the direction perpendicular to the boundary between commensurate domains (Fig. 4e–g) that minimizes the formation energy UD (eqn (8)), should satisfy the Euler–Lagrange equations:

 
image file: c5ra20882c-t9.tif(9)

This means that the dislocation path is the same as the trajectory of a particle on the inverse potential energy surface −V([u with combining right harpoon above (vector)]) with the fixed initial and final positions corresponding to the deepest energy minima of the original potential energy surface. It should be noted that the “particle mass” is anisotropic and equal to G/2 and E/2 in the directions along and across the boundary between commensurate domains, respectively.

For partial dislocations in h-BN layers aligned in the same direction (Fig. 4c and d), the straight line between the adjacent AB minima (Fig. 1a) corresponding to the minimum energy path satisfies eqn (9) and thus gives exactly the dislocation path. For full dislocations in h-BN layers aligned in the opposite directions (Fig. 4a and b), the solution of eqn (9) cannot be found analytically. Moreover, this solution should depend on the angle between the Burgers vector and the boundary between commensurate domains. However, it is clear that strong scattering introduced by rapidly changing potential at slopes of the AB2′ hill (Fig. 1c) should be avoided. Therefore, we assume that the path of full dislocations in h-BN layers aligned in the opposite directions is also roughly described by the minimum energy path AA′–AB1′–AA′ consisting of two straight lines at the angle 2π/3 to each other. More accurate approximations of the dislocation path give a correction to the formation energy within 10%.

Integration of the Euler–Lagrange equations (eqn (9)) provides the analogue of the energy conservation law:

 
image file: c5ra20882c-t10.tif(10)
here we use an integration constant of zero for isolated dislocations63 and denote image file: c5ra20882c-t11.tif, where θ is the angle between the direction [u with combining right harpoon above (vector)]′ of the dislocation path and normal [n with combining right harpoon above (vector)] to the boundary between commensurate domains. Taking into account eqn (10), the formation energy of dislocations per unit length (eqn (8)) can then be expressed as:
 
image file: c5ra20882c-t12.tif(11)

Let us first consider partial dislocations in h-BN layers aligned in the same direction (Fig. 4c and d). In these dislocations, atoms are displaced along straight lines, θ = π + β and u changes from 0 to b = l. According to eqn (11), the formation energy of partial dislocations is:

 
image file: c5ra20882c-t13.tif(12)
where we approximate the interlayer interaction energy V(u) along the minimum energy path on the basis of eqn (2). Fig. 5 demonstrates that the formation energy per unit length is maximal for tensile dislocations (Fig. 4c) and minimal for shear dislocations (Fig. 4d). Such dislocations correspond to the boundaries between commensurate domains in the zigzag and armchair directions, respectively. The estimated formation energies for h-BN layers aligned in the same direction (Fig. 5) are of the order of the value predicted for partial tensile dislocations in bilayer graphene of 0.12 eV Å−1 (ref. 63).


image file: c5ra20882c-f5.tif
Fig. 5 Calculated formation energy of dislocations UD per unit width (in eV Å−1) as a function of angle β (in degrees) between the Burgers vector [b with combining right harpoon above (vector)] and normal [n with combining right harpoon above (vector)] to the boundary between commensurate domains for a full dislocation (solid black line) in h-BN layers aligned in the opposite directions and a partial dislocation (dashed red line) in h-BN layers aligned in the same direction.

Based on the approximation of the interlayer interaction energy along the minimum energy path by the cosine function and eqn (10), it was shown analytically that partial dislocations are described by a soliton with a relatively short incommensurate region separating commensurate domains.63 Although in the present paper we use eqn (2) to approximate the potential energy surface of h-BN layers aligned in the same direction, the calculated path of partial dislocations is very close to the analytical solution (Fig. 4g). Since the slope of the dependence of displacement u on position x across the boundary is nearly constant around the dislocation center and close to the maximum value |[u with combining right harpoon above (vector)]′|max, the dislocation width can be defined in a similar way to the analytical solution:

 
image file: c5ra20882c-t14.tif(13)
where Vmax is the barrier to relative sliding of the layers. It is clear that the dependence of the dislocation width on the angle β between the Burgers vector and normal to the boundary between commensurate domains is the same as that of the formation energy (eqn (12)). In particular, the width of shear dislocations is smaller than the width of tensile ones (Fig. 6). The dependence of the formation energy and width of partial dislocations on the angle between the Burgers vector and boundary between commensurate domains is determined only by the Poisson ratio and should approximately hold for any hexagonal 2D crystals or heterostructures consisting of layers with slightly different lattice constants where partial dislocations are possible. In particular, the calculated angular dependence of the dislocation width is in good agreement with experimentally measured values for graphene7,20,21 ranging from 11 nm for tensile dislocations to 6–7 nm for shear dislocations (Fig. 6). The calculated width of partial tensile dislocations is also of the order of the DFT estimate for graphene63 of 12 nm.


image file: c5ra20882c-f6.tif
Fig. 6 Calculated dislocation width lD (in nm) as a function of angle β (in degrees) between the Burgers vector [b with combining right harpoon above (vector)] and normal [n with combining right harpoon above (vector)] to the boundary between commensurate domains for a full dislocation (solid black line) in h-BN layers aligned in the opposite directions and a partial dislocation (dashed red line) in h-BN layers aligned in the same direction. The experimental data for graphene from ref. 7 are shown by squares (■) with error bars. The DFT result63 for partial tensile dislocations in graphene is indicated by a circle (•).

In the case of full dislocations in h-BN layers aligned in the opposite directions (Fig. 4a and b), two contributions from the straight pieces of the dislocation path at angles β + π/6 and β − π/6 to the normal [n with combining right harpoon above (vector)] to the boundary between commensurate domains need to be taken into account. This can be done by a simple substitution of ϕ(β) by ϕ(β − π/6) + ϕ(β + π/6) in eqn (12) and (13). All the conclusions made above for partial dislocations are also valid for full ones (Fig. 5 and 6) with the only difference being that full tensile and shear dislocations correspond to the boundaries between commensurate domains in the armchair and zigzag directions, respectively (Fig. 4a and b), which is opposite to the case of partial dislocations. It should be especially emphasized that according to Fig. 4e–g, both for partial and full dislocations the slope |[u with combining right harpoon above (vector)]′(x)| is nearly constant for distances from the dislocation center comparable to the dislocation width. This gives the possibility of accurate measurements of the dislocation width and, consequently, of experimental estimation of the barrier to relative sliding of h-BN layers by transmission electron microscopy analysis of stacking at the boundary between commensurate domains in the same way as has been done for graphene.7

4 Conclusions

DFT calculations of the potential surface of interlayer interaction energy at relative in-plane displacement of h-BN layers have been performed using the vdW-DF2 functional.72 According to these results, the ground state corresponds to the AA′ stacking of the layers aligned in the opposite directions. However, the optimal AB stacking for the layers aligned in the same direction is only 0.4 meV per atom less favourable, which is in agreement with the experimental observation42 of both AB and AA′ stackings for bilayer h-BN. The magnitudes of corrugation of the potential energy surfaces are about 18 meV per atom and 16 meV per atom for the layers aligned in the same and opposite directions, respectively. The barriers to relative motion of the layers aligned in the same and opposite directions are 2 meV per atom and 3.6 meV per atom, respectively, which is in close agreement with LMP2 calculations.55 Almost identical results are obtained for bilayer and bulk h-BN, indicating that interactions between non-adjacent layers have a weak effect on the characteristics of the potential energy surface (1% of the magnitude of corrugation). The contribution of vdW interactions into these quantities is also small (within 5% of the magnitude of corrugation) and comparable to that determined in previous studies for h-BN54,55 and graphene.59,60,66

It is shown that the calculated potential surfaces of interlayer interaction energy for h-BN layers aligned in the same and opposite directions can be fitted by simple expressions containing only the first components of Fourier expansions determined by symmetry of the layers. Analogous approximations have been suggested previously for graphene bilayer60,61,66,88 and double-walled carbon nanotubes81,89–93 based both on DFT calculations and semi-empirical potentials. Thus, it can be expected that for other layered materials the potential surface of interlayer interaction energy can also be reproduced closely by the first Fourier components determined by symmetry.

Recently, a concept of the registry index surface was proposed to predict the qualitative features of the potential surface of interlayer interaction energy of h-BN bilayer.54 The expressions introduced here give the possibility to reproduce the shape and quantitative characteristics of the potential energy surface of h-BN bilayer and can be useful for multiscale simulations of such phenomena as atomic-scale slip-stick motion of an h-BN flake attached to STM tip on the h-BN surface, diffusion of an h-BN flake on the h-BN surface and formation of stacking dislocations in h-BN bilayer (analogous to multiscale simulations of the phenomena observed for graphene4–10,18–26).

According to the proposed approximation, a set of physical properties of h-BN materials related to relative motion of the layers are determined by only one or two independent parameters for the layers aligned in the same and opposite directions, respectively. Namely, the shear mode frequency, shear modulus and barrier to relative rotation of the layers have been estimated for bilayer and bulk h-BN in different stackings. The possibility of partial and full dislocations in stacking of the layers is suggested for h-BN layers aligned in the same and opposite directions, respectively. Their width and formation energy are governed by the same parameters of the potential energy surface. In particular, it is shown that a simple link exists between the dislocation width lD and barrier Vmax to relative sliding of the layers Vmax = C(β)/lD2(β), where the coefficient C(β) = ϕ2(β)kl2/4 is determined by the elastic properties of the layers. Experimental measurements of the dislocation width, shear mode frequency or shear modulus can help to refine the parameters of the functions approximating the potential surfaces of interlayer interaction energy in h-BN and to obtain the barriers and other characteristics of these surfaces in the same way as has been done for graphene.7,61

Acknowledgements

The authors acknowledge the Russian Foundation of Basic Research (14-02-00739-a) and computational time on the Supercomputing Center of Lomonosov Moscow State University and the Multipurpose Computing Complex NRC “Kurchatov Institute”. IL acknowledges financial support from the Marie Curie International Incoming Fellowship (Grant Agreement PIIF-GA-2012-326435 RespSpatDisp), EU-H2020 project “MOSTOPHOS” (no. 646259) and Grupos Consolidados del Gobierno Vasco (IT-578-13).

References

  1. Z. Y. Rong and P. Kuiper, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 48, 17427–17431 CrossRef CAS.
  2. Y. Gan, W. Chu and L. Qiao, Surf. Sci., 2003, 539, 120–128 CrossRef CAS.
  3. J. H. Warner, M. H. Rümmeli, T. Gemming, B. Büchner and G. A. D. Briggs, Nano Lett., 2009, 9, 102–106 CrossRef CAS PubMed.
  4. M. Dienwiebel, G. S. Verhoeven, N. Pradeep, J. W. M. Frenken, J. A. Heimberg and H. W. Zandbergen, Phys. Rev. Lett., 2004, 92, 126101 CrossRef PubMed.
  5. M. Dienwiebel, N. Pradeep, G. S. Verhoeven, H. W. Zandbergen and J. W. M. Frenken, Surf. Sci., 2005, 576, 197–211 CrossRef.
  6. A. E. Filippov, M. Dienwiebel, J. W. M. Frenken, J. Klafter and M. Urbakh, Phys. Rev. Lett., 2008, 100, 046102 CrossRef PubMed.
  7. 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. U. S. A., 2013, 110, 11256–11260 CrossRef PubMed.
  8. Q. Zheng, B. Jiang, S. Liu, Y. Weng, L. Lu, Q. Xue, J. Zhu, Q. Jiang, S. Wang and L. Peng, Phys. Rev. Lett., 2008, 100, 067205 CrossRef PubMed.
  9. I. V. Lebedeva, A. A. Knizhnik, A. M. Popov, O. V. Ershova, Y. E. Lozovik and B. V. Potapkin, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 155460 CrossRef.
  10. I. V. Lebedeva, A. A. Knizhnik, A. M. Popov, O. V. Ershova, Y. E. Lozovik and B. V. Potapkin, J. Chem. Phys., 2011, 134, 104505 CrossRef PubMed.
  11. R. Bistritzer and A. H. MacDonald, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 245412 CrossRef.
  12. R. Bistritzer and A. H. MacDonald, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 12233–12237 CrossRef CAS PubMed.
  13. N. A. Poklonski, A. I. Siahlo, S. A. Vyrko, A. M. Popov, Y. E. Lozovik, I. V. Lebedeva and A. A. Knizhnik, J. Comput. Theor. Nanosci., 2013, 10, 141–146 CrossRef CAS.
  14. K.-T. Lam, C. Lee and G. Liang, Appl. Phys. Lett., 2009, 95, 143107 CrossRef.
  15. Y. Qian, K.-T. Lam, C. Lee and G. Liang, Carbon, 2012, 50, 1659–1666 CrossRef CAS.
  16. J. Zheng, P. Guo, Z. Ren, Z. Zhenyi, J. Bai and Z. Zhang, Appl. Phys. Lett., 2012, 101, 083101 CrossRef.
  17. K. K. Paulla and A. A. Farajian, J. Phys.: Condens. Matter, 2013, 25, 115303 CrossRef PubMed.
  18. L. Brown, R. Hovden, P. Huang, M. Wojcik, D. A. Muller and J. Park, Nano Lett., 2012, 12, 1609–1615 CrossRef CAS PubMed.
  19. B. Butz, C. Dolle, F. Niekiel, K. Weber, D. Waldmann, H. B. Weber, B. Meyer and E. Spiecker, Nature, 2013, 505, 533–537 CrossRef PubMed.
  20. J. Lin, W. Fang, W. Zhou, A. R. Lupini, J. C. Idrobo, J. Kong, S. J. Pennycook and S. T. Pantelides, Nano Lett., 2013, 13, 3262–3268 CrossRef CAS PubMed.
  21. M. Yankowitz, J. I.-J. Wang, A. G. Birdwell, Y.-A. Chen, K. Watanabe, T. Taniguchi, P. Jacquod, P. San-Jose, P. Jarillo-Herrero and B. J. LeRoy, Nat. Mater., 2014, 13, 786–789 CrossRef CAS PubMed.
  22. S. Hattendorf, A. Georgi, M. Liebmann and M. Morgenstern, Surf. Sci., 2013, 610, 53–58 CrossRef CAS.
  23. P. San-Jose, R. V. Gorbachev, A. K. Geim, K. S. Novoselov and F. Guinea, Nano Lett., 2014, 14, 2052–2057 CrossRef CAS PubMed.
  24. B. Lalmi, J. C. Girard, E. Pallecchi, M. Silly, C. David, S. Latil, F. Sirotti and A. Ouerghi, Sci. Rep., 2014, 4, 4066 CAS.
  25. M. M. Benameur, F. Gargiulo, S. Manzeli, G. Autès, M. Tosun, O. V. Yazyev and A. Kis, Nat. Commun., 2015, 6, 8582 CrossRef CAS PubMed.
  26. L. Gong, R. J. Young, I. A. Kinloch, S. J. Haigh, J. H. Warner, J. A. Hinks, Z. Xu, L. Li, F. Ding, I. Riaz, R. Jalil and K. S. Novoselov, ACS Nano, 2013, 7, 7287–7294 CrossRef PubMed.
  27. M. Koshino, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 115409 CrossRef.
  28. R. S. Pease, Nature, 1950, 165, 722–723 CrossRef PubMed.
  29. R. S. Pease, Acta Crystallogr., 1952, 5, 356–361 CrossRef CAS.
  30. R. W. Lynch and H. G. Drickamer, J. Chem. Phys., 1966, 44, 181 CrossRef CAS.
  31. V. L. Solozhenko, G. Will and F. Elf, Solid State Commun., 1995, 96, 1–3 CrossRef CAS.
  32. V. L. Solozhenko and T. Peun, J. Phys. Chem. Solids, 1997, 58, 1321–1323 CrossRef CAS.
  33. V. L. Solozhenko, A. G. Lazarenko, J.-P. Petitet and A. V. Kanaev, J. Phys. Chem. Solids, 2001, 62, 1331–1334 CrossRef CAS.
  34. W. Paszkowicz, J. B. Pelka, M. Knapp, T. Szyszko and S. Podsiadlo, Appl. Phys. A, 2002, 75, 431–435 CrossRef CAS.
  35. A. Bosak, J. Serrano, M. Krisch, K. Watanabe, T. Taniguchi and H. Kanda, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 041402 CrossRef.
  36. K. Fuchizaki, T. Nakamichi, H. Saitoh and Y. Katayama, Solid State Commun., 2008, 148, 390–394 CrossRef CAS.
  37. R. J. Nemanich, S. A. Solin and R. M. Martin, Phys. Rev. B: Condens. Matter Mater. Phys., 1981, 23, 6348–6356 CrossRef CAS.
  38. S. Reich, A. C. Ferrari, R. Arenal, A. Loiseau, I. Bello and J. Robertson, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 71, 205201 CrossRef.
  39. A. Marini, P. García-González and A. Rubio, Phys. Rev. Lett., 2006, 96, 136404 CrossRef PubMed.
  40. C. Li, Y. Bando, C. Zhi, Y. Huang and D. Golberg, Nanotechnology, 2009, 20, 385707 CrossRef PubMed.
  41. C. Zhi, Y. Bando, C. Tang, H. Kuwahara and D. Golberg, Adv. Mater., 2009, 21, 2889–2893 CrossRef CAS.
  42. J. H. Warner, M. H. Rümmeli, A. Bachmatiuk and B. Büchner, ACS Nano, 2010, 4, 1299–1304 CrossRef CAS PubMed.
  43. A. Nag, K. Raidongia, K. P. S. S. Hembram, R. Datta, U. V. Waghmare and C. N. R. Rao, ACS Nano, 2010, 4, 1539–1544 CrossRef CAS PubMed.
  44. K. Albe, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 55, 6203–6210 CrossRef CAS.
  45. G. Kern, G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 8551–8559 CrossRef CAS.
  46. N. Ohba, K. Miwa, N. Nagasako and A. Fukumoto, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 63, 115207 CrossRef.
  47. A. Janotti, S.-H. Wei and D. J. Singh, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 64, 174107 CrossRef.
  48. H. Rydberg, M. Dion, N. Jacobson, E. Schröder, P. Hyldgaard, S. I. Simak, D. C. Langreth and B. I. Lundqvist, Phys. Rev. Lett., 2003, 91, 126402 CrossRef CAS PubMed.
  49. L. Liu, Y. P. Feng and Z. X. Shen, Phys. Rev. B: Condens. Matter Mater. Phys., 2003, 68, 104102 CrossRef.
  50. T. Tohei, A. Kuwabara, F. Oba and I. Tanaka, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 064304 CrossRef.
  51. N. Ooi, V. Rajan, J. Gottlieb, Y. Catherine and J. B. Adams, Modell. Simul. Mater. Sci. Eng., 2006, 14, 515–535 CrossRef CAS.
  52. J. Serrano, A. Bosak, R. Arenal, M. Krisch, K. Watanabe, T. Taniguchi, H. Kanda, A. Rubio and L. Wirtz, Phys. Rev. Lett., 2007, 98, 095503 CrossRef CAS PubMed.
  53. I. Hamdi and N. Meskini, Phys. B, 2010, 405, 2785–2794 CrossRef CAS.
  54. N. Marom, J. Bernstein, J. Garel, A. Tkatchenko, E. Joselevich, L. Kronik and O. Hod, Phys. Rev. Lett., 2010, 105, 046801 CrossRef PubMed.
  55. G. Constantinescu, A. Kuc and T. Heine, Phys. Rev. Lett., 2013, 111, 036104 CrossRef PubMed.
  56. C.-R. Hsing, C. Cheng, J.-P. Chou, C.-M. Chang and C.-M. Wei, New J. Phys., 2014, 16, 113015 CrossRef.
  57. W. Gao and A. Tkatchenko, Phys. Rev. Lett., 2015, 114, 096101 CrossRef PubMed.
  58. G. S. Verhoeven, M. Dienwiebel and J. W. M. Frenken, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 70, 165418 CrossRef.
  59. M. Reguzzoni, A. Fasolino, E. Molinari and M. C. Righi, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 245434 CrossRef.
  60. I. V. Lebedeva, A. A. Knizhnik, A. M. Popov, Y. E. Lozovik and B. V. Potapkin, Phys. Chem. Chem. Phys., 2011, 13, 5687–5695 RSC.
  61. A. M. Popov, I. V. Lebedeva, A. A. Knizhnik, Y. E. Lozovik and B. V. Potapkin, Chem. Phys. Lett., 2012, 536, 82–86 CrossRef CAS.
  62. E. Bichoutskaia, M. I. Heggie, Y. E. Lozovik and A. M. Popov, Fullerenes, Nanotubes, Carbon Nanostruct., 2006, 14, 131–140 CrossRef CAS.
  63. A. M. Popov, I. V. Lebedeva, A. A. Knizhnik, Y. E. Lozovik and B. V. Potapkin, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 045404 CrossRef.
  64. A. N. Kolmogorov and V. H. Crespi, Phys. Rev.B: Condens. Matter Mater. Phys., 2005, 71, 235415 CrossRef.
  65. M. Aoki and H. Amawashi, Solid State Commun., 2007, 142, 123–127 CrossRef CAS.
  66. O. V. Ershova, T. C. Lillestolen and E. Bichoutskaia, Phys. Chem. Chem. Phys., 2010, 12, 6483–6491 RSC.
  67. M. Dion, H. Rydberg, E. Schröder, D. C. Langreth and B. I. Lundqvist, Phys. Rev. Lett., 2004, 92, 246401 CrossRef CAS PubMed.
  68. L. A. Girifalco and R. A. Lad, J. Chem. Phys., 1956, 25, 693–696 CrossRef CAS.
  69. L. X. Benedict, N. G. Chopra, M. L. Cohen, A. Zettl, S. G. Louie and V. H. Crespi, Chem. Phys. Lett., 1998, 286, 490–496 CrossRef CAS.
  70. R. Zacharia, H. Ulbricht and T. Hertel, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 155406 CrossRef.
  71. Z. Liu, J. Z. Liu, Y. Cheng, Z. Li, L. Wang and Q. Zheng, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 205418 CrossRef.
  72. K. Lee, E. D. Murray, L. Kong, B. I. Lundqvist and D. C. Langreth, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 081101 CrossRef.
  73. C. R. Dean, A. F. Young, I. Meric, C. Lee, L. Wang, S. Sorgenfrei, K. Watanabe, T. Taniguchi, P. Kim, K. L. Shepard and J. Hone, Nat. Nanotechnol., 2010, 5, 722–726 CrossRef CAS PubMed.
  74. L. A. Ponomarenko, A. K. Geim, A. A. Zhukov, R. Jalil, S. V. Morozov, K. S. Novoselov, I. V. Grigorieva, E. H. Hill, V. V. Cheianov, V. I. Fal’ko, K. Watanabe, T. Taniguchi and R. V. Gorbachev, Nat. Phys., 2011, 7, 958–961 CrossRef CAS.
  75. L. Britnell, R. V. Gorbachev, R. Jalil, B. D. Belle, F. Schedin, A. Mishchenko, T. Georgiou, M. I. Katsnelson, L. Eaves, S. V. Morozov, N. M. R. Peres, J. Leist, A. K. Geim, K. S. Novoselov and L. A. Ponomarenko, Science, 2012, 335, 947–950 CrossRef CAS PubMed.
  76. C. R. Woods, L. Britnell, A. Eckmann, R. S. Ma, J. C. Lu, H. M. Guo, X. Lin, G. L. Yu, Y. Cao, R. V. Gorbachev, A. V. Kretinin, J. Park, L. A. Ponomarenko, M. I. Katsnelson, Y. N. Gornostyrev, K. Watanabe, T. Taniguchi, C. Casiraghi, H.-J. Gao, A. K. Geim and K. S. Novoselov, Nat. Phys., 2014, 10, 451–456 CrossRef CAS.
  77. A. K. Geim and I. V. Grigorieva, Nature, 2013, 499, 419–425 CrossRef CAS PubMed.
  78. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186 CrossRef CAS.
  79. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775 CrossRef CAS.
  80. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 13, 5188–5192 CrossRef.
  81. A. V. Belikov, Y. E. Lozovik, A. G. Nikolaev and A. M. Popov, Chem. Phys. Lett., 2004, 385, 72–78 CrossRef CAS.
  82. Y. E. Lozovik and A. M. Popov, Chem. Phys. Lett., 2000, 328, 355–362 CrossRef CAS.
  83. Y. E. Lozovik and A. M. Popov, Phys. Solid State, 2002, 44, 186–194 CrossRef CAS.
  84. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799 CrossRef CAS PubMed.
  85. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104 CrossRef PubMed.
  86. A. M. Popov, I. V. Lebedeva, A. A. Knizhnik, Y. E. Lozovik, B. V. Potapkin, N. A. Poklonski, A. I. Siahlo and S. A. Vyrko, J. Chem. Phys., 2013, 139, 154705 CrossRef PubMed.
  87. J. E. Huheey, E. A. Keiter and R. L. Keiter, Inorganic Chemistry : Principles of Structure and Reactivity, HarperCollins, New York, USA, 4th edn, 1993 Search PubMed.
  88. I. V. Lebedeva, A. A. Knizhnik, A. M. Popov, Y. E. Lozovik and B. V. Potapkin, Phys. E, 2012, 44, 949–954 CrossRef CAS.
  89. E. Bichoutskaia, A. M. Popov, A. El-Barbary, M. I. Heggie and Y. E. Lozovik, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 71, 113403 CrossRef.
  90. E. Bichoutskaia, A. M. Popov, Y. E. Lozovik, O. V. Ershova, I. V. Lebedeva and A. A. Knizhnik, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 165427 CrossRef.
  91. A. M. Popov, Y. E. Lozovik, A. S. Sobennikov and A. A. Knizhnik, J. Exp. Theor. Phys., 2009, 108, 621–628 CrossRef CAS.
  92. A. M. Popov, I. V. Lebedeva and A. A. Knizhnik, Appl. Phys. Lett., 2012, 100, 173101 CrossRef.
  93. A. M. Popov, I. V. Lebedeva, A. A. Knizhnik, Y. E. Lozovik and B. V. Potapkin, J. Chem. Phys., 2013, 138, 024703 CrossRef PubMed.
  94. P. Pulay, Chem. Phys. Lett., 1980, 73, 393–398 CrossRef CAS.
  95. Q. Peng, W. Ji and S. De, Comput. Mater. Sci., 2012, 56, 11–17 CrossRef CAS.
  96. K.-A. N. Duerloo, M. T. Ong and E. J. Reed, J. Phys. Chem. Lett., 2012, 3, 2871–2876 CrossRef CAS.
  97. D. Boschetto, L. Malard, C. H. Lui, K. F. Mak, Z. Li, H. Yan and T. F. Heinz, Nano Lett., 2013, 13, 4620–4623 CrossRef CAS PubMed.
  98. P. H. Tan, W. P. Han, W. J. Zhao, Z. H. Wu, K. Chang, H. Wang, Y. F. Wang, N. Bonini, N. Marzari, N. Pugno, G. Savini, A. Lombardo and A. C. Ferrari, Nat. Mater., 2012, 11, 294–300 CrossRef CAS PubMed.
  99. L. Duclaux, B. Nysten, J.-P. Issi and A. W. Moore, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 46, 3362–3367 CrossRef CAS.
  100. Z. Xu, X. Li, B. I. Yakobson and F. Ding, Nanoscale, 2013, 5, 6736–6741 RSC.

This journal is © The Royal Society of Chemistry 2016
Click here to see how this site uses Cookies. View our privacy policy here.