Time-evolution patterns of electrons in twisted bilayer graphene
TTime-evolution patterns of electrons in twisted bilayer graphene
V. Nam Do, ∗ H. Anh Le, and D. Bercioux
2, 3 Phenikaa Institute of Advanced Study (PIAS), Phenikaa University, Hanoi 10000, Vietnam Donostia International Physics Center (DIPC),Manuel de Lardizbal 4, E-20018 San Sebast´ıan, Spain IKERBASQUE, Basque Foundation of Science, 48011 Bilbao, Basque Country, Spain
We characterise the dynamics of electrons in twisted bilayer graphene by analysing the time-evolution of electron waves in the atomic lattice. We perform simulations based on a kernel polyno-mial technique using Chebyshev polynomials; this method does not requires any diagonalisation ofthe system Hamiltonian. Our simulations reveal that the inter-layer electronic coupling induces theexchange of waves between the two graphene layers. This wave transfer manifests as oscillations ofthe layer-integrated probability densities as a function of time. For the bilayer case, it also causesa difference in the wavefront dynamics compared to monolayer graphene. The intra-layer spreadingof electron waves is irregular and progresses as a two-stage process. The first one characterisedby a well-defined wavefront occurs in a short time — a wavefront forms instead during the secondstage. The wavefront takes a hexagon-like shape with the vertices developing faster than the edges.Though the detail spreading form of waves depends on initial states, we observe localisation of wavesin specific regions of the moir´e zone. To characterise the electron dynamics, we also analyse thetime auto-correlation functions. We show that these quantities shall exhibit the beating modulationwhen reducing the interlayer coupling.
I. INTRODUCTION
Stacking two-dimensional (2D) materials is a novelmethod based on the lego -principle for creating newvan der Waals heterostructures with the well-controlledproperties. However, accordingly, to the principle of thismethod, the successive layers are only stacked verticallykeeping the same orientation of one to the other. Animportant step forward comes when allowing a changein the relative orientation of the different stacked lay-ers. The simplest system allowing this new stackingmethod is twisted bilayer graphene (TBG). This systemis composed of two graphene layers stacked within a gen-eral manner and has been receiving large considerationlately. It was predicted that twisting two graphene lay-ers allows a strong tuning its electronic properties.
In-terestingly, a very narrow isolated energy band aroundthe charge neutrality level may appear in the spectrumof TBG configurations with tiny twist angles. Recently,Cao et al. have experimentally demonstrated that thisnarrow flat band is responsible to several strongly corre-lated phases, including an unconventional superconduct-ing and a Mott-like phase.
Theoretically, it was shownby Zou et al. that there are obstructions involving thesymmetries of the TBG lattice in constructing effectivecontinuum and tight-binding models to characterise thedynamics of electrons occupying the flat band.
Generally, stacking two layered materials may resultin a system of reduced symmetry compared to the twoconstituent lattices. The atomic configurations of TBGcan be characterised by an in-plane vector τ and a twistangle θ defining, respectively, the relative shift and ro-tation between the two graphene lattices. However, itis shown that only the twist angle governs the com-mensurability of the stacking, regardless of the twist-ing center. In particular, the lattice alignment is commensurate only when the twist angle takes thevalues given by the formula θ = acos[(3 m + 3 mr + r / / (3 m + 3 mr + r )], in which m, r are positivecoprime integers. When the stacking is com-mensurate, the translational symmetry of the TBG lat-tice is preserved, but it usually defines a large unitcell, especially for small twist angles θ . The elec-tronic calculation for such TBG configurations by bruteforce diagonalization is therefore extremely expensive interms of computational resources. Furthermore, theelectronic calculations based on the time-independentSchr¨odinger/Kohn-Sham equation combined with theBloch theorem are not applicable for incommensurateconfigurations because of the loss of the lattice trans-lational invariance. In this work, we show that methodsbased on the time-dependent Schr¨odinger equation in realspace are a powerful alternative to treat the TBG systemof arbitrary twist angles.Following the time-evolution of wave packets in realspace is a useful technique to simulate the dynamics ofelectrons. This method was used for studying the case ofmonolayer graphene and special TBG configurations. Forinstance, Rusin and Zawadzki and Maksimova et al. used the kicked Gaussian wave packet to analyticallystudy the different features of the zitterbewegung motionof electrons in various carbon-based structures, includingcarbon nanotubes. In these works, the wave packets dy-namics was governed by an effective Dirac Hamiltonian,thus the discrete nature of the atomic lattice was nottaken into account. M´ark et al., however, described theevolution of the kicked Gaussian wave packet in a poten-tial field constructed from an atomistic pseudo-potentialmodel. This approach allows taking into account the dis-tortion of the Dirac cones at high energy, and thus showedthe anisotropic dynamics of electrons in the graphenelattice. In the tight-binding framework, Chaves et al. a r X i v : . [ c ond - m a t . m e s - h a ll ] A p r used the discrete Gaussian form to define a wave packetand showed some quantitative differences in the zitter-bewegung motion of electron described by the effectiveDirac model and the tight-binding description. Xian etal. also used the discrete Gaussian wave packet to sim-ulate the transport of electron in a particular commen-surate TBG. They showed the existence of six-preferabletransport directions along which the wave packets arenot broadened, these are along the direction perpendic-ular to the transport direction. Particularly, they dis-cussed the behaviour of the layer-integrated probabilitydensity in each layer. They interpreted its behaviour sim-ilar to the neutrino-like oscillations where the interlayercoupling plays the role of mixing Dirac fermions in eachlayer: the two neutrino flavours.It is well-known that it is the honeycomb structureof graphene as the chiral interlocking of two trianglesub-lattices responsible for all the peculiar properties ofgraphene and related systems. Accordingly, the elec-tronic properties of graphene can be described by using aformulation in terms of relativistic fermions. This for-mulation is the same one used by Schr¨odinger to showthe zitterbewegung phenomena as the result of the in-terference of states at positive and negative energies.
The two-component spinor structure of the low-energyelectron states in graphene is due to the unit cell of thehoneycomb lattice constituted by only two carbon atoms.Stacking two graphene sheets gives rise to the diversityof the TBG configurations. It is therefore natural to posethe question of how the manifestation of the atomic lat-tice structure on the dynamical behaviour of electrons,particularly in the TBG systems with the lack of trans-lational symmetry.In this work, we address the dynamics of electrons inthe real lattice of generic TBG configurations using thetight-binding approach and try to relate it to the lat-tice symmetries. Though the wave packet method hasbeen successfully used to demonstrate the optical anal-ogy of electrons in graphene, its definition dependsexplicitly on some parameters and therefore not able toprovide a full picture of the electronic properties of a sys-tem. Accordingly, we will analyse the time-evolution oflocalised electrons occupying the 2 p z orbitals of carbonatom instead of Gaussian wave packets, whose definitiondepends on a particular wave vector and an initial po-sition. Within this approach, we can study the changesin the evolution pattern of electron wave functions withrespect to the detail of the lattice structure. By artifi-cially tuning the value of the parameters encoding thehybridisation of the 2 p z orbitals between two graphenelayers, we study the role of the interlayer coupling onthe time-evolution of electron states. For studying thetime-evolution of a state, we use the formalism of thetime-evolution operator ˆ U ( t ), i.e. , | ψ ( t ) (cid:105) = ˆ U ( t ) | ψ (0) (cid:105) ;we employ the kernel polynomials method to approxi-mate ˆ U ( t ). This method is efficient and useful to workdirectly in the lattice space of TBG configurations witharbitrary twist angles. Technically, we use the Chebyshev polynomials of the first kind to approximate the opera-tor ˆ U ( t ). Our implementation scheme is efficient becauseit accounts for the recursive relations of these polynomi-als, and, as a matter of fact, we are never performinga numerical diagonalization of the system Hamiltonian.Within this method, we can incorporate the details ofdiscrete atomic lattice into the dynamical properties ofthe 2 p z electrons of the TBGs. We shall study the intra-layer development of the 2 p z orbitals and the transfer ofthe probability density from one graphene layer to theother. The local information of the dynamics is studiedin the time domain.The outline of this paper is as follows. In Sec II,we present an empirical tight-binding model which al-lows characterising the dynamics of the 2 p z electronsin different levels of hopping approximation, i.e. , thenearest-neighbour (NN), next-nearest-neighbour (NNN)and next-next-nearest-neighbour (NNNN), and we alsopresent the method for the investigation of the time-evolution of the states as well as the calculation forseveral physical quantities characterising the dynamicsof electrons. In Sec. III, we present results for variousgraphene systems: single-layer graphene, TBG in the AAand AB configuration and finally for various TBG withgeneric twist angles. Finally, we present conclusions inSec. IV. II. CALCULATION METHOD
In this section, we present the empirical method fordefining the tight-binding Hamiltonian for TBG. Subse-quently, we present the method for evaluating the time-evolution of a state, and the calculation of the prob-ability density and the density of probability currentbased on the kernel polynomial method. Furthermore,we present also a method for evaluating the time auto-correlation function involving the time-evolution of astate, this quantity provides insight on the electronicstructure of a system under study.
A. The empirical tight-binding Hamiltonian
The Hamiltonian defining the dynamics of the 2 p z elec-trons reads: H TBG = (cid:88) ν =1 (cid:88) i,j t νij ˆ c † νi ˆ c νj + (cid:88) i V νi ˆ c † νi ˆ c νi + (cid:88) ν =1 (cid:88) ij t ν ¯ νij ˆ c † νi ˆ c ¯ νj . (1)In this Hamiltonian, the terms in the square bracket de-fine the hopping of the 2 p z electron in two graphenemonolayers (the index ν denotes the layer) with t νij theintra-layer hopping energies between two lattice nodes i and j , and V νi the onsite energies that is generally in-troduced to include local spatial effects. The creationand annihilation of an electron at a layer “ ν ” and a lat-tice node “ i ” is encoded by the operators ˆ c † νi and ˆ c νi ,respectively. The hopping of electron between two layersis described by the last term of the Hamiltonian charac-terised by the hopping parameters t ν ¯ νij . The notation ¯ ν implies that ¯ ν (cid:54) = ν . The values of the hopping parameters t νij and t ν ¯ νij are obtained via the model: t ij = V ppπ exp (cid:18) − R ij − a cc r (cid:19) . (cid:34) − (cid:18) R ij . e z R ij (cid:19) (cid:35) + V ppσ exp (cid:18) − R ij − dr (cid:19) . (cid:18) R ij . e z R ij (cid:19) . (2)This model for the hopping parameters is constructedthrough two Slater-Koster parameters V ppπ ≈ − . V ppσ ≈ .
48 eV. These parameters characterise thehybridisation of the nearest-neighbour 2 p z orbitals inthe intra-layer and inter-layer graphene sheets, respec-tively. The hopping parameters decay with exponentiallaw as a function of the distance between the lattice nodes R ij = | R ij | ; R ij is the vector connecting two lattice sites i and j ; e z is the unit vector along the z -direction perpen-dicular to the two graphene layers, and d ≈ .
335 nm isthe distance between two graphene layers. Accordingly,when i and j belong to the same layer, R ij is perpen-dicular to e z so that we obtain the intra-layer hopping t νij = V ppπ exp[ − ( R ij − a cc ) /r ], otherwise we get t ν ¯ νij .The other parameters are defined as: r ≈ . √ a cc an empirical parameter characterizing the decay of theelectron hopping, and a cc ≈ .
42 ˚A the distance betweentwo nearest carbon atoms in the graphene lattice. Inthis work, we are interested in the intrinsic properties ofTBG, so we simply set the onsite energies V σi to be zero. B. The formalisms for the time-evolution of a state
Let us start by considering an initial state | ψ (0) (cid:105) atthe time t = 0. This state can evolve in time to | ψ ( t ) (cid:105) by acting on it with the time-evolution operator ˆ U ( t ) =exp( − i ˆ Ht/ (cid:126) ): | ψ ( t ) (cid:105) = ˆ U ( t ) | ψ (0) (cid:105) = exp( − i ˆ Ht/ (cid:126) ) | ψ (0) (cid:105) . (3)This equation is the formal solution of the time-dependent Schr¨odinger equation, where ˆ H denotes theHamiltonian operator. We account for the discrete na-ture of the atomic lattice by describing the system withina tight-binding approximation presented in the Sec. II A.In writing the tight-binding Hamiltonian (1), we use alocalised basis set {| j (cid:105) , j = 1 , ..., N } to specify the rep-resentation. Here, the ket | j (cid:105) denotes the 2 p z orbitallocated at the lattice node j and N is the total numberof lattice nodes of the whole system. We can express a state | ψ ( t ) (cid:105) in this basis set in the following way: | ψ ( t ) (cid:105) = N (cid:88) j =1 g j ( t ) | j (cid:105) , (4)where g j ( t ) determines the probability amplitude of find-ing electron at the node j at time t . The probabilitydensity P j ( t ) = |(cid:104) j | ψ ( t ) (cid:105)| = | g j ( t ) | is the quantitydetermining the dynamics of the electron states. Thevalue of g j ( t ) is obtained by solving the time-dependentSchr¨odinger equation or equivalently by performing thecalculation for Eq. (3).In this work, we evaluate the time-evolution operatorˆ U ( t ) by expanding it in terms of the Chebyshev poly-nomials of the first kind Q m ( x ) = cos[ m arcos( x )]. Asfirst, we rescale the spectrum of the Hamiltonian ˆ H tothe interval [ − , H = W ˆ h + E , wherein W is the half of spectrumbandwidth, E the central point of the spectrum, and ˆ h the rescaled Hamiltonian. Practically, we use the powermethod to estimate W . The time-evolution operator istherefore expanded regarding the Chebyshev polynomialsas follows:ˆ U ( t ) = e iE t/ (cid:126) + ∞ (cid:88) m =0 δ m, + 1 ( − i ) m B m (cid:18) W t (cid:126) (cid:19) Q m (ˆ h ) , (5)where B m is the m -order Bessel function of the first kind,and δ m, is the Kronecker symbol. We define the so-called Chebyshev vectors | φ m (cid:105) = Q m (ˆ h ) | ψ (0) (cid:105) which canbe calculated using the recursive relation | φ m (cid:105) = 2ˆ h | φ m − (cid:105) − | φ m − (cid:105) , (6)with | φ (cid:105) = | ψ (0) (cid:105) and | φ (cid:105) = ˆ h | φ (cid:105) . Thus, the state | ψ ( t ) (cid:105) is formally obtained via: | ψ ( t ) (cid:105) = e iE t/ (cid:126) + ∞ (cid:88) m =0 δ m, + 1 ( − i ) m B m (cid:18) W t (cid:126) (cid:19) | φ m (cid:105) . (7)This equation is exact, but we cannot numerically per-form the summation of an infinite series. We thereforeapproximate | ψ ( t ) (cid:105) by a finite series of M terms. Un-fortunately, this truncation breaks the preservation ofthe norm of | ψ ( t ) (cid:105) . Practically, the number of terms M contributing to the summation in Eq. (7) is chosen toguarantee the norm conservation of | ψ ( t ) (cid:105) in a finite, butsufficiently long, evolution time. For instance, in orderto evolve a state in a square TBG sample with 100 nmsize for an evolution time of 50 fs, M should be about1200. To define the initial condition for the time-dependentSchr¨odinger equation, one usually assumes the wave func-tion at t = 0 of a Gaussian form: ψ k ( r , t = 0) = 1 σ √ π exp (cid:20) − ( r − r ) σ (cid:21) φ k ( r ) . In this Gaussian form, φ k ( r ) can be simply chosen as aplane wave exp( i k r ) or generally as a Bloch functiondefining a propagating electron state. The Gaussianpre-factor modulates the extension of the function φ k ( r )localised around the position r with a width of σ . Theadvantage of this choice is that it allows simulating boththe spreading and the moving of the wave centroid. How-ever, the particular behaviour of these phenomena variesconcerning σ and k , two parameters defining a certaininitial state.In this work, we follow a different strategy: we chose alattice node randomly, then we select the corresponding2 p z orbital to be the initial state. It means that we choosethe coefficients g j ( t = 0) = δ ij e iφ , where φ is a randomreal number, and thus | ψ ( t = 0) (cid:105) = N (cid:88) j =1 δ ij e iφ | j (cid:105) = e iφ | i (cid:105) . (8)This choice, though does not allow to simulate the dis-placement of the wave centroid, allows studying thewhole energy spectrum of the 2 p z electron through thespreading of waves in the various graphene systems.To quantify the electron transport, we calculate theexpectation value of the probability current operator. Inthe tight-binding description, the probability current op-erator reads: ˆ J = i (cid:126) N (cid:88) j,k =1 ( r j − r k ) t jk ˆ c † j ˆ c k . (9)Its expectation value on the state | ψ ( t ) (cid:105) is expressed as (cid:104) ˆ J (cid:105) ( t ) = (cid:80) Nj =1 J j ( t ) where J j ( t ) is interpreted as thedensity of the probability current: J j ( t ) = − (cid:126) (cid:88) i ( r j − r i )Im [ t ij g ∗ i ( t ) g j ( t )] . (10)The study of the time-evolution of a state gives usinformation on the electronic structure of the system.Given an initial state | ψ (0) (cid:105) , the time auto-correlationfunction C ψ ( t ) is defined as the projection of | ψ ( t ) (cid:105) onits initial state | ψ (0) (cid:105) : C ψ ( t ) = (cid:104) ψ (0) | ψ ( t ) (cid:105) . (11)In the tight-binding representation with the initial stateschosen as a localised at a particular lattice node | ψ (0) (cid:105) = | i (cid:105) , the time auto-correlation C i ( t ) = (cid:104) i | ψ ( t ) (cid:105) = g i ( t ), i.e. , equal to the local probability amplitude at thenode i . Its power spectrum, defined as the Fourier trans-form of C i ( t ), is the local density of states of electron inthe considered system: ρ i ( E ) = sπ (cid:126) Ω a Re (cid:20)(cid:90) + ∞ dte iEt/ (cid:126) C i ( t ) (cid:21) (12)where Ω a is the volume assigned for each atom in the lat-tice and s = 2 counts the spin degeneracy. We can obtain the system total density of states from Eq. (12) by replac-ing C i ( t ) by an ensemble average of C i ( t ) over a small setof initial states | i (cid:105) . We implemented this procedure forthe first time in Ref. [36], and results for extremely tinytwist configuration of TBG were in agreement with theapproach of continuum models. III. RESULTS AND DISCUSSION
In this section we present results for the three men-tioned physical quantities introduced in the previous sec-tion: the probability density P j ( t ), the density of proba-bility current J j ( t ), and the time auto-correlation func-tion C i ( t ), to characterise the dynamics of electrons inmonolayer and in bilayer graphene systems. A. Monolayer graphene
To better understand the physics of TBG, we shallstart analysing the more straightforward case of mono-layer graphene. We performed the calculation for thetight-binding Hamiltonians accounting for the NN, NNN,and NNNN hopping terms. As we shall see later, thethree models result in different quantitative behaviourfor the time auto-correlation function but have the samespreading pattern of the electron wave function; thus forsimplicity, we will present results only for the NN case.We present in Fig. 1 the distribution of the probabil-ity densities P j ( t ) = | g j ( t ) | and the probability currentdensities J j ( t ) obtained for the spreading of a 2 p z stateinitially located at a single lattice node. At each latticenode, we use the solid-circles and the arrows to representthe probability densities and the probability current den-sities, respectively. The circle radius is proportional tothe value of P j ( t ), which is normalised at each t to themaximal value of the set { P j ( t ) } , ∀ j ∈ { , ..., N } . Sim-ilarly, the arrow length is proportional to the length of J j ( t ), which is also scaled appropriately. The length andthe direction of the arrows indicate the tendency thatthe probability density at a lattice node to transfer tothe neighbour lattice nodes.The time frames taken at t = 0 . . π and ± π/ i.e. , along the direc-tion of the armchair lines (highlighted in red in Fig. 1,frame with t = 0 . t = 0 . . ± π/ i.e. ,also along the direction of the armchair lines. However,at t = 0 . ± π/ , ± π/ ± π/ i.e. , along the zigzag lines of the honeycomb lattice (c.f.Fig. 1, frame with t = 0 . t ∈ [1 ,
2] fs, we find a continuous spreading of the
FIG. 1. Spreading of the distribution of electron probability densities (in red) in the monolayer graphene taken at several timemoments. The arrows denote the vectors of probability current density (in green). The blue lines denote the three mirror-symmetry planes σ v of the hexagonal lattice. We highlight the direction of the armchair- and zigzag-lines in red in the framewith t = 0 . electron wave function, and we observe the formation ofa wavefront with the hexagonal shape. After 2 fs, thewavefront is well established with the corners headingthe directions of the zigzag lines.To quantify the pattern of the wave spreading we di- rectly inspected the distribution of both the probabilitydensities P j ( t ) and the probability current densities J j ( t )on the lattice nodes. We learnt that the distribution ofthese two quantities obeys the features of the point group D h . These symmetry properties are not identical tothose of the graphene lattice, described by the symmor-phic space group p mm , and thus the point group D h . However, we should remember that D h is a sub-group ofthe one D h , and is the point group of the lattice node.We then conclude that the spreading pattern of electronsdepends on not only the lattice symmetry, but also thatof the initial state.To get quantitative information on the energy spec-trum of the π -bands from the observation of the spread-ing of a 2 p z state, we calculated the time auto-correlationfunction C i ( t ) [c.f. Eq. (11)]. The Fourier transform of C i ( t ) provides the local density of states at the latticenode i . In our calculation, we found that the time auto-correlation function, though being a complex function ingeneral, gets purely real when we consider a model withonly the NN hopping. In Fig. 3 we show the value of C i ( t ) obtained by invoking the three hopping models.The result shows the oscillating behaviour of C i ( t ) as afunction of time with decreasing magnitude. This be-haviour implies the declining of the correlation at longevolution times. For the NN hopping approximation, theHamiltonian has only one parameter t cc = V ppπ whichsets the energy scale. In this case, we find that C i ( t ) isperiodic with the oscillation pattern remarked in Fig. 3.By changing the value of V ppπ and measuring the cor-responding frequency f , we verified that the frequencyis determined by f = V ppπ / (2 π (cid:126) ) = 6 . × Hz. Byintroducing in the Hamiltonian higher order hopping pro-cesses, the behaviour of C i ( t ) becomes complex, and wecannot find a clear dominant frequency associated to anyof the higher-order hopping terms. Fourier transforming C i ( t ) via Eq. (12) results in the local density of stateswhich shows the electron-hole symmetry in the NN andNNNN model, but not in the NNN model. B. AA- and AB-stacking bilayer graphene
We will analyse in this section two particular cases oftwisted bilayer graphene: the AA- and the AB-stacking.One should notice that we generate the TBG configu-rations by starting from the AA-stacking configurationand then twisting the two graphene layers about a verti-cal axis going through a pair of carbon atoms. Accord-ingly, the AA- and AB-stacking configurations are char-acterised by a twist angle of 0 ◦ and 60 ◦ , respectively.The point group symmetry of the AA- and AB-bilayergraphene is related to the one of the monolayer. Pre-cisely, the symmetry of the AA-stacking bilayer grapheneis characterised by the symmorphic space group p mm ,generated by the lattice translation and the point group D h , whereas the AB-stacking system is characterises bythe symmorphic space group symmetry p m gener-ated by the lattice translation and the point group D d . We start by considering the inter-layer transfer of elec-tron wave function: we calculate the layer-integratedprobability densities. This quantity is expressed as the R e [ C i (t)] -0.4-0.200.20.40.60.81 Time (fs) I m [ C i (t)] -0.2-0.100.10.20.3 NN-modelNNN-modelNNNN-model
Time (fs)
10 15 20 25 30 R e [ C i (t)] -0.02-0.0100.010.02 = FIG. 2. The real and imaginary part of the auto-correlationfunction C i ( t ) = (cid:104) i | ψ ( t ) (cid:105) calculated from three hopping mod-els: the NN in red, the NNN in green and the NNNN in blue.The red rectangle box remarks the typical oscillation patternpredicted by the NN model. The inset shows the continuousvariation of the real part of C i ( t ) in a longer evolution time. summation of the probability density in each layer: P α ( t ) = (cid:88) j ∈ ( L α ) P j ( t ) ∀ α ∈ { T,B } (13)In Fig. 3, we present the variation of P T and P B as a func-tion of the time-evolution. The layer-integrated probabil-ity density between the two graphene layers presents anoscillatory pattern as a function of time; this behaviouris similar to the phenomenon discussed by Xian et al. as the neutrino-like oscillation. In the case of the AA-stacking configuration, we observe how the wave on thetop layer quickly penetrates into the bottom one com-pared to the AB-stacking configuration. After almost1 . d = 3 .
35 ˚A isthe same. The hybridisation of the 2 p z orbitals betweenthe two graphene layers is also characterised by the sameenergy value of V ppσ = 0 .
48 eV. In order to analyse therole of the interlayer coupling on the electron dynamicsin the two-layer systems, we investigate the effects of tun-ing the inter-layer coupling parameter V ppσ on the layer-integrated probability densities. In Fig. 3 we presentthese probabilities obtained by setting V ppσ = 0 .
48 eV(the solid curves) and 0.12 eV (the dot-dashed curves).We observe how the reduction of V ppσ yields to an in-crease of the characteristic transfer time, which we defineas the evolution time τ T → B needed to transfer 50% of thewave from the top layer to the bottom one. Calculationfor various values of V ppσ shows that τ T → B ∝ /V ppσ . Time (fs) P T = B Time (fs) P T = B Top-layer, V < = 0.48 eVBottom-layer, V < = 0.48 eVTop-layer, V < = 0.12 eVBottom-layer, V < = 0.12 eVa)b) = T ! B = T ! B FIG. 3. Neutrino-like oscillation of the layer-integrated prob-ability density P T/B ( t ) for the (a) AA- and (b) AB-stackingbilayer graphene. Data plotted for two values of the Slater-Koster parameter V ppσ = 0 .
48 eV (the solid curves) and 0.12eV (the dot-dashed curves). Only data in the interval of(0 ,
30) fs are displayed to zoom-in the oscillations. The verti-cal lines highlight the times τ T → B and t = 1 . At very long evolution time, each graphene layer sup-ports about one-half of the initial waves and the layer-interchange transfer becomes almost stationary with veryweak oscillations in time. It is worthy to note that, for theAB-stacking configuration, we distinguished two cases ofthe initial state | i (cid:105) , one at the A-atom on top of the B-atom in the bottom layer, and the other at the B-atom onthe center of the hexagonal ring underneath. We foundthat the layer-integrated probability densities in the twocases are the same, but the in-plane wave spreading pat-terns are different as discussed in next paragraphs. Wenow analyse the features of the intra-layer spreading pat-terns in the AA- and AB-configurations. In Figs. 4 and 5we present the evolution of a 2 p z state initially localisedat a lattice node in the top layer of the two AA- and AB-stacking configurations of bilayer graphene, respectively.We use colours to represent the probability densities ontwo graphene layers, specifically, the red for the top layerand the black for the bottom one. Similar to the case ofmonolayer graphene, the radius of the solid-circle at eachlattice node is proportional to the value of P j normalisedat the maximum value for each value of time.By comparing the wavefront spreading behaviour ofthe electron wave function in the monolayer grapheneand that in the AA-stacking configuration for the evolu-tion time t < . P j ( t ) on the top and bottom graphene layers are simi-lar to the case of monolayer graphene, but becomes dif-ferent for larger evolution time. It should be noticedfrom Fig. 3(a) that in the duration of (0 , .
3) fs the top layer-integrated probability density P T monotoni-cally decreases. It means that the wave continues trans-ferring to the bottom layer and achieves the maximaltransferring percentage at 1 . t , the part of the wave in the bottom layer transferback to the top one. It results in the oscillation behaviourof P T ( t ) and P B ( t ) similar to a Fabry-P´erot resonator.From Fig. 3(a) we can determine a set of characteristictime intervals, e.g. , (0 , .
3) fs, (2 . , .
6) fs, (6 . , .
3) fs,and so on, in which the wave transfers predominantlyfrom the top to the bottom layer; alternatively, in thecomplementary time intervals the wave transfers in theopposite direction. We found that after 1 . P j ( t ) of monolayer graphene is neither co-incident to that on the top nor the bottom layers ofthe bilayer system. This difference is the result of thecombination of the intra-layer and inter-layer spreadinginduced by the hopping terms in the Hamiltonian (1).The wavefront at long evolution time present a hexag-onal shape similar to the monolayer graphene case. Adirect inspection, however, shows that the form of thespreading pattern obeys only the point group C v . Re-member that the symmetry of a node in the AA latticeis characterised by the point group D h but, the succes-sive interlayer penetration of the electron wave lowers thesymmetry of the distribution of probability densities tothe C v symmetry.For the AB-stacking configuration, we use the sametechnique for displaying data as for the AA configura-tion. From Fig. 3(b), we learn that for this configura-tion, the characteristic time τ T → B ≈ t < . P j ( t ) on thetop layer is identical to that of the monolayer graphene,but a quantitative difference becomes appearing when t ∈ (1 . , .
1) fs. When t > . , .
3) fs is the onein which the wave monotonically transfers from the toplayer to the bottom one. Though the percentage of thewave transfer at t ≈ . , .
3) fs, (2 . , .
1) fs, (3 . , .
6) fs,and so on, in which the wave transfers predominantlyfrom the top to the bottom layer, and in the complemen-tary intervals where the wave transfers in the oppositedirection. At long evolution time, the wave spreading isalso characterised by a wavefront in the hexagonal shapethat, similar to the AA lattice case, reflects the planesymmetries of the lattice nodes in the AB-stacking sys-tem, i.e. , the group C v , a sub-group of the point group D d . FIG. 4. Spreading of the electron probability density for the AA-stacking configuration taken at several time moments. Latticenodes in red/black belong to the top/bottom graphene layer. The blue lines in the frames at t = 1 . t = 3 . σ v of the lattice.FIG. 5. Spreading of the electron probability density for the AB-stacking configuration taken at several time moments. Theinitial state | i (cid:105) is set at the position of an A-atom of the top graphene layer which is on top of a B-atom of the bottom layer.Lattice nodes in red/black belong to the top/bottom graphene layer. The blue lines in the frames at t = 1 fs and t = 3 . σ v of the lattice. Time (fs) R e [ C i (t)] -0.01-0.00500.0050.01 Time (fs) R e [ C i (t)] -0.01-0.00500.0050.01 AB-stackingAA-stackingV < = 0.48 eVV < = 0.12 eV FIG. 6. Time auto-correlation functions C i ( t ) of the AA-(red) and AB-stacking (blue) bilayer graphene calculated us-ing the NN model with V ppσ = 0 .
48 eV and 0.12 eV.
We also calculated the time auto-correlation function C i ( t ) for the AA- and AB-stacking configurations. InFig. 6 we present the data for C i ( t ) as a function of thetime-evolution for the two different parameter models: V ppσ = 0 .
48 eV and 0 .
12 eV. We observe, in general,the intricate behaviour of C i ( t ) in the two cases, whichare different from each other. Interestingly, the figureshows the beating behaviour of the auto-correlation func-tions when considering V ppσ = 0 .
12 eV. Our calculationshows that the beating behaviour does not appear clearlywith V ppσ = 0 .
48 eV, but it does when decreasing thevalue of V ppσ to the values smaller than about 0.3 eV.We also realise the beating oscillation behaviour is simi-lar to the oscillation features of the time auto-correlationfunction of the monolayer graphene. It is expected be-cause we should obtain a picture of the two independentgraphene layers in the limit of V ppσ = 0. This observa-tion reflects the fact that the interlayer coupling plays therole of modulating the electronic states between the twographene layers. When a wave is spreading in one layer,it penetrates partially into the other and thus creates twowaves spreading in the two layers. The coupling betweenthe two layers induces the exchange of wave between thetwo layers and forms the wave-interference pattern in thespace limited by the two layers. The interference is sen-sitive to the alignment of the two atomic lattices. Itthus explains the typical evolution features of electronicstates in particular atomic configurations. Though thebehaviour of the auto-correlation functions versus thetime is complicated, its Fourier transform results in thedensity of states of these two configurations. Time (fs) P T = B Time (fs) P T = B = 5.0 ° , top-layer = 5.0 ° , bottom-layer = 2.5 ° , top-layer = 2.5 ° , bottom-layer = 1.0 ° , top-layer = 1.0 ° , bottom-layera)1.5 fs b) FIG. 7. Layer-integrated probability densities P T/B ( t ) in thetop (solid lines) and bottom (dot-dashed lines) layers of threeTBG configurations with the twist angles of 5 ◦ (green), 2 . ◦ (blue), and 1 ◦ (red). Panels (a) and (b) are for the casesthat the initial state | p z (cid:105) locates at the central point of theAA -like and AB-like region, respectively. The parameter V ppσ = 0 .
48 eV. The vertical lines highlight the time t = 1 . C. Twisted bilayer graphene
Regarding the atomic structure, twisted bilayergraphene is a generalisation of the AA- and AB-stackingbilayer systems with a generic rotation angle betweenthe two graphene layers. In general, the alignment be-tween the two graphene lattice in the TBG systems is notcommensurate, i.e. , not defined by a unit cell, and thusthe lattice has very low symmetry. In the case of com-mensurate stacking, the space group characterising theTBG lattice is determined to be either p m p mm depending on both the twist origin and the twist an-gles. . Interestingly, the generic TBG lattice shows aspecial moir´e structure of the hexagonal form. In eachmoir´e zone, we can find regions in which the atomic ar-rangement is close to the AA- and AB- or BA-stackingconfigurations. We illustrate the moir´e zone in Fig. 8with the blue hexagon where we marked the AA- andthe AB-like regions [c.f. the frame with t = 0 . . ◦ and 5 ◦ . The inter-layer0 FIG. 8. Spreading of the electron probability density in the TBG configuration with θ = 2 . ◦ taken at several time moments.The initial state | ψ (0) (cid:105) = | i (cid:105) located in the centre of the AA-like region in the moir´e zone. Lattice nodes in red/black belong tothe top/bottom layer. The blue hexagon denotes the moir´e zone. The characters AA i ( i = 1 , ..., i and BA i ( i = 1 , , θ = 5 ◦ (left panel) and θ = 2 . ◦ (right panel) atlarge evolution times t = 27 fs and t = 30 fs. FIG. 10. Similar to Fig. 9 but with the initial state | ψ (0) (cid:105) = | i (cid:105) located in the centre of the AB-like region in the moir´e zone.Lattice nodes in red/black belong to the top/bottom layer. coupling always induces the wave transfer between thetwo graphene layers. Figure 7 shows that similar to thecase of the AA-configuration, the transmission of elec-tron wave function from the top layer to the bottom onereaches a maximal value in about 1.5 fs. To illustratethe wave transfer between two graphene layers, we studythe variation of the layer-integrated probability densi-ties on time. Figure 7 shows the oscillation behaviourof the layer-integrated probability densities for three in-commensurate TBG configurations with the twist angleof 5 ◦ (green), 2 . ◦ (blue) and 1 ◦ (red). (The last oneis close to the first magic angle θ ≈ . ◦ . ) Further-more, we investigated how the layer-integrated probabil-ity densities change by changing the initial position: thepanels (a) and (b) are for the cases that the initial statelocalised in the centre of the AA -like and AB-like re-gions, respectively. We observe that the percentage ofthe wave transmitted from the top layer to the bottom one depends on the twist angle. For a short time of evo-lution, t < ◦ . However, after5 fs, there is about 60% of the wave propagating on thetop layer and about 40% doing in the bottom one. Theminimal oscillation of the green curves P T/B ( t ) implies avery weak transfer of wave between the two layers. Thisdynamical observation supplements to the explanationof the effective decoupling of the two graphene layers inthe TBG configurations with large twist angles. Inother words, the two parts of the wave become nearlyindependently propagating on the two graphene layersafter a long time-evolution. For the TBG configurationswith much smaller twist angles, e.g. . ◦ and 1 ◦ , afterabout 10 fs, the fluctuation of the blue and red curves P T/B ( t ) is always significant around the value of 50%.It implies the strong interaction between the two wavecomponents when propagating in the TBG lattices with2 Time (fs) R e [ C (t)] -0.03-0.02-0.0100.010.020.03 Time (fs) I m [ C (t)] -0.03-0.02-0.0100.010.020.03 = 5.0 ° = 2.5 ° = 1.0 ° a)b) FIG. 11. The real (a) and imaginary (b) parts of the timeauto-correlation function C ( t ) for three TBG configurationswith the twist angles of θ = 5 ◦ (green), 2 . ◦ (blue) and 1 ◦ (red). small twist angles.We now present in Fig. 8 the intra-layer spreading ofa 2 p z state initially located at the central site of the AA-like region through the distribution of the probabilitydensities P j ( t ). The time-evolution is illustrated simi-larly to the case of the AA- and AB-stacking configura-tions. From the figure we observe that during the timeinterval (0 ,
2) fs, the initial state spreads similarly to caseof the AA-stacking lattice, i.e. , extending along the sixpreferable directions, then a hexagonal wavefront is es-tablished. In this time interval, the wavefront takes thetypical hexagonal shape, and it is still within the AA-like region. For time larger than t = 3 fs, we observe howthe wavefront corners enter the AB-like regions, and thewavefront edges reach the transition regions between theAA j and AA regions. At this time, the probability den-sities start to be redistributed: they become concentratedin the AB-like regions as the six clusters seen in the timeframe at t = 4 fs. These clusters then move in the tran-sition regions between the AA i and AA i +1 regions, i.e. ,along the zigzag lines connecting the AB-like regions inthe first moir´e zone to the other AB-like regions in thenext moir´e zones, while the probability densities on theedges of the hexagonal wavefront become scattered intothe AA i +1 ( i (cid:54) = 0) regions, see the time frames at t = 5 , C symmetry. The wavefronts on the two graphene layersare not coincident due to the misalignment of the latticestacking.Interestingly, for long evolution time, we observe the higher intensity of the probability densities in the AA-like regions ( t >
15 fs), particularly, in the TBGs oftiny twist angles, this higher intensity is observed sig-nificantly in the AA i region only, see Fig. 9. This ob-servation might reflect the “localisation” of low-energyBloch wave functions in the AA-like regions as depictedin Refs. [6, 7, 18, 19, and 36]. Notice that, at the evolu-tion time t , we would expect the dominance of electronstates of energies about (cid:126) /t . It therefore implies that, atlong observation time, the localised signature shown inFig. 9 of the electron wave function in the AA-like re-gions is the behaviour of the states associated with thenarrow energy band around the charge neutrality level.This localisation feature might also be related to topo-logical properties of the wavefront as recently pointedout in Ref. [8, 45, and 46]. Quantitatively, this associ-ation is consolidated by Fourier transforming the timeauto-correlation C ( t ) defined by Eq. (11) to obtain thedensisty of states. The resulted DOS of the TBG con-figurations with the twist angles θ < . ◦ shows a small,but significant peak, around the charge neutrality levelas reported in Refs. [7 and 36].In the following, once again, we will show how thedynamics of wave spreading in TBG strongly dependson the symmetry of an initial localised state. We nowconsider an initial 2 p z state at the central node of oneAB-like region in the moir´e zone. We note that, con-trary to the central node in the AA-like regions, for thischoice, there is no exact symmetry elements containingthe central node of the AB-like regions. For short time-evolution ( t < i.e. , alongthe zigzag lines separating the AA-like regions (c.f. panel t = 2 fs in Fig. 10). Along the opposite directions, thewave spreads into the AA-like regions, and the proba-bility densities become concentrated at the centre ratherthan scattered, (c.f. panels for t = 3 and 4 fs in Fig. 10).Following the distribution of the probability densities atlarger times, the probability densities propagate alongthe zigzag lines in the transition regions between theAA i - and AA i +1 -like regions and concentrated in the cen-tre of the AA-like regions. Due to the “approximated”symmetries about the initial position of the 2 p z state, thewavefront is formed and has an almost hexagonal shape.The six corners of the wavefront orient the preferablyevolved directions. The distribution of P j ( t ) on the twolayers satisfies the “approximated” point group symme-try C .To complete our discussion of the wave evolution inthe TBG lattices, we present in Fig. 11 the time auto-correlation function C ( t ). We remind again that thechoice of the initial condition affects crucially the wavespreading. In Figs. 8 and 10, we have presented thedata for two particular initial conditions which result intypical spreading patterns of the 2 p z state. In order toextract quantitative information on observables, for in-3stance the density of states, from the time-evolution ofelectronic states we have to account for all possible ini-tial conditions. According to Eq. (11), to calculate thetime auto-correlation function C ( t ), we need to calcu-late a set of functions C i ( t ) = (cid:104) i | ψ ( t ) (cid:105) with the initialstates | i (cid:105) = | p z (cid:105) chosen at every lattice node in a suffi-ciently large TBG sample. Though a TBG lattice is notalways defined by a unit cell with translational symme-try, the moir´e zone can be seen as an approximated unitcell. It suggests that we need to consider only the latticenodes in a moir´e zone. However, since the typical length L M defining the size of the moir´e zone is related to thetwist angle θ via the expression L M = √ a cc / θ/ C ( t ). The results are shown in Fig. 11for three TBG configurations with θ = 5 ◦ , . ◦ and 1 ◦ .The figure shows the complex behaviour of C ( t ) as a func-tion of time. Despite that, the Fourier transform of C ( t ),see Eq. (12), results in the density of states with typicalvan Hove peaks are shown in Ref. [7 and 36]. IV. CONCLUSION
We have presented a study of the time-evolution char-acteristics of electrons in the bilayer graphene latticeswith arbitrary twist angles. We used the Chebyshevpolynomials of the first kind to approximate the time-evolution operator for a sufficiently long time-evolutionto calculate time-correlation functions reliably. We haveshown that the inter-layer electronic coupling induces theinterchange transfer of waves between the two graphenelayers, resulting in the oscillating behaviour of the layer-integrated probability densities as a function of time,similar to complex Fabry-P´erot oscillations. This be-haviour can be also interpreted as the precession of elec-trons when describing the moir´e-induced spatial modu-lation in the interlayer coupling in terms of non-Abeliangauge fields. The percentage of the wave transmittedfrom one layer to the other depends on the twist angle, i.e. , smaller than 50% and weak oscillation for large twist angles, θ > . ◦ , and larger than 50% and strong oscilla-tion otherwise. This dynamical observation supplementsthe understanding of the effective decoupling between thetwo graphene layers in the TBG configurations with largetwist angles. For the wave spreading in each graphenelayer, we have indicated that the spreading shape of elec-tron waves is dictated by the dominant hopping mecha-nism of the honeycomb pattern of the monolayer latticeand by the plane symmetries of the bilayer lattices. Thewave spreading is irregular and takes place in two stages:The first one occurs within a very short time-evolution, inwhich the wave spreads to the three nearest neighboursand then develop to the lattice nodes along the direc-tions of the armchair-lines of the honeycomb lattice. Thesecond stage is characterised by the formation of a well-defined wavefront of hexagonal shape with the cornersdeveloping faster the edges. For tiny twist TBG config-urations, we have observed the signature of the electronlocalisation in the AA-like regions inside the TBG’s moir´ezone at long time-evolution. This would associate withthe formation of a narrow energy band around the chargeneutrality level. We have shown the interchange transferof wave between the two graphene layers resulting in thedifference of the distribution of the probability densitieson the TBG lattices from that on the monolayer. Wehave also observed the appearance of a beating patternin the autocorrelation functions for a reduced intra-layercoupling — it is possible to achieve this reduction exper-imentally. It might suggest a way for engineering theelectronic properties of the bilayer systems. This studyprovides a complementary intuitive understand of theelectron behaviours in the twisted bilayer graphene. Thecalculation method implemented here represents an al-ternative paradigm for future studies of exotic electronicproperties of layered materials, including twisted bilayergraphene but also other van der Waals heterostructures.
ACKNOWLEDGEMENTS
The work of VND and HAL is supported by the Na-tional Foundation for Science and Technology Develop-ment (NAFOSTED) under Project No. 103.01-2016.62.The work of DB is supported by Spanish Ministerio deCiencia, Innovation y Universidades (MICINN) under theproject FIS2017-82804-P, and by the Transnational Com-mon Laboratory
Quantum-ChemPhys . ∗ [email protected] M. Xu, T. Liang, M. Shi, and H. Chen, “Graphene-liketwo-dimensional materials,” Chem. Rev. , 3766 (2013). A. K. Geim and I. V. Grigorieva, “Van der waals het-erostructures,” Nature , 419 (2013). A. V. Rozhkov, A. O. Sboychakov, A. L. Rakhamanov,and F. Nori, “Electronic properties of graphene-based bi-layer systems,” Phys. Rep. , 1 (2016). J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Cas-tro Neto, “Graphene bilayer with a twist: Electronic struc-ture,” Phys. Rev. Lett. , 256802 (2009). J. M. Lopes dos Santos, N. M. R. Peres, and A. H. CastroNeto, “Continuum model of the twisted graphene bilayer,”Phys. Rev. B , 155449 (2012). R. Bistritzer and A. H. MacDonald, “Moir bands in twisteddouble-layer graphene,” Proc. Natl. Acad. Sci. U.S.A. , D. Weckbecker, S. Shallcross, M. Fleischmann, N. Ray,S. Sharma, and O. Pankratov, “Low-energy theory for thegraphene twist bilayer,” Phys. Rev. B , 035452 (2016). M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi,K. Kuroki, and L. Fu, “Maximally localized wannier or-bitals and the extended hubbard,” Phys. Rev. X , 031087(2018). Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi,E. Kaxiras, and P. Jarillo-Herrero, “Unconventional super-conductivity in magic-angle graphene superlattices,” Na-ture , 43 (2018). Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken,J. Y. Lou, J. D. Sanchez-Yamagishi, K. Watanable,T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, “Correlated insulator behaviour at half-fillingin magic-angle graphene superlattices,” Nature , 80(2018). Liujun Zou, Hoi Chun Po, Ashvin Vishwanath, andT. Senthil, “Band structure of twisted bilayer graphene:Emergent symmetries, commensurate approximants, andwannier obstructions,” Phys. Rev. B , 085435 (2018). M. Angeli, D. Mandelli, A. Valli, A. Amaricci, M. Capone,E. Tosatti, and M. Fabrizio, “Emergent D symmetry infully relaxed magic-angle twisted bilayer graphene,” Phys.Rev. B , 235137 (2018). S. Shallcross, S. Sharma, and O. A. Pankratov, “Quan-tum interference at the twist boundary in graphene,” Phys.Rev. Lett. , 056803 (2008). E. J. Mele, “Interlayer coupling in rotationally faulted mul-tilayer graphenes,” J. Phys. D: Appl. Phys. , 154004(2012). A. V. Rozhkov, A. O. Sboychakov, Rakhmanov A. L., andF. Nori, “Single-electron gap in the spectrum of twistedbilayer graphene,” Phys. Rev. B , 045119 (2017). J.C. Rode, D. Smirnov, C. Belke, H. Schmidt, and R. J.Haug, “Twisted bilayer graphene: interlayer configurationand magnetotransport signatures,” Ann. Phys. , 1700025(2017). E. Su´arez Morell, J. D. Correa, P. Vargas, M. Pacheco,and Z. Barticevic, “Flat bands in slightly twisted bilayergraphene: Tight-binding calculations,” Phys. Rev. B ,121407 (2010). G. Trambly de Laissardiere, D. Mayou, and L. Magaud,“Local- ization of dirac electrons in rotated graphene bi-layers,” Nano Lett. , 804 (2010). G. Trambly de Laissardiere, D. Mayou, and L. Magaud,“Numer- ical studies of confined states in rotated bilayersof graphene,” Phys. Rev. B 86 , 125413 (2012). Kazuyuki Uchida, Shinnosuke Furuya, Jun-Ichi Iwata,and Atsushi Oshiyama, “Atomic corrugation and elec-tron localization due to moir´e patterns in twisted bilayergraphenes,” Phys. Rev. B , 155451 (2014). P. Lucignano, D. Alfe, V. Cataudella, D. Ninno, andG. Cantele, “The crucial role of atomic corrugation on theflat bands and energy gaps of twisted bilayer graphene atthe “magic angle” θ . ◦ ,” (2019), arXiv:1902.02690. W. Zawadzki and T. M. Rusin, “Zitterbewegung (trem-bling motion) of electrons in semiconductors: a review,” J.Phys.: Condens. Matter , 143201 (2011). G. M. Maksimova, V. Ya. Demikhovskii, and E. V. Frol-ova, “Wave packect dynamics in a monolayer graphene,”Phys. Rev. B , 235321 (2008). G. M´ark, P. Vancs´o, C. Hwang, P. Lambin, and L. P. Bir´o,“Anisotropic dynamics of charge carriers in graphene,”Phys. Rev. B , 125443 (2012). A. Chaves, L. Covaci, Kh. Yu. Rakhimov, G. A. Farias,and F. M. Peeters, “Wave-packed dynamics and valley fil-ter in strained graphene,” Phys. Rev. B , 205430 (2010). L. Xian, Z. F. Wang, and M. Y. Chou, “Coupled diracfermions and neutrino-like oscillation in twisted bilayergraphene,” Nano Lett. , 5159 (2013). A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S.Novoselov, and A. K. Geim, “The electronic propertiesof graphene,” Rev. Mod. Phys. , 109–162 (2009). E. Schr¨odinger, “ ¨Uber die kr¨aftefreie Bewegung in derrelativistischen Quantenmechanik,” Sitzungsber. Preuss.Akad. Wiss. Phys. Math. Kl. , 418 (1930). M. I. Katsnelson, “Zitterbewegung, chirality, and minimalconductivity in graphene,” Eur. Phys. J. B , 157 (2006). M. S. Jang, H. Kim, H. A. Atwater, and W. A. GoddardIII, “Time dependent behavior of a localised electron in aheterojunction boundary of graphene,” Appl. Phys. Lett. , 043504 (2010). V. Ya. Demikhovskii, G. M. Maksimova, A. A. Perov, andE. V. Frolova, “Space-time evolution of dirac wave pack-ets,” Phys. Rev. A , 052115 (2010). Kh. Y. Rakhimov, A. Chaves, G. A. Farias, and F. M.Peeters, “Wavepacket scattering of dirac and schr¨rodingerparticles on topical and magnetic barriers,” J. Phys.: Con-dens. Matter , 275801 (2011). S. T. Park, “Propagation of a relativistic electron wavepacket in the dirac equation,” Phys. Rev. A , 062105(2012). D. E. Fernandes M. Rodrigues and G. Falcao, “Time evo-lution of electron waves in graphene superlattices,” AIPAdvances , 075109 (2016). A. Weibe, G. Wellein, A. Alvermann, and H. Fehske,“The kernel polynomial method,” Rev. Mod. Phys. ,275 (2006). H. A. Le and V. N. Do, “Electronic structure and opticalproperties of twisted bilayer graphene calculated via timeevolution of states in real space,” Phys. Rev. B , 125136(2018). P. Moon and M. Koshino, “Optical absorption in twistedbilayer graphene,” Phys. Rev. B , 205404 (2013). M. Koshino, “Interlayer interaction in general incommen-surate atomic layers,” New. J. Phys. , 015014 (2015). H. N. Nazareno, P. E. de Brito, and E. S. Rodrigues,“Dynamics of wave packets in two-dimensional crystals un-der external magnetic and electric fields: formation of vor-tices,” Phys. Rev. B , 125405 (2007). E. Kogan, V. U. Nazarov, V. M. Silkin, and M. Kaveh,“Energy bands in graphene: comparison between the tight-binding model and ab initio calculations,” Phys. Rev. B ,165430 (2014). Sylvain Latil and Luc Henrard, “Charge carriers in few-layer graphene films,” Phys. Rev. Lett. , 036803 (2006). Mikito Koshino and Edward McCann, “Parity and val-ley degeneracy in multilayer graphene,” Phys. Rev. B ,115315 (2010). A. Luican, G. Li, A. Reina, J. Kong, R. R. Nair, K. S.Novoselov, A. K. Geim, and E. Y. Andrei, “Single-layerbehavior and its breakdown in twisted graphene layers,”Phys. Rev. Lett , 126802 (2011). I. Brihuega, P. Mallet, H. Gonz´alez-Herrero, G. Tramblyde Laissardi`ere, M. M. Ugeda, L. Magaud, J. M. G´omez-Rodrguez, F. Yndur´ain, and J.-Y. Veuillen, “Unravelingthe intrinsic and robust nature of van hove singularities intwisted bilayer graphene by scanning tunneling microscopyand theoretical analysis,” Phys. Rev. Lett. , 196802(2012). Jian Kang and Oskar Vafek, “Symmetry, maximally local-ized wannier states, and a low-energy model for twistedbilayer graphene narrow bands,” Phys. Rev. X , 031088(2018). Hoi Chun Po, Liujun Zou, Ashvin Vishwanath, andT. Senthil, “Origin of mott insulating behavior and su-perconductivity in twisted bilayer graphene,” Phys. Rev.X , 031089 (2018). P. San-Jose, J. Gonz´alez, and F. Guinea, “Non-abeliangauge potentials in graphene bilayers,” Phys. Rev. Lett. , 216802 (2012). J. W. Jeon, H. Kim, H. Kim, S. Choi, and B. H. Kim, “Ex-perimental evidence for interlayer decoupling distance oftwisted bilayer graphene,” AIP Advances8