Designing Flat Band by Strain

TL;DR

This paper investigates how heterostrain influences moiré band structures in twisted bilayer graphene and TMDs, revealing significant band separation changes and tunable flat bands, which help explain experimental observations.

Contribution

It demonstrates that heterostrain can significantly modify moiré bands in bilayer systems, providing new insights into band engineering and experimental phenomena.

Findings

Heterostrain increases band separation in bilayer graphene.

Heterostrain enables tunable flat bands in bilayer TMDs.

Resolves puzzles in spectroscopy and quantum oscillation experiments.

Abstract

We study the effects of heterostrain on moir\'e bands in twisted bilayer graphene and bilayer transition metal dichalcogenide (TMD) systems. For bilayer graphene with twist angle near $1^\circ$, we show that heterostrain significantly increases the energy separation between conduction and valence bands as well as the Dirac velocity at charge neutrality, which resolves several puzzles in scanning tunneling spectroscopy and quantum oscillation experiments at once. For bilayer TMD, we show that applying small heterostrain generally leads to flat moir\'e bands that are highly tunable.

Equations

$$r^{\prime }=(\mathbb{I}+\mathcal{E})r+d_0,$$

$$\mathcal{E}\cong \left(\begin{array}{cc}{ϵ}_{xx}& {ϵ}_{xy}-\theta \\ {ϵ}_{xy}+\theta & {ϵ}_{yy}\end{array}\right)=\mathcal{S}(ϵ)+\mathcal{T}(\theta ),$$

$${\mathbb{A}}_i^{\prime }=(\mathbb{I}+\mathcal{E}){\mathbb{A}}_i,{\mathbb{G}}_i^{\prime }\cong (\mathbb{I}-{\mathcal{E}}^T){\mathbb{G}}_i.$$

$$\mathbb{A}=\frac{\sqrt{3}}{2a}\beta ({ϵ}_{xx}-{ϵ}_{yy},-2{ϵ}_{xy})$$

$${\mathbb{D}}_{\xi }=(\mathbb{I}-{\mathcal{E}}^T){\mathbb{K}}_{\xi }-\xi \mathbb{A},$$

$${\mathbb{G}}_i^{\prime }\cong (\mathbb{I}-{\mathcal{E}}_1^T){\mathbb{G}}_i,{\mathbb{G}}_i^{\prime \prime }\cong (\mathbb{I}-{\mathcal{E}}_2^T){\mathbb{G}}_i.$$

$$g_i\cong {\mathbb{G}}_i^{\prime }-{\mathbb{G}}_i^{\prime \prime }={\mathcal{E}}^T{\mathbb{G}}_i,a_i\cong {\mathcal{E}}^{-1}{\mathbb{A}}_i.$$

$$H=\left(\begin{array}{cc}{h}_1(\mathbb{k})+V_1(r)& T(r)\\ T^{†}(r)& h_2(\mathbb{k})+V_2(r)\end{array}\right),$$

$$h_l(\mathbb{k})=\sum _{\xi =\pm }-\hslash v_F/a[(\mathbb{I}+{\mathcal{E}}_l^T)(\mathbb{k}-{\mathbb{D}}_{l,\xi })]\cdot (\xi {\sigma }^x,{\sigma }^y),$$

$${\mathcal{S}}_{ua}$$

$$T(r)=$$

$$h_l(\mathbb{k})=\sum _{\xi =\pm }-\hslash v_F/a[(\mathbb{I}+{\mathcal{E}}_l^T)(\mathbb{k}-{\mathbb{D}}_{l,\xi })]\cdot (\xi {\sigma }^x,{\sigma }^y)+\frac{m}{2}({\sigma }^z+1),$$

$$S_v=ϵ\left(\begin{array}{cc}\cos \phi & \sin \phi \\ \sin \phi & -\cos \phi \end{array}\right).$$

$$V_l(r)=\sum _{i=1,2,3}\left(\begin{array}{cc}{V}_c{e}^{i(g_i\cdot r+(-1{)}^l{\varphi }_c)}& 0\\ 0& V_v{e}^{i(g_i\cdot r+(-1{)}^l{\varphi }_v)}\end{array}\right)+h.c.,$$

$$H_{eff}=$$

Full text

I Introduction

Recent experimental discoveries of correlated insulator and superconducting states in two-dimensional moiré materials including twisted bilayer graphene (TBG)Cao et al. (2018a, b); Yankowitz et al. (2019) and graphene-boron nitride heterostructures Chen et al. (2018) have stimulated tremendous interest in engineering flat or narrow bands to realize correlated electron phenomena. For this purpose, moiré superlattices provide a unique and highly tunable material platform. With a vast variety of $2d$ materials and heterostructures available Novoselov et al. (2016), moiré superlattice systems may offer unprecedented advantages for studying many-body physics and realizing exotic states of matter.

Moiré patterns appear ubiquitously in layered $2d$ materials with a slight mismatch in the lattice orientation and/or lattice constant of the layers. Aside from twist angles between layers, heterostrain—which refers to relative strains between layers—provides an alternative way to create and modify moiré patternsHuder et al. (2018a). Using strain to tune moiré bands may have practical advantages than using twist angle. Strain can be controlled in situ via piezoelectric substrate. The possibility of tuning band structure and achieving partially flat band with strain has recently been studied in graphene (see for example Guinea (2012); Frisenda et al. (2017)) and surface states of topological crystalline insulator Tang and Fu (2014).

In this work, we systematically study the effect of heterostrain on moiré band structures for homobilayer systems including bilayer graphene and transition metal dichalcogenides (TMD). Our motivation is twofold. First, (unintentional) heterostrain is ubiquitous in TBG samples Kerelsky et al. (2018), likely due to the interaction with the substrate. However, the effect of heterostrain on the flat bands in magic-angle TBG is not well understood. We find a small amount of strain dramatically changes the low-energy band structure, which resolves several puzzles in transport and spectroscopy experiments all at once. Second, engineering a tunable flat band system by strain provides unprecedented opportunities for future studies. To that end, we focus on bilayer TMDs and show the conditions for flat moiré bands induced by heterostrain without any twist.

This work consists of two parts. In the first part, we show that heterostrain in TBG significantly increases the energy separation between moiré conduction and valence bands, for instance, to $\sim 30$meV with $0.5\%$ uniaxial heterostrain as commonly observed in STM experiments Kerelsky et al. (2018). The energy separation of the van Hove singularities saturates over a range of twist angles $\theta$ around $1^{\circ}$, rather than being extremely sensitive to small deviation from the magic angle in the unstrained case. These findings explain the unexpectedly large separation of van Hove singularities observed in the recent STM experiments on TBG with $\theta\sim 1^{\circ}$ Kerelsky et al. (2018).

Moreover, in the presence of heterostrain, the conduction and valence bands of the same valley remain connected by two Dirac points near charge neutrality. However, due to the lowered symmetry, these two Dirac points are away from mini-Brillouin zone corners and are no longer degenerate in energy. This effect may explain the 4 (instead of 8) fold Landau level degeneracy111The possibility of the lifted degeneracy by heterostrain is brought to our attention by Matthew Yankowitz and Cory Dean. observed in the transport experimentsCao et al. (2018a, b); Yankowitz et al. (2019). We also show that heterostrain sets a lower bound for the Dirac velocity, preventing it from vanishing at the magic angle. The Dirac velocity is anisotropic and on the order of $0.14v_{F}$ for a small heterostrain of $0.3\%$-$0.6\%$ ($v_{F}$ being the bare Dirac velocity of monolayer graphene). This value is comparable to the one inferred from quantum oscillation and capacitance measurementsCao et al. (2018a).

In the second part, we study bilayer TMDs with volume preserving heterostrain. In contrast to the twisted bilayer graphene, we find that nearly flat moiré bands are generically present near the top of valence band in heterostrained bilayer TMD without fine tuning and without the need of twist. In addition, the moiré band gap and band structure are highly tunable by strain, pressure and displacement field, thus providing an ideal platform for correlation-driven phenomena.

II 2-dimensional strain and general continuum models

In this section, we consider the general geometrical description of $2$-dimensional strain. Mathematically, the coordinate transformation in $2d$ can be written as

[TABLE]

where $\mathbb{d}_{0}$ is a $2$-dimensional vector that parametrizes the displacement and $\mathcal{E}$ is an arbitrary $2$-dimensional matrix that contains the strain and rotation. In the small deformation limit, $\mathcal{E}$ can be written as

[TABLE]

where the symmetric part, labeled by $\mathcal{S}(\epsilon)$, represents the strain, and the anti-symmetric part, labeled by $\mathcal{T}(\theta)$, represents the rotation. We emphasis that strain for $2d$ materials, different from the rotation angle which is a single parameter, is characterized by a 2 by 2 symmetric matrix, which has 3 free parameters. Therefore, it provides more possibilities to engineer the superlattice structures for layered $2d$ materials.

Let us review the effects of geometric deformation on the properties of monolayer graphene. For free monolayer graphene, we define $\mathbb{A}_{1}=a(1,0)$ and $\mathbb{A}_{2}=a(1/2,\sqrt{3}/2)$ to be the primitive lattice vectors. Correspondingly, the reciprocal lattice vectors are $\mathbb{G}_{1}=\frac{2\pi}{a}(1,-1/\sqrt{3})$ and $\mathbb{G}_{2}=\frac{2\pi}{a}(0,2/\sqrt{3})$. $\mathbb{K}_{\pm}=\mp(2\mathbb{G}_{1}+\mathbb{G}_{2})/3$ are referred as two valley points. The low energy description of monolayer graphene contains two massless Dirac fermions at $\mathbb{K}_{+}$ and $\mathbb{K}_{-}$ points with spin degeneracy. Geometrically, a deformation $\mathcal{E}$ changes the shape of the unit cell as well as the Brillouin zone. Mathematically the rescaled primitive and reciprocal lattice vectors are

[TABLE]

A generic $\mathcal{E}$, including both nonzero strain $\mathcal{S}$ and rotation $\mathcal{T}$, breaks almost all the point group symmetry of the lattice except the inversion $C_{2z}$. In addition to the geometric effect, the strain $\mathcal{S}$ adjusts the distances between atoms, which leads to differences in hopping matrix elements for nearest carbon atoms along different directions. As a result, the locations of the low energy Dirac fermions are shifted away from the rescaled valley points $\mathbb{K}^{\prime}_{\pm}\cong(\mathbb{I}-\mathcal{E}^{T})\mathbb{K}_{\pm}$. Within a simple two center approximation $t(r)\sim t_{0}e^{\beta(r/a_{0}-1)}$, the shift is proportional to the strain in the small strain limit and can be described by an effective gauge connection for the low energy Dirac fermionsNam and Koshino (2017):

[TABLE]

The two Dirac fermions carry opposite charges under this fictitious gauge field. Therefore, the positions of the two Dirac fermions in momentum space are given by

[TABLE]

where $\xi=\pm$ labels the two valleys. The hopping modulus factor $\beta$ is a dimensionless parameter determined by the intrinsic properties of the material. The approximate value of $\beta$ is estimated by first principle calculation for graphene, $\beta_{g}\cong 3.14$Nam and Koshino (2017).

The effects of geometric deformation for monolayer TMDs are similar as for monolayer graphene. For unstrained monolayer TMD, the low energy theory can be modeled by two massive Dirac fermionsLiu et al. (2013); Xiao et al. (2012) at $\mathbb{K}_{\pm}$ points. As a result of the large spin-orbital coupling, the valence bands near two valleys carry opposite spins. Generic deformations reshape the unit cell and the Brillouin zone and break all the point group symmetries. It also shifts the Dirac fermions away from the rescaled valley points $\mathbb{K}^{\prime}_{\pm}$. The shift is again described by Eq. 4 and 5. First principle calculations suggest for WSe2, the hopping modulus factor is $\beta_{WSe_{2}}\cong 2.30$Fang et al. (2018a).

In the rest of the paper, we consider homobilayers systems starting with AA stacking and then apply small twist and heterostrain. The coordinate transformations of the two layers can be described by two deformation matrices $\mathcal{E}_{1}$ and $\mathcal{E}_{2}$. The rescaled reciprocal lattice vectors for the two layers are

[TABLE]

Such geometrical deformation generates a moiré superlattice whose reciprocal lattice vectors $\mathbb{g}_{i}$’s and primitive lattice vectors $\mathbb{a}_{i}$’s are given by

[TABLE]

respectively, where $\mathcal{E}=\mathcal{E}_{2}-\mathcal{E}_{1}$ is the relative deformation matrix. Physical properties only depend on the relative deformation $\mathcal{E}$ in the small twist and strain limit. For all the calculations in this paper, we assume that $\mathcal{E}_{2}=-\mathcal{E}_{1}=\frac{1}{2}\mathcal{E}$, namely the two layers are rotated and strained oppositely with the same magnitude.

In the limit where the moiré superlattice constant is much larger than the atomic scale, the low energy electronic structure can be effectively captured by the continuum modelBistritzer and MacDonald (2011); Lopes dos Santos et al. (2012, 2007). The spirit of the continuum model is the same as the nearly free electron approximation for band structures in solid state physics. Essentially, one should take the low energy bare dispersions of the $2d$ materials and perturb them by the periodic moiré superlattice. A general continuum model for bilayer systems can be schematically written asBistritzer and MacDonald (2011); Yuan and Fu (2018); Koshino et al. (2018)

[TABLE]

where $h_{1}(\mathbb{k})$ and $h_{2}(\mathbb{k})$ are the bare dispersions of layer 1 and 2 respectively. The $T(\mathbb{r})$ describes the spatial dependent interlayer tunneling and the $V_{1}(\mathbb{r})$ and $V_{2}(\mathbb{r})$ describe the intralayer potential induced by the moiré superlattice. The $T(\mathbb{r})$ and $V(\mathbb{r})$ vary between different materials and are also dependent on extrinsic conditions such as pressureYankowitz et al. (2019). In the following sections, we will focus our discussions on two cases $\mathit{1}^{\circ}$ twisted bilayer graphene with uniaxial heterostrain and $\mathit{2}^{\circ}$ bilayer WSe2 with volume preserving heterostrain.

III Twisted bilayer graphene with uniaxial heterostrain

III.1 The continuum model

In bilayer graphene, the bare low energy dispersion for each layer contains two massless Dirac fermions, namely

[TABLE]

where $l=1,2$ labels the two layers, $\xi=\pm$ labels the two valleys. The $\sigma$ matrices act on the pseudospin or sublattice degrees of freedom. Each valley also has two-fold spin degeneracy. $\mathbb{D}_{l,\xi}$ denotes the location of the Dirac fermion as given in Eq. 5. The fermi velocity in monolayer graphene is estimated $v_{F}\cong 10^{6}$ m/sCastro Neto et al. (2009), which gives the kinetic energy scale $\hbar v_{F}/a\cong 2.68$eV.

In this section, we restrict ourselves to twisted bilayer graphene with uniaxial heterostrain. The uniaxial heterostrain refers to a class of strain where the bilayer system is relatively stressed along one direction and unstressed on the perpendicular direction. Scanning tunneling microscope (STM) experiments have indicated such type of heterostrain in twisted bilayer graphene samplesKerelsky et al. (2018); Huder et al. (2018b); Qiao et al. (2018). We anticipate that many generic features of the uniaxial strain should also apply for more general form of strain. Geometrically, uniaxial strain can be described by two parameters, namely the strain magnitude $\epsilon$ and the strain direction $\varphi$. The strain tensor can be written as the followingKerelsky et al. (2018)

[TABLE]

where $\nu=0.16$ is the Poisson ratio for graphene. The combination of twist and strain is described by the relative deformation matrics $\mathcal{E}=\mathcal{T}(\theta)+\mathcal{S}_{ua}(\epsilon,\varphi)$.

For graphene system, $V_{l}(\mathbb{r})$ is parametrically smaller than the interlayer tunneling, therefore we set it to zero. The interlayer tunneling has the following form (following the convention in, for example, Koshino et al. (2018)),

[TABLE]

where $\omega=e^{i2\pi/3}$. The effect of lattice corrugation in TBG on moiré band structure can be included in the continuum model by choosing renormalized tunneling amplitudes $u,u^{\prime}$ and $u^{\prime\prime}$ Koshino et al. (2018), which depend on heterostrain and twist angle. For simplicity, in our calculation we assume the two graphene layers are unrelaxed and use the tunneling parameters $u=u^{\prime}=u^{\prime\prime}=110$meVBistritzer and MacDonald (2011). With this set of parameters, the first magic angle without strain locates around $\theta_{c}\cong 0.95^{\circ}$.

An apparent geometric effect of the combination of twist and heterostrain is that the moiré superlattice is no longer a regular triangular lattice . An example of moiré pattern generated with $\theta=1.05^{\circ}$, $\epsilon=0.7\%$ and $\varphi=0^{\circ}$ is shown in Fig. 1(a). The uniaxial strain makes the moiré dots elliptical, which is visible in local measurements such as STMKerelsky et al. (2018). As a consequence, the Brillouin zone of the moiré superlattice is a distorted hexagonal as shown in Fig. 1(b). In order to compare results for different twists and strains, we will always stretch the irregular hexagonal Brillouin zone to be regular.

A typical band structure of twisted bilayer graphene near magic angle with small uniaxial heterostrain is shown in Fig. 2. Three generic features are worth noticing in these plots. First, the energy separation between conduction and valence bands, as seen from the separation of van Hove singularities, is significantly enlarged compared to the case with no strain. Interestingly, these two bands remain rather flat for most parts of the Brillouin zone. This finding agrees with the large splitting of van Hove singularities of the conduction and valence bands observed in the scanning tunneling spectroscopy (STS) Kerelsky et al. (2018). We will study in detail of the evolution of the bandwidth with respect to the geometrical parameters $\theta$, $\epsilon$ and $\varphi$.

Second, the conduction and valence bands remain connected by two Dirac crossings in each valley, which is due to the $C_{2z}\mathcal{T}$ symmetry. However, these two Dirac crossings are now located at generic points away from Brillouin zone corners because of the lack of three-fold rotational symmetries in the presence of strain. In addition, the energies of the two Dirac points within a valley are shifted by the heterostrain and are no longer degenerate. This may explain the observed 4-fold (instead of 8-fold) Landau level degeneracy in the experimentsCao et al. (2018a, b); Yankowitz et al. (2019) near charge neutrality. We note that the energy offset of the two Dirac points leads to a finite density of electrons and holes at charge neutrality. Nonetheless, the residual electron/hole density is found to be rather small: around $\sim 3\times 10^{10}$cm*-2* (including valley and spin degeneracy) for $0.6\%$ strain. We shall also investigate the dependence of the Dirac fermion energy splitting on the system parameters.

Finally, compared to the unstrained case, the Dirac velocity is found to be greatly enhanced—reaching $0.14v_{F}$ for a small strain of $0.6\%$, and is anistropic due to the lowered symmetry. Moreover, we find that heterostrain sets a lower bound for the Dirac velocity, preventing it from vanishing around magic angle. The enhanced Dirac velocity found here is comparable to the one inferred from quantum oscillation and capacitance experiments Cao et al. (2018a).

III.2 The bandwidth and the Dirac point shift

Now we study $\mathit{1}^{\circ}$ the total bandwidth of the conduction and valence bands and $\mathit{2}^{\circ}$ the energy shift of the two Dirac crossings within a valley as a function of the parameters $\theta$, $\epsilon$ and $\varphi$. We are interested in the regime of twist angle $\theta$ around $1^{\circ}$. Experimental data indicate that the uniaxial heterostrain in twisted bilayer graphene samples can vary from $0.1\%$ to $0.7\%$Kerelsky et al. (2018). In addition, it is easy to convince ourselves that the system is periodic for $\varphi\rightarrow\varphi+\pi/3$ because of the symmetry of the unstrained system. Therefore, we can restrict our studies to $\varphi\in[0,\pi/3)$.

Firstly, let us fix $\theta=1.05^{\circ}$ and $\epsilon=0.3\%$ and plot the density of states the conduction and valence bands for several values of $\varphi$ as shown in Fig. 3. We observe many subtle changes in the structure of van Hove singularities within each bands as the direction of the strain is varied. For instance, $\varphi=0^{\circ}$ has three density of state peaks in each band. However, for $\varphi=20^{\circ}\sim 30^{\circ}$, there is only one prominent peak in each band, which is consistent with the experimental resultsKerelsky et al. (2018). Nonetheless, the total bandwidth stays approximately constant.

Secondly, we study the total bandwidth measured by the separation of the van Hove singularities in the conduction and valence bands, labeled by $\Delta$. Fig. 4 shows the evolution of $\Delta$ as a function of $\theta$ for several different $\epsilon$’s at fixed strain direction $\varphi=25^{\circ}$. An interesting feature is that the bandwidth is insensitive to the twist angle for angles near the magic angle. We also observe that the saturation value of the bandwidth is linearly proportional to the magnitude of the strain as shown in Fig. 5, which suggests the strain effect is dominant near the magic angle. In certain sense, the heterostrain provides an intrinsic kinetic energy cutoff, $\Lambda\sim\epsilon\hbar v_{F}$, for the system and stabilizes the bandwidth.

Finally, we study the energy shift, labeled by $\delta$, of two Dirac points within a valley as a function of $\epsilon$ and $\varphi$ at fixed twist angle. As shown in Fig. 6, the Dirac point shifts highly depend on the direction of the strain. We also find that the maximal value of the Dirac point shift is proportional to the magnitude of the strain. To our surprise, a quite small heterostrain such as $\epsilon=0.1\%$ can create a considerable splitting of the dirac fermions (on the order of $5$meV).

Although a complete analytical understanding of heterostrain effects found in our calculation is currently lacking, we provide some intuition for the observed features in the band structures. In the unstrained twisted bilayer near the magic angle, the nearly flat energy dispersions arise from subtle interference effects between the mismatched energy dispersion within each layer and electron hopping between layers. For example, the state at $\mathbb{K}_{-}$ from layer-1 would be coupled through the moiré interlayer tunneling to three states from layer-2 with the same energy. Similarly, the state at $\mathbb{K}_{-}$ from layer-2 is also coupled to three equal energy states in layer-1. We can crudely think of the magic angle is tuning the bilayer system to a critical situation where the kinetic energy scale vanishes. However, such “interference condition” is explicitly violated by the heterostrain, which introduces a natural kinetic energy cutoff, $\Lambda\sim\epsilon\hbar v_{F}$, to the system. It is easy to see the energies of the three states from layer-2 that can couple to $\mathbb{K}_{-}$ state in layer-1 are now shifted precisely by $\Lambda\sim\epsilon\hbar v_{F}$ because of the heterostrain. This lifting of perfect interference is presumably the cause for the enlarged bandwidth and the Dirac fermion energy shift.

III.3 Higher order van Hove singularities

As shown in Ref. Yuan et al. (2019), with a single tuning parameter such as the twisted angle, one can achieve higher order van Hove singularities in bilayer graphene, which are perfect playgrounds for correlation-driven physics. Here we show that heterostrain is another effective way to generate higher order van Hove singularities. To demonstrate this, we show two examples in Fig. 7 and 8. We find multiple ordinary and higher order van Hove singularities in both conduction and valence bands in these examples. The energies of the van Hove singularities are not required to be the same in general because of the low symmetry of the heterostrained system. Therefore, we may observe multiple peaks in the density of states for certain conditions as shown in Fig. 8.

IV Heterostrain in bilayer transition metal dichalcogenide

IV.1 Heterostrain induced flat bands

For monolayer TMD, the model hamiltonian for one valley can be written as the followingLiu et al. (2013); Xiao et al. (2012); Fang et al. (2018a)

[TABLE]

where $\sigma$ matrices act on sublattice space and there is no spin degeneracy because of the large spin-orbital coupling. This hamiltonian generally applies for a large class of TMD materialsLiu et al. (2013); Xiao et al. (2012); Fang et al. (2018a). In this paper, we take WSe2 as a representative for TMD materials. For WSe2, it is estimated that $\hbar v_{F}/a\cong 1.1$eV and $m\cong 1.2$eV. Similar as the graphene case, with strains, the Dirac points are shifted and their locations in momentum space are again given by Eq. 5 and 4. It is found from first principle that the drift parameter $\beta$ is around $2.3$ for WSe2Rostami et al. (2015); Fang et al. (2018b). Of course, it is hard to pin down the exact value of the $\beta$ and it also varies for different TMD materials. Therefore, we treat $\beta$ as a potential tuning parameter to extract general features of the moiré band structure.

In this section, we consider bilayer transition metal dichalcogenide system with only heterostrains for simplicity. It turns out that the bilayer TMD system is a desirable platform for engineering flat bands with strains. We restrict ourselves to a class of heterostrain which is called volume preserving strain. The general form is given by,

[TABLE]

Physically, $S_{v}$ describes that the material is strained by $\epsilon$ along $\varphi$ direction and by $-\epsilon$ along the normal direction. We consider the bilayer system starting from AA stacking configuration, where the inversion symmetry is not present. The AB stacking configuration is an interesting case which we will address in a separate paper. Applying a generic $S_{v}$ explicitly breaks all the point group symmetries of the system.

For bilayer WSe2 moiré superlattices starting from AA stacking configuration, we need to include both the interlayer tunneling and the intralayer potential in Eq. 8. The interlayer tunneling $T(\mathbb{r})$ has the same form as Eq. 11. It is estimated in Ref.Wu et al. (2018) that, for bilayer WSe2, the tunneling parameters are $u\cong 1.1$meV, $u^{\prime\prime}\cong 9.7$meV and $u^{\prime}$ approximately zero because of the large band gap. The intralayer potential has the following formWu et al. (2018)

[TABLE]

where $\mathbb{g}_{3}=-(\mathbb{g}_{1}+\mathbb{g}_{2})$. For WSe2, the parameters are $V_{c}\cong 6.8$meV, $V_{v}\cong 8.9$meV, $\phi_{c}\cong 89.7^{\circ}$ and $\phi_{v}\cong 91^{\circ}$. These data vary slightly for other TMD materials. Thus, we will focus on generic features of the electron band structure which do not rely much on the precise values of these parameters.

We show the band structure for bilayer WSe2 with heterostrain calculated from the continuum model in Fig. 9. We find generic flat bands near the top of the valance band for small strain magnitude. The top two valence bands are close in energy and have very small bandwidth compared to the large gap separating them to the rest of the spectrum. The third band is also quite flat and energetically separated from other bands. Increasing strain can enlarge the bandwidth and band gap.

The strain angle $\varphi$ provides an additional nob for engineering the band structure. Now we focus on the top two valence bands. In Fig. 10 (c), we show the phase digram as a function of $\beta$ and $\varphi$ (the phase diagram is not sensitive to the magnitude of the strain $\epsilon$ for small $\epsilon$). Interestingly, we find an alternating pattern of topological and trivial bands as a function of $\varphi$. The threefold periodicity of the phase diagram in $\varphi$ is due to the emergent threefold rotational symmetry of interlayer coupling (11) and intralayer potential (14). We show the band structure of a critical point at $\beta_{c}\cong 2.10$ and $\varphi=0$. The topological phase transition closes the gap with a Dirac dispersion at one of the $M$ points in the Brillouin zone.

It is worth to mention the displacement field is a good way to control the relative position of the first three bands. With considerable interlayer bias, we can make one band well separated from other bands on the top of the spectrum. In addition, this band is topologically trivial. An example of such band structure is shown in Fig. 10 (d). As expected, the wavefunction of the top band is mostly polarized in one layer.

IV.2 Pressure tunned topological transition

Pressure on the bilayer systems changes the interlayer distance and in return adjusts the relative ratio between the strength of the interlayer tunneling and intralayer potential energy, which is another way to design the moiré band structure. To demonstrate this effect, we consider twisted or heterostrained bilayer TMD with pressure. Pressure can be modeled by introducing a phenomenological parameter $p$ in the model Eq. 8, $T(\mathbb{r})\rightarrow pT(\mathbb{r})$, which represents the relative ratio of interlayer tunneling and intralayer potential. For twisted bilayer TMD at twist angle $\theta\cong 1^{\circ}$, we find multiple topological phase transitions between the top three moiré bands as we increase $p$ from $1$ to $1.5$ as shown in Fig. 11. At $p\cong 1.195$, there is a quadratic band touching at $\Gamma$ in moiré Brillouin zone between the second and third band which shifts the Chern numbers of the two bands by $\pm 2$. At $p\cong 1.41$, the top two bands touch at $\Gamma$ with a Dirac dispersion. The top band becomes topologically trivial for $p>1.41$. In this system, we can get topological flat bands with chern number $\pm 2$. Similar effects also appear in heterostrained cases.

IV.3 Effective model for the flat bands

With small heterostrains, for instance $\epsilon\cong 1.5\%$ in Fig. 9, for each valley the two top valence bands which are close in energy and very flat (total bandwidth $W\cong 1.4$meV), and they are well separated from the rest of the bands (gap $\Delta\cong 8$meV). One can estimate the Coulomb interaction scale on the moiré superlattice $V=e^{2}/(4\pi\epsilon a_{M})\cong 4.2$meV, with static dielectric constant of WSe2 $\epsilon\cong 15.3$Laturia et al. (2018) and moiré superlattice constant $a_{M}\cong 67a\cong 23.4$nm for $\epsilon_{xy}=1.5\%$. The parameters of the system are in the limit $\Delta\gg V\gg W$. In addition, the total Chern number for the top two valence bands is trivial. Therefore, we can propose an effective model just for the top two bands. Taking into account of the valley/spin degeneracy, the simplest guess is a 2-band Hubbard model on a triangular lattice, which can be schematically written as the following

[TABLE]

where $n_{i,\alpha}=\sum_{v=\pm}c_{i,\alpha,v}^{\dagger}c_{i,\alpha,v}$. The two valleys are found to be degenerate in energy, therefore, have an approximate $SU(2)$ symmetry in the small intervalley-scattering limit. The hopping terms are spatially anisotropic depending on the heterostrain. Half filling the two bands may lead to interesting correlated states, such as valley/spin polarized state, which spontaneously breaks time reversal symmetry.

Aside from the top two valence bands, we can also focus on the third band which is also quite flat and energetically separated from other bands. The Chern number is again trivial for this band with current parameters. The effective model for this band will be a single band triangular lattice $SU(2)$ Hubbard model with anisotropic hopping terms. Similar model could also describe the physics of the displacement field separated flat band in Fig. 10 (d). This model is previously suggested to be the effective model of the quasi-two-dimensional organic charge transfer salts, which potentially host exotic spin liquid statesKanoda and Kato (2011); Powell and McKenzie (2011).

Finally, we emphasis that the Coulomb interaction scales with the strain as $V\sim 1/a_{M}\sim\epsilon$ while the bandwidth of the moiré bands scales with the strain as $W\sim(1/a_{M})^{2}\sim\epsilon^{2}$. Therefore, different strains can produce Hubbard models with different $t/U$ ratio.

V Discussion

In this paper, we explore the effects of uniaxial heterostrain in twisted bilayer graphene. It is clear that uniaxial heterostrain generically broadens the bandwidth of the nearly flat bands at magic angle, which could be the reason for the observed large energy separation of van Hove singularities in the STM experimentKerelsky et al. (2018). Large bandwidth may not be desirable for correlation physics. Therefore, future experiments should try to minimize the heterostrain between layers. On the other hand, heterostrain could be helpful to generate higher order van Hove singularities, which is a good platform to enhance correlation effects. A systematic study of the evolution of the van Hove singularities is an interesting but involved subject. Another natural question to ask is whether there are situations in which heterostrain helps to flatten the band dispersions rather than to broaden the bands. To answer these questions, we need more comprehensive studies where optimization methods such as machine learning could be helpful. On the other hand, a full analytical understanding of the strain effect is also demanded for future study.

The other effects of heterostrain are shifting the energies of the two Dirac crossings within one valley and greatly enhancing the Dirac velocity. This could potentially explain the observed 4-fold Landau level degeneracy in magneto transport experimentsCao et al. (2018a, b); Yankowitz et al. (2019). One natural consequence of such scenario is that there will be finite electron/hole fermi surfaces at the charge neutrality point. Nonetheless, due to the enhanced Dirac velocity, the fermi surfaces are estimated within our continuum model to be really small and could be difficult to resolve with current experimental precision. To get a more accurate estimation theoretically, one needs to include various subtle effects such as lattice relaxations, which is beyond the scope of the current paper but a good subject for future investigations.

We also show heterostrain is a good way of creating and tuning flat bands in bilayer TMD systems. There is no requirement of a “magic” strain here. The band structure is easily controlled by the strain direction, pressure and displacement field.

Acknowledgement

We thank Abhay Pasupathy, Cory Dean, Matt Yankowitz and Brian LeRoy for stimulating discussions and communications. ZB is supported through Pappalardo fellowship at MIT. NY and LF are supported by DOE Office of Basic Energy Sciences, Division of Materials Sciences and Engineering under Award DE-SC0018945. LF is partly supported by the David and Lucile Packard Foundation.