Commensurate and incommensurate double moiré interference in twisted trilayer graphene

When stacking two or more two-dimensional van der Waals materials with a lattice mismatch or a relative twist angle, a moiré superlatice is formedGeim and Grigorieva (2013). For some moiré superlattices, a distinguishing feature is the appearance of flat bands at charge neutrality for which the renormalized Fermi velocity is zero, resulting in the Coulomb interaction strength to significantly exceed the kinetic energy of electrons in the flat band, favouring electron-electron correlationsBistritzer and MacDonald (2011). Flat bands appear in twisted bilayer graphene (TBG) with the twist angle equals 1.05^∘, termed the first magic angle, which exhibits a very interesting range of exotic phenomena including Mott insulatingCao et al. (2018a), superconductivityCao et al. (2018b), ferromagnetismSharpe et al. (2019), Chern insulatorsNuckolls et al. (2020), quantum anomalous Hall effect (QAHE)Serlin et al. (2020) and ferroelectricityKlein et al. (2022). These exciting discoveries have inspired a vast theoretical and experimental search to extend the family of moiré superlattices that exhibit the correlation-driven behaviours, including twisted monolayer-bilayer graphenePolshyn et al. (2020); Xu et al. (2021), twisted bilayer-bilayer grapheneCao et al. (2020); He et al. (2021), trilayer graphene on hexagonal boron nitrideChen et al. (2019a, b), twisted multilayer graphenePark et al. (2022); Burg et al. (2022) and transition metal dichalcogenidesTang et al. (2020); Wang et al. (2020). These moiré superlattices share both similarities and differences with the TBG in symmetries, band topology and interaction strength, which could help us have a deeper understanding of the correlated behaviours in TBG. Compared to the TBG, some of the moiré superlattices may have practical advantages in fabrication or tunability of the physical properties.

Recently, twisted trilayer graphene (TTG) has gained extensive attention due to the presence of unconventional correlated statesPark et al. (2021); Cao et al. (2021). In twisted trilayer graphene, new tunable degrees of freedom are introduced by the addition of an extra third layer on the bilayer graphene. For instance, in the TTG with two consecutive twist angles $\theta_{12}$ and $\theta_{23}$, the beatings of two bilayer moiré patterns may lead to a more complex supermoiré patternZhu et al. (2020a); Zhang et al. (2021). In fact, the experimental technique to realize the TTG is readily available. Different from the TBG, electronic structures of TTG are highly dependent on the original stacking arrangements and on which layer is twistedPolshyn et al. (2020); Wu et al. (2021a), and are more sensitive to external perturbationsWu et al. (2021b); Lopez-Bezanilla and Lado (2020); Chen et al. (2021). In mirror symmetric TTG, a set of dispersive bands coexists with flat bands at charge neutrality, and the magic angle is $\sqrt{2}$ times larger than that of the TBGCarr et al. (2020). Robust and highly tunable superconductivity has been observed in magic angle mirror symmetric TTGPark et al. (2021); Cao et al. (2021). It has been reported that TTG without mirror symmetry hosts a variety of correlated metallic and insulating states, and topological magnetic statesPolshyn et al. (2020); Chen et al. (2021); Xu et al. (2021). The magic angle of such low symmetry TTG approximates that of the TBG, and possesses correlated states that are asymmetric with respective to the external electric fieldXu et al. (2021); Ma et al. (2021).

Compared to the TBG, there are many new challenges in the theoretical calculations of the electronic structures of supermoiré TTG: the large size of the moiré pattern and the lack of commensurate supercell for general twist angles. Firstly, the system size becomes much larger. For example, in the mirror symmetric TTG, the number of atoms in a particular twist angle is half more than that of TBG with the same twist angle. When the twist angles are small, the length of the supermoiré period can be extremely large. The case of the mirror-asymmetric TTG with two independent twist angles is even worse. For example, the moiré length of the TTG in Fig. 1 is 2.6 times larger than that of TBG with the same angle. Secondly, because the two twist angles are independent, the supercell description is no longer valid in some cases since the TTG samples become incommensurate.

In theoretical calculations, electronic structures of twisted multilayer graphene are commonly described by effective continuum modelsZhu et al. (2020b); Li et al. (2019); Lei et al. (2021); Carr et al. (2020). These continuum models with momentum space basis are efficient for computation and require no constraints on the twist angles to construct commensurate supercells. The continuum results show that the van Hove singularities (VHSs) of TTG are significantly dependent on the tuning parameters, for instance, twist anglesZhu et al. (2020a) and original stacking arrangementsLi et al. (2019). Similar to the TBG, the lattice relaxations significantly influence the electronic behaviours of TTG with tiny twist angles and is an important factor to obtain realistic description of the systemShi et al. (2020). Up to now, the lattice relaxation effect is only considered via a continuum model to consider the in-plane distortions and a generalized stacking fault energy to account for the interlayer couplingCarr et al. (2020); Zhu et al. (2020c). An atomistic simulation of the lattice relaxation effects on the electronic structures of supermoiré TTG is still missing. In practice, although the control of the global twist angle with a precision of about 0.1 degrees has been achieved, twist-angle variation across different parts of the device with the order of $\pm 0.01^{\circ}$ are still existsUri et al. (2020); Kazmierczak et al. (2021). How will such angle disorder affect the electronic structures of TTG?

In this paper, we systematically investigate the moiré interference in both commensurate and incommensurate TTG. In particularly, the lattice relaxation, twist angle and angle disorder effects on the electronic structures of TTG with two independent twist angles are studied. To address the challenges of the lack of periodicity and the large system size, we adopt a round disk method to construct the TTG samples with arbitrary twist angles and calculate the electronic properties via the tight-binding propagation method (TBPM) implemented in our home-made package TBPLaSLi et al. (2022). The TBPM is based on the numerical solution of time-dependent Schrödinger equation and requires no diagonalization processes. Importantly, both memory and CPU costs scale linearly with the system size. So the TBPM is an efficient method to calculate electronic properties of large-scale and complex quantum systemsYu et al. (2019); Shi et al. (2020); Kuang et al. (2021). The lattice relaxation is considered via the semi-classical molecular dynamics simulation. We find that the electronic behaviours of TTG are quite sensitive to both the two independent twist angles and angle disorders, in particular for the TTG with mirror symmetry. The systems can have either a constructive or destructive moiré interference depending on the two twist angles. Interestingly, by modulating the angles, novel states, for instance the kagome-like states can be constructed in TTG. The VHSs and local density of states (LDOS) mapping can be utilized as quantities to identify the twist angles of TTG in case that the unit cell size corresponds to a range of different possible sets of twist angle pairs.

This paper is organized as follows. In Sec. II, we introduce the notation and the geometry of the twisted trilayer graphene, the TB model and the computational methods. In Sec. III, for different twist angle configurations, we show the low energy density of states as a function of twist angle for TTG with and without lattice relaxation. The real space distribution of the electron states of energies at VHS near charge neutral point for different twist angle configurations are also investigated. In Sec. IV, we discuss the effect induced by the twist angle disorder on the VHS. Finally, we give a summary of our work.

In this part, we investigate the twist angle and lattice relaxation effects on the electronic properties of TTG. As we mentioned before, the structures of TTG are strongly dependent on the original stacking arrangements and on which layer is twisted. Here we mainly focus on two different cases: the case I is the samples with $-\theta_{12}=\theta_{23}$, which have mirror symmetry; the case II is a stack of graphene layers where each layer is rotated by a constant amount with respect to the previous one. Each sample in case II has two consecutive twist angles with $\theta_{12}=\theta_{23}$. In case I, the moiré length of the TTG is equal to that of TBG with the same twist angle, whereas the moiré length of TTG in case II is much larger than both of the bilayer moiré lengths due to the interference between the two bilayer moiré patterns. According to the Eq. (5), if $\theta_{12}=\theta_{23}=4.4^{\circ}$, the moiré length of the TTG is around 41.7 nm. For case I with $-\theta_{12}=\theta_{23}=4.4^{\circ}$, the moiré length is only around 3.2 nm. The density of states as a function of twist angle $\theta$ to explore the evolution of the van Hove singularities and the lattice relaxation effects are calculated. Due to the incommensurate feature in both cases, we use the round disk model with radius of 172 nm in all calculations. Previous results have been shown that the radius of 172 nm is large enough to get rid of the influence of the edge statesYu et al. (2019). In the round disk method, commensurate or incommensurate systems with any twist angles can be constructed. Therefore, an open boundary is adopted in the calculations.

We first discuss the common features emerging in these two different cases. As shown in Fig. 2, the twist angles significantly modulate the energy position of VHS. The red regions represent VHSs. As the twist angle decreases, the VHS gap decreases first to reach a minimum where the magic angle appears, and then increases in some cases. Moreover, the lattice relaxation obviously modifies the DOS of TTG with tiny twist angles. For samples in case one without lattice relaxation (rigid), shown in Fig. 2(a), two VHSs near charge neutral point (CNP) have their gap narrowing as twist angle $\theta$ decreases, and merge at CNP when angle $\theta=1.57^{\circ}$ where the DOS reaches the maximum magnitude. We name this angle the magic angle. In fact, a sharp DOS peak located at the CNP appears when $1.2^{\circ}\leq\theta\leq 1.6^{\circ}$. With the angle continually decreasing, the VHS gap first increases and then decreases with VHSs merge again at $\theta=0.3^{\circ}$. When consider the lattice relaxation (relaxed) in case one, shown in Fig. 2(b), the evolution of the VHS gap shows similar tendency as the rigid case. The VHSs merge at CNP and reaches maximum (red area) at the magic angle $1.35^{\circ}$. The VHS then splits and diminishes as $\theta$ continues to decrease. Compared to the rigid case, the relaxed samples exhibit a narrower range of magic angles, and have VHSs with reduced magnitudes in case of tiny twist angles.

Compared with ($-\theta$,$\theta$) structures, VHSs evolve differently in ($\theta$,$\theta$) structures as shown in Fig. 2(c) and (d). The DOS in case two is orders of magnitude lower than that of the case one. For samples in case two, adjacent layers form two identical bilayer moiré periods with a relative rotation of $\theta$ forming a supermoiré. The length of the supermoiré period inversely proportional to $1/\theta^{2}$. For the rigid structures shown in Fig. 2(c), as $\theta$ decreases, the VHS gap narrows and reaches minimum ($20$ meV) at $2.1^{\circ}$ but the two VHSs never merge at CNP. This result shows agreement with previous study with a continuum model approachZhu et al. (2020b). For the relaxed structures in Fig. 2(d), the VHSs merge at CNP at $1.57^{\circ}$, which is different from the rigid case.

For general twist angles, Fig. 3 shows the DOS of relaxed samples as a function of $\theta_{23}$ and fixed $\theta_{12}$. In case one with $\theta_{12}=-1.35^{\circ}$, when $\theta_{23}=1.8^{\circ}$, the VHS gap has a value of 30 meV. With $\theta_{23}$ decreasing to the magic angle $1.35^{\circ}$, the first and second VHSs approach the CNP and the first VHS merge at the CNP when $\theta_{23}=1.35^{\circ}$. When the $\theta_{23}$ continually decreases, a sharp DOS peak still appears at the CNP but with a lower magnitude due to the destructive interference effects from the superlattice between 2-3 layer pair. For the sample with $\theta_{23}=-\theta_{12}$ where the top and bottom layers are perfectly aligned, there is a strong increases in the DOS due to the constructive interference effects. If we use the presence of VHS as a proxy for electronic correlations, such correlation-driven behaviour is extremely sensitive to a slight change in twist angles, which is similar to previous resultsZhang et al. (2021). Figure 3(b) shows DOS of relaxed structures in case two with different $\theta_{23}$ and fixed $\theta_{12}=1.57^{\circ}$. As $\theta$ increases from magic angle $1.57^{\circ}$, VHSs split with the gap becomes wider, and the magnitude of the VHS becomes smaller. As $\theta$ decreases from $1.57^{\circ}$, the VHS has a decrease of its magnitude and vanishes as $\theta$ approaches $0^{\circ}$. Similar to the case one, for $\theta_{23}$ equals to the magic angle, the VHS suffers a constructive interference effect from these two bilayer moiré patterns. For $\theta_{23}$ away from the magic angle, the electronic structures can be understood as a magic angle TBG with a destructive interference from the second bilayer moiré structure. However, the TTG in case two is less sensitive to the angle disorder than the case one. For TBG and TTG with mirror symmetry, the twist angle can be identify from the moiré length since the unit cell size corresponds to a unique set of twist angles. On the contrary, for TTG without mirror symmetry, a given unit cell size corresponds to a range of different possible sets of twist angle pairs. In this case, the twist angle can be identified by the above VHS evolution with twist angles in experiment.

The LDOS mapping is another quantity that can be used to identify the twist angles in experiment. The calculated LDOS mapping of relaxed structures in case two with ($1.57^{\circ}$,$\theta$), as shown in Fig. 4, demonstrate the effect of the periodic potential of the moiré pattern on the localization of states in real space. We focus on the energies at the van Hove peaks near the charge neutral point in the DOS and show the mapping of the three graphene layers separately. In TTG, the length of the bilayer moiré period formed by adjacent graphene layers is $a/2\sin{(\theta_{ij}/2)}$. Smaller relative twist angle gives larger moiré period. In the configuration of ($1.57^{\circ}$,$\theta$), the bottom and middle layers with twist angle of $1.57^{\circ}$ forms a moiré period of length $\lambda_{1}=8.98$ nm, the middle and top layers with twist angle of $\theta_{23}$ forms a moiré period of length $\lambda_{2}=a/2\sin{(\theta_{23}/2)}$. Since the ratio of the two moiré lengths determines the features of the LDOS mapping, we categorize the mapping results into four types.

For type $A$ illustrated in Fig. 4 with $\theta_{23}=0.5^{\circ}$, $\lambda_{2}$ is much larger than $\lambda_{1}$, the VHS states are mainly localized on the middle and bottom layers. It is obvious that in the limit $\theta_{12}\gg\theta_{23}$, TTG decomposes into a decoupled TBG moiré supercell and a graphene monolayer. The top layer does not contribute significantly to the electronic structures in the low-energy range. That is, the TTG can be considered as a magic angle TBG modified by an effective potential from the third layer. In the middle layers, features of both moiré patterns are shown, and the localization is enhanced where the $AA$ areas of both moiré patterns coincide. For type $B$ with the ratio of the moiré length $\lambda_{2}/\lambda_{1}$ around 1.6, for instance $\theta_{23}=0.98^{\circ}$, the bottom and top layers show features of bilayer (2L) moiré pattern slightly affected by the other outer layer, the middle layer gives a new complex supermoiré pattern with period $\lambda_{3}$ larger than both 2L moiré periods. For type $C$ with $\theta_{23}=1.89^{\circ}$, the mapping features are similar to that of type $B$. The middle layers also clearly show supermoiré patterns much larger than both 2L moiré patterns. For type $D$ with $\theta_{23}=\theta_{12}$ the moiré length ratio value $\lambda_{1}/\lambda_{2}$ is 1. The LDOS mappings only show bilayer magic angle features, and the VHS states are mainly localized in the middle layer. Obviously, the LDOS mapping is highly dependent on twist angles. By tuning the twist angles, new novel superstructures could constructed. For instance, in case two with $\theta_{12}=1.57^{\circ}$ and $\theta_{23}=0.98^{\circ}$, a Kagome-like lattice based on moiré pattern are constructed on the top layer. All in all, for the case $\theta_{12}\approx\theta_{23}$, the systems suffer a strong moiré interference, and the VHS states are mainly localized in the middle layer. For $\theta_{12}\gg\theta_{23}$, the systems have a weak moiré interference, and can be decomposed into a TBG with $\theta_{12}$ and a graphene monolayer.

Stacking two 2D materials to an arbitrary angle with a precision of $0.1^{\circ}$ is still challenging in experiment. Even using the ’tear and stack’ technique to fabricate the moiré superlattices, twist-angle disorders are still unavoidable from the nonuniformity of the twist angle across the large-scale sample in experiments. In fact, in a high-quality graphene moiré pattern, the main source of disorder is the variations of twist angles across the sample. Previous results have been demonstrated that the correlated phases are extremely sensitive to a slight change in twist anglesZhang et al. (2021); Uri et al. (2020). The angle disorder may explain why two different samples with identical twist angles manifest quite different electronic properties in experiments. Compared with TBG, the more tunable TTG has more chance suffering the angle disorder since one has to precisely control more than one angle in a sample during the fabrication process. In TTG, the moiré length is determined by the two independent twist angles. A small variation of one twist angle could results in a huge increase of the moiré length. Then, how will the angle disorder affects the electronic structures of TTG? As we mentioned before, we can construct TTG with arbitrary twist angles by utilizing the round disk method. That is, the twist angles can be continuously tuned no matter whether the system is commensurate or incommensurate. By combining the round disk method with the TBPM, it is convenient to investigate the angle disorder effects. In this part, we discuss the electronic structures of TTG in the presence of angle disorder. In particular, we focus on the angle disorder effects on the VHSs.

In TTG structures of case one with ($-\theta$,$\theta$) where the two moiré pattern align, we introduce misalignment by an extra twist angle deviation $\Delta\theta=0.02^{\circ}$ for $\theta_{12}$, i.e. ($-\theta\pm 0.02^{\circ}$,$\theta$). Figure 5 shows the effect of twist angle disorders on the DOS of TTG with three different angles $\theta=1.57^{\circ},\;0.98^{\circ},\;0.1^{\circ}$ where $\theta=1.57^{\circ}$ is the magic angle. Let us first focus on the magic angle case. In rigid magic angle TTG, with a variation of $\Delta\theta$ in only $\theta_{12}$, the positions as well as the width of the VHS are unaffected, whereas the peaks of the VHSs are remarkably suppressed. If we quantify the angle disorder effects by a ”BCS superconducting transition temperature”, which is approximately described as $T_{c}\propto exp(-\frac{1}{g\rho(E_{vH})})$, our results suggest that the angle disorder strongly suppresses $T_{c}$. In the $T_{c}$ expression, $E_{vH}$ is the energy of the VHS peak, $g$ is the electron-phonon coupling. On the contrary, in relaxed magic angle TTG, the positions, width and amplitude of the VHSs are unaffected by the angle disorder. The results in Fig. 5(b) shows that the TTG systems with mirror symmetry are protected against twist angle disorders, which is consistent with previous resultsCarr et al. (2020). In TTG with twist angles smaller than magic angles, the VHS is quite sensitive to the twist angle disorder. For TTG with $\theta=0.98^{\circ}$, the angle disorders smear some peaks located round 0.3 eV in the DOS. In TTG with $\theta=0.1^{\circ}$, both the width and the amplitude of the first and second VHS are modulated by the angle disorder.

In summary, we adopt a real-space and atomistic approach for the simulation of the electronic behaviour and the relaxation effect of twisted trilayer graphene for general twist angle pairs. We demonstrate how the position and strength of the VHSs evolve with the twist angle in different situations. We mainly focus on two different structures: one has a mirror symmetry with angle pair (-$\theta$,$\theta$) and the other has two consecutive angles with ($\theta$,$\theta$). Due to the moiré interference, the moiré length of the systems in case two are far more larger than that of case one with the same angle $\theta$. Our results show that the atomic relaxations have significant effects on the VHS properties and the emergence of magic angle. The position of the VHSs are highly dependent on the twist angle. We find that for mirror symmetric ($-\theta$,$\theta$) structures, the magic angles are $1.57^{\circ}$ in rigid case and $1.35^{\circ}$ in relaxed case. For ($\theta$,$\theta$) structures, the magic angles are $2.1^{\circ}$ and $1.57^{\circ}$ in the absence and presence of lattice relaxations, respectively. In TTG with two independent twist angles, when one of the twist angles deviates from the magic angle within a range, the VHS evolution shows a destructive moiré interference effect from the angle changes, and the VHSs follow only with the magic angle. We then show that the LDOS mappings provide an intuitive real-space representation of the double moiré interference and the resulting large complex supermoiré pattern. We find that the mapping features are closely related to the ratio of the two moiré length. For the case that the two angles are identical, the system suffers a strong moiré interference effect and the VHSs states are mainly localized in the middle layer. For one angle far away from the other, the system has a weak interference and can be decoupled into a TBG and a monolayer graphene. By modulating the twist angles, Kagome-like states are constructed due to the interplay between different layers. More importantly, we could use the LDOS and LDOS mapping to identify the twist angles of the TTG in case that the unit cell size corresponds to a range of different possible sets of twist angle pairs.

For mirror symmetric ($-\theta$,$\theta$) TTG that detected a superconductivity, we find that the twist angle disorder which breaks the mirror symmetry strongly affects the VHS property. This may explain why two different samples with identical twist angles manifest different electronic properties in experiments. We also find that the relaxation weakens the disorder effect. In fact, for such small disorders, the relaxation can restore the system to the favorable AA stacking for outer layers. Our results could provide a guide for the experiment to explore the flat band behaviours in supermoiré TTG.