Moiré superlattice on the surface of a topological insulator

Twisting van der Waals heterostructures to induce correlated many body states provides a novel tuning mechanism in solid state physics. In this work, we theoretically investigate the fate of the surface Dirac cone of a three-dimensional topological insulator subject to a superlattice potential. Using a combination of diagrammatic perturbation theory, lattice model simulations, and ab initio calculations we elucidate the unique aspects of twisting a single Dirac cone with an induced moiré potential and the role of the bulk topology on the reconstructed surface band structure. We report a dramatic renormalization of the surface Dirac cone velocity as well as demonstrate a topological obstruction to the formation of isolated minibands. Due to the topological nature of the bulk, surface band gaps cannot open; instead additional satellite Dirac cones emerge, which can be highly anisotropic and made quite flat. We discuss the implications of our findings for future experiments.

Figure 1

Schematic depictions of experimental realizations of a surface-moiré potential. (a) Surface potential induced by gating a periodic patterned dielectric. (b) Twisting a gapped 2D material on the surface of the 3D TI.

Figure 2

Momentum exchanges on the surface Dirac cone, satellite Dirac cone labels, and density of states. (a) The lowest SDC ($X_1$) is a result of the superlattice potential $V_{{\mathbf{Q}}_x}$ coupling the states at momenta $\pm {\mathbf{Q}}_x/2$, shown by the white line connecting two black dots. (A SDC at the same energy occurs from coupling the states at $\pm {\mathbf{Q}}_y/2$, not shown.) On the square lattice, the next lowest energy SDC ($M_1$) results from the superlattice potential coupling the four degenerate states at $\pm {\mathbf{Q}}_x/2\pm {\mathbf{Q}}_y/2$, indicated by the white square connecting four black dots. (b) The resulting downfolded cone and satellites computed with Eq. (45) at $W=0.2t$. (c) The resulting density of states in arbitrary units. The cone ${\mathrm{\Gamma }}_1$ is visible and isolated from other features in the spectrum.

Figure 3

Surface state properties of the 3D TI lattice model in Eq. (37). (a) A depiction of the full surface Brillouin zone and the region where surface state solutions exist is marked in red for the model parameters considered here: $t=1,m_0=-1,m_2=1$. (b, left) The band structure depicting bulk states (gray) and the surface Dirac cones (red) along the high-symmetry cut shown in (a). (b, right) The corresponding surface density of states in the layer $z=L$ for periodic and open boundary conditions (PBCs and OBCs, respectively) along the $z$ direction.

Figure 4

Momentum exchange of the surface Dirac cone in the first and second BZ. (a) The gamma point (black circle) in the first and second Brillouin zone (BZ) is scattered by a momentum exchange $\mathbf{Q}$ to the red circle in the second BZ. (b) The scattered momentum point (red) is shifted back into the first BZ and its angular rotation about the $C_4$ axis is neglected, as described in the main text.

Figure 5

Surface dispersion and density of states due to the surface potential on the 3D TI lattice model in Eq. (37). This data is obtained on a cubic system of linear size $L=89$ with the moiré surface potential in Eq. (45) with a twist parameterized by the wave vector $Q=2\pi 8/89$ [see Eq. (43)]. We use Lanczos with twisted boundary conditions in the $x\text{−}y$ plane to compute the surface dispersions whereas the surface DOS is computed using KPM with an expansion order $N_C=2^{11}$. A selection of relevant SDCs are labeled by ${\mathrm{\Gamma }}_{0,1,2}, X_{1,2}$, and $M_{1,2}$ depending on where in the mini-Brillouin zone they appear. (a) At $W=0.1t$, the Dirac cone is only slightly perturbed. The dispersion reveals the formation of SDCs but the surface DOS is relatively unaffected. (b) At $W=0.2t$, the SDC $M_2$ moves below ${\mathrm{\Gamma }}_1$ in energy, revealing it in the surface DOS. (c) At $W=0.28t$, both ${\mathrm{\Gamma }}_0$ and ${\mathrm{\Gamma }}_1$ are visible in the surface DOS. (d) At $W=0.3t, {\mathrm{\Gamma }}_0$ becomes obscured by higher bands and the band associated with ${\mathrm{\Gamma }}_2$ has become flat, causing a large increase in the density of states. In Appendix pp7 we demonstrate that the states probed in the average surface DOS shown here are not Anderson localized.

Figure 6

Comparison of surface perturbation theory and numerically exact results for the model in Eq. (37) with a surface moiré potential from Eq. (45) characterized by the wave vector $Q\approx 4\pi /(11+5\sqrt{5})$ for various SDCs at the $\mathrm{\Gamma }$ point. Dirac point location (a) and corresponding velocity (b). Dots indicate values extracted from dispersion curves on the surface of a TI as a function of the interlayer tunneling $W$. The darker lines labeled “5th order” indicate perturbative results for that Dirac point computed to fifth order (see Sec. 4d and Appendix pp6-s3), which are in good agreement with our numerically exact results. Perturbatively, ${\mathrm{\Gamma }}_2$ has a magic angle condition, however this feature is rounded out and the velocity remains finite in the exact calculation. Notice that Dirac points cross without any level repulsion, as explained in the text. When the energies cross in (a) the velocity calculation becomes unreliable (taken as finite differences on the dispersion curve), so we have omitted those points from the curves. Results are obtained for a linear system size $L=11$ ($Q=2\pi /11$) using Lanczos. The points indicated with $★$ are extracted from the surface DOS from fitting the data to Eq. (41) with a linear system size $L=89$ ($Q=2\pi 8/89$) and agree well with the results obtained from the dispersion.

Figure 7

Effect of the induced potential as a function of the twist angle on the surface DOS ${\rho }_S\left(E\right)$ computed from the lattice model in Eq. (37). Results for an interlayer tunneling $W_0=W_1=0.3t$ across various values of the twist value $\theta$ parameterized by $Q=2\pi m/L=4\pi sin(\theta /2)$ for the integers $m=1,\cdots ,L$ with a KPM expansion order $N_C=2^{12}$ and a system size $L=55$ averaged over 100 realizations of the potential. The density of states as a function of energy for various twist angles demonstrates the original Dirac cone ${\mathrm{\Gamma }}_0$ is shifted to the close proximity of ${\mu }_V=0.45t$ and experiences its most dramatic velocity renormalization at small angles. The figures also reveal the satellite Dirac cone at ${\mathrm{\Gamma }}_1$ where an additional “V shape” appears in the DOS (i.e., a Dirac semimetal scaling $\rho \left(E\right)\sim |E-E_{{\mathrm{\Gamma }}_1}|$) at lower energy. The location of the Dirac cones are marked with blue (${\mathrm{\Gamma }}_0)$ and green (${\mathrm{\Gamma }}_1$) lines.

Figure 8

The evolution of the original surface Dirac point at $\mathrm{\Gamma }$ (top) and the SDC at ${\mathrm{\Gamma }}_1$ (bottom) in the lattice model in Eq. (37) as a function of the twist via $Q=4\pi sin(\theta /2)$. This data is obtained from the low energy scaling of the surface density of states in Eq. (41) on a cubic system of size $L=55$ and a KPM expansion order $N_C=2^{12}$. (a) The location in energy and (b) the velocity of the Dirac point at $\mathrm{\Gamma }$. (c) The location in energy and (d) the velocity of the SDC at ${\mathrm{\Gamma }}_1$. Near $Q=\pi$ distinct behavior occurs due to a change in the most relevant momentum exchange processes in perturbation theory, which needs to be reformulated about $\pi -Q$. We performed this reformulation only for the velocity in (b), which shows excellent agreement. In each case, the perturbation theory qualitatively describes the numerical results.

Figure 9

DFT-to-KPM pipeline for moiré calculations. Top left: Side view of a three quintuple layer (QL) ${\mathrm{Bi}}_2{\mathrm{Se}}_3$ crystal. The larger, red atoms represent bismuth atoms while the smaller, gray atoms are selenium, and one QL is labeled. We use DFT to compute the Kohn-Sham Bloch states and then find the Wannier functions and hopping matrix elements (pictured in bottom right). Truncating this matrix (see Fig. 10), we perform simulations on the full supercell with the KPM.

Figure 10

Truncating the Wannierized model to efficiently model ${\mathrm{Bi}}_2{\mathrm{Se}}_3$. (a) Comparing the electronic structure from the full DFT calculation (gray) to the reconstructed Wannier tight-binding Hamiltonian with spin-orbit coupling (black) for a slab with 5 QLs. A basis with Bi and Se $p$ orbitals effectively captures the low-energy DFT band structure. (b) The bands from the truncated tight-binding Hamiltonian (red) compared to the full Wannierization (gray), with the zoomed in view of the Dirac cone in (e). (c) The 2D surface Brillouin zone of the ${\mathrm{Bi}}_2{\mathrm{Se}}_3$ slab with the cut pictured in (a),(b) drawn as the gray line. The small red zone indicates crystal momenta that support surface states (as calculated by a cut $250\mathrm{meV}$ above the Dirac cone energy). (d) The dependence of the surface Dirac cone gap (from hybridization between the two surfaces) on the number of QLs. (f) The scatter plot for hopping strength versus the bond pair distance in the tight-binding Hamiltonian, as well as the cutoffs used to construct the truncated Hamiltonian (the Hamiltonian terms retained are colored in red).

Figure 11

Effect of gating with the patterned dielectric on the surface density of states ${\rho }_S\left(E\right)$ of ${\mathrm{Bi}}_2{\mathrm{Se}}_3$. Within the bulk gap, the density of states indicates a Dirac cone (far left, $W=0$). As we tune the gating on the patterned dielectric (with length scale $\approx 9\mathrm{nm}$), features consistent with SDCs emerge (red, blue, green, purple lines track SDCs at the gamma point). For energies less than ${\mathrm{\Gamma }}_0$, there is an enhancement of the DOS for moderate values of $W$.

Figure 12

Tracking SDCs energies and velocities for ${\mathrm{Bi}}_2{\mathrm{Se}}_3$. The original Dirac cone ${\mathrm{\Gamma }}_1$ and the three SDCs originating from the gamma point of the moiré Brillouin zone (${\mathrm{\Gamma }}_{2,3,4}$). (a) Shows the energies of the Dirac cones and (b) shows the velocities of the Dirac cones. Notice that ${\mathrm{\Gamma }}_0$ has unavoided crossings with ${\mathrm{\Gamma }}_{1,2}$. All cones have renormalized velocities except ${\mathrm{\Gamma }}_2$ which, within error, has zero velocity. Error bars are dominated by the energy grid of the density of states from which the energies of the states are extracted.

Figure 13

Lack of surface localization. The typical density of states (red) as a function of energy for the $W$ values in Fig. 5. In all cases, the typical density of states remains a finite fraction of the total surface density of states (black), indicating that the states stay extended and do not localize. The surface density of states is taken at $N_C=2^{12}$ while the typical surface density of states is shown for $N_C=2^9, 2^{10}, 2^{11}$, and $2^{12}$ (with the red shading only shown for $N_C=2^{12}$). All horizontal (density of states) scales are the same.