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

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 phenomena^1–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,⁸ diffusion and drift of a graphene flake on a graphene surface via rotation to incommensurate states^9,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 electronic^22–25,27 and optical²⁶ 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 displacement^11–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,³⁵ shear mode frequency^37,38 (E_2g mode with adjacent layers sliding rigidly in the opposite in-plane directions) and frequency of relative out-of-plane vibrations (B_1g ZO mode with adjacent layers sliding rigidly towards and away from each other)³⁹ 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.⁴² Rotational stacking faults have been observed for few-layer h-BN.⁴² 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,⁵³ shear mode frequency^46,52,53 and frequency of out-of-plane relative vibrations^{39,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.⁴³ The dependence of the interlayer binding energy on the layer stacking has been analyzed both for bulk^49,51,55 and bilayer h-BN^54,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 approximation⁵⁸ 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 constant⁶² 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.⁶³ 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 functional⁶⁷). 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 approach⁷² 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),⁵⁵ 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 discovered^73–76 graphene-h-BN heterostructures (see also ref. 77 for a review) showing the commensurate–incommensurate phase transition.⁷⁶ The calculations performed here represent a necessary intermediate stage (along with the previous calculations for bilayer graphene^60,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 code⁷⁸ with the non-local vdW-DF2 functional⁷² 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).⁷⁹ 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 method⁸⁰ 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.⁶⁰ 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⁻⁶ eV.

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

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 walls⁸¹ 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 calculations⁵⁵ 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 surface^59,60,64 and the barriers to relative motion of the layers⁶⁰ 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.⁵⁵ 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 graphene^{7,59–61,64–66} and close to the LMP2 value ∼2.5 meV per atom for h-BN.⁵⁵ 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 graphene^{7,59–61,64–66} and for bilayer⁵⁵ and bulk^49,51,55 h-BN as well as the LMP2 values of about 20 meV per atom⁵⁵ for bilayer and bulk h-BN.

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-BN²⁸ as well as with the LMP2 calculations⁵⁵ and DFT calculations with account of nonlocal many-body dispersion (MBD)⁵⁷ 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 results⁵⁵ (Tables 1 and 2) and the experimental data⁴² 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 exchange⁷²) 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 bilayer^59,60,66 and h-BN bilayer^54,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 as⁵⁸:

U = 0.5(U₀ + U₁(2 cos(k₁x)cos(k₂y) + cos(2k₁x) + 1.5)), (1)

where k₁ = 2π/(3l), 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 U_1NB = U_1BN and U_0NB = U_0BN. Then for h-BN layers aligned in the same direction we arrive at:

U_I = U_0,tot + 4.5U_1NB + (U_1NN + U_1BB − U_1NB)×(2 cos (k₁x) cos (k₂y) + cos (2k₁x) + 1.5), (2)

where the point x = 0 and y = 0 corresponds to the AA stacking and U_0,tot = 2U_0NB + U_0NN + U_0BB. 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 (U_1NN + U_1BB − U_1NB), 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,

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 (2U_1NB − 0.5U_1NN − 0.5U_1BB) and (U_1NN − U_1BB) 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

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

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 E_AA − E_AB = 4.5(U_1NN + U_1BB − U_1NB) = 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 E_SP − E_AB = 0.5(U_1NN + U_1BB − U_1NB) = 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 layers^60,61,66,88 and between carbon nanotube walls, both for infinite commensurate walls^81,89–92 and in the case where corrugations of the potential surface are determined by the contribution of edges⁹³ or defects.⁸¹ 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 method⁹⁴ up to a residual force of 6 × 10⁻³ eV Å⁻¹. 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⁻² 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⁻², ν = 0.22 (ref. 95) and 291 J m⁻², ν = 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 E_2g 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 state⁶¹:where μ = p(m_N + m_B)/4 for n = 2p and μ = p(p + 1)(m_N + m_B)/2(2p + 1) for n = 2p + 1, m_B and m_N are masses of boron and nitrogen atoms, respectively, p is an integer, l is the bond length and we denote U_eff = (l/2π)²∂²U/∂x². Curvatures of the potential energy surfaces in the states corresponding to the AA′, AB1′ and AB stackings are described by the parameters U_eff(AA′) = (U_1NN + U_1BB − 4U_1NB)/3 = 1.40 meV per atom, U_eff(AB1′) = (U_1BB + 2U_1NB − 2U_1NN)/3 = 0.71 meV per atom and U_eff(AB) = (U_1NN + U_1BB − U_1NB)/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 data^37,38 and the results of previous DFT calculations^46,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

Fig. 2 Calculated shear mode frequencies f (in cm⁻¹) 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:

where 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 experimental^35,99 and theoretical values⁵³ for bulk h-BN.

Fig. 3 Calculated shear moduli C₄₄ (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 calculations⁵³ 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:

U_rot = 〈U_I〉_x,y = 〈U_II〉_x,y = U_0,tot + 1.5(U_1NN + U_1BB + 2U_1NB) (7)

Thus, the barriers to relative rotation of h-BN layers from metastable states corresponding to the AA′, AB1′ and AB stackings are E_rot(AA′) = U_rot − U_II(AA′) = 1.5(U_1NN + U_1BB − 4U_1NB) = 6.29 meV per atom, E_rot(AB1′) = −1.5U_1NN + U_1BB + 3U_1NB = 2.29 meV per atom and E_rot(AB) = 1.5(U_1NN + U_1BB) − U_1NB = 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 E_rot ∼ 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.¹⁰⁰

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 electronic^22–25,27 and optical²⁶ 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.⁷ 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.⁶³

To consider dislocations in bilayer h-BN we follow the formalism of the two-chain Frenkel–Kontorova model⁶² which was applied previously to study dislocations in double-walled carbon nanotubes^62,91 and graphene.⁶³ Elementary dislocations in graphene correspond to partial dislocations with the Burgers vector b⃑ 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 ( 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⃑ normal to the boundary between commensurate domains,⁶³ 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⃑ at an arbitrary angle β to the normal n⃑ 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.

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⃑ normal to the boundary between commensurate domains (black dashed line) and Burgers vector b⃑ are indicated. The magnitudes of the vectors are scaled up for clarity. (a and b) The directions of the Burgers vectors b⃑′₁ and b⃑′₂ corresponding to two straight pieces of the path of full dislocations (b⃑ = b⃑′₁ + b⃑′₂) 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⃑ 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 l_D 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⃑ of atoms of the layers can be formulated as u⃑ = 0 at x → −∞ and u⃑ = −b⃑ at x → +∞, where x is the coordinate along the normal n⃑ to the boundary between commensurate domains and b⃑ 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,

where the axis y is chosen along the boundary between commensurate domains, u⃑′ = du⃑/dx describes the strain across the boundary, V(u⃑) 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 − ν²) 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 − ν²)⁻¹ 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.⁶³ 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 − ν²)⁻¹ = 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⃑ 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 U_D (eqn (8)), should satisfy the Euler–Lagrange equations:

This means that the dislocation path is the same as the trajectory of a particle on the inverse potential energy surface −V(u⃑) 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:

here we use an integration constant of zero for isolated dislocations⁶³ and denote , where θ is the angle between the direction u⃑′ of the dislocation path and normal n⃑ 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:

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:

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 Å⁻¹ (ref. 63).

Fig. 5 Calculated formation energy of dislocations U_D per unit width (in eV Å⁻¹) as a function of angle β (in degrees) between the Burgers vector b⃑ and normal n⃑ 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.⁶³ 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⃑′|_max, the dislocation width can be defined in a similar way to the analytical solution:

where V_max 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 graphene^7,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 graphene⁶³ of 12 nm.

Fig. 6 Calculated dislocation width l_D (in nm) as a function of angle β (in degrees) between the Burgers vector b⃑ and normal n⃑ 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 result⁶³ 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⃑ 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⃑′(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.⁷

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.⁷² 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 observation⁴² 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.⁵⁵ 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-BN^54,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 bilayer^60,61,66,88 and double-walled carbon nanotubes^81,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.⁵⁴ 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 graphene^{4–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 l_D and barrier V_max to relative sliding of the layers V_max = C(β)/l_D²(β), where the coefficient C(β) = ϕ²(β)kl²/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.
2. Y. Gan, W. Chu and L. Qiao, Surf. Sci., 2003, 539, 120–128.
3. J. H. Warner, M. H. Rümmeli, T. Gemming, B. Büchner and G. A. D. Briggs, Nano Lett., 2009, 9, 102–106.
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.
5. M. Dienwiebel, N. Pradeep, G. S. Verhoeven, H. W. Zandbergen and J. W. M. Frenken, Surf. Sci., 2005, 576, 197–211.
6. A. E. Filippov, M. Dienwiebel, J. W. M. Frenken, J. Klafter and M. Urbakh, Phys. Rev. Lett., 2008, 100, 046102.
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.
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.
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.
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.
11. R. Bistritzer and A. H. MacDonald, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 245412.
12. R. Bistritzer and A. H. MacDonald, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 12233–12237.
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.
14. K.-T. Lam, C. Lee and G. Liang, Appl. Phys. Lett., 2009, 95, 143107.
15. Y. Qian, K.-T. Lam, C. Lee and G. Liang, Carbon, 2012, 50, 1659–1666.
16. J. Zheng, P. Guo, Z. Ren, Z. Zhenyi, J. Bai and Z. Zhang, Appl. Phys. Lett., 2012, 101, 083101.
17. K. K. Paulla and A. A. Farajian, J. Phys.: Condens. Matter, 2013, 25, 115303.
18. L. Brown, R. Hovden, P. Huang, M. Wojcik, D. A. Muller and J. Park, Nano Lett., 2012, 12, 1609–1615.
19. B. Butz, C. Dolle, F. Niekiel, K. Weber, D. Waldmann, H. B. Weber, B. Meyer and E. Spiecker, Nature, 2013, 505, 533–537.
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.
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.
22. S. Hattendorf, A. Georgi, M. Liebmann and M. Morgenstern, Surf. Sci., 2013, 610, 53–58.
23. P. San-Jose, R. V. Gorbachev, A. K. Geim, K. S. Novoselov and F. Guinea, Nano Lett., 2014, 14, 2052–2057.
24. B. Lalmi, J. C. Girard, E. Pallecchi, M. Silly, C. David, S. Latil, F. Sirotti and A. Ouerghi, Sci. Rep., 2014, 4, 4066.
25. M. M. Benameur, F. Gargiulo, S. Manzeli, G. Autès, M. Tosun, O. V. Yazyev and A. Kis, Nat. Commun., 2015, 6, 8582.
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.
27. M. Koshino, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 115409.
28. R. S. Pease, Nature, 1950, 165, 722–723.
29. R. S. Pease, Acta Crystallogr., 1952, 5, 356–361.
30. R. W. Lynch and H. G. Drickamer, J. Chem. Phys., 1966, 44, 181.
31. V. L. Solozhenko, G. Will and F. Elf, Solid State Commun., 1995, 96, 1–3.
32. V. L. Solozhenko and T. Peun, J. Phys. Chem. Solids, 1997, 58, 1321–1323.
33. V. L. Solozhenko, A. G. Lazarenko, J.-P. Petitet and A. V. Kanaev, J. Phys. Chem. Solids, 2001, 62, 1331–1334.
34. W. Paszkowicz, J. B. Pelka, M. Knapp, T. Szyszko and S. Podsiadlo, Appl. Phys. A, 2002, 75, 431–435.
35. A. Bosak, J. Serrano, M. Krisch, K. Watanabe, T. Taniguchi and H. Kanda, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 041402.
36. K. Fuchizaki, T. Nakamichi, H. Saitoh and Y. Katayama, Solid State Commun., 2008, 148, 390–394.
37. R. J. Nemanich, S. A. Solin and R. M. Martin, Phys. Rev. B: Condens. Matter Mater. Phys., 1981, 23, 6348–6356.
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.
39. A. Marini, P. García-González and A. Rubio, Phys. Rev. Lett., 2006, 96, 136404.
40. C. Li, Y. Bando, C. Zhi, Y. Huang and D. Golberg, Nanotechnology, 2009, 20, 385707.
41. C. Zhi, Y. Bando, C. Tang, H. Kuwahara and D. Golberg, Adv. Mater., 2009, 21, 2889–2893.
42. J. H. Warner, M. H. Rümmeli, A. Bachmatiuk and B. Büchner, ACS Nano, 2010, 4, 1299–1304.
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.
44. K. Albe, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 55, 6203–6210.
45. G. Kern, G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 8551–8559.
46. N. Ohba, K. Miwa, N. Nagasako and A. Fukumoto, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 63, 115207.
47. A. Janotti, S.-H. Wei and D. J. Singh, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 64, 174107.
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.
49. L. Liu, Y. P. Feng and Z. X. Shen, Phys. Rev. B: Condens. Matter Mater. Phys., 2003, 68, 104102.
50. T. Tohei, A. Kuwabara, F. Oba and I. Tanaka, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 064304.
51. N. Ooi, V. Rajan, J. Gottlieb, Y. Catherine and J. B. Adams, Modell. Simul. Mater. Sci. Eng., 2006, 14, 515–535.
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.
53. I. Hamdi and N. Meskini, Phys. B, 2010, 405, 2785–2794.
54. N. Marom, J. Bernstein, J. Garel, A. Tkatchenko, E. Joselevich, L. Kronik and O. Hod, Phys. Rev. Lett., 2010, 105, 046801.
55. G. Constantinescu, A. Kuc and T. Heine, Phys. Rev. Lett., 2013, 111, 036104.
56. C.-R. Hsing, C. Cheng, J.-P. Chou, C.-M. Chang and C.-M. Wei, New J. Phys., 2014, 16, 113015.
57. W. Gao and A. Tkatchenko, Phys. Rev. Lett., 2015, 114, 096101.
58. G. S. Verhoeven, M. Dienwiebel and J. W. M. Frenken, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 70, 165418.
59. M. Reguzzoni, A. Fasolino, E. Molinari and M. C. Righi, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 245434.
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.
61. A. M. Popov, I. V. Lebedeva, A. A. Knizhnik, Y. E. Lozovik and B. V. Potapkin, Chem. Phys. Lett., 2012, 536, 82–86.
62. E. Bichoutskaia, M. I. Heggie, Y. E. Lozovik and A. M. Popov, Fullerenes, Nanotubes, Carbon Nanostruct., 2006, 14, 131–140.
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.
64. A. N. Kolmogorov and V. H. Crespi, Phys. Rev.B: Condens. Matter Mater. Phys., 2005, 71, 235415.
65. M. Aoki and H. Amawashi, Solid State Commun., 2007, 142, 123–127.
66. O. V. Ershova, T. C. Lillestolen and E. Bichoutskaia, Phys. Chem. Chem. Phys., 2010, 12, 6483–6491.
67. M. Dion, H. Rydberg, E. Schröder, D. C. Langreth and B. I. Lundqvist, Phys. Rev. Lett., 2004, 92, 246401.
68. L. A. Girifalco and R. A. Lad, J. Chem. Phys., 1956, 25, 693–696.
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.
70. R. Zacharia, H. Ulbricht and T. Hertel, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 155406.
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.
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.
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.
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.
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.
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.
77. A. K. Geim and I. V. Grigorieva, Nature, 2013, 499, 419–425.
78. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
79. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
80. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 13, 5188–5192.
81. A. V. Belikov, Y. E. Lozovik, A. G. Nikolaev and A. M. Popov, Chem. Phys. Lett., 2004, 385, 72–78.
82. Y. E. Lozovik and A. M. Popov, Chem. Phys. Lett., 2000, 328, 355–362.
83. Y. E. Lozovik and A. M. Popov, Phys. Solid State, 2002, 44, 186–194.
84. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
85. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
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.
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.
88. I. V. Lebedeva, A. A. Knizhnik, A. M. Popov, Y. E. Lozovik and B. V. Potapkin, Phys. E, 2012, 44, 949–954.
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.
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.
91. A. M. Popov, Y. E. Lozovik, A. S. Sobennikov and A. A. Knizhnik, J. Exp. Theor. Phys., 2009, 108, 621–628.
92. A. M. Popov, I. V. Lebedeva and A. A. Knizhnik, Appl. Phys. Lett., 2012, 100, 173101.
93. A. M. Popov, I. V. Lebedeva, A. A. Knizhnik, Y. E. Lozovik and B. V. Potapkin, J. Chem. Phys., 2013, 138, 024703.
94. P. Pulay, Chem. Phys. Lett., 1980, 73, 393–398.
95. Q. Peng, W. Ji and S. De, Comput. Mater. Sci., 2012, 56, 11–17.
96. K.-A. N. Duerloo, M. T. Ong and E. J. Reed, J. Phys. Chem. Lett., 2012, 3, 2871–2876.
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.
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.
99. L. Duclaux, B. Nysten, J.-P. Issi and A. W. Moore, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 46, 3362–3367.
100. Z. Xu, X. Li, B. I. Yakobson and F. Ding, Nanoscale, 2013, 5, 6736–6741.

---