DFT simulations of inter-graphene-layer coupling with rotationally misaligned hBN tunnel barriers in graphene/hBN/graphene tunnel FETs

The atomically thin two-dimensional (2D) materials translate to excellent electrostatic gate control at nanoscale channel length dimensions,¹ near-ideal two dimensional carrier behavior,^2,3 for conventional⁴ and novel device applications. Examples of the latter include field effect transistors with gate controlled resonant interlayer tunneling,^5–8 which we here refer to as a class as interlayer tunnel FETs (ITFETs), and the proposed Bilayer-pseudospin FET (BiSFET)^9,10 based on the interlayer exciton condensation. Experimental advances in stacking 2D atomic crystals have allowed the fabrication of van der Waal's heterostructures with precisely chosen sequence of atomic layers.^11,12 Stacking 2D materials in vertical van der Waal's heterostructures allows researchers to probe interesting physical phenomena such as interlayer tunneling between channel layers through a tunnel barrier, as demonstrated in a single layer (SL) graphene/hexagonal boron nitride (hBN)/SL-graphene heterostructure.¹³ Aligning the band structures of the opposing graphene layers in momentum space then allows for resonant interlayer tunneling and associated negative differential resistance,^7,8,14 which is the foundation for ITFETs.

Aligning the band structures in momentum space, of course, requires carefully aligning the crystallographic orientation of the two graphene layers. Achieving such alignment experimentally requires substantial care.^6–8,15 During layer-by-layer exfoliation of these van der Waal's materials, rotational misalignment is the norm, and even in such studies with aligned channel layers, the tunnel barrier has not been rotationally aligned. However, misalignment may affect device properties substantially. Prior theoretical studies have calculated the electronic structure of rotationally misaligned graphene bilayers using a continuum approximation¹⁶ and a tight-binding method.¹⁷ The opening of a band gap in the electronic structure of rotationally misaligned bilayer graphene has been demonstrated due to the formation of interlayer-bonded diamond nanodomains.¹⁸ Variations in the electronic structure of graphene deposited on the rotationally misaligned hBN have been even considered in a study of hydrogen desorption.¹⁹ Density functional theory (DFT) calculations have been performed on short-period crystalline structures to derive effective Hamiltonians to describe the influence of moiré pattern superlattices on the electronic properties of graphene/graphene and graphene/hBN heterostructures.²⁰ However, prior theoretical studies of resonant interlayer tunneling have not considered the effect of rotational misalignment of tunnel barrier layers with respect to the channel layers on the tunneling current.^5,14,21

In this work, we consider the effects on the interlayer current flow that rotational misalignment of the tunnel barrier with respect to rotationally aligned channel layers may have. Consider for the sake of illustration that the resonant tunneling currents I of interest here could be approximated in the form I ≈ ηe^−κt at resonance, where η and κ are constants and t is the tunnel barrier thickness. Misalignment between the graphene layers and the hBN barrier structures, in principle, could alter the evanescence coefficient κ. However, extrapolation of experimental results⁷ consistent with the apparent value of κ (corresponding to about an order of magnitude reduction in current per hBN layer) in a bilayer graphene/multilayer-hBN/bilayer graphene system to a few and even zero layers still leads to a prediction of weak interlayer tunneling even with the graphene bilayers aligned to each other (confirmed by resonant interlayer tunneling) but not to the hBN tunnel barrier. This limiting result cannot be explained through changes in evanescence alone. Therefore, we must also consider the changes in the lead coefficient η due to rotational misalignment, such as would be due to weakened coupling across the rotationally misaligned graphene-hBN interface. To consider these effects, we employ the ab-initio atomistic DFT approach. Direct calculation of current flow in this way, however, is impractical with actual ITFET device and current-flow geometries. However, DFT calculations allow us to track the effects of the hBN on the (nominal) graphene band structure as a measure of the effective coupling across the interface, as a function of rotation angle. Moreover, the inter-channel layer coupling that mediates this current flow also leads to degeneracy breaking in the band structure of the channel layers, so that knowledge of one provides information on the other. We consider SL-graphene/hBN/SL-graphene systems, although with multi-layer (ML) hBN tunnel barriers in some simulations, largely to minimize the atom count in our DFT simulations with large in-plane unit cells in rotationally misaligned systems. However, the results may be relevant to all such van der Waal's heterostructures. We show that the rotational misalignment between the graphene and hBN layers can significantly degrade the interlayer-graphene-layer coupling, and that this reduction is consistent with the reduction in the coupling across the graphene-hBN interface. In contrast, we find that there is little apparent effect on the evanescence for the coupling within the hBN barrier. That is, our results are consistent with a reduction in the lead coefficient η and not apparently with increases in the evanescence coefficient κ in the above approximation for current. This interface-localized effect would be qualitatively consistent with the weak tunneling in the few hBN layer limit extrapolated from the experimental results as noted above.

To translate knowledge of degeneracy breaking to information on current flow, consider that in the weak coupling limit—e.g., short of saturating the interlayer current flow to the lead-limited value—the tunneling current I varies as

within a first-order Bardeen transfer Hamiltonian approximation.^22,23 Here, $Mμν$ represents the interlayer coupling and is obtained from the overlap of the opposite channel layer wavefunctions across the interlayer tunnel barrier. Here, $gV$ is the valley degeneracy of graphene, and the summation in Eq. (1) runs over all states μ, ν of the top and bottom graphene layers, respectively. However, we do not need numerical values for I in this work; it is merely the quadratic dependence of the interlayer coupling that is of interest. This interlayer coupling also produces a symmetric/antisymmetric splitting of the otherwise degenerate band structures in the opposite layers under flat band conditions, which is linearly proportional to $Mμν$ in the weak coupling limit. Thus, by directly calculating the energy splitting due to coupling through a given tunnel barrier, we also can infer the greater effect on interlayer tunneling currents, with the relative change in current being proportional to the square of the relative change in degeneracy breaking in the weak coupling limit.

For aligned and misaligned systems, and for SL and ML hBN barriers, DFT simulations were performed to calculate both the absolute band structures and the band splittings due to inter-graphene-layer coupling in the rotationally aligned and misaligned systems. However, first, the atomic structures were relaxed using the projector-augmented wave method with a plane-wave basis set as implemented in the Vienna ab initio simulation package (VASP).^24,25 Graphene/hBN/graphene supercells were created using the Virtual NanoLab software.²⁶ The lattice constants of graphene and hBN are well matched (to within 1.6%), and the considered rotational misalignments between the hBN and the mutually aligned graphene layers are necessarily commensurate for our calculations. Therefore, lattice strain is minimal for the considered supercell structures. Commensurate rotation angles for the graphene/hBN interface are given by¹⁶ 

where n is a non-negative integer (n = 0,1,2,…). The number of atoms in the supercell is given by $3×(6n2+6n+2)$⁠. Figures 1(a) and 1(b) show the commensurately rotated graphene/SL-hBN/graphene supercells with rotation angles, 21.79° and 13.17°, respectively, corresponding to two smallest supercells (⁠$n=1,2$ from Eq. (2)) and visualized here using XCrySDen.²⁷ Hexagonal supercells with varying number of aligned hBN layers are shown in Figure 2. In the case of the aligned system, the C atoms of graphene are placed on top of the B atoms of the hBN layer, which is the energetically most stable configuration.²⁸ With approximately 200 atoms per supercell, use of more computationally intensive hybrid functionals or GW methods is not practical. The local density approximation (LDA)²⁹ was employed for the exchange-correlation (XC) functional as LDA can model materials well with the same orbital character at conduction and valence band edges. In a previous work on 2D materials on high-k dielectrics, we had used LDA, generalized gradient approximation (GGA),³⁰ and hybrid HSE06³¹ for XC functional in simulations of heterostructures with no qualitative differences in the results obtained.^32,33 Here, we also have re-checked some of the key results using the GGA functional. Van der Waal's non-local dispersive forces were modeled using the DFT-D2 scheme, wherein a semi-empirical correction is added to the conventional Kohn-Sham DFT theory.³⁴ A k-mesh grid of 7 × 7 × 3 for the sampling of the first Brillouin zone of the supercell was selected according to Monkhorst-Pack type meshes with the origin being at the Γ point for all calculations except the band structure calculation (where a fine k-space resolution along major axes was sampled). Atomistic relaxations of supercells were allowed to converge when the Hellmann-Feynman forces on the atoms were less than 0.001 eV/Å. After relaxation, inversion symmetry was enforced on the crystal structure so that the induced degeneracy splitting is only due to coupling between the graphene layers across the hBN layer(s) and not any level of incomplete relaxation. The interlayer degeneracy splitting calculated in this way is only weakly k-dependent within a few hundred meV of the Dirac point. For specificity, we calculated the splitting at ≈100 meV above the Dirac (K) point, avoiding the intersecting band structures about the Dirac point, as illustrated in Figure 3.

The interlayer energy splittings with the hBN layer number for the fully aligned graphene/hBN/graphene system calculated by the above outlined procedure were 69.1, 20.4, 2.2, and 0.49 meV for one to four layers of hBN, respectively (which is reasonably in line with the known zero hBN layer energy splitting limit of 370 meV represented by a Bernal stacked graphene bilayer³⁵). For one, two, or even more layers, such coupling essentially would short the graphene layers together in the large area systems, in contrast to the experimental results for rotationally misaligned systems.^7,8

To investigate this discrepancy between theory and experimental results, we first simulated channel/tunnel-barrier/channel heterostructures with a rotationally misaligned tunnel barrier layer between aligned channel layers. For the SL-graphene/SL-hBN/SL-graphene system with the hBN layer misaligned by commensurate rotation with mutually aligned graphene layers, such as for a resonant tunneling device, Figure 4 and Table I exhibit a substantial reduction in the energy splitting with layer misalignment. As also shown in Table I, such a reduction in the interlayer coupling due to this interface-localized effect would drop the estimated interlayer current per unit area by factors of approximately 9, 25, and 53 for the three twist angles of 21.79°, 13.17°, and 9.43°, respectively, within the 1st order perturbation limits for energy splitting and tunneling current (Eq. (1)). Thus, despite the growth in primitive unit cell size by factors of 7, 19, and 37, respectively, the expected interlayer current per primitive unit cell for these three rotation angles would remain roughly constant and actually decrease somewhat to 0.80, 0.76, and 0.70 times that for the aligned system. Although we are not in the perturbational limit for the strongest energy splittings, these results, nevertheless, suggest a strong effect of rotational misalignment on current flow.

We then considered SL-graphene/ML-hBN/SL-graphene heterostructures for graphene/hBN interfaces of fixed rotation angle $θ=21.79°$⁠, corresponding to the smallest non-rotationally aligned supercell from Eq. (2), and varied the number of hBN layers from 1 to 3. After relaxation of the supercell structures, the energy degeneracy splittings were computed and compared with the corresponding values for the aligned systems, and the effect on current in the perturbational limit was estimated, as shown in Table II. Although the reduction in energy degeneracy with increasing number of hBN layers is not entirely logarithmic for the aligned or misaligned systems (Figure 5)—perhaps not unexpectedly with potentially complex tunneling paths among atomic orbitals varying with the number of hBN layers within these atomistic systems—there is no evidence that increasing the number of layers would reduce the tunneling current more rapidly for a misaligned system. For this limited set of data, the effects of increasing tunnel barrier thickness are comparable to or, indeed, even weaker for the misaligned system. That is, to the extent that one may approximate interlayer current flow in the form I ≈ ηe^−κt for the sake of discussion, the results of these simulations, Tables I and II and Figures 4 and 5, are consistent with significant reductions in the lead coefficient η and not with reduction in the evanescence coefficient κ with the increasing unit cell size.

In Figure 6, we also show the band-structure near the K point (nominal Dirac point) for both the rotationally aligned SL-graphene/three-layer-hBN/SL-graphene system and the rotationally misaligned (⁠$θ=21.79°$⁠) SL-graphene/three-layer-hBN/SL-graphene system. For the aligned system, in addition to the degeneracy breaking due to the inter-graphene-layer coupling, a substantial effect on the band structure is observed near the K-point due to the interaction of the graphene layer with the hBN layers as predicted previously.²⁸ However, for the misaligned system, there is little apparent effect on band structure near Dirac point, again suggesting weakened coupling across the graphene-hBN interface in the rotationally misaligned case.

For graphene/SL-hBN/graphene system with misaligned hBN layer between mutually aligned graphene layers, we repeated the simulations with GGA exchange-correlation functional for comparison with the LDA results. The same relaxed crystal structure was used, but separate self-consistent field and band structure simulations were performed using the GGA. Figure 7 shows the band structure and interlayer energy splitting near the Dirac point for graphene/SL-hBN/graphene with hBN layer rotated at an angle $θ=21.79°$⁠, obtained using both the GGA and LDA functionals. As can be seen, the degree of band splitting matches closely. The interlayer energy splitting as a function of twist angle for graphene/SL-hBN/graphene calculated using LDA and GGA has been listed in Table III for comparison. The energy splitting obtained using the GGA functional is a few meV smaller than the corresponding LDA ones for all three cases.

Finally, we point out that as the rotation angle decreases, the interlayer coupling and the expected tunneling current do not tend to the value of aligned heterostructure in these calculations, but decreases instead as the unit cell size increases. In practice, for small rotation angles or, more generally, larger unit cell sizes beyond our capabilities to simulate, loss of coherence on unit cell dimensions or other causes would likely lead to saturation of the interlayer current density, which we expect is the experimentally observed condition. Our results only suggest that rotation leads to lower than otherwise expected interlayer tunneling currents, without addressing the limit of this reduction. Moreover, for sufficiently small rotation angles, the physical structure may self-align, at least locally,³⁶ leading back to a strong interlayer coupling and current flow.

In summary, our simulation results suggest that the interlayer tunneling current for graphene/hBN/graphene heterostructures will be affected strongly by the rotational alignment of the hBN interface with respect to the graphene layers. Indeed, for the three commensurate rotation angles considered in our study (21.79°, 13.17°, and 9.43°), despite the growth in primitive unit cell size, the expected interlayer current per primitive unit cell for these three rotation angles remained roughly constant and actually decreased somewhat. At some point, as the rotation angle decreases or, more generally, the primitive unit cell size increases beyond our ability to simulate, this trend obviously must fail, perhaps due to loss of coherence on unit cell dimensions and, ultimately, perhaps local formation of commensurately aligned regions separated by non-aligned regions.³⁶ However, the essential behavior of interest, reduction of current flow with interlayer rotation should remain. Moreover, we did not see a greater effect on coupling and expected current reduction with the increasing number of hBN layers in a misaligned system (21.79°) as compared to the aligned system. That is, to the extent that one may approximate interlayer current flow in the form I ≈ ηe^−κt for the sake of discussion, the results of these simulations are consistent with significant reductions in the lead coefficient η and not with reduction in the evanescence coefficient κ with the increasing unit cell size. Accordingly, the reduction in the expected current appears attributable to, specifically, the coupling across the interface between the channel layer and the tunnel barrier in this atomistic system, consistent with the weak tunneling current observed experimentally for such graphene/hBN/graphene systems.^7,8 Finally, it is perhaps reasonable to speculate that such effects of rotational misalignment also may be present in other van der Waal's heterostructures not directly simulated here.

See supplementary material for the complete electronic structures of the studied rotationally misaligned graphene/hBN/graphene heterostructures.

This work was supported, in part, through the Southwest Academy of Nanoelectronics (SWAN) by the Nanoelectronic Research Initiative (NRI), a Semiconductor Research Corporation (SRC) program sponsored by Nano Electronics Research Corporation (NERC) and National Institute of Standards and Technology (NIST). The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing high performance computing (HPC) resources that have contributed to the research results reported within this paper.