Pressure induced compression of flatbands in twisted bilayer graphene

The studies on the electronic properties of twisted bilayer graphene (tBLG) has received recently a renewed boost of interest thanks to groundbreaking discoveries of flatband superconductivity [1] and Mott gaps [2, 3] for certain magic angles close to  ∼$1^{\circ}$ where a high density of states is generated near the charge neutral point. The possibility of tailoring narrow flatbands in systems with such remarkable simplicity in composition as graphene consisting of only carbon atoms makes it an attractive testbed to understand its microscopic mechanisms of electron correlations and its coupling with lattice vibrations as fully as possible [6–10]. Other van der Waals systems including trilayer graphene on hBN, transition metal dichalcogenide layers manifesting flatbands and topological superlattice bands have been proposed in recent literature [11–16].

In twisted bilayer graphene the graphene layers' electronic structure can be considered effectively decoupled when the twist angles are large enough, typically above $\theta \gtrsim 10^{\circ}$ as the twisted multilayer graphene stacks grown on SiC substrates [4, 5]. The progressive increase in coupling is manifested in the gradual decrease of the Fermi velocity as the twist angle is reduced starting from the Fermi velocity of graphene for large twist angles [17–21]. It is only in the regime of small twist angles on the order of  ∼$1^{\circ}$ that the moiré patterns are long enough to reduce the Fermi velocity band dispersion to the point of flattening the bands almost completely due to non-perturbative coupling between the electronic states separated by short moiré reciprocal lattice ${\bf G}$-vectors in momentum space. In this limit a series of magic twist angles of $\theta = 1.05^{\circ}, 0.5^{\circ}, 0.35^{\circ}, 0.24^{\circ}$, and $0.2^{\circ}$ have been predicted numerically [20] that approximately follow the $\theta \sim (1/n){}^{\circ}$ sequence for $n=1, 2, \ldots$, where the band slopes assume zero values at the $\tilde{\Gamma}$ point of the moiré Brilouin zone, and the bandwidths achieve a series of minima. One important hurdle for the realization of the flatbands in twisted bilayer graphene lies at the high sensitivity of electronic structure to twist near the magic angles because the bandwidths can undergo variations on the order of  ∼10 meV for small variations in twist angle $\Delta \theta \sim 0.1^{\circ}$ leading to drastic changes in the electronic properties of tBLG.

In this work we show that application of appropriate pressure to the system is an additional control knob that can control the flatband bandwidth of the system. Contrary to the twist angles whose values are determined once for all for each fabricated device, the vertical pressure is a continuously variable system parameter controllable with a hydraulic press [22, 23] that should be applicable also in twisted bilayer graphene devices. Our work presents a roadmap for tailoring flatbands in twisted bilayer graphene even when the twist angle control is not strictly precise. We will begin by briefly introducing the model Hamiltonian used in our calculations, then move on to discuss the bandwidth phase diagram of the system as a function of pressure and twist angle, and further discuss the effects of an interlayer potential difference in the electronic structure of the flatbands before closing the paper with the summary and discussions.

The Hamiltonian of twisted bilayer graphene at valley K is based on the continuum model of the graphene layer Hamiltonian [20] perturbed by stacking-dependent interlayer tunneling and intralayer potential variations [21]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h (\theta)= (\upsilon_{F} \hat{R}_{-\theta}{\bf p} + {\bf A}({\bf r})) \cdot \sigma_{xy} + V({\bf r}) \mathbb{1} + \Delta({\bf r}) \sigma_z \nonumber \end{align} \tag{ 1 }$$

where $\sigma_{xy} = (\sigma_x, \sigma_y)$ and $\sigma_z$ are the graphene sublattice pseudospin Pauli matrices, the momentum is defined in the xy-plane ${\bf p} = (\,p_x, p_y)$ and $\hat{R}_{-\theta}$ introduces a phase shift in the off-diagonal term in the Dirac Hamiltonian ${\rm e}^{\pm {\rm i} \theta_{\bf p}} \rightarrow {\rm e}^{\pm {\rm i} (\theta_{\bf p} - \theta)}$ to account for the rotation of the layers, where $\theta_{\bf p}$ is measured with respect to the x-axis and θ is the rotation of the graphene layer with respect to x. The graphene layers can be coupled through a stacking-dependent interlayer tunneling $T({\bf r})$ [20] resulting in the Hamiltonian

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle H_{\rm tBLG} = \left(\begin{array}{@{}cc@{}} h(-\theta/2) & T({\bf r}) \nonumber \\ T^{\dagger}({\bf r}) & h(\theta/2) \end{array}\right). \nonumber \end{align} \tag{ 2 }$$

The symmetric and opposite rotations of the top and bottom graphene layers allows to preserve the moiré Brillouin zone (mBZ) orientation and therefore of the stacking dependent moiré patterns $V({\bf r})$, $\Delta({\bf r})$, ${\bf A}({\bf r})$ for intralayer potential variations, local mass term and virtual strains due to pseudomagnetic field vector potentials, in addition to the spatially varying interlayer tunneling $T(\bf {r})$ for every position ${\bf r}$. Using the first harmonic approximation we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M({\bf r}) = \sum_{m=1}^{6} {\rm e}^{{\rm i} {\bf G}_m {\bf r}} M_{m} \label{moirepattern} \nonumber \end{align} \tag{ 3 }$$

where the moiré patterns $M({\bf r})$ can be modeled through its Fourier coefficients M_m. More explicitly, for a triangular lattice we can write the scalar and vector moiré patterns as [21, 24]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle V({\bf r}) = 2 C_V \, {\rm Re}\left[ {\rm e}^{{\rm i} \phi_V} f({\bf r}) \right], \,\,\, \Delta ({\bf r}) = 2 C_{\Delta} \, {\rm Re}\left[ {\rm e}^{{\rm i} \phi_{\Delta}} f({\bf r}) \right] \nonumber \end{align} \tag{ 4 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\bf A}({\bf r}) = 2 C_{AB} \, \hat{z} \times \nabla {\rm Re}\left[ {\rm e}^{{\rm i} \phi_{xy}} f({\bf r}) \right] \nonumber \end{align} \tag{ 5 }$$

where we have used the auxiliary function $f({\bf r}) = \sum_{m=1}^6 {\rm e}^{{\rm i} {\bf G}_m \cdot {\bf r}} (1 + (-1){}^m)/2$, where the six first shell mBZ reciprocal lattice vectors are ${\bf G}_m = \hat{R}_{2\pi(m-1)/3} {\bf G}_1$ for m indices running from 1 to 6 and are generated through successive rotation by $2\pi/3$ of the vector ${\bf G}_1 \simeq (0, 4\pi \theta / \sqrt{3} a)$, where $a=2.46~\mathring{\rm A}$ is the lattice constant of graphene. The moiré pattern Hamiltonian parameters obtained from ab initio calculations in sublattice representation for rigid bilayer graphene are [21] $C_{AA} = C_{B'B'} = 1.10 \, {\rm meV}$, $\varphi_{AA} = \varphi_{B'B'} = 82.54^{\circ}$, $C_{BB} = C_{A'A'} = C_{AA}$, $\varphi_{BB} = \varphi_{A'A'} = -\varphi_{AA}$, $C_{AB} = 2.235~{\rm meV}$, $\varphi_{AB} = 0^{\circ}$, and with LDA out of plane relaxation we have $C_{AA} = 2.3~{\rm meV}$, $\varphi_{AA} = 27.5^{\circ}$, $C_{BB} = C_{AA}$, $\varphi_{BB} = -\varphi_{AA}$, $C_{AB} = 2.08 ~{\rm meV}$, $\varphi_{AB} = 0^{\circ}$, where we use the notation $\phi_{xy} = \pi/6 - \varphi_{AB}$. In our modeling of the rigid twisted bilayers the intralayer moiré patterns have a small effect in the electronic structure and they can be neglected, whereas the relaxed moiré pattern parameters lead to particle-hole symmetry breaking as we will show later on. The momentum conservation condition in twisted bilayer graphene

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\bf k}' = {\bf k} + {\bf G} \nonumber \end{align} \tag{ 6 }$$

implies that a Bloch state with momentum ${\bf k}$ from one layer scatters to ${\bf k}'$ at the other layer through a moiré reciprocal lattice vector ${\bf G}$ [21]. In the small angle approximation we have ${\bf G} \simeq - \theta \hat{z} \times {\bf g}$ where the reciprocal lattice vector of graphene is represented through ${\bf g}$. If we consider the ${\bf q} = {\bf k} - {\bf K}$ and ${\bf q}' = {\bf k}' - {\bf K}'$ momenta measured with respect to the Dirac points of each layer relatively displaced by $\Delta {\bf K} = {\bf K}' - {\bf K}= 2K \sin(\theta/2)$ we have the relationship ${\bf q}' = {\bf q} + {\bf K} - {\bf K}' + {\bf G} = {\bf q} + {\bf Q}$, where the three ${\bf Q}_j$ vectors are given by ${\bf Q}_0 = K \theta (0, -1)$ and ${\bf Q}_{\pm} = K \theta (\pm \sqrt{3}/2, 1/2)$ in the small angle approximation to represent the interlayer coupling Hamiltonian

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T({\bf r}) = \sum_{j} {\rm e}^{-{\rm i} {\bf Q}_j {\bf r}} T^{\,j}_{s, s'}, \label{interlayercoupling} \nonumber \end{align} \tag{ 7 }$$

where the interlayer coupling matrices are given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T^0 = \omega \left(\begin{array}{@{}cc@{}} 1 & 1 \nonumber \\ 1 & 1 \end{array}\right), \,\, T^{\pm} = \omega \left(\begin{array}{@{}cc@{}} 1 & {\rm e}^{\mp {\rm i} 2\pi/3} \nonumber \\ {\rm e}^{\pm {\rm i} 2\pi/3} & 1 \end{array}\right). \nonumber \end{align} \tag{ 8 }$$

We note that these T^j matrices result when the twist is applied to a bilayer with $\tau = (0, 0)$ initial AA stacking configuration and differ by an additional phase of ${\rm e}^{-{\rm i} {\bf G}_j \tau}$ with respect to the AB stacking case where $\tau = (0, a/\sqrt{3})$. The smooth variation of $T({\bf r})$ in equation (7) can be traced back to the relatively larger interlayer distance of $c \sim 3.35~\mathring{\rm A}$ when compared to inter-carbon distances of $a_{\rm CC} \sim 1.42~\mathring{\rm A}$ [20] and effectively implies that the bilayer graphene interlayer coupling strength can be described by a single parameter $\omega = t_1/3 \sim 0.113$ eV when $c=3.35~\mathring{\rm A}$, and somewhat weaker $\omega \sim 0.1$ eV when out of plane LDA relaxations are allowed between the layers [21]. We note that ω is proportional to the interlayer coupling term $t_1 \sim 3 \omega$ of commensurate bilayer graphene evaluated at the Dirac point when only the three ${\bf G}$-vector contributions of interlayer coupling nearest to the Dirac point K are considered. The interlayer coupling can in principle be modeled more accurately by including the effects of the Fourier components for larger ${\bf G}$-vectors in momentum space that become more relevant in the presence of out of plane corrugations and in-plane strains [21, 25]. Given the high accuracy of the single parameter interlayer coupling model for describing the electronic properties of rigid twisted layers, we will focus on the role of ω in our precisely defined continuum model and defer the discussions about the moiré strains and larger momenta Fourier components for future work.

In the following we obtain the phase diagram map of the bandwidth as a function of pressure to show how increase in pressure by a few GPa can trigger the appearance of flatbands in tBLG for twist angles that are larger by a fraction of a degree than the magic angles at zero pressure. We have calculated the phase diagram map of the moiré bands bandwidth of twisted bilayer graphene as a function of the three main parameters in the model Hamiltonian, namely the Fermi velocity of each graphene layer $\upsilon_{\rm F}$, the twist angle θ, and the pressure dependent interlayer coupling strength ω. The Fermi velocity $\upsilon_{\rm F}$ can change with the dielectric environment due to many-body effects but otherwise can be assumed to be constant. The twist angle θ can be controlled at will during device fabrication but is difficult to modify afterwards in a controlled manner, whereas the pressure that controls the interlayer coupling ω can in principle be varied continuously in a hydraulic press up to values of 2.5 GPa [22]. Our main results are summarized in figure 1 where we represent the flatband bandwidth colormap obtained as a function of interlayer coupling parameter ω and twist angle θ. Note that the bandwidth for the low energy conduction and valence bands are electron–hole symmetric in the absence of intralayer moiré pattern terms in the Hamiltonian. The role of intralayer moire patterns in the electron–hole asymmetry of the bandwidth is illustrated in figure 2. From the analysis of the numerically calculated bandwidth phase diagrams for different Fermi velocities we find a fitting formula for the three visible flatband magic angle lines

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \theta^{\circ}_{n} = C_n \frac{\omega}{\left| t_0 \right|} \quad ({\rm deg}) \label{mangle} \nonumber \end{align} \tag{ 9 }$$

where C₁  =  27.8, C₂  =  11.5, and C₃  =  5.46, that lead to $\theta_1 = 1.07^{\circ}$, $\theta_2 = 0.44^{\circ}$ and $\theta_3 = 0.21^{\circ}$ with t₀  =  −2.6 eV and $\omega=0.1$ eV. These results are in fair agreement with the numerical magic angles in [20] that approximately follow the $\theta \sim (1/n){}^{\circ}$ sequence for t₀  =  −2.7 eV and $\omega=0.11$ eV Hamiltonian parameters, except for the angles of $\theta = 0.35^{\circ}, 0.24^{\circ}$ which are not clearly resolved in our bandwidth phase diagram. Our flatband magic angles are linearly proportional to the interlayer coupling ω, and are inversely proportional to the Fermi velocity $\upsilon_{\rm F} = \sqrt{3} \left| t_0 \right| a /2\hbar$. We note that in our case the magic angles were obtained from bandwidth minima in the ω and θ parameter space, while previous calculations identified the magic twist angles from the zero slope in the band dispersion at the $\tilde{\Gamma}$-point in the mBZ. One practical implication of our findings is that it should be possible to access the flatbands by increasing ω with pressure by a few GPa when the twist angles in bilayer graphene are a fraction of a degree larger than the zero pressure magic angles.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The Fermi velocity that we use as one of the free parameters in our model is an intrinsic property of graphene that can be enhanced by Coulomb interactions and is therefore subject to specific environment and device quality. Strictly speaking, a logarithmic divergence is expected for the Fermi velocity at close proximity of the Dirac point due to the long range of the Coulomb tail which introduces k-dependent dispersion slope changes [26–28]. Yet a constant enhanced Fermi velocity often gives an excellent fit to experimental data, with the lower $\upsilon_F\sim 1 \times 10^6$ or  ∼$1.05 \times 10^6$ fitting well experiments of graphene on SiC or SiO₂ substrates [4, 29] and CVD grown twisted bilayer graphene [30], while higher $\upsilon_F\sim 1.1 \times 10^6$ m s⁻¹ are better for fitting the experimental data in high quality graphene devices with hexagonal boron nitride barrier materials [31]. In this work, we used ab initio LDA calculation values of t₀  =  −2.6 eV for the intralayer hopping term [32] that is in the lower end of the spectrum with $\upsilon_{\rm F} \sim 0.84 \times 10^6$, and close to $t_0 \sim -2.7$ eV used in [20]. Our tight-binding Fermi velocity choice is more appropriate for band theories that intend to introduce the many-body corrections explicitly on top of the non-interacting model.

The twist angles range examined in the phase diagram lie between 0.05° and 2.5°. Convergence of the eigenvalues and eigenvectors of our continuum model can be expected when truncation of the moiré reciprocal lattice vector is of the order of $k \sim 2 \omega / (\sqrt{3} a \left| t_0 \right|)$ [20], requiring a larger cutoff when ω is larger. We have used a cutoff in momentum space for a radius of about six moiré reciprocal lattice vectors $6 G_1=24\pi \theta/ (\sqrt{3} a)$ using fixed Hamiltonian matrix sizes of 676  ×  676 to obtain the phase diagram in figure 1, which should be valid for sufficiently large $\theta \gtrsim \omega / ( 12 \pi \left| t_0 \right|)$ or in the limit of small ω. Our model assumes rigid in-plane lattices although for systems with small twist angles $\theta \lesssim 0.5^{\circ}$ one expects structural instabilities associated with the commensuration moiré strains due to the reduction of the elastic energies that scale with the twist angle as $\propto \theta^2$ [25].

Pressure can be varied in a continuous manner even after device fabrication to modify the magnitude of interlayer tunneling. Recent experimental progress that make use of a hydraulic press allowed to achieve continuously variable pressures of up to 2.5 GPa for hBN encapsulated graphene devices while values as large as 5 GPa might be achievable by optimizing the design of the press, and even larger pressure would be achievable with diamond anvil cells. These values are still safely below pressures where structural phase transitions from graphite to diamond can start taking place [33, 34]. The relationship between the interlayer coupling parameter ω and pressure P is obtained by combining Fourier transformed interlayer hopping data between maximally localized Wannier functions obtained from LDA ab initio calculations evaluated at three fixed interlayer separation distances of $c = 3.35, \, 3.2, \, 3.1~\mathring{\rm A}$ for different interlayer stacking, see [21, 22, 25] for related studies. The pressures for different interlayer distances at AB and BA stacking configurations within the LDA are given by $P=0, \, 2.01$ and 3.45 GPa, and larger values for AA stacking of $P= 2.09, 4.80, 7.73$ GPa within LDA consistent with the calculations in [35], which can be fitted with the Murnaghan equation of state $P(V) = (K_0/K'_0) \left[ (V/V_0){}^{-K'_0} - 1 \right]$ [36] or the third order Birch–Murnaghan equation of state [37] using the parameters listed in table 1. By assuming approximately equal distribution of AB, BA and AA stacking areas [22, 25] and averaging the values of interlayer tunneling at the Dirac point for different stacking configurations we obtain a polynomial fit for the relationship between pressure and averaged interlayer tunneling ω through

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P = A\, \omega^2 + B\, \omega + C, \label{polynomial} \nonumber \end{align} \tag{ 10 }$$

whose numerical parameters are A  =  455.5 GPa (eV)⁻², B  =  −71.05 GPa eV⁻¹, C  =  3.281 GPa and $\omega_0=0.098$ eV is the tunneling for P  =  0. Because the tunneling is weaker for larger interlayer distances our $\omega_0$ consistent with relaxed LDA results is smaller than $\omega = 0.113$ eV obtained from the average of rigid graphene bilayers separated at fixed $c=3.35~\mathring{\rm A}$ [21]. For sake of definiteness we will use for our zero pressure continuum model Hamiltonian the parameters t₀  =  −2.6 eV and $\omega=0.1$ eV that lead to electronic structure results closely similar to those obtained in [20].

Table 1. Bulk moduli K₀ and $K_0^\prime$ for Murnaghan and B₀ and $B_0^\prime$ in kbar units for the third-order Birch–Murnaghan equations of state consistent with the interlayer potentials obtained in [35] from ab initio calculations for different local stacking configurations.

	AA		AB		BA
	LDA	RPA	LDA	RPA	LDA	RPA
K0	322	353	308	358	308	358
$K_0^\prime$	10.88	12.59	12.51	12.02	12.51	12.02
B0	239	331	302	354	302	354
$B_0^\prime$	14.27	15.00	14.25	13.48	14.25	13.48

The weaker interlayer tunneling regime $\omega < 0.1$ eV in the phase diagram might be achievable in systems whose average interlayer distances can be reduced either through intercalation of ions or addition of a barrier hBN. The use of an intercalation hBN spacer layer between the top and bottom graphene layers to reduce interlayer tunneling can prevent the structural instabilities for small twist angles that are present when both graphene layers are in direct contact.

The sharp increase in the density of states (DOS) and its width at the flatbands relative to Coulomb interaction strength determines how prone the system is towards the development of broken symmetry phases. Here we analyze the impact that the interlayer coupling strength ω has in the density of states (DOS) when its value is increased along the flatband line equation given in equation (9) for a few select twist angles and interlayer tunneling ω. The increase of the interlayer tunneling parameter ω expands the energy range around which a given Bloch state in one layer can scatter to the other layer and affects the eigenvectors associated with the flatbands. An immediate effect of having larger magic twist angles is that we can expect an enhancement of the flatband DOS due to the increase of mBZ area proportionally to $\theta^2$, making more states accessible in momentum space, because more electrons per unit area are required to fill the same number of the moiré flatbands. The reduction in the moiré pattern length implies a relative enhancement of the effective Coulomb interaction scale $\newcommand{\e}{{\rm e}} U_{\rm ef\,\!f}= e^2 / (4 \pi \varepsilon_0 \varepsilon_r \ell_M) \sim e^2 \theta / ( 4 \pi \varepsilon_0 \varepsilon_r a)$ where $\newcommand{\e}{{\rm e}} \ell_M \sim a/ \theta$, and therefore an increase from $\theta \sim 1^{\circ}$ to $1.5^{\circ}$ leads to already a 50% increase. We show in figure 3 the density of states and the local density of states associated with the charge neutrality flatbands for different values of twist angle θ and interlayer tunneling ω where we can observe a steady increase in the flatband DOS either for its peak height and width for progressively larger twist angle sizes. In the sub-panel figure 3(a), we show the density of states (DOS) of the flatbands corresponding to the first, second and third magic angles calculated for different values of interlayer tunneling ω. The progressive increase in the DOS maxima and width variations for larger twist angle flatbands can be related to the increase of the moiré Brillouin zone area. The local density of states plots in figure 3(b) clearly indicates the broadening in energy of the density of states for higher order magic angles of the electrons that localize at the AA stacking region. This broadening of the flatbands for smaller magic angles reflects the relative increase of interlayer coherence. The increase of interlayer coupling strength ω in the band structure shown in figure 3(c) introduces a mild widening of the flatbands as pressure increases along the first magic angle line.

Zoom In Zoom Out Reset image size

In summary we have investigated the phase diagram map of the low energy bandwidth evolution in twisted bilayer graphene as a function of parameters in the continuum model Hamiltonian including the Fermi velocity $\upsilon_{\rm F}$, the twist angle θ and the interlayer tunneling parameter ω in search of the phase space where we can achieve bandwidth minima, in particular when the interlayer tunneling is enhanced by means of external pressure. The flatbands for the continuum Hamiltonian can be summarized in a single line equation relating the minimum bandwidth magic twist angle with interlayer tunneling, and is inversely proportional to the Fermi velocity or intralayer hopping. Our calculations indicate that by applying pressures on the order of  ∼2 GPa achieved in recent experiments through a hydraulic press [22] it should be possible to access the first magic angle around $\theta \sim 1.5^{\circ}$ which is considerably greater than  ∼$1^{\circ}$ and therefore should have considerably greater structural stability against the moiré strains and have enhanced effective Coulomb interaction energy scales $\newcommand{\e}{{\rm e}} U_{\rm ef\,\!f} \propto \ell^{-1}_M \propto \theta$ thanks to the reduction in the moiré pattern size. Hence, we can envision that application of pressure in twisted bilayer graphene nanodevices to achieve larger magic angles should considerably facilitate access to flatbands and electron-electron interaction driven ordered phases.

Support from the Korean National Research Foundation is acknowledged for BLC through NRF-2017R1D1A1B03035932, for NL through NRF-2018R1C1B6004437 and for JJ by the 2018 Research Fund of the University of Seoul. NL acknowledges support by the Korea Research Fellowship Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science and ICT (KRF-2016H1D3A1023826).