{'fraction_non_alphanumeric': 0.031313131313131314, 'fraction_numerical': 0.01164983164983165, 'mean_word_length': 4.897934868943606, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stacking and Registry Effects in Layered Materials: The Case of Hexagonal Boron Nitride

The interlayer sliding energy landscape of hexagonal boron nitride (h-BN) is investigated via a van der Waals corrected density functional theory approach. It is found that the main role of the van der Waals forces is to"anchor"the layers at a fixed distance, whereas the electrostatic forces dictate the optimal stacking mode and the interlayer sliding energy. A nearly free-sliding path is identified, along which bandgap modulations of ~0.6 eV are obtained. We propose a simple geometrical model that quantifies the registry matching between the layers and captures the essence of the corrugated h-BN interlayer energy landscape. The simplicity of this phenomenological model opens the way to the modeling of complex layered structures, such as carbon and boron nitride nanotubes.

The interlayer potential landscape in layered materials is essential for understanding their mechanical and electromechanical behavior. For example, nanoelectromechanical systems (NEMS) based on low dimensional structures of such materials often rely on mechanical deformations such as twisting [1][2][3] and bending [4][5][6]. These processes involve relative sliding of the layers, which exhibits a corrugated energy landscape even in atomically flat systems, such as graphite and hexagonal boron nitride (h-BN). This corrugation arises from the non-uniform charge density distribution around the atomic positions within each layer [7][8][9]. It is well accepted that the most important factors that govern corrugation in such system are electrostatic and dispersion interactions. However, a clear picture of how the complex interplay between these factors determines the corrugated energy landscape and its manifestation in unique material properties has not emerged yet.
Previous efforts towards the understanding of these phenomena have utilized density functional theory (DFT) with (semi-)local approximations [9][10][11][12]. However, these methods do not provide an appropriate description of dispersive interactions. Several approaches within DFT have been developed to overcome this problem, demonstrating successful treatment of dispersion effects in graphite [13][14][15] and h-BN [13,14,16]. However, to the best of our knowledge, such studies have not addressed the question of corrugation.
In this letter, we present a theoretical study of the intricate interplay between dispersion and electrostatic interactions in layered materials. As an illustrative example we choose the case of h-BN, where both types of interactions are expected to have a considerable influence on the landscape of the interlayer potential. Using a first-principles van der Waals (vdW) corrected DFT approach, we show that dispersion interactions play the role of fixing the interlayer distance, while the electrostatic forces determine the optimal stacking mode and the in-terlayer sliding corrugation. In addition, we predict the existence of a nearly free-sliding path along which considerable bandgap modulations are obtained. Finally, we propose a simple geometrical model that quantifies the registry matching between the layers and captures all the important physical features appearing in the corrugation energy landscape of h-BN. We start by demonstrating the importance of appropriate treatment of dispersion interactions when modeling layered materials in general and h-BN in particular. To this end, we calculate the interlayer distance dependence of the binding energy (BE) of h-BN using the generalized gradient approximation (GGA) of Perdew, Burke, and Ernzerhof (PBE) [17] augmented with the Tkatchenko-Scheffler vdW (TS-vdW) correction [18,19]. Within this approach, the leading order pairwise C 6 /R 6 dispersion correction is added to the internuclear energy term, where the C 6 coefficients and the vdW radii are determined directly from the DFT ground state electron density and reference values for the free atoms.
The geometry of a single BN sheet was optimized with the PBE functional using two-dimensional periodic boundary conditions, [20] as implemented in the Gaussian suite of programs [21] using the large quadruplezeta Weigend-Ahlrichs basis set (Def2QZVP) [22]. Unit cells of bilayer and bulk h-BN were then constructed by stacking single BN sheets in the experimentally observed AA' stacking mode (see Fig. 2b). The BE of these structures was calculated as a function of the interlayer distance at fixed intralayer coordinates with and without the TS-vdW correction, as implemented in the FHI-aims code [23]. We used the tier-2 numerical atomiccentered orbitals basis set, known to yield binding energies converged to the meV/atom level [19], verified here by comparison with tier-3 basis set calculations for selected stacking modes.
As shown in Fig. 1, the PBE+vdW equilibrium interlayer distance calculated for bulk h-BN (solid blue line) is 3.33Å, in perfect agreement with the experimental value [24]. This is in stark contrast to the result obtained from the pure PBE functional: an overestimated equilibrium interlayer distance of 4.17Å accompanied by a very modest BE of 2.05 meV/atom (dashed purple line) [25,26]. The lattice constant obtained with the PBE+vdW approach is consistent with the value of 3.31 A [16] obtained previously using the adiabatic-connection fluctuation-dissipation theorem, but here the computational cost is significantly reduced. The PBE+vdW binding energy for bulk h-BN is 85.9 meV/atom, which is somewhat higher then the BE of 56 meV/atom recently reported for graphite [15].
Similarly, in the case of bilayer h-BN, PBE yields a small BE of 1.0 meV/atom and an overestimated interlayer distance of 4.22Å, whereas PBE+vdW predicts a much larger BE of 38.1 meV/atom and an interlayer distance of 3.37Å [27]. We note that the interlayer binding in bilayer h-BN is weaker than that obtained for bulk h-BN. This is due to the fact that in the bulk system each layer interacts with two adjacent layers, instead of one, as well as with additional layers that are farther away. We note that Rydberg et al. [13], reported a nearly indistinguishable behavior of bulk and bilayer h-BN. However, they obtained an overestimated interlayer distance of 3.63Å.
Having established the validity of PBE+vdW approach for the description of the interlayer coupling in h-BN, we now apply it to the study of corrugation in the bilayer system. Starting from the AA' stacking mode with an interlayer distance of 3.37Å, we perform a set of lateral shifts of one h-BN layer parallel to the basal plane of of the other. At each shifted configuration we calculate the total energy of the bilayer system. The resulting sliding energy landscape is presented in Fig. 2a.
In order to quantify the role of vdW interactions for the sliding process, we compare the changes in total en- ergies with the changes in the vdW contribution upon interlayer sliding. We find that the maximal total energy change obtained is 26 meV/unit-cell while the largest change in the vdW correction is 5 meV/unit-cell (not shown). Comparing with the vdW contribution to the BE, which is two orders of magnitude larger (see Fig. 1), we conclude that the main role of the vdW interactions is to "anchor" the layers at the appropriate interlayer distance. Once stacking is established, the sliding energy profile at fixed interlayer distance is governed by electro-static interactions resulting from the polar nature of the B-N bond [7,9]. These conclusions are further supported by the fact that the general shape of the sliding energy surface remains qualitatively unchanged when performing the calculation without the vdW correction, provided that the layers are kept at the same interlayer distance.
Apart from the stable AA' configuration, two other high symmetry configuration are obtained. We mark them as AB 1 (See Fig. 2b) and AB 2 in accordance with the standard nomenclature of graphite stacking. In the AB 1 mode, the boron atoms of the top layer lie exactly above those of the bottom layer, while the nitrogen atoms are located above the center of a BN hexagon. The AB 2 mode is similar but with interchanged positions of the boron and nitrogen atoms. Despite the similar geometry, the nature of the two stacking modes is surprisingly different [9]. The AB 1 mode is stable with a total energy comparable to that of the AA' mode, whereas the AB 2 mode is unstable with a total energy 26 meV/unit-cell higher than the AA' configuration. We note that with the pure PBE functional the AA' stacking mode is the global minimum of the sliding energy landscape, lower by only ∼ 3.5 meV/unit-cell than the AB 1 mode. When applying the vdW correction, the AB 1 stacking mode turns out to be lower by ∼ 1.4 meV/unit-cell than the AA' configuration. Because this difference is within the accuracy limits of our approximations we do not attribute any physical importance to it.
An interesting feature that appears in the sliding energy surface is the existence of a nearly free-sliding path, indicated by the yellow line appearing in Fig. 2a. Fig. 2b presents the total energy differences with respect to the AA' stacking mode along this path, showing total energy variations of ∼ 2 meV/unit cell. This suggests the h-BN will exhibit highly anisotropic mechanical and electromechanical properties.
The calculated corrugation energies presented above indicate that different stacking modes exhibit considerably different energetic stability. Due to the inverse correlation between stability and bandgap, one may expect that upon interlayer sliding the electronic properties of layered h-BN will vary. In order to quantify this electromechanical effect we calculate the bandgap of bilayer h-BN as a function of interlayer sliding. We perform these calculations using the screened-exchange hybrid approximation of Heyd, Scuseria, and Ernzerhof (HSE06) [21,28], which is expected to reproduce experimental optical bandgaps of bulk semi-conductors [29], and the doubleζ polarized 6-31G** Gaussian basis set [30]. This approach has been recently shown to describe the physical properties of graphene-based materials with exceptional success [31].
The sliding induced bandgap variations are presented in Fig. 3. A steep 0.6 eV bandgap decrease is obtained as the layers are laterally shifted away from the AA' configuration. This change is about 10% of the calculated bandgap value of 6.05 eV for the AA' stacking mode, indicating a strong electromechanical response of bilayer h-BN towards interlayer sliding. Furthermore, we find that in addition to the large variations in magnitude, the bandgap changes its nature from direct to indirect during the interlayer sliding process, in agreement with the calculations of Liu et al. [9]. It is interesting to note that large bandgap variations are obtained along the nearly free-sliding path, meaning that the electronic properties of h-BN will be sensitive to stress applied along this direction.
From what we have presented thus far, it is clear that there is an intimate relation between the registry matching of the different 2D crystalline sheets and the energetic stability of the various stacking modes in layered materials. Previous studies have addressed this issue only in qualitative terms regarding different configurations as having "good" or "bad" stacking [9].
Here, we present a simple and intuitive model that quantifies the registry matching using basic geometric considerations. As stated above, the main contribution to the corrugation energy in h-BN comes from the electrostatic interactions between the atomic sites, which carry a partial charge due to the polar nature of the B-N bond [7,9]. In order to mimic this effect, we ascribe to each atom in the unit cell a circle centered around its position. Considering the projection on a plane parallel to the layers, as illustrated in Fig. 2c, we mark by S ij the overlaps between the circle centered around the i atom belonging to the top layer and the j atom of the bottom layer (i and j being either N or B). We now define the registry index as: where the factors S AA ′ ij and S AA ij are the respective overlaps at the AA' and AA stacking modes, introduced for normalization purposes. Here, AA denotes the case were the two layers are completely eclipsed and in perfect registry. With this definition, the registry index, R, is lim-ited to the interval [0, 1] where the minimal (maximal) value is obtained for the AA' (AA) mode.
It is now possible to plot the registry index as a function of the lateral interlayer shifts while using the ratio between the radii of the N (r N ) and B (r B ) circles as a single fitting parameter to obtain good agreement between the registry index and the corrugation energy surfaces. In Fig. 2d we plot the registry index surface for r N = 0.50R BN and r B = 0.15R BN , where R BN = 1.45Å is the equilibrium BN bond length in h-BN. This ratio between the two radii takes into account the non-uniform charge distributions around the B and N atomic sites where the nitrogen (boron) atom has a larger (smaller) electron cloud around it. Clearly, a good agreement is obtained between the registry index calculated via the simple geometric model and the corrugation energy surface obtained from state-of-the art DFT calculations. Therefore, we conclude that it is possible to capture all the important physical parameters that govern the corrugation process using simple geometric considerations [32].
In summary, we have studied the corrugated sliding energy landscape in layered h-BN via a van der Waals corrected DFT approach. The delicate interplay between different contributions was investigated, revealing that the main role of the van der Waals forces is to anchor the layers at a fixed distance, whereas the electrostatic forces dictate the optimal stacking mode and the interlayer sliding corrugation. A nearly free-sliding path has been identified, along which bandgap modulations of ∼ 0.6 eV are predicted to occur. A simple geometrical model that quantifies the registry matching between the layers and captures the essence of h-BN interlayer sliding was presented. The model is not limited to the case of h-BN and can be easily extended to treat other layered materials, such as graphite, and generalized to describe more complex structures such as multilayered nanotubes of different types. Following the guidelines presented in this letter, it is possible to identify stable configurations and quantify the corrugation in complex layered systems based on simple geometrical arguments with a negligible computational effort. This work was supported by the Israel Science Foundation, the Israeli Ministry of Defense, the Minerva Foundation, and the Center for Nanoscience and Nanotechnology at Tel-Aviv University.