Topological phase and chiral edge states of bilayer graphene with staggered sublattice potentials and Hubbard interaction

Gated heterostructures containing bilayer graphene with staggered sublattice potentials are investigated by tight binding model with Rashba spin-orbital coupling and Hubbard interaction. The topological phase diagrams depend on the combinations of substrates and the Hubbard interaction. The presence of the staggered sublattice potential favor the topological phase transition with small Rashba spin-orbital coupling strength. The presence of the Hubbard interaction modified the topological phase boundaries, increasing the minimal spin-orbital coupling strength for topological phase transition. A phase space of topological semi-metal with indirect band gap is identified in the non-interacting systems. For the bilayer graphene with different staggered sublattice potentials in the two layers, the conditions for the zigzag nanoribbons to host edge polarized chiral edge states are discussed. The conditions require moderate or vanishing Rashba spin-orbital coupling strength, as well as proper range of the gate voltage. The conditions for the systems with and without the Hubbard interaction are compared. The edge polarization can be controlled by the gate voltage.


Introduction
Graphene consisted of a single layer or few layers of carbon atoms [1,2] is a promising material for advanced electronic [3], spintronic and valleytronic devices [4]. The Weyl spinor like carriers near to the Fermi level exhibit high electronic mobility [2], long spin relaxation length and time [5][6][7][8][9][10][11][12][13]. Breaking the A-B sublattice symmetry will turn the Weyl spinors to the massive Dirac Fermion like particles. It can be realized by substrates, such as h-BN [14,15] and SiC [16], which induce staggered sublattice potentials Δ being equal to 28 meV and 130 meV, respectively. The localized edge states in the zigzag nanoribbons [17][18][19] have been extensively investigated, which are the candidate for information carrier in integrated spintronic and valleytronic devices. The localization of the edge state can be tuned by strain [20]. We explore the graphene nanoribbons that host edge polarized chiral edge states, whose localization can be controlled by the gated voltage. For this purpose, we firstly study the topological properties of the bulk bilayer graphene with staggered sublattice potentials and Hubbard interaction.
In the AB-stacked bilayer graphene(BLG) [1], an interlayer potential difference V 2 induced by a gated voltage opens a tunable band gap [21]. In the presence of the Rashba spin-orbital coupling(SOC), the gated BLG become topological insulator(TI) with topological invariant Z 2 =1 when the Rashba SOC strength becomes sufficiently large for a given gated voltage [22,23]. Meanwhile, the TI phase has valley Chern number being C V =2, implying the quantum spin Hall(QSH) phase. With the Rashba SOC strength smaller than the critical value of topological phase transition, the BLGs is in the quantum valley Hall(QVH) phase with Z 2 =0 and C V =4. The Rashba SOC strength induced by external electric fields [5,24,25] is far to be sufficient for the TI phase transition. The Rashba SOC could be enhanced by constructing a curve surface [26,27], substrate proximity effect [28][29][30][31] or adatom doping [32,33]. In the vicinity of substrate consisted of heavy metal [34][35][36], the crystal field is largely enhanced, which gives rise to sizeable Rashba SOC. The intercalation of heavy metal into the graphene hollow sites can also enhance the Rashba SOC [37]. Adding staggered sublattice potentials to both layers of the BLG change the phase diagram, which allow topological phase transition with infinitesimally small Rashba SOC strength at V≈Δ [38].
In the single-layer graphene(SLG), the presence of SOC induces topologically nontrivial phase, such as QSH phase [39]. The presence of Rashba SOC modifies the band structure and spin texture, which induces novel optical and electrical properties [40,41]. In the absence of the SOC, the one dimensional armchair nanoribbons exhibit topological band gap as well [42,43]. In comparison, the BLG have more parameters, such as the gated voltage and staggered sublattice potentials of each layer, to control the topological phase. In general, the topological phases of BLG and SLG have different valley Chern numbers.
On the other hand, the electron-electron interaction changes the physical properties of the two dimensional materials, including the topological properties. The Hubbard model has been added to the Kane-Mele (KM) model to incorporate the effect of electron correlation, and it has been found that the phase diagrams depend on the strength of the Hubbard interaction [44][45][46][47][48]. The topological invariant for the interacting systems can be calculated by integral of the Berry phase that is defined by the eigenstates of the topological Hamiltonian [49][50][51][52][53][54]. The topological Hamiltonian is the inverse of the Green's function at zero frequency. The Green's function of the interacting system can be calculated by cluster perturbation theory(CPT) [55][56][57][58][59].
In this article, the topological phase diagrams of BLGs with different types of staggered sublattice potentials are studied. The effects of electron correlation described by the Hubbard model are calculated by the CPT. The minimal requirement of Rashba SOC strength for topological phase transition is not infinitesimally small, but finite. The effects of next nearest neighbor(NNN) hopping in the tight binding model is also discussed. The conditions for the BLGs that support the edge polarized chiral edge states within the bulk gap are investigated. The paper is organized as follows. In section 2, the tight binding model of the BLGs, the CPT method and the calculation of the topological invariant are described. In section 3, the numerical result of the phase diagrams of Z 2 and C V are discussed. The phase boundaries in the presence and absence of the electron correlation are compared. In section 4, the edge polarized chiral edge states of the zigzag nanoribbons are discussed. In section 5, the conclusion is given.

Theoretical
-- Hc . . respectively. The summation with index i j , áá ññ^runs through the inter-layer NNN sites. The parameter δ N equates 1 or 0, for the presence or absence of the NNN hopping, respectively. Applying the periodic boundary condition with Bloch phase, the Hamiltonian in the absence of the Hubbard interaction becomes an eight-by-eight matrix. Eigenstates of the matrix equation give the band structure as well as wave functions and Berry phase, which are used to calculate the topological invariant.

Cluster perturbation theory
In the presence of the Hubbard interaction, the single particle eigenstate is not well defined. As a replacement, the single particle Green's function in a unit cell described the dynamic properties of the system. One of the efficient methods to calculate the Green's function with interaction is the CPT. The CPT method has four steps.
(1) One defines an isolated cluster in the lattice. The cluster that usually consists of multiple primitive cells must be able to tile the extended lattice. (2) The Green's function of the cluster with the Hubbard interaction is calculated by exact diagonalization. The Hamiltonian of the isolated cluster is given by equation (1) with all summations being restricted to the lattice sites within the cluster. Because the Hamiltonian conserves the particle number, the basis states of the many-particle interacting Hamiltonian are the Fock states. For the halffilling N-site cluster, the dimension of the Hilbert space equates to the combination C N N 2 . In this Hilbert space, the Hamiltonian is expressed as a sparse matrix, whose ground state Wñ | is found by the Lanczos algorithm [64]. The cluster Green's function is obtained by the operation where operator c m(n) is the annihilation operator with composite index m(n) for a lattice site within a primitive cell and a spin, E 0 is the ground state energy, z is the frequency of the Green's function. (3) The lattice of the BLG is covered by isolated clusters, which form a superlattice. The Green's function of the superlattice is given by where Q is the wave vector in the first Brillouin zone of the superlattice, V(Q) is the reciprocal superlattice representation of the hopping matrix between adjacent clusters. Finally, the Green's function of the original lattice is obtained as [56,57] where k is the wave vector in the first Brillouin zone of the original lattice, L is the number of primitive cells in the cluster, r i is the center location of the i-th primitive cell. For the BLGs, the number of lattice sites in a cluster must be integral multiple of four, because a primitive cell contains four lattice sites. In our numerical calculation, the isolated cluster contains eight lattice sites(two primitive cells).

Topological invariant
The topological invariant Z 2 of the band structures is defined as is the Berry connection with the summation index n covering all occupied band and u n (k) being the periodic part of the Bloch state, ) is the z component of the Berry curvature. For the non-interacting case with U=0, the Bloch states are given by diagonalization of the Hamiltonian defined in equation (1) with Bloch periodic condition. For the interacting case with U 0 ¹ , the eigenstates of the topological Hamiltonian that equates to inverse of the zero frequency Green's function, h G k, 0 [53], are used in equation (6). Although the eigenstates of h topo are not physical state, they preserve the topological properties of the interacting systems. The numerical integration in equation (6) is performed by Fukui and Hatsugai's procedure [65][66][67][68].
Integrating the Berry curvature through the whole Brillouin zone gives Chern number that is zero, because the Berry curvature is odd under time reversal. The Berry curvature has large magnitude in the vicinity of the K and K′ point, so that on can define the continuous Dirac fermion models for each valley. The Berry curvature is defined by the Dirac spinor in momentum space. By integrating the Berry curvatures of the occupied bands through the whole momentum space, k , x Î -¥ ¥ ( ) and k , y Î -¥ ¥ ( ) , the Chern number of the continuum Dirac Fermion model corresponding to K and K′ valley being denoted as C K and C K¢ , respectively, are defined as with the Berry curvatures of the n-th band in K(K′) valley, k k , , being given as The summation index n covers the four occupied valence bands. v x and v y are the velocity operator of x and y direction in the continuum Dirac Fermion model. ε n and n K K k y ñ ¢ | ( ) are the energy level and wave function of the n-th band in K(K′) valley with wave vector being k. Finally, the valley Chern number is defines as [22,23,38,69]. The valley Chern number is only calculated in the non-interacting model.

Numerical results for the phase diagrams
In this section, the numerical results for the phase diagrams of the BLGs with different options of substrates are discussed. In the absence of the Hubbard interaction, the phase boundary is determined by the gap closing condition, which can be obtained by solving the non-interacting Hamiltonian with Bloch periodic condition at K point. Two eigenvalues are −Δ +1 −V and Δ −1 +V. The other six eigenvalues are roots of two cubic equations, The gap closing conditions for different combinations of (Δ +1 , Δ −1 ) are obtained by the condition that the two cubic equations have one common root. The NNN hopping terms vanish at K points, so that the presence of the NNN hopping does not change the phase diagram of the non-interacting systems. In the presence of Hubbard interaction, the phase boundary is numerically calculated. The Hubbard interaction coefficient takes the value U=1.6t, which have been proven to accurately described the correlation effect in graphene [70,71]. The phase diagram of the suspended BLGs is plotted in figure 1. For the non-interacting systems, the phase boundary between the QSH and QVH phase is g t V R 1 3 2 2 = + [22]. The presence of the Hubbard interaction only modifies the phase boundary slightly at large gated voltage. The Hubbard interaction expand the regime of QSH at moderate gate voltage(V≈50∼300 meV). The phase boundary is driven above the x axis by the interaction, which implies that the phase transition to QSH phase requires finite gated voltage for all g R . For the interacting system, the presence of the NNN hopping changes the phase boundary for gate voltage being smaller than 50 meV.
which is plotted as solid lines in figures 2(a) and (b). At V=Δ, the phase boundary approaches the y axis of the phase diagram with g R =0, implying topological phase transition at infinitesimally small Rashba SOC strength. However, the presence of the Hubbard interaction modifies the phase boundary, as shown in figures 2(c) and (d), and increases the minimal Rashba SOC strength for topological phase transition to a finite value. Because the valley Chern number is only calculated in the non-interacting model, the phase diagram with interaction only present the phase boundary between topological trivial and nontrivial phase with Z 2 =0 and 1, respectively. Additional presence of the NNN hopping significantly modifies the phase boundaries for gate voltage with small magnitude.
For the most realistic model that includes both Hubbard interaction and NNN hopping, the phase diagrams show that the minimal Rashba SOC strength along the phase boundaries decrease as the staggered sublattice potentials are increasing. Specifically, for the BLGs with SiC substrates, the minimal Rashba SOC strength along the phase boundary is g R =81.8 meV at gate voltage being V=85.5 meV. Removing the NNN hopping at the same gate voltage allow the phase transition at smaller Rashba SOC, g R =60 meV. For the corresponding non-  interacting model, the band structures exhibit linear band crossing at the topological phase transition point, signaling band inversion, as shown in figure 3(a). The presence of the NNN hopping bring trigonal warping to the band structure, indicated by the local minimal of the first conduction and valence bands beyond the K point, as shown in figure 3(b). In the presence of the Hubbard interaction and the absence of the NNN hopping, the linear band crossing is clearly visible in figure 3(c). The band crossing occur beyond the K point. In the present of the NNN hopping, the two conduction bands and the two valence bands near to the Fermi level strongly couple together, which is indicated by the non-zero spectral function value in the interval among these bands, as shown in figure 3(d). The impact of the presence of the interaction is exhibited through the effective Hamiltonian including the self-energy, H eff =H+Σ(ω, k). The self-energy depends on the frequency and Bloch wave vector. The diagonal terms of the self-energy are equivalent to local potentials for the corresponding lattice sites and spin indexes. Among the non-diagonal terms of the self-energy, those corresponding to the first NN hopping terms has the largest magnitude. These terms are equivalent to changing the wave vector in the continuous Dirac Fermion model, so that the band crossing is beyond the K point. The presence of these effective potential and hopping terms change the global properties of the systems. When the gate voltage has large magnitude, the self-energy becomes negligible, so that the phase boundaries of the interacting and noninteracting systems are nearly the same. Flipping the sign of the staggered sublattice potential of the l=−1 layer gives the phase diagrams in figure 4. For the non-interacting systems, the phase diagrams are shown in figures 4(a) and (b) for BN and SiC substrates, respectively. The topological semi-metal(TSM) phase with Z 2 =1, C V =0 and zero indirect band gap appears. Another TSM phase was previously identified in suspended BLGs with both intrinsic and Rashba SOC in [72]. At the phase boundary between the TSM and BI phases, the gap closing occurs beyond the K point, so that the phase boundary has to be numerically calculated. In the presence of the Hubbard interaction and the NNN hopping, the phase diagrams are shown in figures 4(c) and (d). In the regime with V<Δ, the presence of the NNN hopping significantly change the phase boundaries.
The numerical result shows that the BLGs in SiC substrate with Δ +1 =Δ −1 =Δ=130 meV host the topological phase transition with the strength of the Rashba SOC as small as 81.5 meV. The relatively large SOC could be obtain by intercalation of heavy metal element [37] between the graphene and SiC substrate or between the two graphene layers. The SOC could be further enhance by applying pressure on the heterostructures that reduce the distance between the graphene and the substrate. An experimental sample of BLG that is encapsulated between two h-BN(SiC) substrates could randomly be BLGs with Δ +1 =Δ −1 or Δ +1 =−Δ −1 . The two structures with the same g R have different band gaps, so that they can be distinguished by measuring the band gap of the sample.
For the BLGs with only one substrate, for example, Δ +1 =Δ=−28 meV and Δ −1 =0, the phase diagram is similar with that in figure 2(a). For the non-interacting systems, the phase boundary between the QSH and | | | |, the phase diagrams are more complicated. Since these phase diagrams do not contains additional novel phase, the results are not shown in the paper.
The experimental implementations of the BLGs with only one substrate could be obtained by deposition of a suspended BLG on the surface of the SiC or h-BN crystal. Intercalation of heavy metal element in the suspended BLG could be performed before the deposition. After the deposition of the BLG on the substrate, direct growing of another substrate on the top surface of the BLG could be performed. Optionally, the top substrate could also be implemented by a modified scanning tunneling microscope(STM) instrument. Instead of the metallic sharp tip, the tip consisted of SiC or h-BN crystal with atomically flat surface is manufactured. Unlike the regular STM, the insulating tip does not have conducting current for measurement of distance between the tip and the BLG. Instead, the parallelity and distance between the STM flat tip and the BLG are measured by optical reflection pattern and absorption spectral, respectively. When the STM flat tip and the BLG are parallel, the Fabry-Perot interfere of optical reflection could generate strong signal for spatial adjustment. When the STM flat tip approach the BLG, 1 D + change from zero to finite value due to the Van der Waals interaction, so that the band gap of the BLG is changed. The longitudinal distance between the STM flat tip and the BLG could be monitored by measuring the optical absorption spectral. Moving the STM flat tip along the transversal direction changed the displacement between the lattices of BLG and STM flat tip, which change Δ +1 between positive and negative values. Thus, the BLG with Δ +1 =Δ −1 and Δ +1 =−Δ −1 could be obtained on demand.

Edge polarized chiral edge states of zigzag nanoribbon
In this section, the edge polarized chiral edge states within the bulk gap is studied. These states exist in the zigzag nanoribbons of BLGs that break the particle-hole symmetric and are in the QVH phase. In the QSH phase, the helical edge states always appear in both two edges. In order to construct the BLGs that support the edge polarized chiral edge states, the non-symmetric staggered sublattice potential(Δ +1 ¹ Δ −1 ), finite gate voltage, and moderate or vanishing Rashba SOC are required. We use the BLGs with Δ +1 =−Δ −1 =Δ=130 meV and g R =28 meV as example. With the paremters being away from the phase boundaries between QVH and QSH phases, the NNN hopping slightly change the band structure or spectral function, so that it is neglected in the discussion. For the interacting systems, the spectral function of a nanoribbon is obtained by the CPT method. Assuming that the nanoribbon is laid on the x-y plane and the zigzag edges extend along the y axis, the unit cell of the nanoribbon is consisted of N rectangular clusters arranging along the armchair direction along x axis. Each cluster contains eight lattices. The Green's function of the isolated unit cell is obtained by the inversion of a block tri-diagonal matrix with the diagonal block being the inverse of the Green's function of each cluster and the non-diagonal block being the hopping matrix between the neighboring block. The Green's function of the nanoribbon is obtained by equation (4).
The edge states are originated from the connection between the flat edge bands and the bulk bands. The flat edge bands are localized states at one of the lattice site at the zigzag edge, so that their energy levels are determined by the potential of the corresponding lattice sites. In the absence of the Hubbard interaction, the potentials of the lattice sites are determined by H l, l D and H l,V in the Hamiltonian. In the following discussion, the flat edge band that is localized at the left(right) zigzag edge of the l-th layer is denoted as F L,l (F R,l ), whose energy level is Δ l +lV(−Δ l +lV ), as plotted in figure 5(a). The band edges of the valence and conduction bulk bands are near to −Δ and Δ+V. Meanwhile, F L,l (F R,l ) connects to the valence(conduction) bulk bands in the band structures. Assuming positive gate voltage and Δ +1 =−Δ −1 =Δ, the energy level of F L,−1 (F R,+1 ) remains in the valence(conduction) bulk bands, forming none edge state within the bulk gap, as shown by the thin lines in figure 5(a). The energy level of F L,+1 is within the conduction bulk bands, forming a gapless edge band for arbitrary gated voltage. The weak Rashba SOC split the edge band into two bands. Tuning the gate voltage into the ranges with Δ<V, Δ<V<2Δ or 2Δ<V drives the energy level of F R,−1 into the conduction bulk band, the bulk gap or the valence bulk band, respectively. Thus, the right zigzag edge hosts none edge states, conductive edge state or gapless edge state, depending on the gate voltage. With Δ<V, the nanoribbon host only edge polarized chiral edge states within the bulk gap, as shown in figure 6(a); with Δ<V<2Δ, the nanoribbon hosts both of edge polarized chiral edge states and conductive edge states, as shown in figure 6(b). With 2Δ<V, the edge states within the bulk gap is not edge polarized. Reversing the gate voltage to negative value exchange the properties of F L,l and F R,l , so that the localization of the chiral edge states is flipped into the opposite side of the nanoribbon.
In the presence of the Hubbard interaction, the self-energy effectively change the on-site potential of each lattice site, so that the energy levels of the flat edge band and the bulk band are changed. The energy levels of the flat edge bands and the band edges of the conduction and valence bulk band are extracted from the spectral functions of the nanoribbons, which are plotted in figure 5(b). The critical values of the gated voltage are obtained from the crossing point between the energy level of F R,−1 and the band edge of the bulk bands. With V0.38Δ, the energy level of F R,−1 remains in the conduction band, so that the nanoribbon host only edge polarized chiral edge states within the bulk gap. The spectral function of the nanoribbon with V=0.38Δ is plotted in figure 6(c). The flat edge bands become weakly dispersive. With 0.38Δ<V<1.54Δ, the nanoribbon hosts both of edge polarized chiral edge states and conductive edge states within the bulk gap. The spectral function with V=0.65Δ is plotted in figure 6(d). Overall, the presence of the Hubbard interaction shrinks the range of the gate voltage that makes the nanoribbon hosting edge polarized chiral edge states, and reduce the effective bulk gap.

Conclusion
The topological phase diagrams of BLG heterostructures with optional combinations of h-BN and(or) SiC substrates and the presence of the Hubbard interaction are studied. The Green's function of the interacting systems are calculated by the CPT method. The topological invariants are calculated by employing the topological Hamiltonian that is the inverse of the Green's function at zero frequency. Comparing to the noninteracting systems, the presence of the Hubbard interaction modifies the phase boundaries between topological trivial and non-trivial phases. Specifically, in non-interacting systems, Z 2 type of topological phase transition can happen at infinitesimally small Rashba SOC strength. The presence of the Hubbard interaction in the same BLG increases the minimal Rashba SOC strength for topological phase transition to a finite value. For BLGs with Δ +1 =−Δ −1 in the absence of the Hubbard interaction, topological semi-metal(TSM) phase with zero indirect band gap and non-trivial topological invariants, i.e. Z 2 =1 and C V =0, is founded.
The conditions that the zigzag nanoribbon of the BLGs can host edge polarized chiral edge states are studied. The staggered sublattice potentials of the two layers are required to be asymmetric. For the typically choice with Δ +1 =−Δ −1 =Δ, in the absence of the Hubbard interaction, the BLGs with 0<V<Δ host pure edge polarized chiral edge states, and the BLGs with Δ<V<2Δ host the mixing of edge polarized chiral edge states and conductive edge states. In the presence of the Hubbard interaction, the range of the gate voltage for these two phases shrink. The localization of the edge polarized chiral edge states can be controlled by the sign of the gated voltage. These features would provide more feasible systems for the design of spintronic and valleytronic devices.