Electronic structure of biased alternating-twist multilayer graphene

We consider a model of vertically stacked $N$ graphene layers with the $\ell$th layer alternatingly twisted by an angle $\theta_{\ell}=(-1)^{\ell}\theta/2$, as shown in Fig. 1. For a perpendicular electric field, we assume that the interlayer potential difference $U$ is the same between the two adjacent layers. Following Leconte et al. Leconte2022 , the Hamiltonian of ATMG in the presence of an interlayer potential difference can be expressed as

where the diagonal blocks $H_{\boldsymbol{k}}^{(\ell)}=\hbar v_{0}(\boldsymbol{k}\cdot\boldsymbol{\sigma}_{\theta_{\ell}})$ with $\boldsymbol{\sigma}_{\theta_{\ell}}=e^{\frac{i}{2}\theta_{\ell}\sigma_{z}}\boldsymbol{\sigma}e^{-\frac{i}{2}\theta_{\ell}\sigma_{z}}$ contain the Dirac cones of twisted graphene layers, $T(\boldsymbol{r})$ is the interlayer tunneling matrix, and $V={\rm{diag}}(V^{(1)}\mathbb{I}_{2},V^{(2)}\mathbb{I}_{2},...,V^{(N)}\mathbb{I}_{2})$ is a diagonal matrix that captures the interlayer potential differences. The Fermi velocity of the monolayer graphene is set as $v_{0}=\sqrt{3}a\lvert t\rvert/2\hbar\simeq 10^{6}$ m/s with the lattice constant $a=2.46$ $\rm\AA$ and the nearest-neighbor intralayer hopping parameter $t=-3.1$ eV.

In our model, the electric potential sequence $V^{(\ell)}$ is defined to satisfy both $V^{(\ell+1)}-V^{(\ell)}=U$ and $V^{(1)}=-V^{(N)}$, and the interlayer tunneling at the twisted interface takes the form

where $\boldsymbol{q}_{0}$, $\boldsymbol{q}_{\pm}$ are given by $\boldsymbol{q}_{0}=2k_{\rm D}\sin(\theta/2)(0,-1)$, $\boldsymbol{q}_{\pm}=2k_{\rm D}\sin(\theta/2)(\pm\sqrt{3}/2,1/2)$ with the Dirac momentum for monolayer graphene $k_{\rm D}=4\pi/3a$. Our model also considers the corrugated lattice structure due to the effect of out-of-plane relaxation, leading to the larger interlayer spacing in AA-stacking region than the AB/BA-stacking region Uchida2014 ; Wijk2015 , thus resulting in unequal intra/inter sublattice hopping terms $w$ and $w^{\prime}$, respectively. Following the convention of an initial AA stacking configuration Jung2014 , the interlayer tunneling matrices are given by

where $w^{\prime}=0.0939$ eV and $w=0.12$ eV Chebrolu2019 . For the full numerical calculations, we include the lattice corrugation $(w\neq w^{\prime})$, whereas we assume the rigid model of equal hopping terms ($w=w^{\prime}=0.12$ eV) for simplicity when we study the Hamiltonian analytically. We choose representative twisted angles $3^{\circ}$$-$$5^{\circ}$ for the numerical calculations above the typical magic angle values that lie between $1^{\circ}$$-$$2^{\circ}$. For these large angles the interlayer coupling substantially reduces the Fermi velocity of the dispersive Dirac cones near the two moiré Brillouin zone (mBZ) corners $\bar{K}$ and $\bar{K}^{\prime}$ but have not completely flattened them.

In the absence of the interlayer potential difference, the effective Hamiltonian of the ATMG at $\bar{K}$ and $\bar{K}^{\prime}$ can be described as a set of TBG models at different angles with an additional decoupled monolayer graphene model at $\bar{K}$ (or at $\bar{K}^{\prime}$ depending on the continuum model we start with) for an odd number of layers Khalaf2019 . The electronic structure of ATMG has a close analogy with Bernal-stacked multilayer graphene where the effective Hamiltonian is described by a set of bilayer graphene models with different effective masses with an additional decoupled monolayer graphene model for an odd number of layers min2008a ; min2008b . This analogy between Bernal stacked multilayer graphene and ATMG can be expanded to their wave functions.

We now construct the wave function of ATMG at $\bar{K}$ or $\bar{K}^{\prime}$ using the first shell model, in which the momentum-space lattice is truncated at the nearest-neighbor shell of the moiré reciprocal lattice $\bm{G}$ vectors, assuming the rigid model ($w=w^{\prime}$). For the bilayer case ($N=2$), the zero-energy eigenstates near $\bar{K}$ and $\bar{K}^{\prime}$ consist of four two-component spinors as follows:

where $\alpha=w/\left[2v_{0}k_{\rm D}\sin(\theta/2)\right]$ is a dimensionless parameter. Following Bistritzer and MacDonald Bistritzer2011 , we define $a_{\lambda}$ as a normalized eigenstate of $\hat{\boldsymbol{k}}\cdot\boldsymbol{\sigma}_{\theta_{\ell}}$ corresponding to the eigenvalue $\lambda=\pm 1$, and $b_{\lambda}=(b_{\boldsymbol{q}_{0},\lambda},b_{\boldsymbol{q}_{+},\lambda},b_{\boldsymbol{q}_{-},\lambda})^{\rm T}$ determined by the equation $b_{\boldsymbol{q}_{j},\lambda}=-h_{j}^{-1}T_{j}^{\dagger}a_{\lambda}$ with $h_{j}=\hbar v_{0}(\boldsymbol{k}+\boldsymbol{q}_{j})\cdot\boldsymbol{\sigma}_{\theta_{\ell}}$. In a similar way to Bernal stacked multilayer graphene, the eigenfunctions of ATMG have the form of the solution of a one-dimensional chain problem, thus we can construct the eigenfunctions of ATMG in the following manner Khalaf2019 :

where $\tau=2-\delta_{r,n+1}$, $\theta_{r}=r\pi/(N+1)$ with $r=1,2,\ldots,n$ for even $N=2n$ or for odd $N=2n+1$ with additional $r=(n+1)$th mode near $\bar{K}$, and $\psi^{\rm{TBG}}_{r,\lambda}$ can be obtained from $\psi^{\rm{TBG}}_{\lambda}$ in Eq. (4) by letting $\alpha\rightarrow t_{r}\alpha$ and $b_{\lambda}\rightarrow t_{r}b_{\lambda}$ with $t_{r}=2\,\cos{\theta_{r}}$. Here, $\Psi_{r,\lambda}=(\Psi_{r,\lambda}^{(1)},\Psi_{r,\lambda}^{(2)},\ldots,\Psi_{r,\lambda}^{(N)})^{\rm T}$ is a normalized eigenstate of the effective Hamiltonian $H_{\rm eff}=\hbar v_{r}(\boldsymbol{k}\cdot\boldsymbol{\sigma})$ of ATMG with $\lvert\Psi_{r,\lambda}\rvert^{2}=1$, where $v_{r}$ is a Fermi velocity of the Dirac cone $\Psi_{r,\lambda}$ given by

Inserting $r=n+1$ in Eq. (6), one finds that the $(n+1)$th mode for the odd number of layers corresponds to an eigenstate of the decoupled monolayer Hamiltonian.

In the presence of the interlayer potential difference $U$, we obtain analytically the low-energy effective Hamiltonian using first-order degenerate-state perturbation theory based on the minimal size Hamiltonian including only the first shell of the moiré $\bm{G}$ vectors. Due to the electric field, the Dirac cones near $\bar{K}$ or $\bar{K}^{\prime}$ are hybridized and split from one another, so the effective Hamiltonian of each Dirac node would be altered in the following form:

where $\Delta(\alpha,U)$ and $v^{*}$ are the energy shift and modified Fermi velocity of the effective Hamiltonian, respectively. The energy shift due to the interlayer potential $U$ can be expressed as $\Delta(\alpha,U)=C(\alpha)U$, and the Fermi velocity $v^{*}$ can be expressed as a linear combination of $v_{r}$. In the following Secs. II.B and II.C, we illustrate the effective Hamiltonian of $N=3$ and $4$ ATMG as examples. We leave the discussions of the analytical results for the $N=5$ case to Appendix A.

Here we derive the effective Hamiltonian of alternating twist trilayer graphene (AT3G) at the $\bar{K}$ and $\bar{K}^{\prime}$ points of the mBZ. There are two Dirac cones centered at $\bar{K}$ with $v_{0}$ and $v_{1}$ Fermi velocities, as shown in Fig. 2(a), thus the size of the perturbation matrix $V_{\bar{K}}$ would be $2\times 2$. Note that $v_{0}$ represents the Fermi velocity of monolayer graphene. Using Eq. (4), we obtain the following normalized wave functions $\Psi_{r,\lambda}$ near $\bar{K}$ in our first shell model:

At $\bar{K}$ in AT3G, the perturbation $\hat{V}$ in the first shell model is given by $\hat{V}=\rm{diag}$($-U\mathbb{I}_{2},\boldsymbol{0}_{6},U\mathbb{I}_{2}$). Then, in the basis of the wave functions in Eq. (8), we obtain the perturbation matrix $V_{\bar{K}}$ using $V_{11}=V_{22}=0$ and $V_{12}=V_{21}=-U/\sqrt{1+12\alpha^{2}}$. By diagonalizing $V_{\bar{K}}$, we obtain the effective Hamiltonian of biased AT3G near $\bar{K}$ as

where $C(\alpha)=1/\sqrt{1+12\alpha^{2}}$ and $v^{*}=\left(v_{0}+v_{1}\right)/2$. Comparing the left inset of Figs. 2(a) and 2(b), we can deduce that the two Dirac bands near $\bar{K}$ are hybridized equally and split by $2C(\alpha)U$ acquiring the average Fermi velocity $v^{*}$ from the unbiased values.

On the other hand, near $\bar{K^{\prime}}$, only one Dirac cone with $v_{1}$ exists, whose wave function is

Then, the perturbation matrix $V_{\bar{K}^{\prime}}$ vanishes and the Dirac cone at $\bar{K}^{\prime}$ remains unaltered to leading order in $U$, resulting in the effective Hamiltonian

In Figs. 2(c) and 2(d), we illustrate the result of the leading-order energy splitting coefficient $C(\alpha)$ and the modified Fermi velocity $v^{*}$ obtained from the analytic model and numerical method as a function of twist angle.

In the following, we derive the effective Hamiltonian of alternating twist tetralayer graphene (AT4G) at $\bar{K}$ and $\bar{K}^{\prime}$. At $\bar{K}$, there are two Dirac cones with the velocities $v_{1}$ and $v_{2}$ as shown in Fig. 3(a), and the corresponding wave functions are given by

with $r=1,2$. The perturbation $\hat{V}$ at $\bar{K}$ in this case is given by $\hat{V}=\rm{diag}$($-\frac{3U}{2}\mathbb{I}_{2},-\frac{U}{2}\mathbb{I}_{6},\frac{U}{2}\mathbb{I}_{2},\frac{3U}{2}\mathbb{I}_{6}$). Since there are two Dirac cones at $\bar{K}$, the size of the perturbation matrix $V_{\bar{K}}$ would be $2\times 2$, and its elements $V_{rr^{\prime}}$ are expressed as

By diagonalizing the matrix $V_{\bar{K}}$, we obtain the effective Hamiltonian of biased AT4G near $\bar{K}$ as

where $C_{\pm}(\alpha)$ and $v_{\pm}^{*}$ corresponding to upward and downward shifted Dirac cones are given by

and

Here, $A_{\pm}$ and $B$ are unnormalized mixing coefficients of the two Dirac cones given, respectively, by

From the above result, we find that the two Dirac cones with the velocities $v_{1}$ and $v_{2}$ are hybridized with the ratio of $A_{\pm}$ and $B$, and shifted by $C_{\pm}(\alpha)U$, as schematically illustrated in Fig. 3(b).

On the other hand, near $\bar{K}^{\prime}$, the wave functions for two Dirac cones with the velocity $v_{r}$ $(r=1,2)$ are given by

Since $\sin{\ell\theta_{r}}=(-1)^{r}\sin{(N+1-\ell)\theta_{r}}$, the wave function at $\bar{K}^{\prime}$ can be obtained by reversing the components of the wave function at $\bar{K}$. Therefore, the effective Hamiltonian of biased AT4G near $\bar{K}^{\prime}$ can be obtained by reversing the sign of the interlayer potential difference $U$ in Eq. (14) as

where $C_{\pm}(\alpha)$ and $v_{\pm}^{*}$ are the same as those at $\bar{K}$. In detail, our model Hamiltonian [Eq. (1)] for $N=4$ has a combined symmetry expressed as

where

We note that $\hat{\Sigma}$ changes only the valley index ($K\leftrightarrow K^{\prime}$) keeping the same mBZ corner points ($\bar{K}\rightarrow\bar{K}$, $\bar{K}^{\prime}\rightarrow\bar{K}^{\prime}$) Moon2013 , whereas the time-reversal operator $\hat{\mathcal{T}}$ changes both the valley index ($K\leftrightarrow K^{\prime}$) and the mBZ corner points ($\bar{K}\leftrightarrow\bar{K}^{\prime}$). This combined symmetry is preserved even in the presence of an interlayer potential difference, thus the effective Hamiltonians between $\bar{K}$ and $\bar{K}^{\prime}$, which are respectively described in Eqs. (14) and (19), are also related as Eq. (20).

As the layer number $N$ is increased, the size of the perturbation matrix is also proportionally increased and it becomes progressively cumbersome to obtain analytically the effective Hamiltonian of ATMG for a large number of layers in the presence of an applied field even if we use the first shell model. Instead, here want to provide the general behavior patterns of the effective Hamiltonian of biased ATMG for arbitrary $N$. Tables 1 and 2 show the summary of the effective Hamiltonian for $N=2$$-$$8$ ATMG in the presence of the interlayer potential difference.

Firstly, for ATMG with an odd number of layers, there are $(N-1)/2$ TBG Dirac cones labeled by $v_{r}$ $\left[r=1,2,\ldots,(N-1)/2\right]$ near the two mBZ corners $\bar{K}$ and $\bar{K}^{\prime}$, and one decomposed monolayer Dirac cone with $v_{0}$ at $\bar{K}$. Regardless of the layer number $N$ and mBZ symmetry points, the form of the perturbation matrix $V$ is solely determined by the number of Dirac cones at $\bar{K}$ or $\bar{K}^{\prime}$, even though its elements depend on $N$. One can thus find the same form of the effective Hamiltonian when the number of Dirac cones is the same, as seen in Table 1. In detail, if there are $m$ Dirac cones at $\bar{K}$ or $\bar{K}^{\prime}$, we have $m/2$ pairs of Dirac-cones shifted by $\pm\Delta_{i}$ ($i=1,2,\ldots,m/2$) for even $m$, whereas we have $(m-1)/2$ pairs of Dirac cones plus one Dirac cone without energy shift for odd $m$. Each pair of Dirac cones has the same effective Fermi velocity.

For ATMG with an even number of layers, there are $N/2$ TBG Dirac cones labeled by $v_{r}$ $\left(r=1,2,\ldots,N/2\right)$ near $\bar{K}$ and $\bar{K}^{\prime}$ when $U=0$. If an external field is applied ($U\neq 0$), the Dirac cones are split with different energy shift $\Delta$ and Fermi velocity $v^{*}$. As already mentioned in Sec. II.C, the effect of an applied electric field at $\bar{K}^{\prime}$ ($\bar{K}$) can be effectively described by flipping its direction ($\hat{z}\rightarrow-\hat{z}$) at $\bar{K}$ ($\bar{K}^{\prime}$). Moreover, we can generalize Eq. (20) by expanding $\hat{\Sigma}$ symmetry in Eq. (21) by conveniently alternating $\sigma_{x}$ and $-\sigma_{x}$. Similarly to the $N=4$ case, the combined $\hat{\Sigma}\hat{\mathcal{T}}$ symmetry is still preserved for ATMG with an even number of layers in the presence of an interlayer potential difference, relating the effective Hamiltonians between $\bar{K}$ and $\bar{K}^{\prime}$ with flipped energy shifts, as seen in Table 2.

Lastly, let us consider the effective Hamiltonian of biased ATMG in the asymptotic limit ($\alpha\rightarrow 0$) where the twist angle $\theta$ becomes much larger than the first magic angle of ATMG. For the first shell model, $\lvert b_{\lambda}\rvert$ becomes proportional to $\alpha$, so only monolayer terms $a_{\lambda}$ of $\Psi_{r,\lambda}$ survive in this limit. Thus, the energy splitting coefficient $C(\alpha)$ of ATMG with arbitrary $N$ at $\bar{K}$ ($\bar{K}^{\prime}$) can be obtained as odd (even) layer components of $\hat{V}$, as schematically shown in Fig. 4.

On the other hand, the modified Fermi velocity $v^{*}$ converges to $v_{0}$ as $\alpha\rightarrow 0$, since all eigenstates of biased ATMG in this limit have just a single monolayer term, giving the monolayer graphene Dirac cone with the velocity $v_{0}$. Therefore, the effective Hamiltonian would be described by a set of monolayer graphene Hamiltonian with the energy shift described by $C(\alpha\rightarrow 0)U$, which can be obtained by the pattern presented in Fig. 4. Figure 5 shows the band structure of $N=5-8$ ATMG in the presence of the interlayer potential difference $U=0.1$ eV at $\theta=5^{\circ}$ along with the analytical result obtained in the asymptotic limit, which agrees closely with the full numerical result except for small deviations in the Fermi velocity of the Dirac cones.