Graphene-hexagonal boron nitride van der Waals heterostructures: an examination of the effects of different van der Waals corrections

The structural and electronic properties of graphene on hexagonal boron nitride (hBN) as 2D van der Waals heterostructure were calculated using density functional theory method with van der Waals corrections. Four van der Waals corrections, along with conventional PBE-DFT, were utilized: the inter-atomic potentials-based DFT-D2, DFT-D3, Tkatchenko-Sheffler (TS), and the ab-initio, non-local correlation terms-based vdW-DF2-B86R. Results show that the structural properties of the 2D heterostructure, especially the inter-monolayer spacing, are consistent with previous theoretical works. In terms of energetics, the conventional PBE-DFT functional resulted to no binding between graphene and hBN, while utilizing the TS correction resulted to graphene-hBN adhesion energy value that is consistent with previous theoretical and experimental works. Electronic structure wise, the conventional PBE-DFT essentially predicted a zero-gap graphene on hBN, while all calculations involving van der Waals corrections resulted to band gaps that are consistent with previous studies. However, with the exception of TS, all van der Waals corrections predicted a Dirac cone that is shifted upward in energy from the Fermi level, making graphene artificially p-doped. As such, TS is recommended as one of the most appropriate van der Waals corrections for graphene-hBN 2D heterostructure. This work demonstrated the variations in graphene-hBN electronic properties as a result of the different implementations of the van der Waals corrections, but could be as useful as the more expensive theoretical methods such as GW.


Introduction
Since its mechanical exfoliation in 2004, graphene has found numerous applications particularly in the field of nanotechnology innovations, by virtue of its outstanding mechanical and electronic properties [1][2][3]. Moreover, graphene is widely utilized in a variety of other applications, from inert gas sensors and ballistic electronic devices to high temperature superconductor [4]. Unfortunately, graphene does not possess a band gap that is important in many applications requiring semiconducting properties. In this respect, graphene is known to interact well with a variety of substrates, particularly SiC, O-terminated quartz and alumina, in which graphene's electronic structure and properties can be perturbed and a band gap opens at its Dirac cone [5]. Although opening a gap in graphene is necessary for it to be useful in digital-related applications, maintaining graphene's inherent electronic properties such as its high carrier mobility is also important. As such, a graphene-substrate heterostructure with relatively weak interaction is usually preferred [6][7][8][9]. With this, graphene's intrinsic properties are retained and at the same time, a band gap opening is still possible. To circumvent these, graphene along with other 2D materials are utilized in building next generation of materials known as 2D van der Waals (vdW) heterorstructures [10][11][12].
In recent years, important surfaces and two-dimensional materials (2D) become a center of attention for researchers due to their excellent properties and potential applications in numerous diverse areas [13][14][15][16][17][18]. In Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. this regard, building a vdW heterostructure using 2D materials, in particular graphene and hexagonal boron nitride (hBN), is increasingly gaining attention [19][20][21][22][23][24][25][26][27][28]. An insulating 2D material with similar structure as that of graphene, hBN possesses a wide band gap (5.97 eV) and can be fabricated by exfoliation from its parent layered bulk BN [29][30][31]. As each hBN monolayer is stacked via weak vdW forces to form bulk BN and due to the absence of dangling bonds (similar to other 2D materials such as phosphorene), it is expected that monolayer hBN would induce weak interaction with graphene thereby retaining the latter's important electronic properties [19,[32][33][34]. In addition, graphene-hBN heterostructure has found many applications including, but not limited to, chemical and pressure sensors [35,36], optical field confinement [37], and spintronics [38]. Moreover, graphene-hBN heterostructure provides a test bed in observing and analyzing unique physics in 2D materials such as highly confined low-loss plasmons [37], diverging thermal conductivity [39], enhanced thermoelectic properties [40], spin tunnel barrier [38], vdW-induced Moiré pattern [41], and photoinduced doping [42].
Recently, a number of methods have arose to properly account for the vdW or dispersion forces within DFT, in an accurate but computationally less expensive ways. This is important especially when dealing with vdW heterostructures, as these could have profound effect on their electronic structure and properties. There are currently two major methodologies of these corrections. One is where the dispersion forces correction based on the interatomic potentials of the form C 6 ·R −6 is added to the total Kohn-Sham energy from DFT. The way the C 6 ·R −6 terms are calculated varies. These could be based on fixed semi-empirical results such as in DFT-D [43], environment-dependent, on-the-fly calculations such as in DFT-D3 [44] or based on self-consistent, electronic density-based calculations such as in Tkatchenko-Sheffler (TS) method [45]. The other major one is through the development of exchange-correlation functionals wherein a proper description of dispersion forces is incorporated (non-local correlation term). Such methods are known as vdW functionals (vdW-DFT). Several flavors of this implementation are available including the vdW-DF methods [46,47] and the so-called opt methods (optPBE-vdW [48], optB86b-vdW [49], among others).
In this work, we showed that the different methodologies of van der Waals corrections affect the energetics and electronic properties of graphene on hBN as a 2D van der Waals heterostructure. It is found that the TS correction gave reasonable adhesion energy between graphene and hBN, consistent with previous theoretical and experimental works, while other corrections tend to give values that are effectively non-binding similar to conventional PBE-DFT. In terms of electronic structure, all van der Waals corrections resulted to a band gap opening in graphene, which is experimentally confirmed. However, with the exception of TS, all van der Waals corrections made the Dirac cone to shift up from the Fermi level, making graphene artificially p-doped. As such, the TS correction is recommended to be an appropriate van der Waals correction for graphene-hBN 2D heterostructure, and that it could provide accurate predictions comparable to more expensive methods such as GW.
A 2×2 graphene-hBN 2D heterostructure was constructed with a k-point mesh of 8×8×1. A vacuum space of around 13 Å is employed to prevent inter-slab interactions. The experimental lattice constants 2.46 Å and 2.504 Å were used for graphene and hBN respectively, and as such, there is a minimal lattice mismatch at 1.7%. Regarding graphene and hBN stacking, the AA stacking was utilized, wherein carbon atoms sit directly on either boron or nitrogen atom in hBN. Although previous works point to the AB stacking (one carbon atom at the center of the hexagonal hollow in hBN) as the energetically most stable stacking of graphene-hBN [20,[55][56][57], an experimental work remarked that in most fabrication/deposition cases, this stacking usually becomes misaligned and the resulting final structure is of random orientation [58]. To simplify matters, we utilized the simplest possible structure which is the AA stacking. Besides, AA-stacked graphene-hBN is just 0.04eV higher in total energy compared to AB-stacked, and therefore, energetics wise, the structures are essentially the same as per DFT level of accuracy, and are both experimentally realizable [36]. Moreover, the electronic structure of AA-and AB-stacked graphene-hBN are basically the same [28]. Applications wise, it was found that in terms of hydrogen evolution reaction (HER) capability, the AA-stacked graphene-hBN exhibited lower Gibbs free energy of hydrogen adsorption (the lone descriptor of HER) than the AB-stacked one, implying that the former could result to better HER performance and thus is more relevant in practical applications [27].
The adhesion energy, E adhesion , of graphene and hBN is calculated using is the total energy of the graphene-hBN heterostructure, E graphene is the energy of pristine graphene and E hBN is the energy of pristine hBN monolayer. The adhesion energy per unit cell (in eV/cell) can be obtained by dividing equation (1) by the size of the simulation cell. It is noted that the total energies of the isolated monolayers correspond to that of the amount of unit cells present in the heterostructure uni cell.

Results and discussions
3.1. Structural properties and energetics of 2D graphene-hBN heterostructure for different vdW corrections We first discuss the physical structure and energetics of the graphene-hBN heterostructure by calculating lattice parameters, inter-monolayer distance and binding energy. Shown in figure 1 is the model structure for the graphene-hBN heterostructure. Tables 1, 2 and 3 show the computed nearest-neighbor distance in graphene under the graphene-hBN heterostructure, the inter-monolayer distance between graphene and hBN and the graphene-hBN adhesion energy for plain DFT, with DFT-D2, with DFT-D3, with TS and with vdW-DF2-B86R corrections respectively. As can be observed from table 1, graphene's nearest-neighbor distance is seemingly not   [60,61]. Experimentally, the same value of 3.32 Åwas reported in [62]. As such, the obtained values using vdW corrections agree with previous reports. It can be stated that in the case of graphene-hBN heterostructure, any flavor of vdW corrections can be used if the structural property under concern is the intermonolayer distance, with DFT-D2 providing the closest value to previous calculations and experimental results. Moving on to the binding energy between graphene and hBN (table 3), it is noted that plain DFT resulted to no adhesion between the two monolayers, while very minimal adhesion energies were obtained for DFT-D2, DFT-D3 and vdW-DF2-B86R corrections. To note, for pure PBE calculation, the obtained inter-monolayer distance between graphene and hBN is actually the minimum distance wherein there is zero interaction between the two monolayers. For shorter distances, the cohesive energy blows up to more positive values indicating instability. These values are comparable with the results of [26], reporting a minimal 6meV/cell adhesion energy using GGA exchange-correlation functional. To note, it is plausible that the obtained minimal binding energy using DFT-D2 might be due to the fixed C 6 and R 6 coefficients. On the other hand, the TS correction resulted to a significant adhesion energy of 0.15 eV/cell, implying a physisorption-type of adhesion between graphene and hBN. In comparison, a previous work [63] reported adhesion energies in the order of 0.1 eV/cell for graphene and hBN with vdW-DF and vdW-DF2 corrections. As experiments showed that graphene-hBN heterostructure can be materialized, these results suggest that the TS correction is one of those vdW corrections that could properly describe the binding between graphene and hBN. Moreover, these results shed light on the role of vdW corrections and the limitation of the conventional GGA exchange-correlation functional in describing the energetics of graphene-hBN heterostructure.

Electronic structure of graphene-hBN heterostructure for different vdW corrections
Experimentally, it was found that a band gap opens in graphene grown on hBN, with reported values ranging from 15-40 meV [22,23,64,65] obtained via transport measurements and magneto-optical spectroscopy. On the other hand, measurements using angle-resolved photoemission spectroscopy (ARPES) obtained band gap as high as 160meV, which was inferred to be due to very small inter-layer distance giving rise to a stronger shortrange interlayer interaction other than van der Waals [66]. As such, experimental band gap values depend on the quality of fabricated heterostucture, and the inter-layer distance. Several theoretical works indicated that the band gap could depend on the stacking alignment and applied electric field, with values ranging from 0 to 130meV [19,59]. One calculation utilizing the vdW-corrected B97-D functional estimated a tremendous gap of 325meV for AA-stacked graphene-hBN. Another study utilized the expensive GW corrections and obtained a gap of 145meV [36,58]. In summary, it can be stated that the value of graphene's band gap on hBN depends on may factors, and that no standard value could be put in place. Furthermore, as observed experimentally, graphene's Dirac cone stays at the Fermi level without any repositioning for grahene-hBN heterostructure [67].  Analyzing further the obtained band gaps, at is inferred that dispersion corrections could replicate the experimentally observed band gap of graphene on hBN. Such methods are far less computationally expensive than GW for instance, but could predict graphene's band gap more accurately. Like most theoretical methods, it   is quite hard to choose which is the best vdW corrections, as different implementations result to a different predictions. Nevertheless, from our results, it can be stated that the TS correction could be one of the most appropriate vdW correction scheme for graphene on hBN. Aside from the predicted gap of around 73meV, which is halfway between the values obtained experimentally via different methods, the Dirac point did not shift away from the Fermi level, something that other vdW corrections failed to achieve. It is suggested that pure abinitio based vdW corrections such as vdW-DF2-B86R might have overestimated the position of the Dirac cone. For TS, the combination of ab-initio based vdW coefficients and semi-empirical implementation of the vdW energy correction could have a balanced effect that retained the position of the Dirac cone in graphene.

Conclusions
In this work, the structural and electronic properties of graphene-hBN heterostructure were assessed using vdW-corrected density functional theory calculations. Four different flavors of vdW corrections were implemented: DFT-D2, DFT-D3, TS and vdW-DF2-B86R, along with the conventional PBE-DFT method. It is found that in terms of structural properties, all implementations of vdW corrections gave values consistent with previous calculations, especially the inter-monolayer distance. The calculated adhesion energy between graphene and hBN with TS correction was found to be consistent with previous theoretical and experimental results, making it one of the most appropriate vdW corrections in describing graphene-hBN binding. With regard to electronic structure, all van der Waals correction schemes resulted to an opening of band gap in graphene, comparable to previous theoretical works and experiments. Interestingly, all corrections with the exception of TS shifted the Dirac cone up the Fermi level by around 100meV on the average, inducing an artificial p-doping to graphene. It is inferred that TS is one of the most appropriate van der Waals corrections in dealing with graphene-hBN 2D heterostructure, especially in predicting graphene's electronic structure. This work demonstrated the variations in prediction capability of different van der Waals corrections with respect to  graphene-hBN 2D heterostructure. Moreover, this work also showed that certain van der Waals corrections, such as TS, could be as accurate as the more expensive methods such as GW.