Electronic properties of twisted bilayer graphene in high-energy k · p-Hamiltonian approximation

A model of twisted bilayer graphene has been oﬀered on the base of developed quasi- relativistic approach with high energy k · p -Hamiltonian. Monolayer-graphene twist is accounted as a perturbation of monolayer-graphene Hamiltonian in such a way that at a given point of the Brillouin zone there exists an external non-Abelian gauge ﬁeld of another monolayer. Majorana-like resonances have been revealed in the band structure of model at a magic rotation angle θ M = 1 . 05 ◦ . The simulations have also shown that a superlattice energy gap existing at a rotation angle 1 . 08 ◦ vanishes at a rotation angle 1 . 0 ◦ .


Introduction
Experimentally it has been shown an appearance of unconventional superconductivity at a magic rotation angle θ M = 1.05 • . 1 A feature of the unconventional superconductivity is accompanying insulator states such as flat bands. There is a "superlattice gap" at a rotation angle θ G = 1.1 • . 2 A tuneable band gap of the H. V. Grushevskaya et al. superlattice makes it possible to design graphene electronics elements. The interference Moire pattern of electron density appears at magic rotation angles due to coherence of the phase states of commensurately placed monolayers rotated relative to each other.Commensurate arrangement of monolayers relative to each other provides a support of exceptional flatness which has been demonstrated to reduce strain fluctuations, the main source of scattering and also the culprit for density inhomogeneities. At small angles θ (θ < 2.0 • ) of rotation of one monolayer relative to another one, interaction of Dirac cones from different monolayers is strongly decreased due to only tunnel interlayer coupling. 3 In this case of electron-hole symmetry, the nontrivial topology of the practically undistorted monolayer structure experimentally manifests itself as van Hove singularities of DOS. Hybridization of electronic states from different monolayers is experimentally observed at energies above the theoretically predicted energies of the van Hove singularities. 3 Theoretically, these features are described by accounting of a third-nearest-neighbor coupling t in a single monolayer graphene and are produced by a trigonal warping of the equi-energy model lines for small t or hexagonal miniband for large t (extra six Dirac points) around the Dirac point K(K ). 4,5 But, due to the fact that the hexagonal miniband is near Γ-point, flattening of the band occurs without reduced to zero Fermi velocity in the vicinity of the K-point that is in contradiction with experiments. Theoretically, the reduced Fermi velocity can be obtained in the presence of a Moire pattern provided interlayer interaction in a model of two monolayers graphene is "switched on." Such a model of bilayer graphene is a massive chiral model with touching bands. Electron density moire pattern for mutually deformed relaxed monolayers at small θ ≈ 1 • and less has very small reduced AA stacking region due to an energy "superlattice gap" closing at small θ ≈ 1 • . 6 However, the Fermi velocity in mutually deformed relaxed monolayers does not trends to zero, also at value θ M , although nonrelaxed deformed monolayers have a flat band at θ → 0 • . The drawback of the model of nonrelaxed deformed twisted monolayers is the high value of the first magic angle (θ M = 1.25 • ) for zero value of the Fermi velocity. The phenomenological tight-binding model of the graphene superlattice without strain, but with interlayer interaction of the graphite type in twisted layers predicts flat bands due to zero Fermi velocity at θ M = 1.05 • . 7 The DFT-consistent construction of an ab initio tight-binding graphene-superlattice model with universal form of interlayer couplings gives an agreement with the phenomenological tight-binding model through the inclusion of maximally localized Wannier orbitals, taking into account weak dipole interactions. 8 However, the predicted value of θ M (1.08 • ) in this construction is quite different from the experimental one. Besides, the results of the ab initio tight-binding diverge strongly in respect to folding approximation. Thus, modern theoretical representations on twisted bilayer graphene do allow to explain neither the appearance of flat zones at experimentally observed rotation angles θ, nor experimentally observed van Hove singularities of DOS at small θ.
In papers 9-16 it has been proposed a superlattice model of graphene on a support, for example, on hexagonal boron nitride. It has been shown that a Hamiltonian of the superlattice represents itself a Hamiltonian of monolayer graphene with corrections to pseudo-electric and pseudo-magnetic fields which strain the graphene monolayer. It can be naturally assumed that for the bilayer twisted graphene the similar approach could be of interest.
We propose to use a high-energy k · p quasi-relativistic monolayer graphene model to describe graphene superlattices. Its advantages are the proven existence of quasiparticle massless Majorana-like excitations, the presence of a dynamic mass leading to reducing of the Fermi velocity with simultaneous flattening of bands, the possibility of the existence of a hexagonal miniband in the vicinity of the Dirac point as well as the presence of a nontrivial non-Abelian Zak phase.

Theoretical Background
Graphene is a 2D semimetal hexagonal monolayer, which is comprised of two trigonal sublattices A, B. Semi-metallicity of graphene is provided by delocalization of π(p z )-electron orbitals on a hexagonal crystal cell. Since the energies of relativistic terms π * (D 3/2 ) and π(P 3/2 ) of a hydrogen-like atom are equal each other 17 there is an indirect exchange through d-electron states to break a dimer. Therefore, a quasirelativistic model monolayer graphene, besides the configuration with three dimers per the cell, also has a configuration with two dimers and one broken conjugate double bond per the cell. The high-energy k · p Hamiltonian of a quasiparticle in the sublattice, for example, A reads where ψ † −σ A |0, −σ is a spinor wave function (vector in the Hilbert space), σ = {σ x , σ y } is the 2D vector of the Pauli matrixes, p = {p x , p y } is the 2D momentum operator, Σ AB , Σ BA are relativistic exchange operators for sublattices A, B respectively; i 2 Σ AB Σ BA is an unconventional Majorana-like mass term for a quasiparticle in the sublattice A, |ψ * BA is a spinor wave function of quasiparticle in the sublattice (1) is a spin-valley-current coupling. One can see that the term with conventional mass in (1) is absent.
The exchange interaction term Σ x rel is determined as 18 Here interaction (2 × 2)-matrices ∆ AB and ∆ BA are gauge fields (or components of a gauge field). Vector-potentials for these gauge fields are determined by the phases α 0 and α ±,k , k = 1, 2, 3 of π(p z )-electron wave functions ψ pz (r) and ψ pz,±δ k (r), k = 1, 2, 3 respectively that the exchange interaction Σ x rel (3) in accounting of the nearest lattice neighbors for a tight-binding approximation reads 18-20 where the origin of the reference frame is located at a given site on the sublattice A(B), V (r) is the three-dimensional (3D) Coulomb potential, designations ψ pz, ±δi (r), ψ pz, ±δi (r 2D ) ≡ ψ pz (r ± δ i ), i = 1, 2, 3 refer to atomic orbitals of p z -electrons with 3D radius-vectors r ± δ i in the neighbor lattice sites δ i , nearest to the reference site; r ± δ i is the p z -electron 3D-radius-vector. Elements of the matrices Σ AB and Σ BA include bilinear combinations of the wave functions so that their phases α 0 and α ±,k , k = 1, 2, 3 enter into ∆ AB and ∆ BA from (4 and 5) in the form Therefore, an effective number N of flavors in our gauge field theory is equal to 3. Then owing to translational symmetry we determine the gauge fields ∆ ±,i in 2040055-4 Mod. Phys. Lett. B Downloaded from www.worldscientific.com by PROEKTNO-KONSTRUKTORSKIJ on 07/08/20. Re-use and distribution is strictly not permitted, except for Open Access articles.
Equation (1) can be rewritten as where BA . The equation similar to (15), can be also written for the sublattice B. As a result, one gets the equations of motion for a Majorana bispinor (|ψ AB , |ψ * BA ) T : 20,21

Accordingly to (9) eigenvalues E
(1) i , i = 1, 2 of (15) and, accordingly, eigenvalues E i , i = 1, . . . , 4 of the 4 × 4 Hamiltonian (16), (17) are functionals of c ± . To eliminate arbitrariness in the choice of phase factors c ± one needs a gauge condition for the gauge fields. The eigenvalues E i , i = 1, . . . , 4 are real because the system of equations (16), (17) is transformed to Klein-Gordon-Fock equation. 21 Therefore we impose the gauge condition as a requirement on the absence of imaginary parts in the eigenvalues E i , i = 1, . . . , 4 of the Hamiltonian (16) and (17): To satisfy the condition (18) in the momentum space we minimize a function | m E i | absolute minimum of which coincides with the solution of the system (18). For the mass case band structures for the sublattice Hamiltonians are the same. Therefore neglecting the mass term the cost function For the non-zero mass case, we assume the same form of the function f due to smallness of the mass correction.
Topological defect pushes out a charge carrier from its location. The operator of this non-zero displacement presents a projected position operator PrP with the projection operator P = N n=1 |ψ n,k ψ n,k | for the occupied subspace of states ψ n,k (r). Here N is a number of occupied bands, k is a momentum. Eigenvalues of PrP are called Zak phase. 22 The Zak phase coincides with a phase γ mn = i C(k) ψ m,k |∇ k | ψ n,k · dk, n, m = 1, . . . , N of a Wilson loop W mn = T exp(iγ mn ) being a path-ordered (T) exponential with the integral over a closed contour C(k). 23

High-Energy k · p-Hamiltonian for Twisted Bilayer Graphene
Let us write down a non-Abelian connection A which determines the topologically non-trivial non-Abelian Zak phase (19). Matrix elements A mn of A in a momentum space are given by the the integrand of the entering (19): where u m (k) are eigenfunctions of charge carriers in m-th hole or electron band, n, m = 1, 2. A hopping at non-zero non-Abelian gauge fields has been considered at neglecting of the Majorana-like mass term M AB (M BA ). The eigen problem has been solved for the following Hamiltonian in momentum representation: where origin for the wavevectors q is the Dirac point K(K ), the operators of relativistic exchange Σ BA depend explicitly on q and upon the introduced gauge fields α ± providing the realness of eigenvectors of the Hamiltonian.
Since the van-der-Waals interaction between the monolayers is small in comparison with an impact of topological charge density defects of the charge carriers, it has been assumed that a bilayer-graphene Hamiltonian can be represented as an unperturbed Hamiltonian of two graphene sheets with an addition term V t . The last is an operator of potential energy of charge carriers of one monolayer in an external field of topological defects of another monolayer.
A change of the flux at parallel transport of a graphene charge carrier gives a topological phase of a wave function for the carrier. For small angles of twist a Hamiltonian of bilayer graphene is the unperturbed Hamiltonian of two monolayers with correction on a perturbation of charge carrier density of one monolayer in an external non-Abelian gauge field acting from the another monolayer.
Let us choose a center of cell of hexagon in the graphene hexagonal lattice as a center of rotation. At rotation on an angle θ points in the Brillouin zone are rotated by the same angle relative to the Γ point. Therefore, the transformation of an arbitrary wave vector q has the following form: where M is an orthogonal matrix of two-dimensional rotation on angle θ, K A is the wavevector of the Dirac point. Then, the hopping V t of second-monolayer charge carriers in the non-Abelian field of first monolayer is presented though minimal coupling of gauge fields as where κ is a coupling constant for the non-Abelian gauge field. The perturbed Hamiltonian H b of our bilayer graphene model: is similar to a Hamiltonian describing a production of electron-holes pairs in a rotating electrical field. We solve the eigen problem for the Hamiltonian H b of the monolayer graphene in external field for the values of wave numbers near of the Dirac touching. Since q → 0, matrix elements of (22), constructed on eigenfunctions of the unperturbed graphene Hamiltonian are matrix elements of the first-order perturbation-theory corrections to non-perturbed energy bands of bilayer graphene.

Numerical Results and Comparison with Experiment for Majorana Quasi-Particle Excitations in Bilayer Graphene
The energy-band contribution of the changes in the external gauge-field flux has been simulated numerically for three values of the angles θ = 1.0 • , 1.05 • , 1.08 • . The corrections (22) calculated and presented in Fig. 1 are added to positive eigenvalues of the non-perturbed Hamiltonian.
The rotation of the gauge field produces three Majorana resonances as Figs. 1(a, c, (Figs. 1(a, b, e, f)). The degree of narrowing for the Majorana peaks is larger at 1.0 • than at 1.08 • . Besides, the height of the peaks decreases sharply at 1.0 • . The Dirac touching and the Dirac band at small q remain that they are nested as subgap states in the perturbed bilayer-graphene bands (Figs. 1(b, d, f)). This gap is negligibly small for θ = 1.0, 1.08 • because the resonances in the vicinity of Dirac point are absent (Figs. 1(b, f)).
Wide Majorana resonances hold in the neighborhood of Dirac point at θ = 1.05 • (Fig. 1(c, d)). The presence of these wide Majorana resonances leads to very large energy band gap signifying that there are Cooper pairs and accordingly the bilayer is in a superconducting state. Since a statistics of the charge carriers is non-Abelian this superconductivity phenomenon is a unconventional superconductivity one at θ = 1.05 • .
The physical meaning of the resonances found is a creation of additional dipole moment in such a way that the resulting dipole moment of the second monolayer is aligned along the dipole moment of the first monolayer. Majorana resonances remain sufficiently intensive ones at high θ (∼ 1.08 • ) and form an energy "superlattice gap". Closing superlattice gap stems from vanishingly small Majorana peaks at θ (∼ 1.0 • ).

Conclusion
To summarize our findings. An approach based on high-energy k · p Hamiltonian for a quasi-relativistic graphene model has been extended for the bilayer graphene description. To account of a monolayer-graphene twist, an unperturbed Hamiltonian of every monolayer gains a perturbation in such a way that at a given point of the Brillouin zone there exists an external non-Abelian gauge field of another monolayer. Corrections to the energy bands stipulated by the interactions with non-abelian gauge fields have been simulated numerically for several twist angles. Majorana-like resonances in energy bands corrections have been revealed and their appearance is discussed in the context of unconventional superconductivity in magic-angle graphene superlattices.