An accurate description of the structural and electronic properties of twisted bilayer graphene-boron nitride heterostructures

We first consider TBG on hBN where the geometrical information, including twist angle and stacking, are described in Fig. 1. The displacements of the atomic sites, the magnitude of the pseudo-magnetic field, the local deformation potentials and the effect of the hBN on the band structure are shown in Fig. 2 for $\theta_{tbg}=1.05^{\circ}$ (see Supplemental Material  [40] Sec. C) for additional results). Some interesting effects are in order here: for a twist angle $\theta_{bot}=0.53^{\circ}$, at the corners of the unit cell the atomic displacements in $G_{bot}$ do not have a circular “vortex” shape as in TBG (see Supplemental Material  [40] Sec. B) [42, 43, 44], instead, they have a triangular-like shape. As shown in Fig. 2(a), the displacements have a counterclockwise behaviour at the corners (AAA stacking configuration in Fig. 1) and clockwise displacement in the inner regions (ABC and ACB). The magnitude of the in-plane twist of the atoms with respect to their original position, $\Delta\theta$, has different values along the moiré because the corresponding atomic binding energies are different. This is due to the different stackings found in each region [31, 33]. Figure 2(e) displays the distribution of the deformation. As we can see, in $G_{bot}$ its maximum value is about $60$ meV at the ABC region and a minimum of $-40$ meV is found at the AAA region. Notice that in the ACB region the on-site energies are small, as in the AAA centers. This behaviour is reminiscent of that of monolayer graphene on hBN [31, 45]. In addition, these magnitudes are quite close to those of graphene on a hBN substrate [24]. In $G_{top}$ the on-site potential is one order of magnitude smaller. Figure 2(f) shows the vector potential $\mathbf{A}$ (red arrows) and the pseudomagnetic field (in color) induced in both $G_{bot}$ and $G_{top}$. The pseudomagnetic field in $G_{bot}$ has a complex shape similar to a “fidget spinner”, and is completely different to that of monolayer graphene on hBN [33] or free-standing TBG [37]. The resulting non-uniform pseudo-magnetic field in $G_{bot}$ has a maximum value of about $18$ T (which is almost twice as in pristine TBG [37]). In $G_{top}$, the pseudomagnetic field is smaller with a maximum of $9$ T, as shown in Fig. 2(a).

In recent studies, the effect that a hBN substrate has on TBG is introduced by means of an effective periodic potential acting only on the nearest graphene layer [46, 22, 47, 48] with parameters obtained from calculations of a graphene monolayer on hBN. As shown in Fig. 2 our results indicate that this approach is correct, however, a renormalized set of substrate parameters should be considered while describing the periodic potentials. We will discuss this in the following sections.

By increasing $\theta_{bot}$ the effect of the hBN on the TBG band structure is reduced. Our results indicate that, for the considered stacking, the effect of the substrate survives for angles larger than $3^{\circ}$. Nevertheless, it should be pointed out that small variations in the value of $\theta_{bot}$ for angles below that threshold, have a huge impact in the electronic properties of the system. As $\theta_{bot}$ increases, shown in Fig. 2(a), the period of the moiré length between $G_{bot}$ and $hBN_{bot}$ is reduced [49]. For example, as detailed in Supplemental Material [40] Sec. A, for $\theta_{bot}=1.59^{\circ}$ we have $q=3$ resulting in $L_{TBG}=3L_{hBN}$. The effect on the band structure is shown in Fig. 2(b)-(d), where the narrow bands are separated by a gap resulting from the breaking of inversion symmetry ($\mathcal{C}_{2}$). In the nearly aligned situation, $\theta_{bot}=0.53^{\circ}$, the gap between the narrow bands is $\sim 25-30$ meV. As the twist angle of the substrate increases, the second moiré length is reduced and the effects of the periodic potentials are suppressed. For $\theta_{bot}=2.65^{\circ}$ the gap is nonzero (see Fig. 2(c)). The persistence of the gap for angles far from alignment is due to the large value of the mass term resulting from relaxation effects and by the large difference of on-site energies between graphene and hBN sites [50, 51]. In addition, as shown in Fig. 2(b) and Fig. 2(c), the density of states (DOS) for two different angles is quite similar. We expect to find similar DOS for smaller values of $\theta_{bot}$. The persistence of the band gap for a window of small angles and the isolated narrow bands may explain the presence of quantum anomalous Hall effects (QAH) and other non-trivial band topology effects found in TBG with a nearly aligned hBN substrate [16, 15, 17]. Our results are consistent with previous studies which argue that the QAH would only appear if one or both $\theta_{bot}$ and $\theta_{top}$ are near alignment [22].

On the other hand, in the trilayer system shown in Fig. 1, we can tune both the graphene twist angle $\theta_{tbg}$ and the substrate $\theta_{bot}$. As mentioned before, we found that, by increasing $\theta_{bot}$, the effect of the substrate on the TBG electronic properties is reduced. By modifying the graphene twist angle $\theta_{tbg}$ while keeping fixed $\theta_{bot}$, the band structure is also modified. Figure 3 shows the band structure for $\theta_{bot}\approx 0.53^{\circ}$ and different twist angles, $\theta_{tbg}=1.05,1.12,1.21\text{ and }2.28^{\circ}$. Here, it is important to note that, in our model, the “magic” angle is around $\theta_{tbg}\sim 1.21^{\circ}$. As it can be seen, at this particular angle, the substrate induces a finite gap of about $\approx 30$ meV. If we decrease the TBG angle, the bands become wider, especially the valence bands. Furthermore, the gap between the valence and conduction bands gap becomes slightly larger. Contrastingly, by increasing the angle, we can see that for $\theta_{tbg}=2.28^{\circ}$, the AA bilayer graphene band structure is recovered [52, 53]. In this particular case, the presence of hBN has the effect of opening a small gap at the Dirac cones.

In free standing TBG, the layer degree of freedom is disentangled from spin and valley, providing eight-fold degeneracy in the low energy states. In the unperturbed system, the Dirac points are at the same energy (dashed black lines in Fig. 3), and at low energies, Dirac cones only interact with their analogous at opposite layers. The presence of the hBN substrate induces a mass gap. However, the hBN has a larger effect on $G_{bot}$ than on $G_{top}$ and the gapped Dirac cones of $G_{top}$ and $G_{bot}$ are no longer at the same energy. Since the $K$ point corresponding to layer $1$ is mapped to a $K^{\prime}$ point of layer $2$ in the opposite valley, a “splitting” in each of the conduction and valence bands is obtained in the full TB model due to the energy difference. Therefore, this splitting is a breaking of the layer degeneracy because the hBN has a larger effect on $G_{bot}$. In fact, as we will show in the following section, by adding a second hBN layer on the top and for certain twist angles, the degeneracy is recovered. A similar phenomena is obtained in TBG in a uniform electric field [54, 55]. Next, to identify the bands corresponding to each valley we define a valley operator of the form

where $i,j$ denotes next-nearest-neighbor sites, $\eta_{ij}=\pm 1$ for clockwise or counterclockwise hopping, respectively and $\sigma_{z}^{ij}$ the Pauli matrix associated with the sublattice degree of freedom. The expectation value of this operator ranges from $-1$ to $+1$, which corresponds to the $K$ and $K^{\prime}$ valley, respectively. Figure 2(b) displays the band structure where the band corresponding to each valley is identified by different colors.

We now consider the case of TBG encapsulated between two layers of hBN. Depending on the orientation and twist angle of each hBN layer with respect to TBG, several configurations can be obtained. Similar to the trilayer case, our starting point is an AAAA stacking configuration at the unit cell origin. Additional configurations can be obtained by modifying the stacking, see for example Refs. [56, 57]. In Fig. 4(a) we show a structure with arbitrary $\theta_{bot}$, $\theta_{top}$ and $\theta_{tbg}$, where the graphene layers are depicted in red and blue and the hBN substrate in green and black. Furthermore, we also show the different stacking configurations that can be found at the high symmetry points along the moiré superstructure.

Due to the three different degrees of freedom to build this system, the number of different heterostructures that can be generated are almost limitless. Here we are going to focus on three different cases where we fix $\theta_{tbg}$ to $1.05^{\circ}$ and modify both $\theta_{bot}$ and $\theta_{top}$. In Figs. 4(b)-(d) we show the corresponding band structure, strain fields, deformation potential and pseudomagnetic fields for the different studied configurations. Particularly, for these structures, the high-symmetry stackings are (AAAA, ABCA, ACBA), (AAAA, ABCB, ACBC) and (AAAA, ACAC, ABAB), respectively.

In Fig. 4(b) we show the case when $\theta_{bot}=\theta_{top}$. It is important to note that, since all the twist angles are measured with respect to $G_{bot}$, the twist angle with respect to $G_{top}$ is $\theta^{\prime}_{top}=-0.53^{\circ}$. Therefore, the lattice structure is such that $hBN_{top}$ and $hBN_{bot}$ are aligned with respect to each other. In this case, a large gap of about $\sim 50$ meV is found in the band structure. The corresponding pseudomagnetic field has a three-lobed structure with a maximum magnitude of about $36$ T, which is larger than in the trilayer case. The magnitude and direction of this field is opposite in each graphene layer and, as shown in the right panel of Figs. 4(b)-(d) is highly sensitive to the hBN twist angle. In the second column of Fig. 4, we can see how the in-plane twist of the atoms with respect to their original position, $\Delta\theta$, has a complex shape which also strongly depends on the hBN angle. Interestingly, in $G_{top}$ the in-plane strain field is oriented in the opposite direction with respect to the field in $G_{bot}$. The deformation potential is proportional to the change in the effective area around each atom with respect to its equilibrium position, Eq. (8). As shown in the third column of Fig. 4, the potential at the corners of the unit cell is negative and has a smaller absolute value than those in the interior regions, where there is a maximum. This behaviour is reminiscent of monolayer graphene on hBN where is know that the moiré pattern in the rigid system smoothly changes between $AA$ type, $AB$ type (carbon on boron) and $BA$ type (carbon on nitrogen). Each of these configurations have a different adhesion energy which is minimum at the AB stacking, while BA and AA are roughly similar [31, 33, 34]. The different energies create in-plane forces which tends to maximize the area of the AB regions and minimizing the other regions. This behaviour is completely capture by our results in the third column in Fig. 4 where in the red regions the local area around each atomic site was maximised and by Eq. (8) the on-site potential is positive. Interestingly, this color difference allow us to distinguish the different stack configurations between each TBG layer and its neighboring hBN layer by simply looking for the maximum or minimum in the corresponding deformation potential.

Interestingly, in the tetralayer structure hBN/tBG/hBN, for twist angles where $\theta_{bot}=\theta_{top}$, Fig. 4(b), the eightfold degeneracy is recovered. As we will elucidate in the following section, this occurs because the periodic potentials acting in each graphene layer are opposite or with the same magnitude. On the contrary, if $\theta_{bot}\neq\theta_{top}$, as shown in Fig. 4(c), the effect of the substrate is different in each graphene layer and the layer degeneracy is broken resulting in a band splitting. Notice that, the pseudomagnetic fields have different maximum values on each layer. The same happens for the in-plane strain. In principle, we can assume that the layer degeneracy is broken in structures with different $\theta_{bot}$ and $\theta_{top}$, however, this is not always the case. Fig. 4(d) shows the TBG band structure with different twist angles $\theta_{bot}$ and $\theta_{top}$ where the layer degeneracy is recovered. By simple geometry we notice that the twist angle of $hBN_{top}$ with respect to $G_{top}$ is $\theta_{top}=0.54^{\circ}$ which means that each hBN layer is twisted in opposite directions with the same angle with respect to its nearest graphene layer and, therefore, the corresponding pseudomagnetic and strain fields, Fig. 4(d), have opposite directions and equal magnitudes. Interestingly, the density of states close to charge neutrality has a single peak. This indicates that even if $\theta_{bot}$ and $\theta_{top}$ are small (or both hBN layers are nearly aligned with its corresponding graphene layer), the band structure can be modified completely. This can have important effects in experimental measurements. The strong dependence of the electronic properties of TBG for small values of $\theta_{bot/top}$ may explain why in some experiments the valley anomalous Hall conductivity in suspended or encapsulated structures of TBG with hBN is not always present [16, 15, 58]. In addition, our results indicate that adding a single hBN layer to TBG breaks the layer degeneracy giving rise to a splitting in the band structure. This effect will produce a double peak in the local density of states similar to the effect produced by a magnetic field [37].

The effect of the substrate can be introduced as an effective potential acting on the closest graphene layer which is periodic in the moiré unit cell [35, 33]:

with amplitudes $V_{\text{SL}}(\bm{G_{j}})$ given by

where

with $M_{j}=(-i\sigma_{2}G_{j}^{x}+i\sigma_{1}G_{j}^{y})/|G_{j}|$. The $2\times 2$ Pauli matrices act on the sublattice index in a single graphene layer. $\Delta$ is a parameter that represents a spatially uniform mass term [68]. Surprisingly, having identified that some of the potential profiles in the TB solutions indicate the existence of additional harmonics (see Supplemental Material  [40] Sec. G), we found that the modulation of $V_{\text{SL}}$ at smaller wavelengths has some effects on the narrow bands and, therefore, Eq. (22) has to be expanded up to two harmonics. As described in Ref. [35, 33], the minimal model for a hBN substrate depends on eight parameters: uniform on-site term $\omega_{0}$, constant mass gap $\Delta$, position-dependent scalar terms $V_{s}^{e}$ and $V_{s}^{o}$ which are even and odd under spatial inversion, respectively. Two position dependent mass terms, $V_{\Delta}^{o}(e)$ and two gauge terms $V_{g}^{o}(e)$, odd and even respectively. In the case of the second harmonic, six additional periodic potential parameters are required. Therefore, the Hamiltonian of TBG/hBN depends on several parameters, four for TBG (Eq. (12) and Table 1) and eight for the substrate, Eq. (22). The exact values and their dependence with the twist angles are known for single layer graphene on hBN [30] but they are still unknown for the combined system TBG/hBN or hBN/TBG/hBN. Obtaining exact values for a large set of parameters is a difficult task since several combinations can give similar results and a complete fitting is out of the scope of this work. In addition, it is important to note that, the induced periodic potentials and mass gap due to the presence of a hBN substrate dominate over $V_{0}(r)$ of the free standing system (Eq. (21)), which, in this particular case, is in fact negligible.

Some qualitative estimations of the substrate parameters can be given by analyzing the band structures and the periodic potentials in Fig. 2. For the considered stacking in Fig. 1, an estimated set of parameters is given by $(w_{0},\Delta,V_{1,s}^{e},V_{1,s}^{o},V_{1,\Delta}^{e},V_{1,\Delta}^{o},V_{1,g}^{e},V_{1,g}^{o})=(3.0,31.62,-0.75,0.68,-0.02,3.4,-5.14,18.6)$ meV, for the first harmonic amplitudes and the only non-zero amplitudes for the second harmonic are the gauge fields $(V_{2,g}^{e},V_{2,g}^{o})=(-5.14,-9.3)$. Fig. 5(b) shows the band structure with the full TB model (black line) and the fitting with a continuum model (red line) using our estimated set of parameters which agree qualitatively with the TB calculations (see Supplemental Material  [40] Sec. H). In particular, in order to achieve this kind of agreement, we found that: i) the mass gap is large and ii) the gauge terms require up to two harmonics. The first feature is in agreement with Ref. [56], where the band structure with a similar stacking has a large mass gap. As shown in Fig. 2, as $\theta_{bot}$ increases, the corresponding moiré length is reduced. This strongly suppresses the effect of the periodic potentials induced by the hBN and only the constant mass terms survives for large angles. As a final note, the strong gauge field in ii) results from the effects of the relaxation of the trilayer system (TBG/hBN).

By fitting the continuum model to the TB band structures we have found that the dominant terms are the gauge fields which give rise to a pseudomagnetic field and the mass term which is responsible for the large gap. In the following, for simplicity, we only consider the effects of the mass term. A similar analysis can be performed with the gauge potential. Fig. 6 shows the band structure of TBG in both valleys with only the constant mass term $\Delta$ of Eq. (22). We define $\Delta_{b}$ and $\Delta_{t}$ as the mass terms acting in $G_{bot}$ and $G_{top}$ respectively. Figs. 6(a) and (d) clearly show the band splitting due to the different mass acting on each graphene layer. On the contrary, Figs. 6(b) and (c) are still degenerated due to the fact that the mass terms for the bottom and top layers are the same. This explains why in some of our TB results the presence of two hBN layers does not split the narrow bands. Our results indicate that in an experimental setup, if a hBN layer is nearly aligned with its closest graphene layer, the layer degeneracy is broken, resulting in a double peak in the local density of states, similar to the effect produced by a magnetic field [37]. We would like to recall that, in the encapsulated case, depending on the orientation of the hBN layers the degeneracy can be recovered.