Large-area, periodic, and tunable intrinsic pseudo-magnetic fields in low-angle twisted bilayer graphene

TL;DR

This study provides the first experimental evidence of large-area, tunable pseudo-magnetic fields in low-angle twisted bilayer graphene, revealing vortex lattices and pseudo-Landau levels, and demonstrates control via rotation and strain.

Contribution

It experimentally confirms the existence of tunable pseudo-magnetic fields in twisted bilayer graphene and shows how to manipulate their properties through geometric and strain modifications.

Findings

Observation of vortex lattices matching moiré patterns

Detection of pseudo-Landau levels via STM/STS

Pseudo-magnetic fields can be tuned by rotation angle and heterostrain

Abstract

A properly strained graphene monolayer or bilayer is expected to harbour periodic pseudo-magnetic fields with high symmetry, yet to date, a convincing demonstration of such pseudo-magnetic fields has been lacking, especially for bilayer graphene. Here, we report the first definitive experimental proof for the existence of large-area, periodic pseudo-magnetic fields, as manifested by vortex lattices in commensurability with the moir\'e patterns of low-angle twisted bilayer graphene. The pseudo-magnetic fields are strong enough to confine the massive Dirac electrons into circularly localized pseudo-Landau levels, as observed by scanning tunneling microscopy/spectroscopy, and also corroborated by tight-binding calculations. We further demonstrate that the geometry, amplitude, and periodicity of the pseudo-magnetic field can be fine-tuned by both the rotation angle and heterostrain applied to the system. Collectively, the present study substantially enriches twisted bilayer graphene as a powerful enabling platform for exploration of new and exotic physical phenomena, including quantum valley Hall effects and quantum anomalous Hall effects.

Equations

$$\Delta d(r)$$

$$\Delta z(r)$$

$$V$$

$$A_x$$

$$A_y$$

$$|\mathbf{B}|\approx \frac{\hslash c}{e}\left(\frac{2\beta }{3}\right)\frac{\bar{u}}{a_{m}d},$$

$$\bar{u}=\frac{\sqrt{3(d_3-d_2{)}^2+(d_2+d_3-2d_1{)}^2}}{2d}$$

$$t_{ij}=n^2{V}_{pp\sigma }(r_{ij})+(1-n^2)V_{pp\pi }(r_{ij}),$$

$$V_{pp\pi }(r_{ij})$$

$$V_{pp\sigma }(r_{ij})$$

Full text

Strain in graphene can be served as an effective tuning parameter in tailoring its exotic properties reflected by different degrees of freedom such as spin and valley Guinea2010 ; Vozmediano2010 ; Katsnelson2015 ; Gargiulo2018 ; Levy2010 ; Klimov2012 . In particular, it has been predicted that a non-uniform strain distribution imposed on graphene arising from either external stress Guinea2010 ; Vozmediano2010 or interfacial coupling Katsnelson2015 ; Gargiulo2018 can induce a strong gauge potential that in turn can act as a pseudo-magnetic field (PMF) on the Dirac electrons. Indeed, it has been demonstrated experimentally that locally strained graphene in the form of nanobubbles can induce PMFs, as manifested by the emergent pseudo-Landau levels observed using the scanning tunneling microscopy/spectroscopy (STM/S) Levy2010 ; Klimov2012 . More recently, existence of larger-area, spatially distributed PMFs has also been reported for single-layer graphene on a black phosphorus substrate, which provides strain textures due to the mismatched symmetries and lattice constants of the heterostructure Liu2018 . Separately, intensive theoretical and experimental studies have demonstrated that low-angle twisted bilayer graphene (TBG) can exhibit a variety of exotic phenomena, including superconductivity and quantum phase transitions Bistritzer2011 ; Laissardiere2012 ; Kim2017 ; Cao_insulator ; Cao_super ; yankowitz2019tuning ; chen2019evidence ; yoo2019atomic ; kerelsky2019maximized ; xie2019spectroscopic ; jiang2019charge ; choi2019electronic . In addition, it has been reported that a small uniaxial heterostrain in low-angle TBG can suppress Dirac cones and lead to the emergence of a zero energy resonance Huder2018 . In practice, low-angle TBG itself possesses large-scale, natural, and periodic strains due to the incommensurate moiré superstructures Guinea2012 ; Katsnelson2015 ; Gargiulo2018 . Such strain has been theoretically studied to significantly affect the electronic properties of the large-scale moiré pattern liu2019pseudo ; Aline2018 . However, yet to date, experimental explorations of such natural and intrinsic strain effects arising from lattice relaxations in the strongly coupled homostructural system, especially the pseudo-magnetic behaviors, are still lacking yan2013strain ; qiao2018twisted ; liu2019pseudo .

In this Letter, by combining STM/S experiments and large-scale tight-binding calculations, we study the electronic properties of TBG with extremely small twist angles ($\theta<1^{\circ}$). We directly detect pseudo-Landau levels generated by the intrinsic PMFs, which originate from the interplay between interlayer interactions and in-plane strain fields. The PMFs form a ring-structured vortex lattice associated with the moiré pattern, and are periodically distributed over the whole sample. The vortex lattice, in which all angular incommensurability is focused at specific singularity points, is chiral for the charge carriers at different valleys in k-space, and can be further tuned, together with the PMFs, by the twist angle and uniaxial heterostrain in the bilayer system. Indeed, the existence of periodic PMFs in such hexagon-structured system provides a possible platform to further explore quantum anomalous Hall effects if a tiny external magnetic field is applied to break the time-reversal symmetry haldane1988model .

Our samples were prepared by chemical vapour deposition on Cu-Ni alloyed substrates. Cu-Ni alloyed substrates have been widely used to grow bilayer graphene Wu2012 . High quality bilayer graphene samples with multiple domains were obtained after appropriate annealing processes. Figure 1a shows an STM topography image with hexagonal symmetry moiré patterns. The atomic structure of the moiré pattern can be determined explicitly by the topography image Hermann2012 ; Wijk2014 ; Huder2018 . Here, we follow the same method to identify the twist angle in Fig. 1a, which gives $\theta=0.48^{\circ}$ (see Supplementary Information Sec. E).

Next, we characterize the electronic properties of the TBG samples by the dI/dV tunneling spectra. Generally, dI/dV conductance is proportional to the local density of states in the vicinity of the STM tip position. To rule out the possible tip effects on the measured dI/dV conductance, we have performed dI/dV measurements on TBG with large twist angles. The obtained spectra show typical angle-dependent van Hove singularities which are highly consistent with previous reports Laissardiere2012PRL (Supplementary Information Sec. B). Figure 1b shows dI/dV spectra taken in both AA and AB stacking regions. The AA region spectra show a series of pronounced peaks (indicated by numbers) with nearly equal energy spacing. Surprisingly, these peaks are typically different from those of TBG in the absence of external magnetic fields Li2010 ; Kim2017 ; Huder2018 ; Laissardiere2012PRL , but rather similar to the Landau levels of a two dimensional electron gas in graphene under strong external magnetic fields Li2007 ; Rutter2011 . In stark contrast, these resonances are hardly found in the AB region. For extremely low-angle TBG, the charge density in the AB region resembles that of ideal Bernal-stacking bilayer graphene Gargiulo2018 .

To interpret our experimental observations, we have performed theoretical calculations of the superstructure shown in Fig. 1a. The as-shown moiré pattern unit cell contains 57964 atoms that is too large even for state-of-the-art first-principles methods. Therefore, we adopted a widely used full tight-binding model of bilayer graphene with an angle-dependent Hamiltonian Laissardiere2012 ; Huder2018 . Figure 1c shows the calculated local density of states by utilizing the Lanczos recursion method in real space recursion_method . Considering the lattice deformation of the bilayers (upper panel), the result reproduces the energy peaks that appear in the measured spectra (Fig. 1b) whereas the appearance of high energy resonances are dramatically different in the rigidly twisted one (bottom panel). The implementation of the lattice distortion in a rigidly twisted bilayer graphene will be discussed later. As a side issue, our calculations (Fig. 1c) also reproduce the observed resonance peak at the Fermi energy. Recent investigations of such system have provided experimental evidence of such strong electronic correlations near the magic angle TBG by means of STM/S kerelsky2019maximized ; xie2019spectroscopic ; jiang2019charge ; choi2019electronic .

For bilayer graphene under external magnetic fields, the massive Dirac electrons are quantized with energies $E_{N}=\pm\hbar\omega_{\mathrm{c}}\sqrt{N(N-1)}$, $\omega_{c}=eB/m^{*}$, N = 0,1,2,…, where $\hbar$ is Planck’s constant and $m^{*}$ is the effective mass of electrons yan2013strain ; Rutter2011 . Here, the resonant peak positions are linear with $\sqrt{N(N-1)}$, indicating the feature of pseudo-Landau levels in bilayer graphene. For the TBG with twist angle of $0.48^{\circ}$, the fitted PMF value is about 9 T, as shown in Fig. 2a. In addition, the PMF in low-angle TBG does not break the time-reversal symmetry and we further calculated the pseudo-Landau levels under external magnetic fields (Fig. 2b). An interesting observation is that, under the strong external magnetic fields, the pseudo-Landau levels are splitted into two peaks with an energy gap of about 15 meV due to the breaking of valley degeneracy. Similar phenomena have been experimentally verified in the strained monolayer graphene Li2015Observation .

To illustrate the delicate physics involved in reaching the agreement between the theory and experiment shown in Fig. 1, we note that, in graphene-based van der Waals structures, compelling evidences have revealed that the lattice relaxation effects are crucial in determining their electronic properties Wijk2014 ; Katsnelson2015 ; Gargiulo2018 ; yoo2019atomic . Especially, in low-angle TBG, the superstructure undergoes lattice distortions and structural transformations, which can be described as a lattice deformation of the moiré superlattice due to the interplay between the interlayer interaction and the in-plane strain fields. That is, on one hand, the in-plane forces move atoms to maximize the area of AB/BA stacking domains, which has the minimum binding energy. On the other hand, the strain induced by the atomic displacements hinders such in-plane atomic rearrangement. Katsnelson2015 ; Gargiulo2018 . By fitting the atomic structures obtained from molecular dynamics simulations of TBG in Refs. Katsnelson2015 ; Gargiulo2018 , the in-plane displacement $\Delta d$ and out-of-plane displacement $\Delta z$ of individual atoms can be approximately expressed as:

[TABLE]

where $r$ is the distance between the individual atoms and the nearest AA point in the x-y plane, $\Theta$ is the Heaviside step function, $\Delta D$ is the maximum in-plane displacement and $\Delta Z$ is the maximum out-of-plane displacement that are dependent on the rotation angle Gargiulo2018 . The deformed samples are constructed using the following procedure: First, the rigidly twisted bilayer graphene is constructed according to the commensurability (see Supplementary Information Sec. E) and the position of each atom is well known. Next, we modify the in-plane and out-of-plane positions of atoms according to the position-dependent expression in Eq. (1). For instance, for the in-plane deformation, individual atom moves $\Delta d$ towards the nearest AA point. In the out-of-plane deformation, each individual atom has a displacement in the $z$ direction with $+$ sign for the top layer and $-$ sign for the bottom layer. The maps of $|\Delta D|$ and $|\Delta Z|$ for TBG with $\theta=0.48^{\circ}$ are plotted in Fig. 3c and d, respectively. The movement of the atoms have the same tendency as that obtained from molecular dynamics simulations in reference Gargiulo2018 . In this part, we assume that the deformed sample keeps the period of the rigid TBG. In the rigid sample, we use the AA point as the rotation center. The moiré supperlattice has point group $D_{3}$ generated by a three-fold rotation about the axis formed by the AA point ($C_{3}$) and a two-fold rotation about the axis perpendicular to the AA point ($C_{2}^{\prime}$). The deformed sample maintains the $D_{3}$ point group symmetry guinea2019continuum ; nam2017lattice . For another sample with uniaxial heterostrain, the supperlattice for both rigidly twisted and deformed cases lose the $D_{3}$ symmetry. In Eq. (1), the parameters $\Delta D$, $l_{D}$, $\sigma_{D}$, $\Delta Z$, $l_{Z}$, and $\sigma_{Z}$ determine the profile of the deformed structure and are fitting parameters to generate proper structures with properties matching the experimental observations. In fact, in our samples, the substrate suppresses the lattice distortion, and the movements of the atoms are weaker than those in a free-standing twisted bilayer graphene walet2019lattice . Moreover, the in-plane deformation has a dominant effect on the electronic properties, which can be seen from the LDOS in Fig. 3e. In the deformed sample in Fig. 1c, the fitting parameters are the sets B (in-plane) and E (out-of-plane) labelled in Fig. 3.

The structure deformation leads to an electrostatic potential $V$, proportional to the local compression/dilatation, and a pseudo-vector potential $\boldsymbol{\mathit{A}}$, associated with the shear deformation at one site. They can be written as Vozmediano2010 ; Katsnelson2015 :

[TABLE]

where $g_{1}\approx 4\;\mathrm{eV}$ is the deformation potential for graphene, $d_{1}$, $d_{2}$, and $d_{3}$ are the first neighbour inter-atomic distances in the deformed lattice, $d=1.42\;{\mathrm{\AA}}$ is the equilibrium carbon-carbon distance, $c\approx 1$ is a numerical factor depending on the detailed model of chemical bonding, $\beta\approx 2$ is the electron Grüneisen parameter, $\gamma_{0}=3.2\;\mathrm{eV}$ is the strength of first-neighbor interaction in the plane, $u_{ij}$ is the deformation tensor. The PMF is given by $\boldsymbol{\mathit{B}}=\nabla\times\boldsymbol{\mathit{A}}$, which can be estimated as Katsnelson2015 :

[TABLE]

where $a_{m}$ is the period of the moiré pattern. And $\bar{u}$ is the magnitude of the shear deformation, which can be calculated from the first neighbor interatomic distances in the deformed sample as:

[TABLE]

The local deformation potential $V$ and the magnitude of the pseudo-magnetic fields $\boldsymbol{B}$ are plotted in Fig. 1d and 1e, respectively. Such structural deformation results in a non-uniform $\boldsymbol{B}$ with a maximum value of about $6$ T. One of the effects of applying a strong external magnetic field in a two-dimensional electronic structure is the Landau quantization of the eigenstates and, subsequently, the quantum Hall effects in the transport properties. In our samples, the PMFs are not uniform, but form a ring-structured vortex lattice around the centers of the AA regions, similar to an Abrikosov vortex in superconductors Katsnelson2015 . Here, these new pseudo-Landau levels are strongly correlated to the space-dependent PMFs, similar to the Landau levels that appear in the presence of a real magnetic fields. It is important to emphasize that the PMFs in low-angle TBG intrinsically arise from the lattice deformation. We expect that, as the displacements of the atoms are more significant in a structure with smaller twist angle Katsnelson2015 , the induced PMFs will also be stronger.

The pseudo-Landau level behaviour with twisted angle of $0.65^{\circ}$ has also been studied by STM/S, in which the observed pseudo magnetic fields are fitted to be 8 T (details can be found in Supplementary Information Sec. G). The calculated PMFs are in qualitative agreement with the fitted one from the $E_{N}$ formula. The fitted PMFs decrease more smoothly in the space comparing to the calculated values, and this difference can be explained in the following: i) The expression of the Landau energies $E_{N}$ is valid for the case of bilayer graphene under an uniform external magnetic field. On the contrary, in our sample, the deformation-induced PMFs are non-uniform. ii) As can be seen from the Landau energies $E_{N}$, the fitted PMFs can be modulated by the effective mass $m^{\ast}$, of which the value is taken from other experimental results yan2013strain ; Rutter2011 , and the effective mass can also be space-dependent due to the local stacking structure. All in all, in this paper, both the theoretical and experimental values of maximum PMFs increase when the twist angle decreases, which is consistent with the fact that the strength of the lattice distortion becomes stronger for a TBG with smaller rotation angle.

The results presented so far have demonstrated that the twisted angle can be an effective tuning parameter in tailoring the lattice deformation and subsequently tuning the electronic and magnetic properties in low dimensional systems. Next, we investigate heterostrain effects on the PMFs and their distributions in the low-angle TBG entity. The heterostrain is uniaxial, which originates from the pinning of the graphene layer at its grain boundaries or step edges during the epitaxial growth (Supplementary Information Sec. D) Huder2018 . Figure 4a shows the topography image with severely distorted moiré pattern. Such distorted moiré structures are due to the existence of uniaxial heterostrain, where the moiré pattern can serve as a magnifying glass Hermann2012 ; Wijk2014 ; Huder2018 ; Chengdong2018 . Namely, when the honeycomb lattice is slightly strained, the resulting moiré pattern is significantly distorted, in which the twist angle and strain can be determined by the STM images. For the sample in Fig. 4a, we obtain the twist angle and heterostrain to be $\theta=0.98^{\circ}$ and $\sigma=0.78\%$, respectively (see Supplementary Information Sec. E).

Figure 4c and 4d show dI/dV spectra and calculated LDOS curves taken along a line from an AA to an AB region. Similar to the findings in the $\theta=0.48^{\circ}$ sample, pseudo-Landau levels are also observed as a series of pronounced peaks. Significantly, the energy spacings between the pseudo-Landau levels vary from the AA to AB region, which indicate the non-uniform characteristics of the PMFs. The splitting of the degenerated $N=1$ Landau level near the Fermi level is due to the heterostrain-induced low-energy bandgap opening in TBG yan2013strain ; Rutter2011 . Figure 4b shows the corresponding linear fitting of pseudo-Landau levels, which gives PMF value of 6 T. By linear extrapolating the LDOS at high energy in Fig. 4c, we deduced the center of the states, which changes with the positions. For our sample with a considerable heterostrain of $\sigma=0.78\%$, it is still hard to visualize the atomic displacements from the high resolution STM image. Despite this, the fast Fourier transform image of the large area moiré patterns shows distorted six-fold symmetry, which is in good agreement with our atomic simulations (Figure. S5).

In accordance with recent reports Huder2018 ; qiao2018twisted , high energy peaks appear in AB regions, as shown in Figs. 4c and 4d. These peaks, associated with a partial band gap opening in a high energy moiré band Huder2018 ; Wong2015 or localized domain wall modes arising from strongly confined moiré potentials qiao2018twisted , suffer minor changes in the presence of heterostrain and atomic deformations (see Supplementary Information Sec. H for the effects of heterostrain and deformation on the electronic properties of TBG). This is expected since these lattice deformations and pseudo-magnetic fields only occur around AA regions, as shown in Fig. 4e and 4f. Due to the considerable heterostrain applied to the bottom layer here, the shape and amplitude of the arising PMFs in each layer are no longer identical. Obviously, the strained layer has been more extended than the unstrained one, resulting in larger PMFs, which go up to 7 T. This tendency is in good agreement with previous findings Katsnelson2015 ; Gargiulo2018 . Besides, the PMFs are stronger along the direction of the uniaxial heterostrain within each circular regions. As a consequence of uniaxial strain, the local deformation potential $V$ on the strained layer is 0.035 eV higher than that of unstrained one. Thereby, both the rotation angle and uniaxial heterostrain can be utilized to fine-tune the geometry, amplitude, and periodicity of the PMFs. In Fig. 5, we plot the twisted angle dependent PMFs and their spatial distributions around the AA regions. The locations in experiments are equally taken with interval of 1 nm or 2 nm. The experimental fitted results shows that the PMFs increase with the decreasing twist angles and the PMF area occurs near the AA regions which reaches its maximum at the AA/AB transitions. The experiments and calculations show good agreements.

Collectively, the central findings achieved here may have important far-reaching implications. First of all, the existence of PMFs in TBG systems may finally offer testing grounds for observing quantum anomalous Hall effects as originally proposed haldane1988model . Secondly, such PMFs are opposite for electrons at the K and K’ valleys, leading to clockwise or anticlockwise valley vortex currents around the AA regions, which can be exploited for establishing quantum valley Hall effects Guinea2010 . Thirdly, the present study has demonstrated that the atomic deformations play a critical role in determining the electronic properties of the low-angle TBG, which in turn should also influence the emergent physics of correlated electrons yankowitz2019tuning ; chen2019evidence ; yoo2019atomic ; kerelsky2019maximized ; xie2019spectroscopic ; jiang2019charge ; choi2019electronic . Furthermore, similar physical phenomena, including the appearance of PMFs, may also be present in other van der Waals systems (for example, transition metal dichalcogenides bilayers or heterostructures tran2019evidence ; seyler2019signatures ; jin2019observation ), and their effects on the electronic properties remain to be fully explored.

Data availability  

The whole datasets are available from the corresponding author on reasonable request.

 

Code availability  

The codes used to construct the twisted bilayer graphene and calculate its electronic properties are available from the corresponding author on reasonable request.