Resonant interaction in chiral, Eshelby-twisted van der Waals atomic layers

We consider an infinitely stacked CTG, where each layer is rotated with respect to the lower layer by a fixed angle $\theta$. Only the structures with $0<\theta\leq 30^{\circ}$ are distinct due to the symmetry of the hexagonal lattice. Such a structure is a helical lattice which has the symmetry of an in-plane rotation by $\theta$ accompanied by a translation to the out-of-plane direction by the interlayer spacing, i.e., a $\theta$ screw rotational symmetry, and has no in-plane periodicity. We set $xy$ coordinates parallel to the graphene layers and $z$ axis perpendicular to the plane. We define the atomic structure of the system by starting from a nonrotated arrangement, where the hexagon center of all the layers share the same in-plane position $(x,y)=(0,0)$, and the $A$-$B$ bonds are parallel to each other. We choose ${\bf a}_{1}=a(1,0)$ and ${\bf a}_{2}=a(1/2,\sqrt{3}/2)$ ($a=0.246~{}\mathrm{nm}$) as the primitive lattice vectors of graphene, and $\bm{\tau}_{A}=-\bm{\tau}_{1}$ and $\bm{\tau}_{B}=\bm{\tau}_{1}$ [$\bm{\tau}_{1}=-(1/3)({\bf a}_{1}-2{\bf a}_{2})$] as the coordinates of the $A$ and $B$ sublattices in the unit cell. Then, we rotate the $l$-th layer by $l\theta$ and get the atomic positions of the $l$-th layer

$${\bf R}_{X}^{(l)}=n_{1}{\bf a}_{1}^{(l)}+n_{2}{\bf a}_{2}^{(l)}+\bm{\tau}_{X}^{(l)}+ld{\bf e}_{z},$$

(1)

where $X=(A,B)$ denotes the sublattice index, $n_{i}$ ($i=1,2$) are integers, ${\bf a}_{i}^{(l)}=\mathcal{R}(l\theta){\bf a}_{i}$ and $\bm{\tau}_{X}^{(l)}=\mathcal{R}(l\theta)\bm{\tau}_{X}$, where $\mathcal{R}(l\theta)$ is a counterclockwise rotation by $l\theta$, $d=0.335~{}\mathrm{nm}$ is the interlayer spacing between two adjacent layers and ${\bf e}_{z}$ is the unit vector normal to the layer. The structure has a nonsymmorphic symmetry around the center of the rotation. We show the layers with $l=0,1,2$ of CTG with $\theta=10^{\circ}$ in Fig. 2(a). We define the reciprocal lattice vectors ${\bf b}_{i}^{(l)}$ of each layer so as to satisfy ${\bf a}_{i^{\prime}}\cdot{\bf b}_{i}=2\pi\delta_{i^{\prime}i}$ and ${\bf b}_{i}^{(l)}=\mathcal{R}(l\theta){\bf b}_{i}$, and plot the reciprocal space of Fig. 2(a) in Fig. 2(b).

Note that, although another similar material which constituent layers stacked with an alternating twist angle, $\theta$, $-\theta$, $\theta$, $-\theta$, $\dots$, will also host an infinite chain of resonant interaction at the middle of the two Dirac points of the neighboring layers (i.e., $\bar{M}$ point of twisted bilayer graphene), such structures are an intuitive expansion of the periodic bilayer moiré superlattices. Their electronic structures can be easily obtained by the conventional effective theory which is based on the moiré periodicity. Thus, we do not consider such structures with obvious in-plane periodicity in this work, except the structure with $\theta=30^{\circ}$ which belongs also to the CTG class of materials.

No CTG has translational symmetry along the in-plane direction, since the moire patterns defined by each pair of layers are not commensurate. In addition, most CTGs are not periodic along the $z$ axis as well (hereafter “incommensurate CTG”). At specific $\theta$, however, $l$-th layer can have an in-plane lattice configuration the same as $(l+iN)$-th layers ($i\in\mathbb{Z}$) for a fixed number $N\in\mathbb{Z}$. Then the system is periodic with respect to the translation by $Nd{\bf e}_{z}$ (hereafter “commensurate CTG”), and we can choose the $N$ successive layers as a primitive cell. The allowed $\theta$ for commensurate CTG is

$$\theta=60^{\circ}\frac{M}{N}$$

(2)

where $M\in\mathbb{Z}^{+}$, $M\leq N/2$, and $\gcd(M,N)=1$. $N$ and $M$ fully define the lattice geometry of commensurate CTG. The structures with $N=2,3,4$, where the primitive cell is bi-, tri-, quad-layer, have only one lattice configuration with $\theta=30^{\circ},20^{\circ},15^{\circ}$, respectively, while that with $N=5$ has two distinct configurations with $\theta=12^{\circ}$ and $24^{\circ}$, which correspond to different $M$ (=1, 2). The primitive cell with $N=2$ has the geometry the same as that of the twisted bilayer graphene quasicrystal [15, 17].

We use a tight-binding model of carbon $p_{z}$ orbitals to describe the electronic structure of the general CTGs. We define the Bloch state of the $l$-th layer with the two-dimensional Bloch wave vectors ${\bf k}$ as

$$|{\bf k},X^{(l)}\rangle=\frac{1}{\sqrt{\tilde{N}}}\sum_{{\bf R}_{X}^{(l)}}e^{i{\bf k}\cdot{\bf R}_{X}^{(l)}}|{\bf R}_{X}^{(l)}\rangle$$

(3)

in incommensurate CTGs and as

	$\displaystyle|{\bf k},X^{(l)}\rangle$
	$\displaystyle=\frac{1}{\sqrt{\tilde{N}}N_{z}}\sum_{i_{z}}\sum_{{\bf R}_{X}^{(l+i_{z}N)}}e^{i{\bf k}\cdot{\bf R}_{X}^{(l+i_{z}N)}}e^{ik_{z}(l+i_{z}N)d}|{\bf R}_{X}^{(l+i_{z}N)}\rangle$		(4)

in commensurate CTGs. Here $|{\bf R}_{X}^{(l)}\rangle$ is the atomic orbital at the site ${\bf R}_{X}^{(l)}$, $\tilde{N}=S_{\mathrm{tot}}/S$ is the number of the graphene unit cells with an area $S=(\sqrt{3}/2)a^{2}$ in the total system area $S_{\mathrm{tot}}$, and $N_{z}$ and $k_{z}$ are the number of the primitive cells and the Bloch wave number along the vertical direction, respectively. The Brillouin zone along ${\bf k}_{z}$ of commensurate CTG with $N$ layers in a primitive cell is

$$-\frac{\pi}{Nd}\leq k_{z}<\frac{\pi}{Nd},$$

(5)

owing to the periodicity along the vertical direction [Eq. (2)]; this is analogous to the fact that the electronic states in a rhombohedral graphite are periodic with respect to $2\pi/(3d)$ shift of $k_{z}$ [22]. We use a two-center Slater-Koster parametrization [23, 24] for the transfer integral between any two $p_{z}$ orbitals,

$$-T({\bf R})=V_{pp\pi}\left[1-\left(\frac{{\bf R}\cdot{\bf e}_{z}}{|{\bf R}|}\right)^{2}\right]+V_{pp\sigma}\left(\frac{{\bf R}\cdot{\bf e}_{z}}{|{\bf R}|}\right)^{2},$$

(6)

where ${\bf R}$ is the relative vector between two atoms, and

	$\displaystyle V_{pp\pi}$		$\displaystyle=V_{pp\pi}^{0}e^{-(|{\bf R}|-a/\sqrt{3})/r_{0}},$
	$\displaystyle V_{pp\sigma}$		$\displaystyle=V_{pp\sigma}^{0}e^{-(|{\bf R}|-d)/r_{0}},$		(7)

$V_{pp\pi}^{0}\approx-3.38~{}\mathrm{eV}$ [25], $V_{pp\sigma}^{0}\approx 0.48~{}\mathrm{eV}$, and $r_{0}\approx 0.0453\,\mathrm{nm}$ [26, 6]. The total tight-binding Hamiltonian of CTG is expressed as

$$\mathcal{H}=\mathcal{H}_{\mathrm{G}}+\mathcal{U}$$

(8)

where $\mathcal{H}_{\mathrm{G}}$ and $\mathcal{U}$ represent the Hamiltonian for the intralayer and interlayer interaction, respectively. The intralayer interaction in each layer is given by

	$\displaystyle\langle{\bf k}^{\prime},X^{\prime(l)}|\mathcal{H}_{G}|{\bf k},X^{(l)}\rangle=h_{X^{\prime}X}^{(l)}({\bf k})\delta_{{\bf k}^{\prime},{\bf k}},$
	$\displaystyle h_{X^{\prime}X}^{(l)}({\bf k})=\sum_{{\bf L}^{(l)}}-T({\bf L}^{(l)}+\bm{\tau}_{X^{\prime}X}^{(l)})e^{-i{\bf k}\cdot({\bf L}^{(l)}+\bm{\tau}_{X^{\prime}X}^{(l)})},$		(9)

where ${\bf L}^{(l)}$ is the lattice vectors of the $l$-th layer and $\bm{\tau}_{X^{\prime}X}^{(l)}=\bm{\tau}_{X^{\prime}}^{(l)}-\bm{\tau}_{X}^{(l)}$. And the interlayer matrix element between the layer $l$ and $l^{\prime}$ is written as [2, 3, 4, 5, 6, 7]

	$\displaystyle\langle{\bf k}^{\prime},X^{\prime(l^{\prime})}|\mathcal{U}|{\bf k},X^{(l)}\rangle=u_{X^{\prime}X}^{(l^{\prime},l)}e^{ik_{z}(l-l^{\prime})d}\delta_{{\bf k}+{\bf G}^{(l)},{\bf k}^{\prime}+{\bf G}^{(l^{\prime})}},$
	$\displaystyle u_{X^{\prime}X}^{(l^{\prime},l)}=-\sum_{{\bf G}^{(l)}}\sum_{{\bf G}^{(l^{\prime})}}{t}({\bf k}+{\bf G}^{(l)})e^{-i{\bf G}^{(l)}\cdot\mbox{\boldmath\scriptsize$\tau$}_{X}^{(l)}+i{\bf G}^{(l^{\prime})}\cdot\mbox{\boldmath\scriptsize$\tau$}_{X^{\prime}}^{(l^{\prime})}},$		(10)

where ${\bf G}^{(l)}=m_{1}^{(l)}{\bf b}_{1}^{(l)}+m_{2}^{(l)}{\bf b}_{2}^{(l)}$ and ${\bf G}^{(l^{\prime})}=m_{1}^{(l^{\prime})}{\bf b}_{1}^{(l^{\prime})}+m_{2}^{(l^{\prime})}{\bf b}_{2}^{(l^{\prime})}$ ($m_{1}^{(l)},m_{2}^{(l)},m_{1}^{(l^{\prime})},m_{2}^{(l^{\prime})}\in\mathbb{Z}$) run over all the reciprocal points of layer $l$ and $l^{\prime}$, respectively. Here

$\displaystyle{t}({\bf q})=\frac{1}{S}\int T({\bf r}+z_{X^{\prime}X}^{(l^{\prime},l)}{\bf e}_{z})e^{-i{\bf q}\cdot{\bf r}}d{\bf r}$

(11)

is the in-plane Fourier transform of the transfer integral, where $z_{X^{\prime}X}^{(l^{\prime},l)}=(\bm{\tau}_{X^{\prime}}^{(l^{\prime})}-\bm{\tau}_{X}^{(l)})\cdot{\bf e}_{z}$. Since the interaction strength $T({\bf R})$ exponentially decays with the interatomic distance, the interlayer interaction is meaningful only between the adjacent layers ($|l^{\prime}-l|=1$). Likewise, $|t({\bf q})|$ also exponentially decays as $|{\bf q}|$ increases. Note that both $T(\textbf{R})$ and $t(\textbf{q})$ are isotropic along the in-plane direction, i.e., $T({\bf R})=T(|{\bf R}|)$ and $t({\bf q})=t(|{\bf q}|)$, if the two orbitals involved have the same magnetic quantum number, such as $p_{z}$ in this work [17].

Since CTG does not have an in-plane periodicity that is common to the entire system, one needs to find a general bases which does not rely on such periodicity. It is straightforward to show that the Hamiltonian $\mathcal{H}$ spans the subspace

$$\{|{\bf k},X^{(l)}\rangle~{}|~{}{\bf k}={\bf k}_{0}+\sum\limits_{l^{\prime}\in\mathcal{L}^{\backslash(l)}}\sum\limits_{{\bf G}^{(l^{\prime})}}{\bf G}^{(l^{\prime})}\},$$

(12)

for any ${\bf k}_{0}$ in the momentum space, where $\mathcal{L}$ is $\mathbb{Z}$ for incommensurate CTG and $\{0,1,...,N-1\}$ for commensurate CTG, and $\mathcal{L}^{\backslash(l)}=\{0,1,\dots,N-1\}\setminus{\{l\}}$. If we take each Bloch state as a “site”, the whole subspace can be recognized as a tight-binding lattice in the momentum space, which is the dual counterpart of the original tight-binding Hamiltonian in the real space [15]. In this momentum-space tight-binding model, the hopping between different sites (the interlayer interaction $\mathcal{U}$) of van der Waals multilayers is an order of magnitude smaller than the potential landscape (the band energies of the monolayers). Thus, in a similar manner to the Aubry-André model in one dimensional real-space lattice under an incommensurate perturbation [27], the eigenfunctions in our model tend to be localized to a few sites in momentum space. The analysis on the degree of the localization in Ref. [15] shows that most states in van der Waals bilayers are made up of an interaction of 20 or fewer (in most cases, just two or three) monolayer states. For most ${\bf k}_{0}$ in CTG, the length of such an interaction chain does not scale with the number of the constituent atomic layers owing to the mismatch between either the momentum or the monolayer state energies. Thus, although the size of the subspace [Eq. (12)] increases drastically with the size of $\mathcal{L}$, we only need a limited number of bases spanned from $|{\bf k}_{0},X^{(l)}\rangle$ of an arbitrary layer, chosen by applying suitable cut-off to both $|t({\bf q})|$ and the energy difference between the two monolayer states $\Delta E$, to describe the electronic structures near ${\bf k}_{0}$ for any practical calculation. This is one of the largest merits of using the momentum-space model, compared to the real-space model which requires infinitely many atomic orbital bases to represent incommensurate systems. We can, then, obtain the quasiband dispersion of the system by plotting the energy levels against ${\bf k}_{0}$. Here the wave number ${\bf k}_{0}$ works like the crystal momentum for the periodic system, and so it can be called the quasi-momentum for the current structures.

Figures 3(a) and (b) show the band dispersion at $k_{z}=0$ in the extended Brillouin zone of the commensurate CTG with $N=2$, calculated by the momentum-space tight-binding model; (a) shows the dispersion along the line $K^{(0)}-\bar{M}-K^{(1)}$, where $K^{(l)}$ is the Dirac point of $l$-th layer and $\bar{M}$ is the point in the middle of the two Dirac points, while (b) shows the dispersion along the line passes through $\bar{M}$ and perpendicular to $K^{(0)}-K^{(1)}$. Since CTG does not have an in-plane periodicity, it has an infinitesimal distinct Brillouin zone, and the bands are replicated incommensurately into the extended Brillouin zone (thin gray lines). We can reveal the distinct, unfolded band dispersion at wave vector ${\bf k}_{0}$ and energy $\varepsilon$ by the spectral function which is defined as

$$A_{l}({\bf k}_{0},\varepsilon)=\sum_{\alpha,X}|\langle\alpha|{\bf k}_{0},X^{(l)}\rangle|^{2}\delta(\varepsilon-\varepsilon_{\alpha}),$$

(13)

where $|\alpha\rangle$ and $\varepsilon_{\alpha}$ are the eigenstate and the eigenenergy, respectively. The red and blue lines in Figs. 3(a) and (b) show $A_{l}$ for $l=0,1$, respectively. The left panel shows that the Dirac cones at $K^{(0)}$ and $K^{(1)}$ interact with each other at $\bar{M}$ resulting in the complex structure which is not observed in the usual moiré superlattices with in-plane periodicity [6]. The right panel shows that such a rich structure stems from the multiple bands at ${\bf Q}_{4}$ and ${\bf Q}_{12}$. In Sec. III.4, we will show that the band dispersion at ${\bf Q}_{4}$ and ${\bf Q}_{12}$ represents the resonant states with 4- and 12-fold screw rotational symmetry, respectively.

Figure 3(c) shows the band dispersion at $k_{z}=0$ of the commensurate CTG with $N=3$ along the line shown in the left inset; this is the line $K^{(0)}-\bar{\Gamma}-\bar{M}-K^{(-1)}$ in terms of the high symmetry points of the moiré Brillouin zone between the layers with $l=0$ and $-1$. The red, green, blue lines show the $A_{l}$ for $l=0,1,-1$, respectively, and the right inset shows their mixing. Again, although the overall structure looks like the band dispersion of usual moiré superlattices [6], due to the momentum mismatch at most wave vectors in quasicrystalline systems [9, 15], we can see rich structures at ${\bf Q}_{9}$ and ${\bf Q}_{18}$, which arise from the resonant states with 9- and 18-fold screw rotational symmetry, respectively (see Sec. III.4).

In Sec. II.3, we showed that the interlayer interaction at specific ${\bf k}_{0}$ makes a chain of interaction which strongly couples the degenerate monolayer states of every constituent layers [Fig. 1(f)]. The electronic structures of such states are predominantly described by the interaction between these degenerate states. To reveal such resonant interactions in incommensurate CTG, we need to find a set of $\{{\bf C}_{j}\}$ ($j\in\mathbb{Z}$), where (i) all $|{\bf C}_{j},X^{(j)}\rangle$ are degenerate and (ii) $|{\bf C}_{j},X^{(j)}\rangle$ in layer $j$ and $|{\bf C}_{j+1},X^{(j+1)}\rangle$ in layer $j$+1 interact with each other by $\mathcal{U}$. Then, due to the screw rotational symmetry of the system, all $|{\bf C}_{j},X^{(j)}\rangle$ inevitably form an infinite chain of resonant interaction. To satisfy (i) regardless of the band distortion, such as trigonal warping, the relative position of ${\bf C}_{j}$ to the Dirac point of the $j$-th layer ($K^{(j)}$, either $K$ or $K^{\prime}$) must be the same in all the layers. Since there are six Dirac points in each hexagonal Brillouin zone, this condition requires ${\bf C}_{j+1}-K^{(j+1)}=\mathcal{R}(\phi)({\bf C}_{j}-K^{(j)})$, where $\mathcal{R}(\phi)$ is a counterclockwise rotation by

$$\phi=\theta+60^{\circ}r\quad(r=0,1,\dots,5).$$

(14)

This can be further reduced to

$${\bf C}_{j+1}=\mathcal{R}(\phi){\bf C}_{j},$$

(15)

since $K^{(j+1)}=\mathcal{R}(\phi)K^{(j)}$. For (ii), ${\bf C}_{j}$ and ${\bf C}_{j+1}$ should satisfy the generalized Umklapp scattering condition [Eq. (10)], i.e.,

$${\bf C}_{j+1}+\tilde{{\bf G}}^{(j+1)}={\bf C}_{j}+{\bf G}^{(j)},$$

(16)

for the reciprocal lattice vectors of the $j$-th layer ${\bf G}^{(j)}$ and the ($j$+1)-th layer $\tilde{{\bf G}}^{(j+1)}$, respectively. Without loss of generality, however, Eq. (16) can be reduced to

$${\bf C}_{j+1}={\bf C}_{j}+{\bf G}^{(j)}$$

(17)

(see Appendix A). Then, $|{\bf C}_{j},X^{(j)}\rangle$ and $|{\bf C}_{j+1},X^{(j+1)}\rangle$ interact with a magnitude of $|t({\bf q})|$, where ${\bf q}={\bf C}_{j}+{\bf G}^{(j)}$ (= ${\bf C}_{j+1}$). We can obtain the set of $\{{\bf C}_{j}\}$ which satisfies (i) and (ii) by first finding ${\bf C}_{0}$ that satisfies Eqs. (15) and (17),

	$\displaystyle{\bf C}_{0}$	$\displaystyle=[\mathcal{R}(\phi)-I_{2}]^{-1}{\bf G}^{(0)}$
		$\displaystyle=\frac{-1}{2\sin(\phi/2)}\mathcal{R}(90^{\circ}-\frac{\phi}{2}){\bf G}^{(0)}$		(18)

where $I_{2}$ is a $2\times 2$ identity matrix, for a reciprocal lattice vector of the 0-th layer ${\bf G}^{(0)}$. Then, due to the geometry, ${\bf G}^{(j)}=\mathcal{R}(j\phi){\bf G}^{(0)}$ for any $j$ is always a reciprocal lattice vector of the $j$-th layer which satisfies Eq. (17), and all the states $|{\bf C}_{j},X^{(j)}\rangle$ with ${\bf C}_{j}=\mathcal{R}(j\phi){\bf C}_{0}$ are degenerate. Since all of these interactions have the same magnitude of coupling $t_{0}=t(|{\bf q}|)$, where

$$|{\bf q}|=\frac{1}{2\sin(\phi/2)}|{\bf G}^{(0)}|,$$

(19)

as long as $t({\bf q})$ is isotropic (Sec. II.2), the set of the Bloch states $\{|{\bf C}_{j},X^{(j)}\rangle\}$ forms resonant interaction.

CTG becomes commensurate, i.e., gains periodicity along the $z$ axis, at specific $\theta$. In addition to the resonant conditions (i) and (ii) for incommensurate CTG, the $\{{\bf C}_{j}\}$ of commensurate CTG has to satisfy a periodic condition; (iii) $|{\bf C}_{j},X^{(j)}\rangle$ and $|{\bf C}_{n+j},X^{(n+j)}\rangle$ with some $n\in\mathbb{Z}$ are equivalent up to the $k_{z}$ phase. This requires that the lattice configuration of the ($n$+$j$)-th layer is the same as that of the $j$-th layer, and also that ${\bf C}_{n+j}={\bf C}_{j}$. The former requires

$$n=uN\quad(u\in\mathbb{Z}),$$

(20)

while the latter is equivalent to

$$\sum_{j=0}^{n-1}{\bf G}^{(j)}=\left[\sum_{j=0}^{n-1}\mathcal{R}(j\phi)\right]{\bf G}^{(0)}=0,$$

(21)

which requires

$$\phi=360^{\circ}\frac{M^{\prime}}{n},$$

(22)

with $M^{\prime}\in\mathbb{Z}$. Without loss of generality, we choose the smallest positive $n$ which satisfy (i)-(iii), so that we describe the shortest periodic unit of the interaction loop. Then, we get

$$Nr+M=\frac{6}{u}M^{\prime}$$

(23)

($n,u,M^{\prime}\in\mathbb{Z}^{+}$) from Eqs. (2), (14), (20), and (22), and it is straightforward to show that $\gcd(u,M^{\prime})=1$ and $\gcd(N,M^{\prime})=1$ (see Appendix B). Accordingly, $\gcd(n,M^{\prime})=1$ and $u$ can only have a value in $\{1,2,3,6\}$.

As we will see later (Secs. III.3 and III.4), $\phi$ represents the angle of the screw rotational symmetry of each resonant state, and $n$ for commensurate CTG represents the discrete screw rotational symmetry of the state. In commensurate CTG, we can implement any $n$, except 1, 2, 3, 6, with a suitable choice of the geometry ($N$, $M$) and $r$.

Above equations indicate that there are infinitely many resonant states in each CTG. In a given geometry $\theta$, different $r$ give the resonant states with different $\phi$, and in commensurate CTG, the states with different $r$ may have different screw rotational symmetry $n$ [Eqs. (14) and (23)]. We plot the five strongest resonant interactions in CTGs with $N=2,3,5$ at the top, middle, and bottom panels in Fig. 4. We can see that a CTG with $(N,M)=(2,1)$ ($\theta=30^{\circ}$), of which lattice configuration has a 12-fold screw rotational symmetry (top panels), can host not only the resonant states with a 12-fold rotational symmetry ($r=0,2,3,5$) but also the states with a 4-fold rotational symmetry ($r=1,4$). Besides, the system with $N=3$ ($N=4$) can host the resonant states with $n=9,18$ ($n=8,24$), while that with $N=5$ (bottom panels) can have the states with $n=5,10,15,30$, and so on. In addition, even for a fixed $r$, different ${\bf G}^{(0)}$ make the resonant states with distinct $t({\bf q})$ which appear at different set of wave vectors $\{{\bf C}_{j}\}$ in the Brillouin zone [Eqs. (18) and (19)]. Despite the infinitely many resonant states, however, the number of distinct screw rotational symmetry $n$ in each commensurate CTG is finite; this number is determined by the geometry $(N,M)$ and can be up to 4, each of which corresponds to $u=1,2,3,6$, at most.

Since $t({\bf q})$ decays exponentially as $|{\bf q}|$ increases, the interaction with the shorter $|{\bf G}^{(0)}|$ exhibits the stronger interaction between the monolayer Bloch states [Eqs. (18) and (19)]. We, however, do not consider the interaction with ${\bf G}^{(0)}={\bf 0}$ in this work, as such interaction merely represents the interaction between the monolayer states at $\Gamma$ of every layers which is common to multilayer systems with any $\theta$, and the coupled states appear at the band edges of pristine graphene, which are very far from the charge neutrality point. Then, the strongest resonant interaction occurs for the next shortest $|{\bf G}^{(0)}|$, i.e., $|{\bf G}^{(0)}|=|{\bf b}_{i}|$, and $r=3$ (if $\theta=30^{\circ}$, both $r=2,3$ give the strongest interaction ). And without loss of generality, we choose ${\bf G}^{(0)}=-{\bf b}_{2}^{(0)}$; a different choice of ${\bf G}^{(0)}$ having the same $|{\bf G}^{(0)}|=|{\bf b}_{i}|$ will only shift ${\bf C}_{j_{1}}$ to ${\bf C}_{j_{2}}$ in the same set of $\{{\bf C}_{j}\}$, and can be obtained by a similarity transformation.

The monolayer states $|{\bf k}^{(j)},X^{(j)}\rangle$ (${\bf k}^{(j)}={\bf C}_{j}+{\bf k}_{0}$) in each layer $j$ form an infinitely long chain of resonant interaction. Although the interlayer interaction [Eq. (10)] couples these states also to other states, the energies of such states usually are very different from the energy of $|{\bf k}^{(j)},X^{(j)}\rangle$, except for accidental degeneracy. Thus, we can get the details of the resonant states by using these degenerated monolayer states as bases, rather than using the full bases model described in Sec. II.2, since the interaction to non-degenerate states does not break the symmetry and degeneracy of the involved monolayer states and merely shifts the resonant states by some, mostly small, constant energies (see Appendix C). Note that there are more resonant chains which are associated with $|{\bf k}^{(cN)},X^{(0)}\rangle$ ($c\in\mathbb{Z}$) [gray lines in Fig. 5(a)], in addition to the chain from $|{\bf k}^{(0)},X^{(0)}\rangle$ [black line in Fig. 5(a)]. These chains are associated with the ${\bf G}^{(0)}$ different from the ${\bf G}^{(0)}$ of $|{\bf k}^{(0)},X^{(0)}\rangle$ but having the same $|{\bf G}^{(0)}|$. Each interaction chain makes the same band dispersion at different ${\bf k}_{0}$.

To represent the Hamiltonian of the resonant states in this minimal bases model, it is convenient to use new coordinate vectors which are defined by the rotation by $\phi$, i.e.,

	$\displaystyle{\bf a}_{i}^{(j)}$	$\displaystyle=\mathcal{R}(j\phi){\bf a}_{i},$
	$\displaystyle\bm{\tau}_{X}^{(j)}$	$\displaystyle=\mathcal{R}(j\phi)\bm{\tau}_{X},$
	$\displaystyle{\bf b}_{i}^{(j)}$	$\displaystyle=\mathcal{R}(j\phi){\bf b}_{i},$		(24)

($i=1,2$, $X=A,B$), instead of the vectors defined by the rotation by $\theta$ in Sec. II.1. This is a unitary transformation which makes the Hamiltonian a highly symmetric form regardless of the distinct geometry and rotational symmetry of the resonant states [$r$ in Eq. (14)]. Note that $|{\bf k}^{(j)},X^{(j)}\rangle$ and $|{\bf k}^{(j+cN)},X^{(j+cN)}\rangle$ ($c\in\mathbb{Z}$) represent the Bloch states in the atomic layers with the same orientation, but they are expressed by lattice vectors rotated by $60^{\circ}crN$. And the resonant states with different $r$ in the same lattice configuration are expressed by coordinate vectors with different orientations. In this choice of the coordinates, all ${\bf G}^{(j)}=m_{1}^{(j)}{\bf b}_{1}^{(j)}+m_{2}^{(j)}{\bf b}_{2}^{(j)}$ are represented by the same coefficients $m_{1}^{(j)}=m_{1}$, $m_{2}^{(j)}=m_{2}$. Likewise, all ${\bf C}_{j}=c_{1}^{(j)}{\bf b}_{1}^{(j)}+c_{2}^{(j)}{\bf b}_{2}^{(j)}$ are represented by the same coefficients $c_{1}^{(j)}=c_{1}$, $c_{2}^{(j)}=c_{2}$, where

$$\begin{pmatrix}c_{1}\\ c_{2}\end{pmatrix}=\left[\tilde{\mathcal{R}}(\phi)-I\right]^{-1}\begin{pmatrix}m_{1}\\ m_{2}\end{pmatrix}.$$

(25)

Here,

$$\tilde{\mathcal{R}}(\phi)=\frac{2}{\sqrt{3}}\begin{pmatrix}\cos(\phi-30^{\circ})&-\sin\phi\\ \sin\phi&\cos(\phi+30^{\circ})\end{pmatrix}$$

(26)

is a counterclockwise rotation by $\phi$ applied to the coefficients of ${\bf b}_{1}$ and ${\bf b}_{2}$.

Then the Hamiltonian of any resonant interaction near ${\bf k}_{0}=0$ of CTG can be expressed by the Hamiltonian of one-dimensional chain,

		$\displaystyle\hat{\mathcal{H}}_{\mathrm{ring}}({\bf k}_{0})=\sum_{j\in\mathbb{Z}}\left\{\psi_{j}^{\dagger}H^{(j)}\psi_{j}+[\psi_{j+1}^{\dagger}W\psi_{j}+\mathrm{H.c.}]\right\},$
		$\displaystyle H^{(j)}({\bf k}_{0})=\begin{pmatrix}h_{AA}^{(j)}&h_{AB}^{(j)}\\ h_{BA}^{(j)}&h_{BB}^{(j)}\\ \end{pmatrix},$
		$\displaystyle W=\left\{\begin{aligned} &-t_{0}\begin{pmatrix}\tilde{\omega}^{*}&\tilde{\omega}\\ \tilde{\omega}^{*}&\tilde{\omega}\end{pmatrix},&&\text{incommensurate CTGs}\\ &-t_{0}\begin{pmatrix}\tilde{\omega}^{*}&\tilde{\omega}\\ \tilde{\omega}^{*}&\tilde{\omega}\end{pmatrix}e^{-ik_{z}d},&&\text{commensurate CTGs}\end{aligned}\right.$		(27)

where $\psi_{j}^{\dagger}$ ($\psi_{j}$) is a creation (annihilation) operator for $(|{\bf k}^{(j)},A^{(j)}\rangle,|{\bf k}^{(j)},B^{(j)}\rangle)$, $h_{X^{\prime}X}^{(j)}=h_{X^{\prime}X}^{(0)}[\mathcal{R}(-j\phi){\bf k}_{0}+{\bf C}_{0}]$ in the new coordinate system, $t_{0}=t({\bf C}_{0})$, where we neglect the ${\bf k}_{0}$ dependence of the interlayer matrix element $t({\bf q})$, and $\tilde{\omega}=\omega^{m_{1}-2m_{2}}$ ($\omega=e^{2\pi i/3}$). Since $\hat{\mathcal{H}}_{\mathrm{ring}}$ is obviously symmetric under the substitution of $j$ by $j+1$, the resonant states are symmetric under the screw rotation by $\phi$, i.e., in-plane rotation by $\phi$ followed by a translation along ${\bf e}_{z}$ by $d$, regardless of the value of $k_{z}$.

In commensurate CTG, the resonant chain has a periodic unit of $\{{\bf C}_{j}\}$ ($j=0,1,\dots,n-1$). Then the Hamiltonian of resonant states near ${\bf k}_{0}=0$ can be expressed by an one-dimensional ring model,

	$\displaystyle\mathcal{H}_{\rm ring}({\bf k}_{0})$
	$\displaystyle=\begin{pmatrix}H^{(0)}&W^{\dagger}&&&&W\\ W&H^{(1)}&W^{\dagger}\\ &W&H^{(2)}&W^{\dagger}\\ &&\ddots&\ddots&\ddots\\ &&&W&H^{(n-2)}&W^{\dagger}\\ W^{\dagger}&&&&W&H^{(n-1)}\end{pmatrix},$		(28)

in the bases of $(|{\bf k}^{(j)},A^{(j)}\rangle,|{\bf k}^{(j)},B^{(j)}\rangle)$, $(j=0,1,\dots,n-1)$. Again, since $\tilde{\mathcal{H}}_{\mathrm{ring}}$ is symmetric under a rotation by a single span of the ring (i.e., moving ${\bf C}_{j}$ to ${\bf C}_{j+1}$), the resonant states is symmetric with respect to the $n$-fold screw rotation regardless of the value of $k_{z}$.

It should be emphasized that the resonant states in commensurate CTG are $u$ ($=n/N$)-fold degenerate along $k_{z}$ whereas the other, non-resonant states in the same system are not. This becomes clear by a similarity transformation of $\mathcal{H}_{\mathrm{ring}}$

	$\displaystyle\tilde{\mathcal{H}}_{\mathrm{ring}}({\bf k}_{0})=U^{-1}\mathcal{H}_{\mathrm{ring}}({\bf k}_{0})U$
	$\displaystyle=\begin{pmatrix}H^{(0)}&\tilde{W}^{\dagger}&&&&\tilde{W}e^{-ink_{z}d}\\ \tilde{W}&H^{(1)}&\tilde{W}^{\dagger}\\ &\tilde{W}&H^{(2)}&\tilde{W}^{\dagger}\\ &&\ddots&\ddots&\ddots\\ &&&\tilde{W}&H^{(n-2)}&\tilde{W}^{\dagger}\\ \tilde{W}^{\dagger}e^{ink_{z}d}&&&&\tilde{W}&H^{(n-1)}\end{pmatrix},$		(29)
	$\displaystyle\tilde{W}=We^{ik_{z}d}=-t_{0}\begin{pmatrix}\tilde{\omega}^{*}&\tilde{\omega}\\ \tilde{\omega}^{*}&\tilde{\omega}\end{pmatrix},$		(30)

with a transformation matrix

$$U=\mathrm{diag}(1,e^{-ik_{z}d},e^{-2ik_{z}d},\dots,e^{-i(n-1)k_{z}d})\otimes I_{2}.$$

(31)

Equation (29) is periodic with $k_{z}\rightarrow k_{z}+2\pi/(nd)$, so the resonant states with a $n$-fold rotational symmetry are $u$-fold degenerate along the $k_{z}$ direction. Such a repetition unit in the Brillouin zone

$$-\frac{\pi}{nd}\leq k_{z}<\frac{\pi}{nd}$$

(32)

which is $u$ times smaller than that of the general, non-resonant states in CTG [Eq. (5)] can be regarded as a reduced Brillouin zone. Moreover, the infinitely many resonant states in the same system (Sec. III.2) also have different reduced Brillouin zone size if they have different $u$.

The CTG with $N=2$ is, however, special in that the monolayer state $|{\bf C}_{j},X^{(j)}\rangle$ in layer $j$ not only interacts with $|{\bf C}_{j+1},X^{(j+1)}\rangle$ in layer $j+1$, like the CTGs with other $N$ [Fig. 5(a)], but also interacts with $|{\bf C}_{j+1},X^{(j-1)}\rangle$ in layer $j-1$ with the same momentum difference [Fig. 5(b)]. Thus, $W$ should be replaced by $W(1+e^{2ik_{z}d})$ for $N=2$, and such an extra phase changes the periodicity along ${\bf k}_{z}$ of the resonant states from $2\pi/(nd)$ to $\mathrm{lcm}(2\pi/(nd),\pi/d)$ [Eq. (29)]. Since only $n=4,12$ are allowed in $N=2$, the size of the Brillouin zone along ${\bf k}_{z}$ of the resonant states in $N=2$ is the same as the Brillouin zone of the other, non-resonant states, $\pi/d$. Except for the extra phase, $2\cos k_{z}d$, in the interlayer matrix elements, $\mathcal{H}_{\mathrm{ring}}$ of $N=2$ is the same as that of the twisted bilayer graphene quasicrystal in Ref. [15]. Thus, the resonant states of the CTG with $N=2$ are natural generalization of the quasicrystalline states of the twisted bilayer graphene quasicrystal to the three-dimensional structures. In CTG, however, the wave functions of the upper and lower layers [Fig. 5(b)] interfere constructively at $k_{z}=c\pi/d$ ($c\in\mathbb{Z}$) and double the magnitude of the interlayer interaction $W$ of $N=2$ compared to that of the twisted bilayer graphene quasicrystal. At $k_{z}=\pi/(2d)+c\pi/d$, on the other hand, the waves interfere destructively and $W$ completely vanishes. As a result, regardless of the matching of momentum and energies, $|{\bf C}_{j},X^{(j)}\rangle$ at these $k_{z}$ are completely decoupled from each other, and show the energies of the decoupled monolayer states. Such a perfect constructive and destructive interference occur only at $N=2$.

We plot the band structure of the resonant states with 12- and 4-fold screw rotational symmetry in commensurate CTG with $N=2$ against ${\bf k}_{0}$ in Figs. 6(a) and (b), respectively, and the resonant states with 9- and 18-fold screw rotational symmetries in commensurate CTG with $N=3$ in Fig. 6(c) and (d), respectively. The dispersion in these plots are consistent with the rich structure in Fig. 3. As shown in the previous section, $N=2$ gives the dispersion the same as that of the quasicrystalline twisted bilayer graphene, but the interlayer interaction is doubled at $k_{z}=c\pi/d$ ($c\in\mathbb{Z}$) and vanishes in the middle between them.

At ${\bf k}_{0}={\bf 0}$, we can analytically obtain the energies of the resonant states (neglecting the constant energy)

$$E_{m}=-2t_{0}\cos q_{m}\cos q_{G}\pm\sqrt{(2t_{0}\sin q_{m}\sin q_{G})^{2}+|h_{0}-2t_{0}\tilde{\omega}\cos q_{m}|^{2}}$$

(33)

for incommensurate CTGs,

$$E_{m}(k_{z})=-2t_{0}\cos(q_{z}+q_{m})\cos q_{G}\pm\sqrt{(2t_{0}\sin(q_{z}+q_{m})\sin q_{G})^{2}+|h_{0}-2t_{0}\tilde{\omega}\cos(q_{z}+q_{m})|^{2}}$$

(34)

for commensurate CTGs with $N\neq 2$, and

$$E_{m}(k_{z})=-4t_{0}\cos q_{z}\cos q_{m}\cos q_{G}\pm\sqrt{(4t_{0}\cos q_{z}\sin q_{m}\sin q_{G})^{2}+|h_{0}-4t_{0}\tilde{\omega}\cos q_{z}\cos q_{m}|^{2}}$$

(35)

for commensurate CTGs with $N=2$, as well as the eigenvectors, i.e., the coefficients to the Bloch bases,

	$\displaystyle{\bf v}_{m}^{\pm}$	$\displaystyle=(1/\sqrt{N^{\prime}})(\dots,e^{i(j-1)q_{m}},e^{ijq_{m}},e^{i(j+1)q_{m}},\dots)^{\mathrm{T}}$
		$\displaystyle\otimes(c_{m,1}^{\pm},c_{m,2}^{\pm})^{\mathrm{T}}$		(36)

for incommensurate CTG and

	$\displaystyle{\bf v}_{m}^{\pm}$	$\displaystyle=(1/\sqrt{n})(1,e^{iq_{m}},e^{2iq_{m}},\dots,e^{i(n-1)q_{m}})^{\mathrm{T}}$
		$\displaystyle\otimes(c_{m,1}^{\pm},c_{m,2}^{\pm})^{\mathrm{T}}$		(37)

for commensurate CTG, by analytically diagonalizing $\tilde{\mathcal{H}}_{\mathrm{ring}}$. Here $q_{z}=k_{z}d$, $q_{G}=(2\pi/3)(m_{1}-2m_{2})$, $h_{0}=h_{AB}^{(0)}({\bf C}_{0})$, $\pm$ corresponds to the conduction band and valence band, respectively, $N^{\prime}$ is the number of layers, and see Appendix D for the expression of $c_{m,1}^{\pm}$ and $c_{m,2}^{\pm}$. And $q_{m}=\phi m$ ($m\in\mathbb{Z}$ for incommensurate CTG and $m\in\{0,1,\dots,n-1\}$ for commensurate CTG) is a Bloch number of the state which corresponds to the screw rotation by $\phi$. Note that Eqs. (33)-(37) are universal to any CTG and any resonant states within, although the parameters vary. Each resonant interaction gives the eigenstates twice as many as the number of distinct $q_{m}$, i.e., infinite in incommensurate CTG and $2n$ in commensurate CTG.

Equation (35) clearly shows that the resonant states at ${\bf Q}_{12}$ of the CTG with $N=2$ show a larger energy spacing in the valence band than the conduction band, just like the well-known resonant states in quasicrystalline twisted bilayer graphene, while those at ${\bf Q}_{4}$ show larger spacing in the conduction band [Fig. 3(b)]. We plot $E_{m}$ at ${\bf Q}_{12}$ and ${\bf Q}_{4}$ of a CTG with $N=2$ against

$$\tilde{q}_{m}=\frac{2\pi}{n}m$$

(38)

($m\in\{0,1,\dots,n-1\}$) in Figs. 7(a) and (b), respectively. The red and blue lines show the dispersion for $k_{z}=0$ and $k_{z}=\pi/(2d)$, respectively. The set $\{\tilde{q}_{m}\}$ is congruent to the set $\{q_{m}\}$, where

$$q_{m}=\frac{2\pi M^{\prime}}{n}m\ (=\phi m),$$

(39)

modulo $2\pi$, due to $\gcd(n,M^{\prime})=1$. In incommensurate CTGs, likewise, the set $\{\tilde{q}_{m}\}$ [$\tilde{q}_{m}\equiv q_{m}$ (mod $2\pi$)] is congruent to $\{q_{m}\}$, since $\phi$ is incommensurate with $2\pi$. As predicted in Sec. III.3, a CTG with $N=2$ experiences a destructive interlayer interference at $k_{z}=\pi/(2d)+c\pi/d$ ($c\in\mathbb{Z}$), so $E_{m}$ at such $k_{z}$ are simply the monolayer energy at ${\bf C}_{0}$ regardless of $\tilde{q}_{m}$.