Disorder in twisted bilayer graphene

We develop a theory for a type of disorder in condensed matter systems arising from local twist-angle fluctuations in two strongly coupled van der Waals monolayers twisted with respect to each other to create a flat-band moiré superlattice. The paradigm of “twist-angle disorder” arises from the currently ongoing intense research activity in the physics of twisted bilayer graphene. In experimental samples of pristine twisted bilayer graphene, which are nominally free of impurities and defects, the main source of disorder is believed to arise from the unavoidable and uncontrollable nonuniformity of the twist angle across the sample. To address this physics of twist-angle disorder, we develop a real-space, microscopic model of twisted bilayer graphene where the angle enters as a free parameter. In particular, we focus on the size of single-particle energy gaps separating the miniband from the rest of the spectrum, the Van Hove peaks, the renormalized Dirac cone velocity near charge neutrality, and the minibandwidth. We find that the energy gaps and minibandwidth are strongly affected by disorder while the renormalized velocity remains virtually unchanged. We discuss the implications of our results for the ongoing experiments on twisted bilayer graphene. Our theory is readily generalized to future studies of twist-angle disorder effects on all electronic properties of moiré superlattices created by twisting two coupled van der Waals materials with respect to each other.

Figure 1

The density of states $\rho \left(E\right)$ as a function of energy $E$ for the lattice model of twisted bilayer graphene at a twist angle $\theta =1.{05}^{\circ }$, a linear system size $L=569$, a kernel polynomial method [26] expansion order $N_C=2^{17}$, and a weak breaking in the interlayer tunneling between AA and AB sites $[w_0/w_1=0.75$, $w_1=w$ where $w_0$ ($w_1$) is the strength of AA and BB (AB and BA) tunneling], which captures lattice relaxation effects [27, 28] and it opens a hard gap on both sides of the semimetal miniband. We note that at small angles, a single parameter controls the physics: $w/[2v_F{k}_{D}sin(\theta /2)]$, so lowering the angle is equivalent to increasing $w_1$. Therefore, one can read the plots of smaller $w_1$ as at an angle larger than $1.{05}^{\circ }$. This density of states has a number of features relevant to the physics: Van Hove peaks, gaps, and the velocity (as determined by the scaling of the density of states). Dark (light) blue lines give the calculated density of states for finite (zero) values of the parameter $w$ as shown in the inset of the figure.

Figure 2

(a) A schematic of graphene and the notation we use for our model. The A (B) sublattice is represented by the blue (orange) lattice sites. The unit cell for the triangular lattice is shown by the dashed central hexagon. The lattice vectors are ${\mathbf{a}}_1=(\sqrt{3}/2,3/2)$ and ${\mathbf{a}}_2=(-\sqrt{3}/2,3/2)$, and we further define ${\mathbf{a}}_3\equiv {\mathbf{a}}_2-{\mathbf{a}}_1$, ${\mathbf{a}}_4\equiv -{\mathbf{a}}_1$, ${\mathbf{a}}_5\equiv -{\mathbf{a}}_2$, as well as ${\mathbf{a}}_6\equiv -{\mathbf{a}}_3$. (b) A course-grained view of the tunneling between the layers calculated from ${\mathcal{T}}_0$ and ${\mathcal{T}}_1$ in Eq. (5) which defines the energy parameters $w_0$ and $w_1$; the color represents whether AA, AB, or BA hopping is dominant based on the chance for an electron on a site in layer 1 to hop onto sublattice A or B on layer 2, given by $P_X\left(\mathbf{r}\right)=|{\left[{\mathcal{T}}_0\left(\mathbf{r}\right)\right]}_X{|}^2+6{\left|{\left[{\mathcal{T}}_1\left(\mathbf{r}\right)\right]}_X\right|}^2$. Note that $C_3$ is broken and the moiré unit cell is larger than in real TBG. Both of these effects are relatively small. (c) Complementary to the real-space picture, in momentum space the lattice Brillouin zone is effectively downfolded by a factor of 3 from the moiré Brillouin zone after unrotating the two graphene layers; this introduces small gaps in the band structure at these points. (d) In our model, the effect of the twist is entirely contained within interlayer coupling, so we model disorder by changing the continuous twist parameter $\theta$ within different regions of space. In this common example, we break up the system into four equal regions and pick a value of ${\theta }_j$ that are drawn from the box distribution $[(1-W_R/2)\theta ,(1+W_R/2)\theta ]$ with $\theta =1.{05}^{\circ }$.

Figure 3

(a) The calculated density of states $\rho \left(E\right)$ for TBG without disorder as a function of energy $E$ for various interlayer tunneling strengths $w=w_1$ [keeping $w_0/w_1=0.75$ where $w_0$ ($w_1$) denote the strength of AA and BB (AB and BA) tunneling] at a low twist angle of $\theta =1.{05}^{\circ }$ close to the magic angle, a system size of $L=569$, and a KPM expansion order of $N_C=2^{17}$ in the lattice model. The calculated minibandwidth in the magic-angle regime $w=110\mathrm{meV}$ is consistent with other studies of the continuum model and the KPM numerical resolution limits to the extent we can access the low-energy regime near charge neutrality. (b)–(d) Comparisons between our lattice model and the continuum theory near $E=0$ and $\theta =1.{05}^{\circ }$ for $w=80$, 100, and $110\mathrm{meV}$, respectively, we find remarkable agreement. The insets show the details of the miniband. At $\theta =1.{05}^{\circ }$ and $w=100\mathrm{meV}$, inset of (c), we see a splitting of the Van Hove peaks that is missing from the continuum model associated with additional zone folding in this model. This is seen clearly in the right inset; the left inset shows how the gap of the lattice model here and in the continuum model also match rather well. In (d) at the magic angle $\theta =1.{05}^{\circ }$ and $w=110\mathrm{meV}$, we see that the Van Hove peaks never clearly merge as they do in the continuum model. Again, this is clearly seen in the right inset. The continuum model data here include 338 bands and have $N_C=2^{13}$ or $2^{14}$ whereas the lattice model has $L=569$ and $N_C=2^{17}$. Overall, the agreement with the continuum TBG model is quite excellent.

Figure 4

The effects of twist disorder on the low-energy density of states. The density of states $\rho \left(E\right)$ as a function of energy $E$ for a clean twist angle $\theta =1.{05}^{\circ }$, linear system size $L=569$, and a KPM expansion order of $N_C=2^{17}$ starting in the semimetal regime of the the TBG model (top) as well as in the magic-angle regime (bottom), for different twist disorder strengths $W_R$ (that characterizes the width of a box distribution $[(1-W_R/2)\theta ,(1+W_R/2)\theta ]$ with $\theta =1.{05}^{\circ }$ from which we sample the random twist angle in each patch). In each case, the randomness smoothly fills in the gap while also smearing out the Van Hove peaks. The insets in the bottom two figures are a zoom-in of the band gap that clearly fills in with increasing disorder.

Figure 5

Summary of results on the miniband properties in the TBG model with a clean twist angle $\theta =1.{05}^{\circ }$ extracted from system sizes $L=569$ and a KPM expansion order $N_C=2^{17}$. (a), (b) The estimated gap size ${\mathrm{\Delta }}_{\mathrm{MB}}$ as a function of disorder strength in the twist angle $W_R$ and the interlayer tunnelings $w$ (where $w=w_1$ and the ratio of AA and BB tunneling to AB and BA tunneling is $w_0/w_1=0.75$). (c), (d) The velocity $v/v(w=0)$ as calculated from the density of states as a function of disorder $W_R$ remains approximately unchanged in the presence of disorder $W_R$ (for each value of $w$). (e), (f) The minibandwidth $D_{\mathrm{MB}}$ for interlayer tunneling $w$ and disorder $W_R$. Note that for larger disorder strength ($W_R=6\%$ or above) in (e) the bandwidth appears to plateau; this is just an artifact arising from disorder completely filling out the gap at this point. While the gap and bandwidth are strongly affected by disorder, the velocity remains unchanged. The red dashed line in (f) that sets the maximum that the minibandwidth can achieve, is determined from the gaps in (b).

Figure 6

The effects of twist disorder on the properties of the Van Hove peaks for a clean twist angle $\theta =1.{05}^{\circ }$, a linear system size $L=569$, and a varying KPM expansion order ($N_C$) in (a) whereas in (b)–(d) we use $N_C=2^{17}$. (a) As we scale the Chebyshev expansion order, we see that the Van Hove peak is logarithmically divergent (with a fit shown as a black dashed line), but once we add disorder, it rounds out and saturates to a finite value. (b) The energy separation between Van Hove peaks remains stable as disorder increases even though we find (c) that the full width at half-maximum (FWHM) of the Van Hove peaks becomes broader as disorder increases. (d) The estimated BCS critical temperature or the effective coupling constant [see Eq. (8) in the main text] from the density of states at the Van Hove peak as disorder is tuned up for various values of $w$.

Figure 7

(a) The calculated density of states $\rho \left(E\right)$ as a function of energy $E$ for the spin-orbit-coupled model of Dirac points perturbed by a quasiperiodic potential, with a quasiperiodic wave vector $Q=2\pi F_{n-2}/F_n$ with the system size $L=F_n=144$ and a KPM expansion order $N_C=2^{14}$. (b) A depiction of how we break up the SOC square lattice model into regions of different quasiperiodic wave vector $Q_i$ (to simulate disorder), which are taken from a box distribution about a central value. We vary both the number of regions and the size of disorder in each region.

Figure 8

The disorder-free density of states $\rho \left(E\right)$ as a function of energy $E$ obtained from a linear system size $L=144$ and a KPM expansion order $N_C=2^{14}$ starting in the semimetal regime of the model, comparing the case of a fixed random phase across the entire sample (b), (d) and a different random phase in each patch (a), (c) for different strengths of disorder in the wave vector and $n_P=7$ randomly placed patches. Note that the random phase in each patch is disordered even for $W_Q=0$.

Figure 9

Density of states as a function of energy in the semimetallic regime of the SOC model focusing on the miniband at low energy using a linear system size $L=144$ and a KPM expansion order $N_C=2^{14}$. We focus on the effects of the different number of random patches used for various different disorder strengths in the quasiperiodic wave vector $W_Q$ from $W=0.35$. Here, we are taking one global phase across the sample to isolate the effects of randomness in $Q$ and choice of patches alone.

Figure 10

Density of states as a function of energy in the magic-angle regime ($W=0.54$) of the SOC model focusing on the miniband at low energy with a linear system size $L=144$ and a KPM expansion order $N_C=2^{14}$. We are displaying the effects of different number of patches of a random wave vector across the sample.

Figure 11

The estimated critical temperature [or effective coupling, see Eq. (8) in main text] from the Van Hove peaks in the DOS as a function of randomness in the twist vector comparing two choices for the random phase for different number of randomly placed patches. (Top row) One fixed phase, corresponding to a single rotation origin. (Middle row) Random phases ${\varphi }_{\mu }\left(i\right)$ in each block. The left panels are $W=0.35$ in semimetal phase, while the right panels are $W=0.54$ at the magic angle. Random phases in each block produce very strong randomness in the model and smear out the Van Hove peaks more easily. (Bottom left) The critical temperature $T_c$ with random phases ${\varphi }_{\mu }\left(i\right)$ in each block but without randomness in $Q$, as function of number of patches $n_p^2$, and normalized by $T_c$ with only one patch. (Bottom right) The gap size as function of randomness. Comparing to the suppression of $T_c$, the gap is filled in for $W_Q\approx 0.5\%$, which is much smaller than the critical $W_Q$ ($\sim 10\%$) needed for Van Hove peaks to be smeared out. These results are obtained from data using a linear system size $L=144$ and a KPM expansion order $N_C=2^{14}$.

Figure 12

Effects of disorder on the renormalization of the velocity of the Dirac cone and the minibandwidth using a linear system size $L=144$ and a KPM expansion order $N_C=2^{14}$. (a) Effective velocity of the Dirac cone and how it is rounded out due to randomness in the wave vector. The finite velocity in the magic-angle regime for $W_Q=0$ is just a finite-size effect [24]. (b) Minibandwidth as a function of disorder in the quasiperiodic wave vector, which monotonically broadens for increasing disorder until the gap is filled in and the miniband is no longer separated from the rest of the band (marked as dashed lines). We include both $W=0.35$ for semimetallic phase and $W=0.54$ for the magic-angle regime. Note that we have set $t=1$ here.

Figure 13

The calculated velocity renormalization in second-order perturbation theory. The dashed line is the result from Ref. [21] while the upper and lower curves are the maximum and minimum velocities of the lattice model written in Eq. (4) near the $K$ and $K^{\prime }$ points. Notice that the velocity for all states never vanishes as it does in the continuum model even though we have demonstrated in the main text that the density of states matches very well. This is plotted for the situation considered in the main text where ${\alpha }_0=0.75{\alpha }_1$ and $\theta =1.{05}^{\circ }$.