Formation of topological domain walls and quantum transport properties of zero-line modes in commensurate bilayer graphene systems

Abstract

We study theoretically the construction of topological conducting domain walls with a finite width between AB/BA stacking regions via finite element method in bilayer graphene systems with tunable commensurate twisting angles. We find that the smaller is the twisting angle, the more significant the lattice reconstruction would be, so that sharper domain boundaries declare their existence. We subsequently study the quantum transport properties of topological zero-line modes which can exist because of the said domain boundaries via Green’s function method and Landauer−Büttiker formalism, and find that in scattering regions with tri-intersectional conducting channels, topological zero-line modes both exhibit robust behavior exemplified as the saturated total transmissionG_tot ≈ 2e²/h and obey a specific pseudospin-conserving current partition law among the branch transport channels. The former property is unaffected by Aharonov−Bohm effect due to a weak perpendicular magnetic field, but the latter is not. Results from our genuine bilayer hexagonal system suggest a twisting angle aroundθ ≈ 0.1° for those properties to be expected, consistent with the existing experimental reports.

Graphical abstract

Keywords

twistronics / lattice reconstruction / topological domain wall / zero-line mode / quantum transport

1 Introduction

Twistronics is a newly-emerged scientific research subject, which studies the influence on electron bandstructures originating from relative rotations about an out-of-plane axis (aka twisting) between layers of two-dimensional materials. An immediate and generic consequence of twisting is the inevitable local misalignment of lattices, even for those that are completely matched before twistings. A commensurate twisted system still preserves spatial periodicity (usually larger than the original one possessed by untwisted layers) after the twisting, otherwise it is incommensurate. As the typical representative, twisted graphene system has been reported to display both superconductivity similar to that of high-critical-temperature copper oxides [1-5] and Mott-like strong-correlated effect [2,3,6,7]. The former provides clues for the long-outstanding problem of unconventional superconductivity and the latter demonstrates a practical controlling approach able to be used in systems with non-negligible interactions, which are usually difficult to control. Therefore, twisted two-dimensional materials attract much attention. Additionally, studies on interlayer couplings [8], experimental observation on Moiré excitons [9-12], its relation to twisting angles [13], and exciton phase transition [14] have as well been conducted in transition metal dichalcogenide twisted systems. Furthermore, in twisted bilayer graphene, people have also studied topological nontrivial phases [15], including Chern insulator [16-19] and higher-order topological insulator [20-22]. To properly describe the major characteristic in its electron structure, such as flat bands and subgaps, low-energy effective models [23-30] and tight-binding models [27,31-36] have been explored.

Commensurate twisted bilayer graphene forms in real space Moiré lattice, where AB/BA stacking regions are domains and domain walls as the transitional regions between them go through AA stacking regions. If each of the monolayer hexagonal lattice is assumed rigid, then the aforementioned domains are relatively small and at the same time wide domain walls are present. This assumption validates itself for cases with large twisting angles ( \theta \gtrsim {1.2}^{\circ }) [37], but turns out invalid for those minimally twisted, because for the latter atomic relaxation indispensably comes into play [38-41]. To capture the effect from this lattice reconstruction, people have developed empirical potential method based on density functional theory [37,42,43], continuum model with parameters [44], and molecular dynamics method [45]. Additionally, the canonical functional optimization method with a very clear physical picture by the celebrated Euler−Lagrange equation [31,46-50] is of special interest, which transforms the energy functional optimization problem into a boundary value problem of a group of simultaneous partial differential equations.

Basically, lattice reconstruction in commensurate bilayer graphene system with small twisting angles enlarges the areas of AB/BA stacking regions and simultaneously diminishes those of AA stacking regions, as well as domain walls. Subjected to an interlayer electric potential difference, the system then can accommodate topological zero-line modes across the vicinity of the width-decreased domain walls, the understanding of whose existence can be facilitated sufficiently by Jackiw−Rebbi model [51,52], which describes an exponential decaying behavior of the topologically confined states in proportion to the Dirac mass difference between the two domains separated by an abrupt domain wall. In graphene-based materials, the properties of zero-line modes have been extensively studied [53-64]. They feature as gapless dispersion in bandstructures, follow peculiar current partition law as “acute angle first” at the conducting channel intersections, and have quantized total transmission in different transport setups with various channel configurations, such as a single channel without intersection, bi-intersectional channels, as well as tri-intersectional channels [60-66]. The advantage of slightly twisted bilayer systems, in regard to forming domain walls and thus hosting zero-line modes, lies in its requirement of just a single unified gate voltage for the whole system to form a tri-intersectional channel configuration, successfully avoiding the subtleties in gate alignment [59,67].

So in this work, for topological zero-line modes to exist in the domain walls of a genuine bilayer graphene system, we first follow the Euler−Lagrange equation establishing sequence, and solve the lattice reconstruction problem using finite element method [68-72], which we believe is not limitedly only applicable to this commensurate case. Having obtained the in-plane atomic relaxation vector field, we then confirm the existence of the desired zero-line modes, and finally we study their quantum transport properties both in the absence and in the presence of a weak vertical magnetic field. We find that a total transmission/conductance plateau of G_{\text{tot}}\approx 2e^2/h can be pinpointed regardless of Aharonov−Bohm effect, but the pseudospin-conserving current partition law is relatively more sensitive. Our results suggest a twisting angle \theta \approx {0.1}^{\circ } for the described properties to be observed.

2 Atomic reconstruction problem of two-dimensional Moiré lattice solved by finite element method

2.1 Statement of the problem

Starting from an AA-stacked bilayer graphene, the mutual relative in-plane displacement between atoms in each layer consists of two parts: (i) rigid relative twisting, and (ii) lattice reconstruction, which quantitatively can be expressed as \mathit{\delta }(\mathit{x})={\mathit{\delta }}_0(\mathit{x})+{\mathit{u}}^{-}(\mathit{x}), where {\mathit{\delta }}_0(\mathit{x})= [1-R^{-1}(\theta )]\cdot \mathit{x} and R(\theta ) is the rotation matrix with a twisting angle \theta; {\mathit{u}}^{-}(\mathit{x})={\mathit{u}}^{(2)}(\mathit{x})-{\mathit{u}}^{(1)}(\mathit{x}) is the relative relaxation vector field between the two layers, labeled respectively by the parenthesized superscripts l=1,2. Only in-plane displacement is considered in this work, because it suffices for stiff substrate [73] and also is the major factor to constitute sharp domain walls. The total lattice energy contains the intralayer elastic energy and the interlayer binding energy [31,47-50]. In this continuum elastic model, relevant position-independent constants are Lamé coefficients for graphene \lambda \approx 19.67\text{eV}/a^2,\mu =57.91\text{eV}/a^2 [74,75], and lattice stacking characteristic energy density V_0=9.717\times {10}^{-3}\text{eV}/a^2 [76,77], where a\approx 0.246\text{nm} is the lattice constant for rigid monolayer graphene lattice. The above three constants in actuality provide an estimation of the twisting-angle-independent width of domain wall as w_{\text{dw}}= (a/4)\sqrt{(\lambda +\mu )/V_0}\approx 5.5\text{nm}. One of the consequences of this constant width of domain walls is that, large domains are only expectable in minimally-twisted systems. In order to minimize the total energy functional, one usually applies the celebrated Euler−Lagrange equation \frac{\delta \mathcal{L}}{\delta u^{\nu }}={\mathrm{\partial }}_{\mu }\left(\frac{\delta \mathcal{L}}{\delta ({\mathrm{\partial }}_{\mu }{u}^{\nu })}\right), and here the resultant simultaneous equations of the relative relaxation vector field between the two layers are

(1)\{\begin{array}{l}{\mathrm{\partial }}_x^2{u}_x^{-}+{\mathrm{\partial }}_{xy}^2{u}_y^{-}+{\displaystyle \frac{\mu }{\lambda +\mu }{\mathrm{\nabla }}^2{u}_x^{-}+\frac{4V_0}{\lambda +\mu }{\displaystyle \sum _{j=1}^3{g}_j(\mathit{x})b_j^x=0}}\\ {\mathrm{\partial }}_y^2{u}_y^{-}+{\mathrm{\partial }}_{xy}^2{u}_x^{-}+{\displaystyle \frac{\mu }{\lambda +\mu }{\mathrm{\nabla }}^2{u}_y^{-}+\frac{4V_0}{\lambda +\mu }{\displaystyle \sum _{j=1}^3{g}_j(\mathit{x})b_j^y=0}}\end{array},

where, for brevity, the position function g_j(\mathit{x})= \sin ({\mathit{G}}_j^{\mathrm{\top }}\cdot \mathit{x}+{\mathit{b}}_j^{\mathrm{\top }}\cdot {\mathit{u}}^{-}) is used, which contributes the nonlinearity to Eq. (1), and vectors {\mathit{b}}_j and {\mathit{G}}_j are the jth ( j=1,2,3) reciprocal vectors for monolayer graphene hexagonal lattice and commensurate twisting bilayer graphene Moiré lattice, respectively. We also have {\mathit{b}}_3=-({\mathit{b}}_1+{\mathit{b}}_2) and similar relation holds for {\mathit{G}}_js. And finally, the superscript symbol “ \mathrm{\top }” means vector/matrix transposition.

Because it enables us to study both bulk electron bandstructures and quantum transport properties of a tri-intersection-channel setup, we content ourselves to solve Eq. (1) within a Moiré Wigner−Seitz primitive cell (denoted as domain D, and its boundary a regular hexagon as \mathrm{\partial }D) with {\text{C}}_6 symmetry and an AA-stacked region as its center, as well as each Bernal-stacked region centered at its six corners. Then the problem becomes a boundary value problem of Eq. (1) within D, subjected to a proper boundary condition on \mathrm{\partial }D.

2.2 Preparation: Boundary condition, initial trial function, and strategic discretization

Inspired by the one-dimensional Frenkel−Kontorova model [31,76], a simple vanishing boundary condition

(2)\begin{array}{rl}\text{BC I:}{\mathit{u}}^{-}(\mathit{x})& =\mathbf{0},\mathit{x}\in \mathrm{\partial }D\cup \{\mathbf{0}\}\end{array}

initializes the solution seeking attempt. But just with a quick observation of its result, especially in the vicinity of each saddle point, it is recognized as insufficient in minimizing the total energy functional for this two-dimensional problem. So a more intricate boundary condition is in demand. As is stated in Ref. [39], one possible atomic relaxation pattern consists of two kinds of local relative rotations: (i) counterclockwise in the vicinity of AA stacking region (center of the Moiré primitive cell), and (ii) clockwise around the Bernal stacking regions (six corners of the Moiré primitive cell), as shown by the red arc arrows in Fig.1(b). Therefore we can devise a specific boundary condition conveying these key features above. First the envelope function we choose is

(3)

where k_1=\pi /(2\pi -w),k_2=\pi /w and w is a parameter controlling the width of the middle transitional window. Actually one can notice that f_{\pi }(x)=\sin x. Specifically, this envelope function shows curves depicted in Fig.1(a). For simplicity we further assume that {\mathit{u}}^{-}(\mathrm{\partial }D)\perp \mathrm{\partial }D, then finally we can have the boundary condition in this form:

(4)\text{BC II:}\{\begin{array}{ll}{\mathit{u}}^{-}(\mathit{x})={\displaystyle \frac{f_w(\vartheta )}{deno}}(\cos \phi ,\sin \phi {)}^{\mathrm{\top }},& x\in \mathrm{\partial }D\\ {\mathit{u}}^{-}(\mathit{x})=\mathbf{0},& x\in \{\mathbf{0}\}\end{array},

where the two relevant angles are \vartheta =6[\arctan (y/x)-\theta /2], \theta is the commensurate twisting angle, and \phi = \lfloor \vartheta /(\pi /3)\rfloor \pi /3+\pi /6+\theta /2, \lfloor \cdots \rfloor being the floor function; and typically an amplitude controlling factor can be deno=1/0.12. Starting from \mathrm{\partial }D, {\mathit{u}}^{-}(\mathrm{\partial }D) in this boundary condition is shown as the black arrows in Fig.1(b), where for the sake of a better visualization the lengths of arrows are enlarged from its actual values. In this way the boundary condition given by Eq. (4) forms the right local rotational pattern at every Bernal stcaking center, additionally it also satisfies the requirement of spatial periodicity of a triangular lattice of the commensurate Moiré pattern.

Unlike the indispensability of an appropriate boundary condition for a boundary value problem in principle, we employ further another two techniques, not necessary but constructive, to assist the solution to cases with narrow domain walls with respect to the edge length of a Moiré primitive cell. One is the trial solution function and the other is a specific discretization tactic of the target domain, both of which are based on observation about the symmetry property of the system. The initial trial function bears an in-plane full rotation symmetry ( C_{\mathrm{\infty }}) with an anticlockwise rotation to meet the requirement near the AA stcaking center, as shown by the dashed circle with arrow heads in Fig.1(b). Specifically, its quantitative expression is

(5)\begin{array}{rl}{\mathit{u}}_0^{-}(\mathit{x})& ={\displaystyle \frac{(-y,x{)}^{\mathrm{\top }}}{2{\rho }_0}}\exp \left(-{\displaystyle \frac{3}{2}}\left|{\displaystyle \frac{\sqrt{x^2+y^2}}{{\rho }_0}}-1\right|\right),\end{array}

where {\rho }_0 is a tunable constant where the magnitude of the field maximizes, as shown in Fig.2(a). As for the discretization, we notice that some parts of the domain need special attention. We find that if samplings are made more dense in the vicinity of the central AA stcaking region, each line connecting AA stcaking region to saddle point (AA-SP), and the domain boundary ( \mathrm{\partial }D), as shown in Fig.2(b), satisfactory solution can be obtained with less effort.

2.3 Solution: In-plane vector relaxation field

Having gathered all the pieces of the above preparation endeavors, we finally can find the relaxation field within a commensurate twisted bilayer graphene Moiré primitive cell. Exemplary results are shown in Fig.3, where comparison between the two boundary conditions BC I [Eq. (2)] and BC II [Eq. (4)] can be made.Fig.3(a1) shows the spatial dependence within a Moiré primitive cell of magnitude of relaxation field from BC I. As the boundary condition has required, it vanishes on all the borders of the hexagon as well as at the center, starting from which and going outward one can see how it makes the transition from a full rotation symmetry to a sixfold one. It also can be seen that atomic reconstruction occurs with the largest extend around the AA stcaking area, but will not exceed \sim a/2, which means for each monolayer the relaxation displacement is smaller than \sim a/4. The corresponding lattice interlayer binding energy density within the range [-3,6]V_0 (whose lower and upper limits correspond to pure Bernal and AA stcaking cases, respectively) is plotted in Fig.3(a2), where although Bernal stacking areas cover most of the domain, a deficiency can be identified that the energy profile in the vicinity of every saddle point (SP), the midpoint of each hexagon fringe, has not been minimized sufficiently. Fortunately, amendment can be achieved by the utility of BC II, as shown in Fig.3(b1) and (b2). Furthermore because the relaxation field is vector-valued in nature, we also check its direction scenario in Fig.3(c), where the field directions are indeed as what has been expected from the initial relaxation pattern described previously. Armed with these relaxation fields, because of the symmetry between the two layers regarding relaxation so that they two should get an equal share from the total, therefore we can then get the relaxed coordinates of atoms in each layer of twisted bilayer graphene Moiré lattice via

(6)\begin{array}{rl}{\mathit{x}}_{\pm }^{\prime }& ={\mathit{x}}_{\pm }\pm {\displaystyle \frac{1}{2}}{\mathit{u}}^{-}({\mathit{x}}_{\pm }),\end{array}

where {\mathit{x}}_{\pm } is the in-plane coordinates of atoms in the rigid hexagonal lattice before atomic reconstruction in the upper/lower layer, respectively.

3 Effects of atomic reconstruction on bandstructure

With relaxed atomic coordinates given by Eq. (6), we can build up tight-binding model for the system by using overlapping integral from Slater−Koster method [32-35,78]

(7)\begin{array}{rl}-t(\mathit{d})& =V_{\text{pp}\sigma }(d)(\hat{\mathit{d}}\cdot e_z{)}^2+V_{\text{pp}\pi }(d)[1-(\hat{\mathit{d}}\cdot e_z{)}^2],\end{array}

where d is the length of the relative position vector \mathit{d} between two atoms, \hat{\mathit{d}}=\mathit{d}/d is the unit vector, and the exponentially decaying hopping functions are (V_{\text{pp}\sigma }(d),V_{\text{pp}\pi }(d))=(V_{\text{pp}\sigma },V_{\text{pp}\pi })\exp \left[-(d-a_{\text{CC}})/r_0\right], where a_{\text{CC}}=a/\sqrt{3} is the distance between two nearest carbon atoms in a rigid monolayer graphene lattice, and r_0 is a model specific character length constant. This model is suitable for both intralayer and interlayer hoppings. Typical bandstructures are displayed in Fig.4, where the first row of figures recovers important features reported previously that as the twisting angle decreases flat bands and subgaps appear, without applying interlayer potential difference. Furthermore, as we desire to investigate the zero-line modes living in the domain walls under interlayer potential difference, we then try to find them as shown in the second row of figures, where as the interlayer potential difference increases, for the \theta ={1.05}^{\circ } case, the subgaps disappear gradually, and unfortunately no obvious signal of zero-line modes can be identified. As we can see later, although it is called “magic angle”, {1.05}^{\circ } is still too large to furnish a decent arena to accommodate topological zero-line modes.

4 Construction of a six-terminal quantum transport device

Actually as mentioned previously, atom displacements due to relaxation in most part of the Moiré primitive cell are well smaller than a/4, so noticeable difference can be hardly observed from a single monolayer of lattice before and after relaxation. But when the two relaxed hexagonal lattices of monolayers are combined as a whole twisted bilayer system, the effect of relaxation in real space becomes sensible. As shown in Fig.5(a), AA stcaking area and Bernal stacking areas simultaneously have been respectively shrunk and swollen considerably, as a consequence of which the AB−BA transitional regions are also diminished, i.e., sharp domain walls are formed accordingly. Now we attach six semi-infinite leads to a relaxed Moiré primitive cell as described above as the central scattering region, together making it a six-terminal quantum transport device. The bilayer leads are constructed from the relaxed atoms near the borders of the primitive cell, so they can constitute armchair nanoribbon and in-gap zero-line modes can be found when interlayer potential difference is applied [Fig.5(b) and (c)], despite the fact that for a same central scattering region there are two kinds of leads, which also is a smoking gun of the existence of zero-line modes in the central scattering region when the same condition is granted.

5 Quantum transport properties of topological zero-line modes

We commence the quantum transport investigation by exploitation of Green’s function method with Landauer−Büttiker formalism [79-88], as well as making use of a special method of division of the central scattering region [89]:

(8)\begin{array}{rl}{T}_{pq}& =\text{Tr}({\Gamma }_{p}G{\Gamma }_q{G}^{†}),\end{array}

(9)\begin{array}{rl}{A}_p& ={\displaystyle \frac{1}{2\pi }}G{\Gamma }_p{G}^{†},\end{array}

where T_{pq} is the transmission from the qth to the pth leads, and conductance is related to it proportionally by G_{pq}=(e^2/h)T_{pq} so we thereafter use transmission and conductance interchangeably; the bare symbol G=[(E+\text{i}{0}^{+})-H-\Sigma {]}^{-1} is the retarded Green's function for a given energy E, the total retarded self-energy \Sigma ={\mathrm{\Sigma }}_p{\Sigma }_p is the sum of the ones contributed from each lead; {\Gamma }_p=\text{i}({\Sigma }_p-{\Sigma }_p^{†}) is the line width function due to the pth lead, and the diagonal elements of the spectrum function A_p gives the local density of states in real space.

Now we treat lead 1 in Fig.5(a) as the input, from which to other leads the transmissions and real-space local density of states near the charge neutrality point ( E_{\text{cn}}\approx 0.15V_{\text{pp}\pi }) are studied. Without applying a magnetic field, the results of branch transmission coefficients G_{i1}(i=2,3,4,5,6), and the total conductance G_{\text{tot}}=\sum _{i=2}^6{G}_{i1} as well as the corresponding local density of states for different twisting angles are shown in Fig.6. From the transmission curves [Fig.6(a1), (b1), and (c1)], despite the change of twisting angles and fluctuations due to finite-size effect, some common features can be extracted, especially for the energy range E\in [0.12,0.18]V_{\text{pp}\pi }, that: (i) outgoing transmissions from leads 3 and 5 are basically the same and strongly suppressed G_{31}\approx G_{51}\approx 0, (ii) those from leads 2 and 6 are essentially in phase with each other and they oscillate with respect to the incident energy, (iii) in contrast, that from lead 4 changes with respect to incident energy out of phase from those of leads 2 and 6, and (iv) the total transmission is mostly quantized at G_{\text{tot}}\approx 2. And we can reasonably extrapolate that in the limit case with infinitely narrow domain walls ( w_{\text{dw}}\to 0) compared with the Moiré lattice length, or minimally twisted cases, there would be the branch relations G_{31}=G_{51}=0, G_{21}=G_{61}, and exactly quantized total transmission G_{\text{tot}}=2 [64]. However, opposite to the situation presented by a monolayer model in Ref. [64], around the charge neutrality point, G_{41} (blue curves in Fig.6) assumes local maximum while G_{21} and G_{61} (red and magenta curves in Fig.6) take local minimum. Besides supporting the transmission curves, the corresponding density of states [Fig.6(a2), (b2), and (c2)] exhibits lit-up regions surrounding as well as a little faraway from the AA stcaking center, which indicates that propagation from one channel to another does not have to take place at the channel intersection region, and there is at least a portion (maybe small) of conducting occurs in the area originally deemed insulating.

Next we take the effect of a weak magnetic field into consideration by Peierls substitution of the overlapping integrals of atomic orbitals for the central scattering region [90-96]

(10)\begin{array}{rl}{t}_{\mathit{r}{\mathit{r}}^{\prime }}& \to t_{\mathit{r}{\mathit{r}}^{\prime }}{\text{e}}^{\text{i}{\varphi }_{\mathit{r}{\mathit{r}}^{\prime }}},\end{array}

(11)\begin{array}{rl}{\varphi }_{\mathit{r}{\mathit{r}}^{\prime }}& ={\displaystyle \frac{q}{\hslash }}{\int }_{{\mathit{r}}^{\prime }}^{\mathit{r}}\text{d}\mathit{l}\cdot \mathit{A}(\mathit{x}),\end{array}

where \mathit{A}(\mathit{x}) is magnetic vector potential changing slowly over space and its line integral contributes an extra phase {\varphi }_{\mathit{r}{\mathit{r}}^{\prime }} to the wavefunction of an electron. Because it is a finite system without translational symmetry in any direction, the choice of magnetic vector potential is nearly arbitrary as long as its curl provides a proper uniform magnetic induction in the z-direction. If \mathit{A}=Bx{\mathit{e}}_y, it can be found out that {\varphi }_{\mathit{r}{\mathit{r}}^{\prime }}=2\pi {\tilde{\varphi }}_B(y-y^{\prime })(x+x^{\prime })/2,{\tilde{\varphi }}_B=a^{2}B/(h/q). So for {\tilde{\varphi }}_B=1/10000 the corresponding resultant counterpart to Fig.6 is shown in Fig.7, where the effect of the weak magnetic field can be summarized as: (i) magnetic field brings about a stronger fluctuation of the transmission curves, (ii) the output signals of leads 3 and 5 are no longer heavily suppressed, (iii) the in-phase behavior between leads 2 and 6 is weakened, but (iv) the total transmission can be still regarded as quantized at G_{\text{tot}}\approx 2.

Last but also important, in the construct of transport device, we have used a larger characteristic energy density than the realistic one V_0\to 20V_0, which would lead to a narrower domain wall for a given twisting angle, because the elastic model gives a domain wall width without twisting angle dependence but proportional to the reciprocal square root of V_0: w_{\text{dw}}\propto 1/\sqrt{V_0}, and makes the number of atoms in the central scattering region be not very large. On the other hand, the Moiré lattice constant depends on the commensurate twisting angle as L\propto 1/\sin (\theta /2). A smaller w_{\text{dw}}/L ratio means actually a smaller \theta. So correspondingly the twisting angles are up to a transform

(12)\begin{array}{r}\theta \to 2\arcsin \left({\displaystyle \frac{\sin (\theta /2)}{2\sqrt{5}}}\right),\end{array}

and the real twisting angles shown respectively in Fig.6(a1), (b1), and (c1) are actually {0.122}^{\circ },{0.092}^{\circ }, and {0.074}^{\circ }; in Fig.7(a1) and (b1) are {0.122}^{\circ } and {0.046}^{\circ }.

6 Summary and conclusion

To summarize, we investigate the properties of topological zero-line modes when they participate quantum transport processes in bilayer graphene system with commensurate changeable twisting angles. To that end, we must at first confirm their existence in the system under study. Besides the applying an interlayer electric potential difference, the lattice reconstruction is indispensable. Therefore, in the first place we study the lattice relaxation problem with a continuous elastic model through finite element method and successfully find the in-plane vector relaxation fields with appropriate boundary conditions and obtain the relaxed Moiré lattice. This kind of lattice relaxation in the momentum space leads to flat bands and subgaps as expected, and in the real space results in the formation of narrow domain wall, which is constructive to the presence of topological zero-line modes. So armed with all these, we finally study the branch transmissions as well as the total conductance of a six-terminal transport device constructed from a whole Moiré primitive cell and find that the branch transmissions obey consistently the “acute angle first” law as a result of conservation of pseudospin, and the nearly quantized total conductance at G_{\text{tot}}\approx 2e^2/h. The former can be affected be a weak magnetic field but the latter is considerably sturdy. Furthermore, our results recommend a twisting angle at about \theta \approx {0.1}^{\circ } for the interesting behaviors to be detected.

References

[1] Y. Cao , V. Fatemi , S. Fang , K. Watanabe , T. Taniguchi , E. Kaxiras , P. Jarillo-Herrero . Unconventional superconductivity in magic-angle graphene superlattices. Nature, 2018, 556( 7699): 43

[2] H. C. Po , L. Zou , A. Vishwanath , T. Senthil . Origin of Mott insulating behavior and superconductivity in twisted bilayer graphene. Phys. Rev. X, 2018, 8( 3): 031089

[3] E. Codecido , Q. Wang , R. Koester , S. Che , H. Tian , R. Lv , S. Tran , K. Watanabe , T. Taniguchi , F. Zhang , M. Bockrath , C. N. Lau . Correlated insulating and superconducting states in twisted bilayer graphene below the magic angle. Sci. Adv., 2019, 5( 9): eaaw9770

[4] H. S. Arora , R. Polski , Y. Zhang , A. Thomson , Y. Choi , H. Kim , Z. Lin , I. Z. Wilson , X. Xu , J. H. Chu , K. Watanabe , T. Taniguchi , J. Alicea , S. Nadj-Perge . Superconductivity in metallic twisted bilayer graphene stabilized by WSe₂. Nature, 2020, 583( 7816): 379

[5] Z. Hao , A. M. Zimmerman , P. Ledwith , E. Khalaf , D. H. Najafabadi , K. Watanabe , T. Taniguchi , A. Vishwanath , P. Kim . Electric field–tunable superconductivity in alternating-twist magic-angle trilayer graphene. Science, 2021, 371( 6534): 1133

[6] Y. Cao , V. Fatemi , A. Demir , S. Fang , S. L. Tomarken , J. Y. Luo , J. D. Sanchez-Yamagishi , K. Watanabe , T. Taniguchi , E. Kaxiras , R. C. Ashoori , P. Jarillo-Herrero . Correlated insulator behaviour at half-filling in magic-angle graphene superlattices. Nature, 2018, 556( 7699): 80

[7] M. J. Calderon , E. Bascones . Correlated states in magic angle twisted bilayer graphene under the optical conductivity scrutiny. npj Quantum Mater., 2020, 5 : 57

[8] C. Zhang , C. P. Chuu , X. Ren , M. Y. Li , L. J. Li , C. Jin , M. Y. Chou , C. K. Shih , Interlayer couplings . Moiré patterns, and 2D electronic superlattices in MoS₂ /WSe₂ hetero-bilayers. Sci. Adv., 2017, 3( 1): e1601459

[9] C. Jin , E. C. Regan , A. Yan , M. I. B. Utama , D. Wang , S. Zhao , Y. Qin , S. Yang , Z. Zheng , S. Shi , K. Watanabe , T. Taniguchi , S. Tongay , A. Zettl , F. Wang . Observation of Moiré excitons in WSe₂/WS₂ heterostructure superlattices. Nature, 2019, 567( 7746): 76

[10] K. L. Seyler , P. Rivera , H. Yu , N. P. Wilson , E. L. Ray , D. G. Mandrus , J. Yan , W. Yao , X. Xu . Signatures of Moiré-trapped valley excitons in MoSe₂/WSe₂ heterobilayers. Nature, 2019, 567( 7746): 66

[11] K. Tran , G. Moody , F. Wu , X. Lu , J. Choi , K. Kim , A. Rai , D. A. Sanchez , J. Quan , A. Singh , J. Embley , A. Zepeda , M. Campbell , T. Autry , T. Taniguchi , K. Watanabe , N. Lu , S. K. Banerjee , K. L. Silverman , S. Kim , E. Tutuc , L. Yang , A. H. MacDonald , X. Li . Evidence for Moiré excitons in van der Waals heterostructures. Nature, 2019, 567( 7746): 71

[12] H. Liu , Y. Zong , P. Wang , H. Wen , H. Wu , J. Xia , Z. Wei . Excitons in two-dimensional van der Waals heterostructures. J. Phys. D Appl. Phys., 2021, 54( 5): 053001

[13] L. Zhang , Z. Zhang , F. Wu , D. Wang , R. Gogna , S. Hou , K. Watanabe , T. Taniguchi , K. Kulkarni , T. Kuo , S. R. Forrest , H. Deng . Twist-angle dependence of Moiré excitons in WS₂/MoSe₂ heterobilayers. Nat. Commun., 2020, 11( 1): 5888

[14] S. Brem , C. Linderalv , P. Erhart , E. Malic . Tunable phases of Moiré excitons in van der Waals heterostructures. Nano Lett., 2020, 20( 12): 8534

[15] Z. Song , Z. Wang , W. Shi , G. Li , C. Fang , B. A. Bernevig . All magic angles in twisted bilayer graphene are topological. Phys. Rev. Lett., 2019, 123( 3): 036401

[16] K. P. Nuckolls , M. Oh , D. Wong , B. Lian , K. Watanabe , T. Taniguchi , B. A. Bernevig , A. Yazdani . Strongly correlated Chern insulators in magic-angle twisted bilayer graphene. Nature, 2020, 588( 7839): 610

[17] S. Wu , Z. Zhang , K. Watanabe , T. Taniguchi , E. Y. Andrei . Chern insulators, van Hove singularities and topological flat bands in magic-angle twisted bilayer graphene. Nat. Mater., 2021, 20( 4): 488

[18] P. Stepanov M. Xie T. Taniguchi K. Watanabe X. Lu A. H. MacDonald B. A. Bernevig D. K. Efetov, Competing zero-field Chern insulators in superconducting twisted bilayer graphene, arXiv: 2012.15126 (2020)

[19] C. Repellin , T. Senthil . Chern bands of twisted bilayer graphene: Fractional Chern insulators and spin phase transition. Phys. Rev. Res., 2020, 2( 2): 023238

[20] M. J. Park , Y. Kim , G. Y. Cho , S. B. Lee . Higher-order topological insulator in twisted bilayer graphene. Phys. Rev. Lett., 2019, 123( 21): 216803

[21] B. Liu , L. Xian , H. Mu , G. Zhao , Z. Liu , A. Rubio , Z. F. Wang . Higher-order band topology in twisted Moiré superlattice. Phys. Rev. Lett., 2021, 126( 6): 066401

[22] M. J. Park , S. Jeon , S. B. Lee , H. C. Park , Y. Kim . Higher-order topological corner state tunneling in twisted bilayer graphene. Carbon, 2021, 174 : 260

[23] J. M. B. Lopes dos Santos , N. M. R. Peres , A. H. Castro Neto . Graphene bilayer with a twist: Electronic structure. Phys. Rev. Lett., 2007, 99( 25): 256802

[24] R. Bistritzer , A. H. MacDonald . Moiré bands in twisted double-layer graphene. Proc. Natl. Acad. Sci. USA, 2011, 108( 30): 12233

[25] S. Carr , S. Fang , Z. Zhu , E. Kaxiras . Exact continuum model for low-energy electronic states of twisted bilayer graphene. Phys. Rev. Res., 2019, 1( 1): 013001

[26] L. Balents . General continuum model for twisted bilayer graphene and arbitrary smooth deformations. SciPost Phys., 2019, 7 : 048

[27] F. Guinea , N. R. Walet . Continuum models for twisted bilayer graphene: Effect of lattice deformation and hopping parameters. Phys. Rev. B, 2019, 99( 20): 205134

[28] V. T. Phong , E. J. Mele . Obstruction and interference in low-energy models for twisted bilayer graphene. Phys. Rev. Lett., 2020, 125( 17): 176404

[29] M. Koshino , N. N. T. Nam . Effective continuum model for relaxed twisted bilayer graphene and Moiré electron−phonon interaction. Phys. Rev. B, 2020, 101( 19): 195425

[30] E. Y. Andrei , A. H. MacDonald . Graphene bilayers with a twist. Nat. Mater., 2020, 19( 12): 1265

[31] N. N. T. Nam , M. Koshino . Lattice relaxation and energy band modulation in twisted bilayer graphene. Phys. Rev. B, 2017, 96( 7): 075311

[32] V. Perebeinos , J. Tersoff , P. Avouris . Phonon-mediated interlayer conductance in twisted graphene bilayers. Phys. Rev. Lett., 2012, 109( 23): 236604

[33] G. Trambly de Laissardière , D. Mayou , L. Magaud . Numerical studies of confined states in rotated bilayers of graphene. Phys. Rev. B, 2012, 86( 12): 125413

[34] P. Moon , M. Koshino . Energy spectrum and quantum Hall effect in twisted bilayer graphene. Phys. Rev. B, 2012, 85( 19): 195458

[35] M. Anđelković , L. Covaci , F. M. Peeters . DC conductivity of twisted bilayer graphene: Angle-dependent transport properties and effects of disorder. Phys. Rev. Mater., 2018, 2( 3): 034004

[36] H. C. Po , L. Zou , T. Senthil , A. Vishwanath . Faithful tight-binding models and fragile topology of magic-angle bilayer graphene. Phys. Rev. B, 2019, 99( 19): 195455

[37] F. Gargiulo , O. V Yazyev . Structural and electronic transformation in low-angle twisted bilayer grapheme. 2D Mater., 2017, 5 : 015019

[38] W. Yan , W. Y. He , Z. D. Chu , M. Liu , L. Meng , R. F. Dou , Y. Zhang , Z. Liu , J. C. Nie , L. He . Strain and curvature induced evolution of electronic band structures in twisted graphene bilayer. Nat. Commun., 2013, 4( 1): 2159

[39] H. Yoo , R. Engelke , S. Carr , S. Fang , K. Zhang , P. Cazeaux , S. H. Sung , R. Hovden , A. W. Tsen , T. Taniguchi , K. Watanabe , G. C. Yi , M. Kim , M. Luskin , E. B. Tadmor , E. Kaxiras , P. Kim . Atomic and electronic reconstruction at the van der Waals interface in twisted bilayer graphene. Nat. Mater., 2019, 18( 5): 448

[40] H. Shi , Z. Zhan , Z. Qi , K. Huang , E. van Veen , J. Á. Silva-Guillén , R. Zhang , P. Li , K. Xie , H. Ji , M. I. Katsnelson , S. Yuan , S. Qin , Z. Zhang . Large-area, periodic, and tunable intrinsic pseudo-magnetic fields in low-angle twisted bilayer graphene. Nat. Commun., 2020, 11( 1): 371

[41] Y. W. Liu , Y. Su , X. F. Zhou , L. J. Yin , C. Yan , S. Y. Li , W. Yan , S. Han , Z. Q. Fu , Y. Zhang , Q. Yang , Y. N. Ren , L. He . Tunable lattice reconstruction, triangular network of chiral one-dimensional states, and bandwidth of flat bands in magic angle twisted bilayer graphene. Phys. Rev. Lett., 2020, 125( 23): 236102

[42] K. Uchida , S. Furuya , J. I. Iwata , A. Oshiyama . Atomic corrugation and electron localization due to Moiré patterns in twisted bilayer graphenes. Phys. Rev. B, 2014, 90( 15): 155451

[43] S. K Jain , V. Juričić , G. T Barkema . Structure of twisted and buckled bilayer grapheme. 2D Mater., 2016, 4 : 015018

[44] S. Dai , Y. Xiang , D. J. Srolovitz . Twisted bilayer graphene: Moiré with a twist. Nano Lett., 2016, 16( 9): 5923

[45] M. M. van Wijk , A. Schuring , M. I. Katsnelson , A. Fasolino . Relaxation of Moiré patterns for slightly misaligned identical lattices: Graphene on graphite. 2D Mater., 2015, 2 : 034010

[46] H. Goldstein C. Poole J. Safko, Classical Mechanics, 3rd Ed., Higher Education Press, Beijing, 2005

[47] X. Lin , H. Zhu , J. Ni . Pressure-induced gap modulation and topological transitions in twisted bilayer and twisted double bilayer graphene. Phys. Rev. B, 2020, 101( 15): 155405

[48] P. San-Jose , A. Gutierrez-Rubio , M. Sturla , F. Guinea . Spontaneous strains and gap in graphene on boron nitride. Phys. Rev. B, 2014, 90( 7): 075428

[49] P. San-Jose , A. Gutierrez-Rubio , M. Sturla , F. Guinea . Electronic structure of spontaneously strained graphene on hexagonal boron nitride. Phys. Rev. B, 2014, 90( 11): 115152

[50] H. Suzuura , T. Ando . Phonons and electron−phonon scattering in carbon nanotubes. Phys. Rev. B, 2002, 65( 23): 235412

[51] R. Jackiw , C. Rebbi . Solitons with fermion number 1/2. Phys. Rev. D, 1976, 13( 12): 3398

[52] S. Q. Shen, Topological Insulators: Dirac Equation in Condensed Matter, 2nd Ed., Springer-Verlag GmbH, Singapore, 2017

[53] L. J. Yin , H. Jiang , J. B. Qiao , L. He . Direct imaging of topological edge states at a bilayer graphene domain wall. Nat. Commun., 2016, 7( 1): 11760

[54] T. Hou , G. Cheng , W. K. Tse , C. Zeng , Z. Qiao . Topological zero-line modes in folded bilayer graphene. Phys. Rev. B, 2018, 98( 24): 245417

[55] Y. Ren , Z. Qiao , Q. Niu . Topological phases in two-dimensional materials: A review. Rep. Prog. Phys., 2016, 79( 6): 066501

[56] Y. T. Zhang , Z. Qiao , Q. F. Sun . Detecting zero-line mode in bilayer graphene via the quantum Hall effect. Phys. Rev. B, 2013, 87( 23): 235405

[57] L. Ju , Z. Shi , N. Nair , Y. Lv , C. Jin , J. Jr Velasco , C. Ojeda-Aristizabal , H. A. Bechtel , M. C. Martin , A. Zettl , J. Analytis , F. Wang . Topological valley transport at bilayer graphene domain walls. Nature, 2015, 520( 7549): 650

[58] X. Bi , J. Jung , Z. Qiao . Role of geometry and topological defects in the one-dimensional zero-line modes of graphene. Phys. Rev. B, 2015, 92( 23): 235421

[59] J. Li , K. Wang , K. J. McFaul , Z. Zern , Y. Ren , K. Watanabe , T. Taniguchi , Z. Qiao , J. Zhu . Gate-controlled topological conducting channels in bilayer graphene. Nat. Nanotechnol., 2016, 11( 12): 1060

[60] K. Wang , Y. Ren , X. Deng , S. A. Yang , J. Jung , Z. Qiao . Gate-tunable current partition in graphene-based topological zero lines. Phys. Rev. B, 2017, 95( 24): 245420

[61] Y. Ren , J. Zeng , K. Wang , F. Xu , Z. Qiao . Tunable current partition at zero-line intersection of quantum anomalous Hall topologies. Phys. Rev. B, 2017, 96( 15): 155445

[62] K. Wang , T. Hou , Y. Ren , Z. Qiao . Enhanced robustness of zero-line modes in graphene via magnetic field. Front. Phys., 2019, 14( 2): 23501

[63] T. Hou , Y. Ren , Y. Quan , J. Jung , W. Ren , Z. Qiao . Metallic network of topological domain walls. Phys. Rev. B, 2020, 101( 20): 201403

[64] T. Hou , Y. Ren , Y. Quan , J. Jung , W. Ren , Z. Qiao . Valley current splitter in minimally twisted bilayer graphene. Phys. Rev. B, 2020, 102( 8): 085433

[65] Z. Qiao , J. Jung , C. Lin , Y. Ren , A. H. MacDonald , Q. Niu . Current partition at topological channel intersections. Phys. Rev. Lett., 2014, 112( 20): 206601

[66] Z. Yan , T. Hou , Y. Han , X. Xu , Z. Qiao . Electronic properties of zero-line modes in bilayer graphene: An ab initio study. Phys. Rev. B, 2022, 105( 3): 035425

[67] M. Kim , J. H. Choi , S. H. Lee , K. Watanabe , T. Taniguchi , S. H. Jhi , H. J. Lee . Valley-symmetry-preserved transport in ballistic graphene with gate-defined carrier guiding. Nat. Phys., 2016, 12( 11): 1022

[68] M. A. Bhatti, Advanced Topics in Finite Element Analysis of Structures: With Mathematica and MATLAB Computations, John Wiley & Sons, 2006

[69] A. Logg K. A. Mardal G. Wells, Automated Solution of Differential Equations by the Finite Element Method, 1st Ed., Springer-Verlag GmbH, Berlin Heidelberg, 2012

[70] O. C. Zienkiewicz R. L. Taylor J. Z. Zhu, The Finite Element Method: Its Basis and Fundamentals, 7th Ed., Butterworth-Heinemann Elsevier, Amsterdam, 2013

[71] J. N. Reddy, An Introduction to Nonlinear Finite Element Analysis: With Applications to Heat Transfer, Fluid Mechanics, and Solid Mechanics, 1st Ed., Oxford University Press, Oxford, 2014

[72] D. Logan, A First Course in the Finite Element Method, 1st Ed., Cengage Learning, Boston, MA, 2017

[73] S. T. Carr, Moire Patterns in 2D Materials, PHD thesis, Harvard University, 2020

[74] K. V. Zakharchenko , M. I. Katsnelson , A. Fasolino . Finite temperature lattice properties of graphene beyond the quasiharmonic approximation. Phys. Rev. Lett., 2009, 102( 4): 046808

[75] J. Jung , A. M. DaSilva , A. H. MacDonald , S. Adam . Origin of band gaps in graphene on hexagonal boron nitride. Nat. Commun., 2015, 6( 1): 6308

[76] A. M. Popov , I. V. Lebedeva , A. A. Knizhnik , Y. E. Lozovik , B. V. Potapkin . Commensurate-incommensurate phase transition in bilayer graphene. Phys. Rev. B, 2011, 84( 4): 045404

[77] I. V. Lebedeva , A. A. Knizhnik , A. M. Popov , Y. E. Lozovik , B. V. Potapkin . Interlayer interaction and relative vibrations of bilayer graphene. Phys. Chem. Chem. Phys., 2011, 13( 13): 5687

[78] J. C. Slater , G. F. Koster . Simplified LCAO method for the periodic potential problem. Phys. Rev., 1954, 94( 6): 1498

[79] S. Datta, Electronic Transport in Mesoscopic Systems, 1st Ed., World Publishing Corporation, Beijing, 2004

[80] S. Datta, Quantum Transport: Atom to Transistor, 1st Ed., Beijing World Publishing Corporation, Beijing, 2007

[81] Z. Qiao, Charge and spin transport in two-dimensional mesoscopic systems, PHD thesis, University of Hong Kong, 2009

[82] D. K. Ferry S. M. Goodnick J. Bird, Transport in Nanostructures, 2nd Ed., Cambridge University Press, Cambridge, 2009

[83] K. Hirose, Quantum Transport Calculations for Nanosystems, 1st Ed., CRC Press, Boca Raton, 2014

[84] M. D. Ventra, Electrical Transport in Nanoscale Systems, 1st Ed., Cambridge University Press, Cambridge, 2016

[85] D. Ryndyk, Theory of Quantum Transport at Nanoscale, 1st Ed., Springer-Verlag GmbH, Switzerland, 2015

[86] M. P. L. Sancho , J. M. L. Sancho , J. Rubio . Quick iterative scheme for the calculation of transfer matrices: Application to Mo(100). J. Phys. F Met. Phys., 1984, 14( 5): 1205

[87] M. P. L. Sancho , J. M. L. Sancho , J. M. L. Sancho , J. Rubio . Highly convergent schemes for the calculation of bulk and surface Green functions. J. Phys. F Met. Phys., 1985, 15( 4): 851

[88] Z. Z. Yu , G. H. Xiong , L. F. Zhang . A brief review of thermal transport in mesoscopic systems from nonequilibrium Green’s function approach. Front. Phys., 2021, 16( 4): 43201

[89] Z. Qiao , J. Wang . A variant transfer matrix method suitable for transport through multi-probe systems. Nanotechnology, 2007, 18( 43): 435402

[90] R. Peierls . Zur Theorie des Diamagnetismus von Leitungselektronen. Z. Phys., 1933, 80 : 763

[91] G. H. Wannier . Dynamics of band electrons in electric and magnetic fields. Rev. Mod. Phys., 1962, 34( 4): 645

[92] E. I. Blount . Bloch electrons in a magnetic field. Phys. Rev., 1962, 126( 5): 1636

[93] W. Kohn . Theory of Bloch electrons in a magnetic field: The effective Hamiltonian. Phys. Rev., 1959, 115( 6): 1460

[94] J. M. Luttinger . The effect of a magnetic field on electrons in a periodic potential. Phys. Rev., 1951, 84( 4): 814

[95] D. R. Hofstadter . Energy levels and wave functions of Bloch electrons in rational and irrational magnetic fields. Phys. Rev. B, 1976, 14( 6): 2239

[96] R. P. Feynman R. B. Leighton M. Sands, The Feynman Lectures on Physics, The New Millennium Ed., Basic Books, New York, 2011

RIGHTS & PERMISSIONS

Higher Education Press