Van der Waals interaction in a boron nitride bilayer

We have carried out quantum Monte Carlo (QMC) calculations to study the interlayer interaction in a boron nitride (BN) bilayer. The binding energy, 81 meV/2BN after finite-size corrections, was found to be larger than that obtained by density functional theory (DFT) with local density approximation, and smaller than those using van der Waals density functionals, both by considerable amounts. The QMC calculated interaction beyond the equilibrium interlayer separation was found to have a longer-range behavior than all the available DFT schemes.


Introduction
The van der Waals (vdW) interaction is one of the most fundamental physical quantities resulting from quantum fluctuation of charges. However, it remains a challenge to account for this interaction quantitatively in both theory and experiment. The vdW interaction is usually described in terms of the − R 6 attraction between two neutral and polarizable particles separated by a distance ≫ R r, where r is the size of the particles. It is only recently that the direct measurement of the vdW interaction between two Rydberg atoms has been accomplished, which confirmed its − R 6 dependence and determined the strength [1]. By summing up the − R 6 contributions between the composing microscopic subsystems, the vdW interaction between two macroscopic bodies separated by a distance D can be determined. The distant attraction between two thin macroscopic layers is concluded to follow − D 4 according to the standard pairwise summation method. However, it has been pointed out that the interaction for zero-gap layer systems should decay more slowly [2]. This was certified later for the bilayer system of a two-dimensional (2D) homogeneous electron gas by quantum Monte Carlo (QMC) methods [3] and for graphite, using the adiabatic-connection fluctuation-dissipation theorem in the random phase approximation [4], though a different conclusion of − D 4.2 behavior for graphite was also reported using QMC calculations [5]. For insulating layer systems, it is usually presumed that the − D 4 behavior applies.
The vdW interaction is one of the physical properties that the standard local and gradient approximation (LDA [6] and GGA [7]) of density functional theory (DFT) [8] methods fail to describe correctly [9]. In recent years, there have been many proposed vdW density functionals (vdWDF) to overcome this deficiency [10]. However, discrepancies in binding energy among these vdWDF results are usually apparent, and these vdWDF schemes have been used to study behavior beyond the equilibrium distance. An alternative approach to tackle this issue is to perform one of the currently most accurate electronic calculations for materials, i.e., QMC methods [11]. In this article, we present the results of studying the BN bilayer system by QMC methods. We calculated the binding energy and studied the behavior of the vdW interaction beyond the equilibrium distance (but not yet reaching the asymptotic region) in this prototypical insulating bilayer system. These results are also compared with those obtained from the standard DFT as well as those employing vdW density functionals. The outcome of the present QMC study would provide a benchmark for future generation of vdWDF and guidance for prospective experiments.

Calculation methods
The QMC calculations were performed using the CASINO code [12]. The trial wave functions of Slater-Jastrow form are used, and the Slater orbitals are generated from the DFT calculations using the CASTEP [13] plane-wave-basis code. The tabulated relativistic Hartree-Fock pseudopotentials from the CASINO pseudopotential library are used to generate the initial wavefunctions for Slater orbitals [14,15], and Casulaʼs approach [16] is used in DMC calculations to minimize the nonlocal pseudopotential error. In the orbital-generation calculations, the plane-wave energy cutoff is taken as 2000 eV. This is much higher than those of the ordinary DFT calculations in order to reduce the variance in the QMC calculations [17].
The QMC calculations start with the variational quantum Monte Carlo method to obtain the optimized parameters in the Jastrow factor, which include isotropic electron-electron, electron-nucleus, and electron-electron-nucleus terms. That is, J can be written as where N and N ions are the numbers of electrons and ions, respectively, and r ij and r iI denote the distance between electrons i and j and also between electron i and ion I respectively [18]. The u, χ, and f terms are power expansions containing variational parameters and depend on the spins of the electrons. The optimization was performed with variance minimization [19] and followed by energy minimization [20]. The nodal surface of the wave function is constrained to that of the Slater determinant constructed from the DFT-generated orbitals, i.e., the fixed-node approximation [21]. The DFT-generated orbitals are represented numerically using a real-space splines grid in order to improve the system-size scaling of the QMC calculations [22]. The optimized wave functions are subsequently used to generate the initial ensemble of electronic configurations for the diffusion quantum Monte Carlo (DMC) calculation through which the ensemble is propagated by the imaginary-time Schrödinger equation toward the ground state [11]. The model periodic Coulomb interaction is used to reduce finite-size effects related to the long-ranged Coulomb interaction [23]. Throughout this study, the supercells consisting of the AA' stacking of two BN layers are used, with c = 30 Å as the size of the supercell perpendicular to the basal planes. The supercell with the equally spaced BN layers, i.e., D = c/2 = 15 Å, is taken as the reference system when the interlayer energy is evaluated. We expect excellent error cancellations in the interlayer interaction energy with respect to the k-point sampling in the DFT methods by using these corresponding reference systems 4 .
We have carried out QMC calculations using both the lateral (3 × 3) and (4 × 4) supercells in order to estimate the finite-size bias of the QMC methods. The orbital-generation DFT calculations used the shifted Monkhorst-Pack (MP) [24] k-point meshes of (3 × 3 × 1) and (4 × 4 × 1) in the primitive supercell for the two QMC supercells, respectively. The lateral (3 × 3) supercell with c = 20 Å was also considered in order to assess the vacuum-size correction. The intralayer lattice constant is fixed at 2.5 Å, which is the experimental value [25] as well as the equilibrium lattice constant of the LDA calculations. Time steps of 0.003 a.u. were used for the DMC production calculations to ensure an acceptance ratio of more than 99.9%. The target population used in all the DMC calculations is 3840. A total of more than 5 million core-hours are required in order to achieve a standard deviation of less than 2.5 meV/ 2BN for the DMC total energy of all the various supercell sizes considered in this study.
For comparison purposes, we also carried out the LDA, Perdew-Burke-Ernzerhof (PBE) [7] GGA and vdWDF calculations. We used the vdW density functionals proposed by Dion et al (vdW-DF) [26] as well as those with a different exchange functional [27,28] or correlation functional (vdW-DF2) [29], which are all implemented by the algorithm of Roman-Perez and Soler [30] in the plane-wave-based Vienna ab initio simulation program (VASP) [31,32]. The vdW-DF and vdW-DF2 are first-principles vdW functionals, while the generations of vdW-optXXX involve empirical fitting to a benchmark dataset of weakly interacting dimers and water clusters. The interaction between ions and valence electrons is described by the projector augmented wave method [33], and the numbers of valence electrons included for B and N are three and five, respectively. The energy cutoff that determines the number of plane waves is more than 900 eV [34], and the (9 × 9 × 1) MP k-point mesh is used. Tests on higher energy 4 For example, considering the shifted Monkhorst-Pack [24] (MP) k-point meshes of (3 × 3 × 1) and (4 × 4 × 1) in the LDA [6] calculations of the lateral (1 × 1) supercell, the bilayer binding energy differs only by 0.2 meV/2BN for these two k-point meshes, although the total energy at the equilibrium distance (3.25 Å in both cases) differs by 95 meV/2BN. cutoffs, denser MP meshes, and larger vacuum size show that the binding energy is converged to within 1 BN meV 2 , and the equilibrium distance stays the same.

DFT results
The interaction energies of all the different DFT schemes are plotted in figure 1 and their binding energy and equilibrium distances (D 0 ) are presented in table 1. The GGA scheme leads to a tiny binding energy and too-large equilibrium distance, as demonstrated in previous works [29,35]. The LDA approach, performed by using VASP and CASTEP (with different pseudopotentials as stated in the previous section), yields binding energies of around 55 meV/ 2BN and equilibrium distance of 3.3 Å, which is close to the interlayer distance of the The calculated interaction energies of the BN bilayer as a function of interlayer distance by the different DFT schemes, as explained in the text. Also plotted for comparison is one of the calculated DMC interaction curves whose notation is explained in figure 2 and table 2. experimental value of 3.33 Å for the hexagonal BN bulk [25]. The binding energy obtained from both the vdW-DF and vdW-DF2 schemes are about 40 meV/2BN larger than those of the LDA scheme. A similar result of 100 meV/2BN was published previously [35], using the earlier version of the vdW-DF scheme (denoted as vdW-DF-layered) [36]. The vdWDF results based on the vdW-DF version but with different exchange functionals (optPBE, optB88, and optB86b) [27] lead to an extensive range in binding energy, i.e. from 116-127 meV/2BN. However, the largest binding energy of 152 meV/2BN has also been reported previously [37] by using a different vdWDF, i.e., vdW-TS [38]. It will be demonstrated in the following section that the two first-principles vdW functions, i.e., vdW-DF and vdW-Df2, give binding energies closest to the QMC results.

QMC results
The interaction energy calculated by the DMC method is plotted in figure 2. The binding energy close to the equilibrium interlayer distance obtained by the DMC method are listed in table 2. Due to the enormous amount of CPU time required in the QMC calculations, only selected data points at quarters of integer Ds (in unit of Å) are included in the QMC study. The lowest-energy interlayer distance was found to be either 3.25 Å or 3.50 Å. The binding energy of the larger supercell (4 × 4) was found to be around 72 meV/2BN while that of the smaller cell around 66 meV/2BN. Using the similar procedure as used previously for a Ne solid [39], the binding energy extrapolated to the system with an infinite lateral size is 79.3 meV/2BN. Extrapolation for the binding energy with an infinite vacuum between two bilayer systems of neighboring supercells can be estimated by assuming the reference energy at D = c/2 decreases as D  The DFT schemes closest to this value are the two first-principles vdW functionals, i.e., vdW-DF and vdW-DF2. All the results in table 2 used the PBE-generated orbitals except the one denoted with an asterisk, which used LDA-generated orbitals. However, we see that these two different initial many-electron wavefunctions lead to about the same DMC evaluated binding energies.

Behavior beyond the equilibrium interlayer separation
In addition to the binding energy, another desirable property of the BN bilayer system is its asymptotic behavior at large interlayer separations. Nevertheless, a full convergence is unreachable in practice at this moment. As shown earlier for a graphene bilayer [4], the interlayer separation must be larger than 9 Å to see a correct asymptotic behavior. In the present system, the interaction is already smaller than the numerical error at the 7 Å interlayer separation. However, an analysis for the behavior beyond equilibrium distance (bD 0 ) predicted by the DMC and different DFT schemes would still shed light on the long-range behavior of this system. Note that in the conventional two-term description (one repulsive and the other attractive term) of the interaction, the attractive interaction becomes dominant already at interlayer separations around 1.3 D 0 , as discussed in the following. For the 6-12 Lennard-Jones potential, the repulsive energy is about one tenth, in magnitude, of the attractive energy at D ≈ 1.3D 0 . For the 4-10 potential (resulting from the summation of the 6-12 potential for 2D bilayer systems), a similar condition can be reached at D ≈ 1.25D 0 . If the equilibrium distance is taken as 3.50 Å, the long-range attractive behavior is dominant for D beyond 4.55 Å and 4.375 Å, respectively, for the 6-12 and 4-10 potentials. For the three fitted Buckingham-like potentials shown in figure 2, the repulsive energies are already less than one seventeenth of the attractive energy at D = 4.5 Å, in magnitude.
The behavior beyond the equilibrium separation for the different DFT schemes (figure 1) can be sorted into four groups, i.e., PBE, LDA, vdW-DF2 and that including vdW-DF, vdW_optPBE, vdW_optB88, and vdW_optB86b. For the 4th group, we will take the vdW_optPBE as the representative scheme for later comparison with those by the QMC methods. The bD 0 behavior for these four DFT schemes are fitted with both the power law and  6 3.50 exponential form. We found that, overall, the exponential form fitted better for these four DFT schemes. However, the power-law fitting provides us with a direct comparison to the D −4 law. All these interaction energies are plotted on a logarithmic scale in order to better display their power-law behaviors ( figure 3). The bD 0 features for both the PBE and LDA results are extremely short-range with the D −10 dependence. The vdW-DF2 and vdW_optPBE schemes result in the D −6 and D −5.3 dependence, respectively, a considerably longer-range behavior as compared to the LDA scheme. For the DMC obtained results, the behavior beyond equilibrium distance refers to the interaction between 4 Å and 7 Å interlayer separation which includes nine data points for DMC (3 × 3) and five data points for both DMC(4 × 4), and DMC(4 × 4)*. The best fitted exponents (using the least absolute deviation fitting) are −4.50, −4.48, and −4.49 for DMC(3 × 3), DMC (4 × 4), and DMC(4 × 4)*, respectively. The average residues are all smaller than the error of the DMC total energy, i.e. 2.2, 0.9, and 2.0 meV/2BN for DMC(3 × 3), DMC(4 × 4) and DMC (4 × 4)* respectively. The bD 0 DMC data can therefore be described as being best fitted with the D −4.5 dependence. The previous study for a graphene bilayer [4] showed that the exponent can be increased (decreased in magnitude) by as much as 1.2 for interaction beyond 9 Å when compared to those of smaller separations (but beyond equilibrium distance). To what extent this conclusion can be applied to the present system is not apparent, as graphene is a zero-gap layer, while BN has a wide band gap. However, the presently obtained D −4.5 dependence for the behavior beyond equilibrium separation is shown to be longer-range than all the available DFT schemes.

Conclusion
In conclusion, we have carried out QMC calculations to determine both the binding energy and distance-dependent power law of the interaction beyond the equilibrium separation in the BN bilayer. The binding energy obtained by the DMC method, i.e., 81 meV/2BN, was found to be larger than the LDA value and smaller than all the considered and previously published vdWDF ones by sizable quantities. We determined the power law of the BN interlayer interaction beyond equilibrium distance as D −4.5 , while those of the available DFT schemes all fall off faster than this dependence.