Local Optical Conductivity of Bilayer Graphene with Kink Potential
aa r X i v : . [ c ond - m a t . m e s - h a ll ] F e b Local Optical Conductivity of Bilayer Graphene with Kink Potential
Zhen-Bing Dai,
1, 2
Zhiqiang Li, and Yan He ∗ College of Physics, Sichuan University, Chengdu, Sichuan 610064, China Department of Physics, Sichuan Normal University, Chengdu, Sichuan 610066, China
We study the optical response of bilayer graphene with a kink potential composed of a domain wallseparating two AB regions with opposite interlayer bias. The band structure and the local opticalconductivity in real space are investigated in details based on a continuum model. We find that theone-dimensional chiral states localized at the domain wall contribute significantly to the local opticalconductivity, which shows a clear distinction in different regions. The effects of domain wall stateson optical conductivity can be detected by spatially and frequency resolved spectroscopic features.From the spectrum at various Fermi energies, important features in the band structure such as theenergy separation between two chiral states can be directly measured. When the domain wall regionis broad, the spatial distribution of local optical conductivity can provide important information onthe bound states as well as the topological domain wall states.
I. INTRODUCTION
In bilayer graphene(BLG), the chiral counterpropagat-ing boundary modes occur at the interface between theAB and BA stacking configurations [1–4], or at the linejunction of two oppositely gated regions [5–8], becausethe Berry curvatures in these two regions near the do-main wall(DW) have opposite signs [9–12]. These edgestates are topologically protected because the inter-valleyscattering is almost absent [13]. A triangular networkof topologically protected layer stacking DW states hasbeen found in twisted BLG with a tiny twist angle of θ = 0 . ◦ [14], due to the atomic scale lattice recon-struction [15]. The area of the DW regions and AB andBA domains become comparable as the twist angle in-creases to a magic angle( θ c = 1 . ◦ ), in which super-conductivity and Mott insulator-like states [16, 17] havebeen observed [1, 18].Compared to the DW between different layer stack-ings, DW due to a kink potential in BLG could offer arobust and flexible platform to realize topological andvalleytronic operations [13, 19–21], such as long lived 1DDW plasmons [12] , waveguides [22], quantum valley Hallstates [23] and valley filters [24]. In addition, net valleycurrents can be induced by the counter-propagating con-ducting channels of the gated BLG [25–27], which can bedetected directly by the measurements of the valley Halleffect. Moreover, the tunability of the dual-split-gatingconfiguration has been conveniently utilized in real de-vices. Each gate pair can operate independently and in-duce DW states within the split-gate region, which pro-vides many opportunities for potential device applica-tions.Optical methods are widely used to quantitativelystudy the electronic band structure and physical prop-erties of materials. It is reported [23, 28–31] that thesurface plasmon reflection at DW solitons can be stud-ied in detail by the near-field infrared nanoscopy for the ∗ Electronic address: heyan [email protected] reversal stacking configurations in BLG. The DW plas-mon also displays a rich phenomenology including tun-able subwavelength confinement of light, higher prop-agation length, plasmon oscillation over a wide rangeof frequencies as well as supporting the propagation oftransverse electric modes [12, 32, 33]. Meanwhile, opti-cal conductivity measurements can provide a feasible wayto probe excitation spectrum of BLG between the gaplessstates and the continuum bands [34]. Up to now, mostworks focus on the electronic structure and the frequencydependence of the optical conductivity of homogeneousgraphene [32, 35–38]. The local optical conductivity ofBLG is strongly sensitive to the atomic-scale stackingorder [28] and the band structure which associated withelectric-field-induced energy asymmetry [32], leading tothe inhomogeneous distributions across DW. To the bestof our knowledge, a systematic theoretical study of theoptical properties of opposite gated BLG is still lacking.In this paper, we will provide a comprehensive way toinvestigate the band structure and the local optical con-ductivity of opposite gated BLG.We start with the four-band continuum model of BLG,and discuss the calculation method of electronic struc-ture, topological properties and the related local opticalconductivity for the opposite gated BLG. In the numer-ical calculation, we provide a method to avoiding thefermion doubling problem that is encountered in the lat-tice Hamiltonian. We analyze the band structure of thegated BLG and the wave functions of a few selected gap-less modes. The origin of the gapless modes can be dis-tinguished from the profiles of the wave function. In theend, we present the spatial and frequency distributionsof the real part of the local optical conductivity acrossDW with different chemical potentials.
II. THEORETICAL BACKGROUND
In order to derive and discuss the distributions of theoptical conductivity in real space of BLG under the asym-metry electronic potential, we will examine the bandstructure of gated BLG. To this end, we begin with acontinuum model [5, 35] of BLG with a bias potentialapplied to these two layers.The Hamiltonian is given by H ξ = − V τ ⊗ σ + v F k x τ ⊗ σ + ξv F k y τ ⊗ σ + t ⊥ τ ⊗ σ + τ ⊗ σ ) , (1)where ( τ , τ , τ ) and ( σ , σ , σ ) are Pauli matrices act-ing on the layer space and sublattice space, respectively.In addition, τ and σ are the corresponding identity ma-trices, the band velocity of the Dirac cone v F = √ aγ ~ is expressed in terms of the in-plane nearest neighbor in-teraction γ = 2 . a ≈ .
246 nm, and t ⊥ = 0 .
39 eVdenotes the interlayer coupling term between adjacentlayers. The corresponding wave-function is arranged as ψ = ( ψ A , ψ B , ψ A , ψ B ) T where A, B label two sublat-tices and 1 , Kξ valley with ξ = ±
1. Theadvantage of the continuum model is that we can focuson a single Dirac cone at one valley without mixing ef-fects of two inequivalent valleys. For simplicity, we set ~ = 1 hereafter.We first consider the topological properties of gatedBLG case with constant bias voltage. In this case,one can directly diagonalize the Hamiltonian in momen-tum space to find out four eigenvectors ψ m ( k ), where m = 1 , · · · , A mµ ( k ) = i h ψ m | ∂∂k µ | ψ m i and F m ( k ) = ∂A my ∂k x − ∂A mx ∂k y , with µ = x, y .Furthermore, the Chern number for the lower two occu-pied bands is given by C = 12 π X m =1 , Z d kF m ( k ) , (2)where the integration is over the whole infinite 2D k-space. In the numerical calculations, a fairly large mo-mentum cutoff is good enough to guarantee that the re-sulting Chern number is quantized as expected. For themodel of Eq.(1), it is found that the Chern number itsvalley dependent, which reads C Kξ = ξ sign( V ) (3)where sign( x ) = ± x , and itmeans that C K = ± C K its valleydependent, it is also called the valley Chern number [11].For the DW region formed by potential, the electricfield usually depend on spatial position, V ( x ). Fig. 1illustrates schematically the topological chiral modes inBLG. In this paper, we mainly focus on a pair of kink andanti-kink DWs in the BLG for one valley K, as a resultof requirement of the periodic system. For a fixed layer,the electric potential is initially positive and switch to Figure 1. Schematic representation of gated BLG for the cre-ation of a kink potential. The vertical electric field generatedby the gated voltage applied to the upper and lower layerswith opposite sign, is indicated by the black arrows. Red andgreen wavy arrows represent counterpropagating chiral modeslocalized to the domain wall. negative at the kink DW. Then it changes back to posi-tive at the anti-kink DW. This potential can be writtenas V ( x ) = V [tanh( x − x w ) − tanh( x − x w ) − V isthe overall magnitude of bias voltage and w denotes thewidth of region in which the potential switches its signin each layer. Since the translational symmetry is pre-served along the y direction, we can keep the momentum k y . But we have to represent k x as a differential operator − i∂ x in real space along the x direction. For concrete-ness, we assume that the system is located inside a finiteinterval − L/ < x < L/ ψ ( − L/
2) = ψ ( L/ − i∂ x will changesign under complex conjugation, we find that particle-hole symmetry is broken in this real space model.The above 4-band continuum model has particle-holesymmetry. In this form, it is easy to see that U † H ∗ U = − H, (4)with U = τ ⊗ σ . This means that if | u n i sat-isfies H | u n i = E n | u n i , one also has H ( U | u n i ∗ ) = − E n ( U | u n i ∗ ). Therefore the spectra are symmetricabout E = 0.In numerical calculation, one can discretize the x in-terval with N lattice sites. Then the 4-band model canbe written as H = − τ ⊗ σ ⊗ V m − iv F τ ⊗ σ ⊗ h + v F k y τ ⊗ σ ⊗ h + t ⊥ τ ⊗ σ + τ ⊗ σ ) ⊗ h (5)Here V m = V ( x i ) δ ij , h = δ i,j and h = x ( δ i,j − − δ i,j +1 ) with i, j = 1 , · · · , N and ∆ x = L/N . Note thatin Eq.(5) the derivative is approximated by a symmet-ric difference operator as ∂ x → h , thus we can call itsymmetric difference model. As discussed before, thedifferential operator breaks the particle-hole symmetryand one might expect that the spectrum is asymmet-ric about E = 0. But the use of symmetric differenceactually introduce an extra symmetry to the system andmake its spectrum symmetric again. To be more explicit,one can introduce a spatial symmetric transformation as U = ( − i +1 δ i,N − j +1 . It is easy to verify the followingproperties as U † V m U = − V m and U † h U = h . Withthis matrix U , one can show the Hamiltonian of Eq.(5)has a generalized particle-hole symmetry as U † H ∗ U = − H, (6)with U = τ ⊗ σ ⊗ U . Therefore, the Hamiltonianof Eq.(5) also has a spectrum symmetric about E = 0.The numerical results suggest that there appears sev-eral chiral gapless modes localized at the DW. Due tothe above generalized particle-hole symmetry, there willbe two branches of chiral modes propagating along thetwo opposite directions. The reason for the appearingof two branch of fermions with opposite chirality is theso-called fermion doubling problem first encountered inlattice Quantum Chromodynamics [40, 41].To avoid this fermion doubling problem [42–44], wefollow the suggesting of Stacey [45] to replace the sym-metric difference h by the the forward difference h = x ( δ i,j − − δ i,j ). Consequently, the 4-band BLG modelwith a lattice realization along x -direction can be writtenas H = − τ ⊗ σ ⊗ V m − iv F τ ⊗ σ + ⊗ h + iv F τ ⊗ σ − ⊗ h T + v F k y τ ⊗ σ ⊗ h + t ⊥ τ ⊗ σ + τ ⊗ σ ) ⊗ h , (7)where we defined σ ± = ( σ ± iσ ) / σ ab ( x , x ) = ig s π Z d k X m = n f ( E m ) − f ( E n )( E m − E n )( ω + iη − ( E m − E n ))+ X m df ( E m ) dE m ω + iη ! h u m ( x , k y ) | v a | u n ( x , k y ) i×h u n ( x , k y ) | v b | u m ( x , k y ) i , (8)where g s = 2 denotes the spin degeneracy in graphene, f ( E m ) is the Fermi-Dirac distribution f ( E m ) = 1 / (1 + e ( E m − µ ) / ( k B T ) ). The velocity operator can be obtainedfrom v a = ∂H∂k a where the spatial indices a, b = x, y . Forthe model of Eq.(5), we have v x = v F τ ⊗ σ . The localoptical conductivity can be obtained as σ ab ( x ) = Z dx ′ σ ab ( x, x ′ ) (9) III. RESULTS AND DISCUSSIONSA. Band structure
As discussed in the previous section, the valley Chernnumbers of the BLG are ± − k y . The parameters used in this cal-culation is displayed in the figure caption. Note thatthere are two pairs of mid-gap states (the red and bluesolid lines) connecting the valence and conduction band.Since the two pairs of gapless modes propagate in theopposite direction, it seems that there may be a fermiondoubling problem. To investigate the physical origin ofchiral states, we plot the real part of the wave functionsof sublattice A (black dashed curves) and B (red dottedcurves) and the probability density for the states indi-cated by the pink points in Fig. 2(a), corresponding to k y = 0 . − . The wave functions are the eigenvectorsobtained by diagonalization of the Hamiltonian Eq.(7).According to the previous results [5, 6], we notice thatthe wave function have prominently A1 and B2 characterat low energies and constant gate voltage. Obviously, theelectron states are localized at the left DW region with x = x = −
75 nm, as shown in Fig. 2(1a) and (3a). Inother word, the chiral states at the left DW region arerepresented by the blue curves in panel (a). Meanwhile,the red curves also indicate the localized states stay atthe right DW region with x = x = 75 nm, in Fig. 2(2a)and (4a). Therefore, these gapless modes actually are lo-calized at the different DW regions with x = x and x and there is indeed no fermion doubling problem, whichcorresponds to the two blue and red curves, respectively.Since the two DW regions have opposite orientations, thetwo pairs of gapless modes should propagate in the oppo-site directions. In addition, the propagation direction ofthe gapless states can be reversed by reversing the signof electric field V ( x ), are indicated by the red and bluecurves in plane (a).We notice that the peak shape of probability density, | Ψ | = | Ψ A | + | Ψ B | , in panel (1a) and (4a) is slightlydifferent from that of panel (2a) and (3a). The peaksof (1a) and (4a) are wider and almost split into doublepeaks. This is because that the energy of the states of(1a) and (4a) are much more closer to the continuumspectrum than that of (2a) and (3a). As the topologicalstates approach to the bulk states, the confining potentialis weakened and the wave function is less localized thanbefore.Next, we investigate the energy spectrum for a broader -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5-0.15-0.10-0.050.000.050.100.15 E ne r g y ( e V ) k y (nm -1 ) (1a)(2a)(3a)(4a)(a) -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5-0.15-0.10-0.050.000.050.100.15 k y (nm -1 ) E ne r g y ( e V ) (b)(1b) (2b)(3b) (4b) -150 -100 -50 0-0.10.00.1 x (nm)j B |y | j A (1a) x (nm)(2a) -150 -100 -50 00.000.050.100.15 x (nm)(3a) x (nm)(4a) -150 -75 0 75 1500.00.10.2 x (nm)j A j B |y | (1b) -150 -75 0 75 150-0.2-0.10.00.1 x (nm)(2b) -150 -75 0 75 150-0.2-0.10.00.10.2 x (nm)(3b) -150 -75 0 75 150-0.10.00.10.2 x (nm)(4b) Figure 2. Energy spectrum, wave functions and probability densities of a kink-antikink profile with V = 0 . V = 0 . w = 2 nm. The black dashed line indicates k y = 0 . − . (b): Energy spectrum of BLG with a kink-antikink potential for w = 15 nm. The black dashed line indicates k y = 0 .
05 nm − . (1a)-(4a): The real parts of the wave functions and probability densities of the states corresponding to the fourpink points in panel (a). (1b)-(4b): The real parts of the wave functions and probability densities of the states correspondingto the four pink points in panel (b). The blue dashed lines indicate the center of the two DWs. kink potential profiles, i.e. w = 15 nm. Fig. 2(b) showsthe band structure of gated BLG with V = 0 . k y = 0 .
05 nm − in Fig. 2(1b)-(4b). The bound states are also locatedaround x = x or x , but these states are more extendedand their amplitudes cross zero near the center of theDW region. In panel (3b) and (4b), we find that thesebound states are localized at both kink and anti-kinkDW regions since these two states are very close to eachother. While the chiral states are only confined at oneDW region as shown in (1b) and (2b). Meanwhile, theprobability distribution of bound states is quite differentfrom that of topological states since these two types ofstates are quite different topological features. As pointedout in [6], the bound states are no longer unidirectionaland non-chiral. They have a quasi-1D free-electron typeof spectrum and are less localized at the DW region ascompared to the chiral states. Their probability distri-bution appears similar to those of excited states. B. Optical conductivity
We now present the results of the local longitudinal op-tical conductivity σ xx of BLG with a DW like bias, whichis computed by the Kubo formula of Eq.(8) and (9). Fig.3(a) shows the frequency dependence of the real part ofthe local optical conductivity Re σ xx in gated BLG withthe parameters w = 2 nm, V = 0 . µ = 0. Theoptical conductivity are scaled by σ = e / ~ , which isfour times the background optical conductivity of mono-layer graphene. This includes the spin degeneracy factor g s = 2 in graphene since the charged impurity scatteringdoes not mix spins. We notice that the longitude opticalconductivity of the DW region is quite different from ofthat AB region. A strong peak in the optical conductiv-ity away from DW region is visible near ω = 0 . ω > . σ xx of gated BLGin the DW and AB region as a function of frequency andspatial coordinate at a series of fixed chemical potential.Fig. 3(b) shows the spatial dependence of the local op-tical conductivity for several incident photon frequenciesat charge neutral point. As illustrated in panel (a), thereis a prominent peak located at the DW region at ω < . ω increases, the contrast of the optical conductivity be-tween the DW and AB region become much smaller, asindicated by the blue curve in panel (b). At ω = 0 . . σ , which is twice the optical conductivity ofmonolayer graphene. For a non-zero chemical potential µ = 0 .
07 eV, where the bulk states are occupied, the de-pendence of local optical conductivity mentioned aboveis similar as shown in Fig. 3(c).Interestingly, the value of optical conductivity at cen-ter of DW region varies very little. In order to obtainmore information, we plot the optical conductivity at R e s xx / s x (nm)(c)m=0.07 eV R e s xx / s w (eV)(d)DW R e s xx / s w (eV)(e)AB R e s xx / s x (nm)(b)m=0 Figure 3. The spatial and frequency dependence of the lo-cal optical conductivity of BLG with DW width w = 2 nm.(a) The real part of the local optical conductivity in unit of σ = e / ~ as a function of spatial coordinate and frequency.The intensity of optical conductivity is indicated by the color.(b) The spatial distribution of the real part of the optical con-ductivity for ω = 0 . , . , . µ = 0. (c) The sameas (b) but with µ = 0 .
07 eV. (d) The real part of the opticalconductivity at the center of DW region as a function of fre-quency for the different chemical potentials µ = 0 , . , . V = 0 . t ⊥ = 0 . T = 300 K, and the damping rate is taken to be η = 3meV. DW region as a function of incident photon frequency forseveral typical chemical potentials in panel (d). At lowenergy, there are a shoulder and peak at ω = 0 . , . ω = 0 . ω = 0 . µ is increasing but still inside the band gap,the optical conductivity curve will stay almost the same.On the other hand, when µ is outside the band gap, thesetwo peaks of optical conductivity of the DW region disap-pear. Therefore, the frequency dependence of the opticalconductivity can give us more evidence of localized statesat DW region.Fig. 3(e) shows the optical conductivity of BLG at theuniform AB region as a function of frequency for the dif-ferent chemical potentials. We will call Re σ xx versus ω curve as the conductivity spectrum. The absorptionpeak appears at ω = 0 . . ω = 0 .
39 eV, which is consistent withthe inter-layer coupling constant t ⊥ = 0 .
39 eV. Whenchemical potential enters into bulk bands, for example µ = 0 .
07 eV, the absorption peak of optical transitionassociated with the band gap is significantly decreased,but the band gap can still be clearly determined. More-over, the intensity of optical conductivity associated withthe inter-layer coupling is significantly increased, with apeak value of roughly σ . Therefore, we can attributethe huge difference between the conductivity spectrumof the DW and AB regions to the different topologicalproperties of eigenmodes.Next, we increase the width of the DW region andinvestigate how the local optical conductivity spectrumchanges. As mentioned before, for wide DW region, addi-tional 1D bound states appear. These bound states willaffect the spatial distribution of the local optical conduc-tivity because they are less strongly localized at the DWregion as compared to the chiral states. To see this differ-ence more clearly, we also plot the frequency and spatialdependence of the optical conductivity of the gated BLGwith parameters w = 15 nm, V = 0 . µ = 0 inFig. 4. In panel (a), we find that there is a transition froma single peak at ω = 0 .
05 eV to twin peaks line shape atfrequency around 0 .
05 eV < ω < . ω = 0 . µ = 0 .
07 eV, as in Fig. 4 (c). In otherword, we can tell if a bound state is present by the ap-pearance of twin peak line shape from the spatial resolvedspectrum of local optical conductivity at one particularfrequency.Furthermore, we also show the local optical conductiv- m=0 w=0.05eVw=0.08 eVw=0.1 eV R e s xx / s x (nm)(b) m=0.07 eV w=0.05eVw=0.08 eVw=0.1 eV R e s xx / s x (nm)(c) R e s xx / s w (eV)(d)DW R e s xx / s w (eV)(e)AB Figure 4. The spatial and frequency dependence of the lo-cal optical conductivity of BLG with DW width w = 15 nm.(a) The real part of the local optical conductivity in unit of σ = e / ~ as a function of spatial coordinate and frequency.The intensity of optical conductivity is indicated by the color.(b) The spatial distribution of the real part of the optical con-ductivity for ω = 0 . , . , . µ = 0. (c) The sameas (b) but with µ = 0 .
07 eV. (d) The real part of the opticalconductivity at the center of DW region as a function of fre-quency for the different chemical potentials µ = 0 , . , . V = 0 . t ⊥ = 0 . T = 300 K, and the damping rate is taken to be η = 3meV. ity spectrum at DW region (panel (d)) and AB region(panel (e)) for several different chemical potentials. Wenote that there is a small peek near ω = 0 .
09 eV due tothe optical transition between the bound states when thechemical potential stays in the band gap. When µ = 0 . IV. CONCLUSION
In this work, we have investigated the electronic struc-ture and the spatial distribution of optical conductivitynear the DWs formed in the reversal gated BLG. We havedemonstrated that the propagating direction of the chi-ral modes, confined to DW, depend on the sign of electricfield and valley index. The origin of the gapless statescould be deduced from the amplitude profiles of the wavefunctions of the gapless modes. In this way, we have de-termined that the two topological edge modes are locatedat the DW. Moreover, we find that the conductivity spec-trum of the DW and AB regions are distinctly differentas a result of the different topological states. When thephoton energy below the band gap of system, the locallongitudinal optical conductivity of BLG with a reversalbias display a prominent peak located at DW. The en-hancement of optical conductivity at DW is owing to theemergence of the topological states. When the energyincreases and exceeds the band gap, a valley-like pro-file of optical conductivity appears at DW region. For asmooth kink potential, twin peaks profile in local opticalconductivity revealed that bound states are less stronglylocalized at DW as compared to the chiral states. There-fore, important features in the band structure such as energy separation between chiral states or bound states,localization of these states and band gap of system canbe easily estimated by the profile of local optical con-ductivity. From our method, we can also describe thespatial distribution of the conductivities for a given val-ley instead of the combination of the contributions fromthe two different valleys.The properties of optical conductivity in real spacenear DW region can help us to better understand the phe-nomena associated with topology, such as the propaga-tion of topological DW plasmons [12]. Our results presentan analysis of the local optical properties of gapped BLGthat not only can provide useful insights for designinglow-power topological quantum devices [13, 19, 47], butalso might open up a novel way for exploring of the fas-cinating edge-state and domain-wall physics in BLG.
ACKNOWLEDGMENTS
We are grateful to Prof. Hao-Ran Chang for valu-able discussions. This work was supported by the Na-tional Natural Science Foundation of China under GrantNo. 11874271 and 11874272. Y.H. is also supported byScience Specialty Program of Sichuan University underGrant No. 2020SCUNL210. We thank the High Perfor-mance Computing Center at Sichuan Normal University. [1] R. Bistritzer and A. H. MacDonald, Proc. Natl. Acad.Sci. USA , 12233 (2011).[2] A. Vaezi, Y. Liang, D. H. Ngai, L. Yang, and E.-A. Kim,Phys. Rev. X , 021018 (2013).[3] L. Ju, Z. Shi, N. Nair, Y. Lv, C. Jin, J. Velasco,J., C. Ojeda-Aristizabal, H. A. Bechtel, M. C. Martin,A. Zettl, et al., Nature , 650 (2015).[4] N. R. Walet and F. Guinea, 2D Mater. , 015023 (2019).[5] I. Martin, Y. M. Blanter, and A. F. Morpurgo, Phys.Rev. Lett. , 036804 (2008).[6] M. Zarenia, J. M. Pereira, G. A. Farias, and F. M.Peeters, Phys. Rev. B , 125451 (2011).[7] J. Li, K. Wang, K. J. McFaul, Z. Zern, Y. Ren, K. Watan-abe, T. Taniguchi, Z. Qiao, and J. Zhu, Nat. Nanotech-nol. , 1060 (2016).[8] H. Chen, P. Zhou, J. Liu, J. Qiao, B. Oezyilmaz, andJ. Martin, Nat. Commun. , 1202 (2020).[9] W. Yao, S. A. Yang, and Q. Niu, Phys. Rev. Lett. ,096801 (2009).[10] F. Zhang, J. Jung, G. A. Fiete, Q. Niu, and A. H. Mac-Donald, Phys. Rev. Lett. , 156801 (2011).[11] F. Zhang, A. H. MacDonald, and E. J. Mele, Proc. Natl.Acad. Sci. USA , 10546 (2013).[12] E. H. Hasdeo and J. C. W. Song, Nano Lett. , 7252(2017).[13] Z. Qiao, J. Jung, Q. Niu, and A. H. Macdonald, NanoLett. , 3453 (2011).[14] S. S. Sunku, G. X. Ni, B. Y. Jiang, H. Yoo, A. Stern-bach, A. S. McLeod, T. Stauber, L. Xiong, T. Taniguchi, K. Watanabe, et al., Science , 1153 (2018).[15] F. Gargiulo and O. V. Yazyev, 2D Mater. , 015019(2017).[16] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken,J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe,T. Taniguchi, E. Kaxiras, et al., Nature , 80 (2018).[17] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi,E. Kaxiras, and P. Jarillo-Herrero, Nature , 43(2018).[18] K. Kim, A. DaSilva, S. Huang, B. Fallahazad, S. Lar-entis, T. Taniguchi, K. Watanabe, B. J. LeRoy, A. H.MacDonald, and E. Tutuc, Proc. Natl. Acad. Sci. USA , 3364 (2017).[19] Z. Qiao, J. Jung, C. Lin, Y. Ren, A. H. MacDonald, andQ. Niu, Phys. Rev. Lett. , 206601 (2014).[20] J. Jung, F. Zhang, Z. Qiao, and A. H. MacDonald, Phys.Rev. B , 075418 (2011).[21] H. M. Abdullah, A. El Mouhafid, H. Bahlouli, and A. Jel-lal, Mater. Res. Express , 025009 (2017).[22] F.-M. Zhang, Y. He, and X. Chen, Appl. Phys. Lett. ,212105 (2009).[23] L. J. Yin, H. Jiang, J. B. Qiao, and L. He, Nat. Commun. , 11760 (2016).[24] A. Rycerz, J. Tworzyd lo, and C. W. J. Beenakker, Nat.Phys. , 172 (2007).[25] R. V. Gorbachev, J. C. W. Song, G. L. Yu, A. V. Kre-tinin, F. Withers, Y. Cao, A. Mishchenko, I. V. Grig-orieva, K. S. Novoselov, L. S. Levitov, et al., Science , 6208 (2014). [26] Y. Shimazaki, M. Yamamoto, I. V. Borzenets, K. Watan-abe, T. Taniguchi, and S. Tarucha, Nat. Phys. , 1032(2015).[27] J. Lee, K. Watanabe, T. Taniguchi, and H. J. Lee, Sci.Rep. , 6466 (2017).[28] B. Y. Jiang, G. X. Ni, Z. Addison, J. K. Shi, X. Liu,S. Y. F. Zhao, P. Kim, E. J. Mele, D. N. Basov, andM. M. Fogler, Nano Lett. , 7080 (2017).[29] L. Jiang, Z. Shi, B. Zeng, S. Wang, J. H. Kang, T. Joshi,C. Jin, L. Ju, J. Kim, T. Lyu, et al., Nat. Mater. , 840(2016).[30] Z. Fei, A. S. Rodin, W. Gannett, S. Dai, W. Regan,M. Wagner, M. K. Liu, A. S. McLeod, G. Dominguez,M. Thiemens, et al., Nat. Nanotechnol. , 821 (2013).[31] L. Ju, L. Wang, T. Cao, T. Taniguchi, K. Watanabe,S. G. Louie, F. Rana, J. Park, J. Hone, F. Wang, et al.,Science , 907 (2017).[32] S. M. Farzaneh and S. Rakheja, J. Appl. Phys. ,153101 (2017).[33] S. A. Mikhailov and K. Ziegler, Phys. Rev. Lett. ,016803 (2007).[34] D. N. Basov, R. D. Averitt, D. van der Marel, M. Dressel,and K. Haule, Rev. Mod. Phys. , 471 (2011).[35] E. McCann and M. Koshino, Rep. Prog. Phys. , 056503 (2013).[36] P. Moon and M. Koshino, Phys. Rev. B , 205404(2013).[37] M. Koshino, Phys. Rev. B , 115409 (2013).[38] H. Min, D. S. L. Abergel, E. H. Hwang, andS. Das Sarma, Phys. Rev. B , 041406(R) (2011).[39] G. Trambly de Laissardiere, D. Mayou, and L. Magaud,Nano Lett. , 804 (2010).[40] J. B. Kogut, Rev. Mod. Phys. , 775 (1983).[41] M. Creutz, Rev. Mod. Phys. , 119 (2001).[42] J. Tworzyd lo, C. W. Groth, and C. W. J. Beenakker,Phys. Rev. B , 235438 (2008).[43] A. R. Hern´andez and C. H. Lewenkopf, Phys. Rev. B ,155439 (2012).[44] B. Messias de Resende, F. C. de Lima, R. H. Miwa,E. Vernek, and G. J. Ferreira, Phys. Rev. B ,161113(R) (2017).[45] R. Stacey, Phys. Rev. D , 468 (1982).[46] Y. Luo, R. Engelke, M. Mattheakis, M. Tamagnone,S. Carr, K. Watanabe, T. Taniguchi, E. Kaxiras, P. Kim,and W. L. Wilson, Nat. Commun. , 4209 (2020).[47] C. Cardoso, D. Soriano, N. A. Garcia-Martinez, andJ. Fernandez-Rossier, Phys. Rev. Lett.121