Moiré fractals in twisted graphene layers

Twisted bilayer graphene (TBLG) subject to a sequence of commensurate external periodic potentials reveals the formation of moiré fractals (MFs) that share striking similarities with the central place theory of economic geography, thus uncovering a remarkable connection between twistronics and the geometry of economic zones. MFs arise from the self-similarity of the emergent hierarchy of Brillouin zones (BZs), forming a nested subband structure within the bandwidth of the original moiré bands. We derive the fractal generators for TBLG under these external potentials and explore their impact on the hierarchy of the BZ edges and the wave functions at the Dirac point. By examining realistic supermoiré structures and demonstrating their equivalence to MFs with periodic perturbations under specific conditions, we establish MFs as a general description for such systems. Furthermore, we uncover parallels between the modification of the BZ hierarchy and magnetic BZ formation in Hofstadter's butterfly, allowing us to construct an incommensurability measure for MFs versus twist angle. The resulting band structure hierarchy bolsters correlation effects, pushing more bands within the same energy window for both commensurate and incommensurate TBLG.

Figure 1

(a) $L_N,\beta$, and $x$ over triangular coordinates (the axes are ${60}^{\circ }$ with respect to each other). The intersection of $\beta$ (magenta) and $x$ (green) lines identify the intersection points ($x,y=x+\beta$) at which $L_N$ exists. (b) The real-space creation of each iteration: TBLG is created by stacking graphene layers with the top layer being the zeroth iteration $j=0$ and the bottom layer being the next $j=1$. The iterations $j=2,3,\cdots$ are created by applying the $z$-independent external periodic potentials identically to both the graphene layers. (c) The $K=3$ CPT hierarchy (since a hexagon of each layer encloses three hexagons of an adjacent layer) where each layer corresponds to a particular order of settlement [6]. $L_N$ is analogous to $K$.

Figure 2

(a1), (b1) The BZ for $j=1$ (black), 2 (magenta), and 3 (green), where $\{{\mathit{b}}_1^{\left(1\right)},{\mathit{b}}_2^{\left(1\right)}\}$ are the reciprocal lattice primitive vectors for $j=1$. (a2), (b2) The fractal structures at the BZ edges. The dashed red hexagon represents the initiator, i.e., the BZ of SLG rotated by $\theta /2$. The solid blue lines outline the fractal structures' outer boundary. The generators (inside) attach alternately to the initiator, forming the fractal structure. The copies of the BZ at each iteration are added such that they overlap with the BZ of SLG [red solid lines in (a1) and (b1)] and these overlaps lead to IFs. The insets between (a1), (a2) and (b1), (b2) display reciprocal and real-space lattice vectors for $q=3,p=1$ and $q=2,p=1$, respectively. The shift between the Dirac points is equal to the hexagon's side length for $q=3,p=1$, and twice its side for $q=2,p=1$. (a3), (b3) The band structures for $V_0=1.2\mathrm{meV}$ and the DOS (with a Gaussian smearing of $0.002\mathrm{meV}$). (a4), (b4) ${\rho }_{n\mathit{k}}\left(\mathit{r}\right)$ of the lowest conduction band (dashed-dotted line) at the Dirac point (see text, Appendix pp8).

Figure 3

$\frac{\mathrm{\Delta }A}{A_{\text{FBZ}}^{(j-1)}}$ vs $\theta$: The blue rectangles represent commensurate (closest) angles for odd $q$ and $p=1$. The first BZ area for incommensurate angles is determined using moiré vectors [96]. The inset confirms the absence of the apparent multivaluedness in the lowest values of the V regions between the dashed lines, with a separation of $\sim 0.{002}^{\circ }$.

Figure 4

The band structures at ${\theta }_r\sim 21.{79}^{\circ }$ in (a) and ${\theta }_r\sim 13.{17}^{\circ }$ in (b) along $K^{\prime }-M-\mathrm{\Gamma }-K$. The twist between the layers is the first magic angle $\theta \sim 1.{05}^{\circ }$. The interlayer hopping parameters: $t_{\text{AA/BB}}=t_{AB/BA}=110\mathrm{meV}$ (Appendix pp6) and $V_0=1.2\mathrm{meV}$. Left to right: The emergent fractal structure in reciprocal space, 2D band structure: orange (blue) bands without (with) mEP, the DOS, and the spatial profiles of ${\rho }_K\left(\mathit{r}\right)$ of the lowest conduction band with $j=2$. Changing ${\theta }_r$ from (b) to (a) modifies $D_f$, shifting more bands towards $E_F$.

Figure 5

(a) The stepwise creation of the FG for $L_N=7$ with $q=3$ and $p=1$, and (b) $L_N=13$ with $q=2$ and $p=1$. The red dashed line is the side $\mathit{u}$ of $A_0$ as defined in Appendix pp1.

Figure 6

(a) The fractal corresponding to one more entry in Table 1, where $q=11$ and $p=3$, having $L_N=31$ with $\beta =4$ at an angle of ${\theta }_r=17.{89}^{\circ }$ for $j=1,2$. (b) The fractal corresponding to the iteration $j=3$.

Figure 7

(a) The real-space structure of commensurate TBLG at $\theta \sim 21.{79}^{\circ }$ and two levels of periodic potentials for $j=2,3$. The direct-lattice primitive vectors are also shown for each level. (b) The corresponding moiré fractal along the edges of the Wigner-Seitz unit cell at $j=3$. (c) The initiator, the arm of the hexagon corresponding to $j=3$, and the FG.

Figure 8

$p_1$ vs $p_2$ for two maximum values of $q$. Evidently, the points arrange themselves in hexagons and each point is associated with a Löschian number $L_N$.

Figure 9

The inclusion of corrugation effects leads to different interlayer hopping parameters, namely, $t_{\text{AA,BB}}=79.7\mathrm{meV}$ and $t_{\text{AB,BA}}=97.5\mathrm{meV}$ [96]. This does not change the $2e^{2ln\left(n_c\right)/D_f}-2$ bands inserted within the band gap of the two lowest bands at the $\mathrm{\Gamma }$ point (shown by the two arrows) in the absence of mEP and hence the corrugated TBLG also does not affect the emergence of fractality.

Figure 10

(a) The renormalized Fermi velocity $v_F^{*}/v_F$ as a function of ${\alpha }^2$ and hence the twist angle $\theta$. (b) The Fermi energy ($E_F$) as a function of ${\alpha }^2$. Similar to the case of a pristine TBLG [11] or an unrotated mEP to TBLG [83], the renormalized Fermi velocity in the presence of a rotated mEP remains unaffected both in the absence and the presence of the corrugation effect.

Figure 11

(a) Trilayer-, (b) four-layer-twisted graphene, and (c) a trilayer boron nitride-graphene-boron nitride system. The above $•$ denote carbon, $•$ boron and $•$ nitrogen atoms, respectively. The angle ${\theta }_r$ is one of the commensurate angles and $d_0\sim 0.335\mathrm{nm}$ is the interlayer separation in TBLG. The number beside each layer represents the layer index used in the text in Appendix pp8.

Figure 12

(a) The spatial variation of the inversion symmetry preserving part $V_{{11}^{\prime }}^{\text{eff}}\left(\mathit{r}\right)$, where the color bar is in the units of $t_{{11}^{\prime }}^2/V_{\text{STM}}$, and (b) the spatial variation of the mass-dependent term $M_{{11}^{\prime }}^{\text{eff}}\left(\mathit{r}\right)$ that breaks the inversion symmetry. Here, the color bar is in the units of $\sqrt{3}{t}_{{11}^{\prime }}^2/V_{\text{STM}}$.

Figure 13

The band structure corresponding to the trilayer graphene system with increasing values of the interlayer bias voltage (a) $V_{\text{STM}}=10\mathrm{meV}$ and (b) $V_{\text{STM}}=40\mathrm{meV}$. The solid orange lines show the band structure in the absence of the potential while the blue lines represent the eigenvalues of $H_2\left(\mathit{r}\right)$ without the perturbation, while the cyan lines show the band structure of $H_2\left(\mathit{r}\right)+H_{\text{pert}}\left(\mathit{r}\right)$ calculated up to the first order in perturbation. The dotted cyan lines show $4e^{2ln\left(n_c\right)/D_f}-2$ bands within the band gap of lowest two bands at the $\mathrm{\Gamma }$ point.

Figure 14

The band structure of the tetralayer graphene system as shown in Fig. 11 for $V_{\text{STM}}=20\mathrm{meV}$. The solid orange lines show the band structure in the absence of the potential while the blue lines represent the eigenvalues of $H_2\left(\mathit{r}\right)$ without the perturbation, while the cyan lines show the band structure of $H_2\left(\mathit{r}\right)+H_{\text{pert}}\left(\mathit{r}\right)$ calculated up to the first order in perturbation. The dotted cyan lines show $4e^{2ln\left(n_c\right)/D_f}-2$ bands within the band gap of lowest two bands at the $\mathrm{\Gamma }$ point.

Figure 15

(a) The two blue and dashed blue hexagons show the moiré BZ of the top and bottom moiré interfaces for the trilayer hBN-G-hBN shown in Fig. 11. The inner red hexagon is the super-moiré BZ. (b) The band structures for two different values of $q,p$ along the high-symmetry path $X-Y-K-X$ [73] as shown in (b). The solid blue lines show the band structures of the graphene-hBN system with $V_2\left(\mathit{r}\right)=0$ while the solid red lines represent graphene sandwiched between two hBN layers as in Fig. 11. The two dark-gray arrows show the band gap between two lowest bands at the $Y$ point where the $2e^{2ln\left(n_c\right)/D_f}-1$ bands are inserted.

Figure 16

The bandwidth of the conduction and valence band in the presence of mEP for $q=3,p=1$ for (a) $t_{AA,BB}=t_{AB,BA}=110\mathrm{meV}$ [11], (b) $t_{AA,BB}=79.7\mathrm{meV}$, and $t_{AB,BA}=97.5\mathrm{meV}$ [24, 96]. The bandwidth $t_w$ is in the unit of $t_{w0}$, which is the bandwidth of pristine TBLG at $\theta \sim 1.{05}^{\circ }$.