Zero-energy modes and valley asymmetry in the Hofstadter spectrum of bilayer graphene van der Waals heterostructures with hBN
We investigate the magnetic minibands of a heterostructure consisting of bilayer graphene (BLG) and hexagonal boron nitride (hBN) by numerically diagonalizing a two-band Hamiltonian that describes electrons in BLG in the presence of a moiré potential. Due to inversion-symmetry breaking characteristic for the moiré potential, the valley symmetry of the spectrum is broken, but despite this, the zero-energy Landau level in BLG survives, albeit with reduced degeneracy. In addition, we derive effective models for the low-energy features in the magnetic minibands and demonstrate the appearance of secondary Dirac points in the valence band, which we confirm by numerical simulations. Then, we analyze how single-particle gaps in the fractal energy spectrum produce a sequence of incompressible states observable under a variation of carrier density and magnetic field.
Long-period moiré patterns are characteristic of well-aligned graphene heterostructures with hBN. They offer an experimentally-viable route (1); (2); (3) to observe fractal magnetic miniband spectra generic for electrons in two-dimensional superlattices in strong magnetic fields (4); (5); (6). These fractal miniband spectra, also known as Hofstadter’s butterflies (7), reflect the fact that the translational symmetry of a lattice, suppressed by the presence of a magnetic vector potential in the electron’s Schrödinger equation, is restricted to magnetic field values corresponding to a rational fraction of magnetic flux quantum per unit cell area of the lattice (4); (5); (6). The generic features of fractal magnetic miniband spectra have been observed (1); (3); (8); (9) and modeled (10); (11); (8); (12); (13); (14); (15); (16) in monolayer graphene heterostructures with hBN, where the properties of Dirac electrons at the conduction/valence band edge allow for a straightforward interpretation of experimental observations, and have also been studied in slightly misaligned pairs of graphene flakes (17); (18). Modeling of the fractal magnetic miniband spectra of a bilayer graphene/hBN heterostructure is much more limited (13), even though one of the first observations of this phenomena was in this system (2).
In this Article, we study fractal magnetic minibands in bilayer graphene (BLG) subject to a moiré superlattice perturbation due to an hBN underlay. The low-energy Hofstadter butterfly spectra of a BLG-hBN heterostructure is dominated by bands related to the degenerate “zero-energy” Landau level (LL) states and of unperturbed BLG (19). We use a symmetry-based approach to model the influence of hBN on electrons in BLG, and we show that the spectra in BLG-hBN is a superposition of two very different miniband spectra associated with electrons in opposite valleys (Brillouin zone corners) of BLG’s band structure: the miniband spectrum in one valley is only weakly broadened by the superlattice potential and incorporates one unperturbed LL, while the miniband spectrum in the other valley is widely broadened. This valley asymmetry arises from spatial inversion symmetry breaking produced by the fact that the moiré perturbation only directly affects one of the two BLG layers. Since the zero energy LL states reside on different layers in opposite valleys, in one valley they are strongly influenced by the moiré perturbation, but in the other valley they are not.
These qualitative features agree with the results of the tight-binding model of Moon and Koshino (13). In addition, we characterize different types of moiré superlattice perturbations including the possibilities that the perturbation creates potential asymmetry between the carbon atoms in BLG or introduces a spatial modulation of the nearest-neighbor carbon-carbon hopping amplitude. Hence, we find that the Dirac point at the conduction-valence band edge can be either gapless or gapped, and that a secondary Dirac point can appear in the valence band of BLG-hBN. Our model Hamiltonian is presented in Section II and its numerical diagonalization in the presence of a magnetic field, including methodology, results and discussion, is described in Section III, where we also derive simple effective Hamiltonians to describe low-energy features in the magnetic minibands and, also, an effective Hamiltonian to describe the secondary Dirac point. In Section IV, we show how gaps in the fractal energy spectrum are manifest as observable incompressible states under a variation of carrier density and magnetic field.
Ii Moiré superlattice Hamiltonian
ii.1 Four band model
To describe the sublattice () and bottom/top () layer composition of electron states in Bernal-stacked BLG on a hBN underlay, we write down their Hamiltonian (19) as a matrix,
Here, we use the basis of Bloch functions for valley (), for valley (), and the Pauli matrices act in the space of sublattice components. The term on the diagonal takes into account the electrons’ Dirac spectrum in each layer where cms (20) is the Dirac velocity, , , , and , while eV (21) is the inter-layer hopping between and sublattices.
The term in Eq. (1) accounts for the moiré superlattice perturbation produced by the hBN substrate (10); (11), and is applied only to the layer of the bilayer which is closest to the hBN (layer ’’) (22). Within the term , the functions , with , are written using the shortest non-zero Bragg vectors of the moiré superlattice , where describes anti-clockwise rotation by angle , takes into account the relative lattice mismatch between graphene and hBN (23), is the misalignment angle between the two hexagonal lattices and Å is the lattice constant of graphene. For , so that the energy scale obtains its minimum value of eV at . The strength of the various terms included in are controlled by separate dimensionless parameters, , which describe in turn an electrostatic potential which does not distinguish between the two carbon sublattices (), a sublattice-asymmetric part of the potential (), and a spatial modulation of the nearest-neighbor carbon-carbon hopping amplitude ().
ii.2 Two band model
One of the most interesting features of electrons in BLG at strong magnetic fields is the degeneracy of two orbital LLs, with and , which appear at , the edge between the valence and conduction band. The mixing of these degenerate LLs by the moiré superlattice potential determines the main features of the lower-energy part of the magnetic miniband spectrum, shown in Fig. 1 for several different choices of moiré perturbation parameters . These low-energy electron states in BLG can be described (19) using a simplified two-band Hamiltonian, which can be obtained from Eq. (1) by a Schrieffer-Wolf projection onto the basis of Bloch states residing on and sublattices. For a BLG-hBN heterostructure, such a projection results in (22),
Hamiltonian is written in the basis of the Bloch states for the valley and for the valley using , and
Examples of the zero magnetic field -valley bandstructure of moiré perturbed BLG with are displayed in the left panels of Fig. 1. These spectra were calculated by numerical diagonalization of the Heisenberg matrix constructed from Hamiltonian (2) in a basis of unperturbed plane wave states. The corresponding dispersion in the -valley is obtained using the relation, , which follows from the time-reversal symmetry of Hamiltonians (1) and (2) for .
Each panel in Fig. 1 corresponds to a different set of moiré perturbation parameters . The choice of parameters used in the top panel corresponds to the predictions of a pair of microscopic models, one based on hopping between graphene carbon atoms and hBN atoms, the other on scattering of graphene electrons by the quadropole electric moments of nitrogen atoms. Interestingly, these two models predict the same relative values for the moiré potential parameters (10),
where depends upon the specific microscopic parameters used to describe the hBN underlay. Here we take , chosen to make the moiré perturbation strong enough that an almost gapless secondary Dirac point (sDP) is produced between the red and blue minibands in the valence band of the upper left panel of Fig. 1. Signatures of this feature were observed experimentally in Ref. (2).
The other panels in Fig. 1 illustrate the influence of each parameter taken individually and they exemplify three additional scenarios: for , the original Dirac point is gapped and there is no sDP (in the two valence bands closest to zero energy); for , the original Dirac point and the sDP are both gapless (and there is a sDP in the conduction band); for , the original Dirac point is gapless and there is no sDP.
Iii The BLG-hBN heterostructure in a magnetic field studied using the two band model
Here we describe the numerical diagonalization of Hamiltonian (2) in the presence of a magnetic field. To simplify this calculation, we take account of the hexagonal symmetry of the moiré pattern and use a non-orthogonal coordinate system (11) , where and are direct moiré lattice vectors and . In this basis the Landau gauge has the form which leads to
where , and is the electron charge. For the wave-vector space, this determines with and . Hence, we employ the basis set of magnetic oscillator functions,
where are the Hermite polynomials. For free electrons in BLG, the LL spectrum contains two degenerate states at , and pairs of conduction/valence band states at , ,
In general, the magnetic vector potential breaks the symmetry of the Hamiltonian with respect to translations of the moiré superlattice. However, for magnetic flux where and are co-prime natural numbers and is the unit cell area of the moiré superlattice, translational symmetry is restored. Because of this, we consider a unit cell of the magnetic superlattice that is times larger than the unit cell of the moiré superlattice in both directions (hence its area is times larger) (4); (11). The magnetic translational group , where describes geometrical translations and and are integers, commutes with Hamiltonian (2) and is isomorphic to the group of geometrical translation, so that its eigenstates form a plane wave basis . Bloch functions exist in the magnetic Brillouin zone which is times smaller than that of the moiré superlattice (4) and in which magnetic minibands are -fold degenerate. For states with momentum within this magnetic Brillouin zone,
where , and indexes the magnetic sub-bands, and indexes the above mentioned -fold degeneracy Without loss of generality, we set and omit it, using the notation , from now on. In order to obtain the energy dispersion, we calculate the matrix elements of the perturbation, , in this basis and diagonalize the resulting matrix (11).
iii.2 Results and discussion
The main panels of Fig. 1 show the magnetic spectrum of the BLG-hBN heterostructure for the four choices of moiré perturbation described above. For small flux, the magnetic miniband spectra can be traced to the sequence of Landau levels for moiré perturbed BLG. Near zero energy, the gap at the original Dirac point is seen in the magnetic spectra of the top two panels. The presence or absence of a secondary Dirac point (depending on the particular parameter set) is also clearly reflected in the magnetic spectra at small flux at the corresponding energy. At a higher flux, the Landau levels broaden and split, forming a fractal pattern, with the most striking features around zero energy.
For all parameter sets, the valley symmetry of the spectrum, preserved in the absence of a magnetic field, is lifted. This is because the moiré perturbation only affects the layer of BLG that is adjacent to the hBN and, thus, it breaks inversion symmetry. In conjunction with time-inversion symmetry breaking by a magnetic field, this allows the energy spectra for electrons in valleys and to be different (13). In particular, note that, in the absence of the perturbation, the distribution of the wave function among the layers is, for a given Landau level, exactly inverted in the two valleys, Eqs. (5,6). The moiré potential affects only the wave function component in the layer adjacent to hBN, hence breaking the layer symmetry and leading to valley-asymmetric spectra with a gap.
Importantly, the spectrum in the valley for which the wave function sits on the layer further from the substrate contains a zero-energy Landau level completely decoupled from the rest of the spectrum (shown as a red line in the -valley in Fig. 1). The spinor structure Eq. (5) of the two zero-energy states, , for unperturbed electrons shows that the states in valley () are localized on the bottom (top) graphene layer only (19). Since the moiré perturbation does not directly influence the top layer, it does not couple the state in the valley to any other state and, thus, the state in this valley always persists. This coexistence of a zero-energy LL and a fractal spectrum of magnetic sub-bands creates a unique opportunity to observe the interplay between electron-electron interaction and Hofstadter’s quantization (24); (9).
Each of the panels in Fig. 1 display further features of interest. For the parameter set (upper panel), the miniband spectrum displays a slightly gapped sDP found in the first valence miniband, which produces a sequence of LL including tilted zeroth LLs, at the corresponding energy and weak magnetic fields ( ). We shall discuss the sDP further in Section III.4 (a similar feature is present in the lower middle panel). For potential modulation taken in isolation, , the band structure obeys six-fold rotation symmetry in each valley (left upper middle panel), in contrast to all other band structure images which only display symmetry under three-fold rotation. For the spatial modulation of the carbon-carbon hopping amplitude, (lower middle panels), the spectrum at both zero magnetic field and finite magnetic field obeys electron-hole symmetry, and the zeroth LL and first LL of the original Dirac point are completely degenerate in both valleys (discussed in Section III.3). For the sublattice-asymmetric potential, (lower panel), the bands at zero magnetic field are symmetric under the operation which combines electron-hole-reflection and a rotation by , which is reflected in the electron-hole symmetry of the magnetic miniband widths in a magnetic field. Additionally for this perturbation, we find that the magnetic miniband structure around can form gapless linear or quadratic Dirac points, which produce sequences of LLs in the magnetic miniband spectrum, best visible around in the -valley and in the -valley (also see Section III.3).
iii.3 Effective model of low energy features in the magnetic minibands
The main features of the broadened magnetic miniband spectrum around zero energy at a simple fraction can be described by truncating the basis of Bloch functions Eq. (7) to the zeroth and first LL only, yielding effective Hamiltonians,
Here is used for the -valley and written using the basis (,). Also , , , , , , and . The factor leads to rapid broadening for which slows down for . In the symbol denotes a direct sum, where the state is decoupled from the rest of the spectrum.
Similarly, for simple fractions , with integer , the spectrum around zero energy may be described by effective Hamiltonians
where the block matrices are
Here we use the basis ,,, and is the zero matrix.
Figure 2 shows the excellent agreement between the result of diagonalizing effective Hamiltonians (8)-(9) (yellow), with the fully numerical diagonalization of Hamiltonian (2) (turquoise), for the various choices of magnetic field, valley, and superlattice parameters indicated on the figure. The computed dispersion surfaces display a wide array of possible forms, including the possibility of gapless linear or quadratic Dirac points for the superlattice perturbation . These features are found either in the -valley whenever in is even (Fig. 2C), or in the -valley whenever in is odd (Fig. 2D).
For the perturbation with two zero-energy LLs of Hamiltonian (2) remain unperturbed (lower middle panel in Fig. 1), which is reflected in the fact that all matrix elements in Hamiltonians (8)-(9) vanish. In this case, Hamiltonian (2) obeys the “electron-hole” symmetry , which implies that its matrix elements must obey . Consequently, the resulting Heisenberg matrix has at least linearly dependent rows, resulting in zero-energy eigenvalues which are -fold degenerate. As the only perturbation can be considered to be a periodic pseudo-magnetic field (10), it could alternatively be created using an Abrikosov lattice, i.e., a vortex structure generated from a type-II superconductor in a magnetic field. Such a system is also expected to possess degenerate zero energy eigenvalues (25).
iii.4 Effective models for the secondary Dirac points
Here we give an analytical description of the almost gapless sDP found at the corner superlattice Brillouin zone in the valence band in the top row of Fig. 1. To do this, we note that zone folding using Bragg vectors brings together three degenerate plane wave states , and to each of the two inequivalent corners of the moiré Brillouin zone where , , and,
with the polar angle of momentum , and the band index. Using theory, the vicinity of each moiré Brillouin zone corner can then be described using an effective Hamiltonian acting on a three component vector of smoothly varying envelope functions, written in a basis of the standing wave states,
Using basis (11) the Hamiltonian for the -valley is diagonal exactly at the Brillouin zone corner,
where and denote the real and imaginary parts. The corresponding Hamiltonian for the -valley is obtained using time reversal symmetry, =.
Figure 3 shows the LL spectra of Hamiltonian (12) using blue(red) dots for the () corner, where the signature of an almost gapless sDP is evident at for the corner. To calculate this spectra, we diagonalized the Heisenberg matrix of Hamiltonian (12) in a basis consisting of the products of magnetic oscillator functions (4) with the standing waves states (11), , where we take , and , where is sufficiently large that the resulting spectrum converges. Any eigenstates with a large weight on high-index magnetic oscillator functions, with , are discarded because they violate the approximation used to construct Hamiltonian (12). Also, to improve the comparison of its LL spectra with that of Hamiltonian (2), we calculate the moiré perturbation correction to the band-edge energies, , using a higher order of perturbation theory than is explicit in Eq. (12).
The result of numerical diagonalization of the two band Hamiltonian (2) is displayed using black dots in Fig. 3. The two spectra agree well, confirming that Eq. (12) is a good description of the sDP. Moreover, we can diagonalize Hamiltonian (12) for arbitrarily small magnetic fields, where the size of the Heisenberg matrix of magnetic Bloch functions (7) needed to diagonalize Hamiltonian (2) becomes prohibitively large.
The energy of the“zeroth” LL originating from the sDP in Fig. 3 depends on the magnetic field, due to a non-zero magnetic momentum generated from the influence of a third band (residing mainly on for Hamiltonian (12)) that mixes with the LLs of the first two bands. For the parameter set , this third band is separated by a large gap from the sDP, which permits the approximate two band description for the sDP,
where the parameters for the energy shift , gap , effective velocity , and magnetic momentums and can be fitted to the result of numerically diagonalizing Eq. (2). Such a fitting is illustrated in Fig. 3 using the grey lines, and produces the best accuracy for low indexed LLs.
Iv Magnetic miniband spectrum and compressibility of 2D electron liquid in BLG-hBN heterostructures
Fig. 4 shows the maps of incompressible states for each parameter set in Fig. 1, created by filling the corresponding magnetic miniband spectrum with electrons. A larger gap is depicted with a darker color, resulting in lines, the gradient of which corresponds to the filling factor. For the regions shown in yellow, the gap is negligibly small as a result of incomplete filling of a magnetic miniband, each of which can accommodate a density of electrons (including spin degeneracy). In contrast, the unperturbed zero-energy Landau level present in the valley accommodates electron density of , creating a large gapless region portrayed in red. Because of its large electron capacity at exactly zero energy, the presence of this zero-energy mode should be clearly visible in capacitance measurements. Electron-hole symmetry is displayed in the plots for (lower left), and (lower right), whereas it is absent in the upper two plots.
Many of the largest gaps in the magnetic spectrum plots (Fig. 1) occur between the zeroth and first LLs of Dirac points (whether original or secondary). In the map of incompressible states, these large gaps are manifest as a solid blue line, which intersect the axis at a density , where for the main Dirac point or for an sDP at the edge of the first conduction/valence miniband. For example, there is a (gapped) sDP in the valence band of the magnetic spectrum for the moiré perturbation (upper panel in Fig. 1). In the upper panel of Fig. 4, the gap between the two zeroth LLs of this sDP is represented as a blue vertical line, which starts from (, ), whereas the two gaps between the zeroth and first LLs are shown as two tilted blue lines, with gradients of .
In summary, we have shown that the presence of the hBN substrate lifts the valley degeneracy of bilayer graphene, producing different magnetic Hofstadter’s butterflies in each valley. The zero-energy Landau level located on the layer furthest from the hBN substrate remains unaffected by the moiré perturbation, which makes the BLG/hBN spectrum unique in comparison to other known magnetic spectra, for which all Landau levels split into sub-bands. We investigated the influence of different possible characteristics of moiré perturbation including an electrostatic potential which does not distinguish between the two carbon sublattices, a sublattice-asymmetric part of the potential, and a spatial modulation of the nearest-neighbor carbon-carbon hopping amplitude, and we identified how they give rise to different features in the electronic spectra including gapped or overlapping bands, or bands connected by a secondary Dirac point. In addition to determining the fractal electronic spectra by numerical diagonalization of a model Hamiltonian, we also derived simple effective Hamiltonians to describe low-energy features in the magnetic minibands and to describe the secondary Dirac point. Finally, we showed how gaps in the fractal energy spectrum lead to the formation of incompressible states that may be observed under a variation of carrier density and magnetic field.
Acknowledgements.We thank A. K. Geim, I. Aleiner and M. Koshino for useful discussions. This work was funded by the EU Flagship Project, EPSRC Science and Innovation Award, EPSRC Grant EP/L013010/1, and ERC Synergy Grant Hetero2D.
- L. A. Ponomarenko, R. V. Gorbachev, G. L. Yu, D. C. Elias, R. Jalil, A. A. Patel, A. Mishchenko, A. S. Mayorov, C. R. Woods, J. R. Wallbank, M. Mucha-Kruczyński, B. A. Piot, M. Potemski, I. V. Grigorieva, K. S. Novoselov, F. Guinea, V. Falâko and A. K. Geim, Nature 497, 594 (2013).
- C. R. Dean, L. Wang, P. Maher, C. Forsythe, F. Ghahari, Y. Gao, J. Katoch, M. Ishigami, P. Moon, M. Koshino, T. Taniguchi, K. Watanabe, K. L. Shepard, J. Hone, P. Kim, Nature 497, 598 (2013).
- B. Hunt, J. D. Sanchez-Yamagishi, A. F. Young, M. Yankowitz, B. J. LeRoy, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero, R. C. Ashoori, Science 340, 1413 (2013).
- E. Brown, Phys. Rev. 133, A1038 (1964).
- J. Zak, Phys. Rev. 134, A1602 134, A1607 (1964).
- J. Zak, Phys. Rev. 134, A1607 134, A1607 (1964).
- D. R. Hofstadter, Phys. Rev. B 14, 2239 (1976).
- G. L. Yu, R. V. Gorbachev, J. S. Tu, A. V. Kretinin, Y. Cao, R. Jalil, F. Withers, L. A. Ponomarenko, B. A. Piot, M. Potemski, D. C. Elias, X. Chen, K. Watanabe, T. Taniguchi, I. V. Grigorieva, K. S. Novoselov, V. I. Fal’ko, A. K. Geim and A. Mishchenko, Nat. Phys 10, 525 (2014)
- L. Wang, Y. Gao, B. Wen, Z. Han, T. Taniguchi, K. Watanabe, M. Koshino, J. Hone, C. R. Dean, Science 350, 1231 (2015).
- J. R. Wallbank, A. A. Patel, M. Mucha-Kruczyński, A. K. Geim, and V. I. Fal’ko, Phys. Rev. B 87, 245408 (2013).
- X. Chen, J. R. Wallbank, A. A. Patel, M. Mucha-Kruczyński, E. McCann, and V. I. Fal’ko, Phys. Rev. B 89, 075401 (2014).
- M. Diez, J. P. Dahlhaus, M. Wimmer, and C. W. J. Beenakker, Phys. Rev. Lett. 112, 196602 (2014).
- P. Moon and M. Koshino, Phys. Rev. B 90, 155406 (2014).
- L. A. Chizhova, F. Libisch, and J. Burgdörfer, Phys. Rev. B 90, 165404 (2014).
- J. R. Wallbank, M. Mucha-Kruczyński, X. Chen, and V. I. Falâko, Ann. Phys. 527, 359 (2015).
- G. J. Slotman, M. M. van Wijk, Pei-Liang Zhao, A. Fasolino, M. I. Katsnelson, and Shengjun Yuan Phys. Rev. Lett. 115, 186801 (2015).
- R. Bistritzer and A. H. MacDonald, Phys. Rev. B, 84, 035440 (2011).
- J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. 99, 256802 (2007).
- E. McCann and V. I. Fal’ko, Phys. Rev. Lett. 96, 086805 (2006).
- R. S. Deacon, K.-C. Chuang, R. J. Nicholas, K. S. Novoselov, and A. K. Geim, Phys. Rev. B 76, 081406 (2007).
- A. B. Kuzmenko, I. Crassee, D. van der Marel, P. Blake, and K. S. Novoselov, Phys. Rev. B 80, 165406 (2009).
- M. Mucha-Kruczyński, J. R. Wallbank and V. I. Fal’ko, Phys. Rev. B 88, 205418 (2013).
- J. Xue, J. Sanchez-Yamagishi, D. Bulmash, P. Jacquod, A. Deshpande, K. Watanabe, T. Taniguchi, P. Jarillo-Herrero, and B. J. LeRoy, Nat. Mater. 10, 282 (2011).
- A. Ghazaryan and T. Chakraborty, Phys. Rev. B 92, 235404 (2015).
- X. Chen and V. I. Fal’ko, Phys. Rev. B 93, 035427 (2016).