# Unconventional topological superconductivity and phase diagram for a two-orbital model of twisted bilayer graphene

###### Abstract

We consider the superconducting and Mott-insulating states for the twisted bilayer graphene, modeled by two narrow bands of electrons or holes with appreciable intraatomic Coulomb interactions. The interaction induces the kinetic exchange which leads to the real-space, either triplet- or singlet-spin pairing, in direct analogy to heavy-fermions and high-temperature superconductors. We construct explicitly the phase diagram as a function of electron concentration for the spin-triplet paired case, as well as determine the topological edge states. The model reproduces principal features observed experimentally in a semi-quantitative manner. The essential role of electronic correlations in driving both the Mott-insulating and superconducting transitions is emphasized.

Introduction—The nature of electron states and, in particular, the microscopic mechanism of unconventional pairing in strongly correlated matter is one of the fundamental problems in condensed matter physics. This is because in the systems such as heavy fermions, Pfleiderer (2009) high-temperature superconductors (SC), Hott et al. (2015) or selected atomic systems in optical lattice,Lewenstein et al. (2012) the interparticle interaction energy can exceed by far the single-particle (kinetic band) energy. In that situation, specific phenomena induced by the interelectronic correlations occur, such as the Mott or the Mott-Hubbard localization,Gebgard (1997) unconventional superconductivity (SC) associated with real-space pairing,Ogata and Fukuyama (2008) as well as specific magnetic behavior, such as metamagnetic transition to localized state and spin-dependent masses of quasiparticles,Spałek (2006) and quantum critical behavior.Carr (2011) In this context, the recent discovery of SC and the Mott insulating behavior Cao et al. (2018a, b); Chen et al. (2018) of twisted bilayer graphene (TBG) provides a model situation for studying phenomena ascribed to high-temperature SC. This is because TBG represents a truly two-dimensional system and the ratio of interaction amplitude to the Fermi energy can be varied experimentally in a controlled manner. What is particularly important here is that the behavior close to the SC-insulator boundary can be studied systematically by changing the carrier concentration without introducing the ubiquitous atomic disorder, as is the case in the high-temperature SC. Also, the whole concentration dependent phase diagram can be sampled by changing the gate voltage and hence, the data represent intrinsic system properties. However, there are also some differences with respect to the high- SC. The main difference is that TBG is inherently multi-orbital,Yuan and Fu (2018) whereas the high- SC may be mapped onto single-orbital models.Zegrodnik and Spałek (2017) This difference forms a basis for an extension of our SGA methodJędrak and Spałek (2011); Zegrodnik et al. (2014); Spałek and Zegrodnik (2013) to the present situation.

Model and pairing— We start with the two-orbital model () on a triangular lattice for which a single site represents one Moiré unit cell. The two orbitals correspond to the two original valleys at the Brillouin zone corners. The starting Hamiltonian isXu and Balents (2018)

(1) |

where , , and denote the hopping, and intraorbital and interorbital Coulomb interactions, respectively. The operator creates an electron with spin on orbital at site and is the orbital particle-number operator (we use the notation for , respectively). Hereafter we assume approximate spin-orbital symmetry by taking . This means that the Hund’s rule coupling is disregarded at this point (see below). Also, the interlayer hopping is neglected, since it can be shown that with such an orbitally independent form of hybridization term can be incorporated into an effective canonical structure without it.Spałek et al. (1985)

The dominant paring channels can be identified by referring to canonical perturbation expansionSpałek et al. (1985); Spałek and Chao (1980) in the manner analogous to that for the one-band Hubbard model. In the simplest case of , the interaction takes the form , where is the total particle number operator on lattice site . The configurations with large local density of electrons (say ) are thus disfavored by interactions, and eliminated by means of the canonical transformation. This procedure, after employing standard approximations, leads to the kinetic exchange taking a general functional form

(2) |

where and are Pauli matrices acting on the spin- and orbital- indices, respectively (the summation is performed over with ). We have used the compact notation , denotes transposition, and sets the effective kinetic exchange scale . Note now that, for , an additional minus sign is generated from the term by performing the transposition, rendering the interaction attractive in some pairing channels. This occurs for spin-singlet, orbital-triplet (, and ) and for spin-triplet, orbital-singlet (, and ) cases. All other pairing symmetries are disfavored. Hereafter we adopt the point of view that an additional Hund’s rule interorbital interaction that breaks full symmetry, which is not explicitly included in the original Hamiltonian (Unconventional topological superconductivity and phase diagram for a two-orbital model of twisted bilayer graphene), would tip the balance towards the spin-triplet, orbital-singlet counterpart.Xu and Balents (2018) The approximate symmetry implies, however, that those states should be energetically close to each other. In general, the spin-singlet pairing should appear as the system approaches the half-filling, since then the intraorbital antiferromagnetic kinetic exchange becomes dominant.Chao et al. (1977); Spałek and Chao (1980) However, it can be shown that upon changing the sign on , the formalism for the spin-triplet case with formally coincides with that for the spin-singlet situation. Zegrodnik and Spałek (2014)

In either case, the parity of the SC order parameter is even. Previous studies of the triangular-lattice Hubbard model Chen et al. (2013) suggest that, in this geometry, unconventional symmetry might be realized in order to optimize the condensation energy. This is due to the fact that, contrary to the usual -wave pairing, entire Fermi surface becomes then gapped upon the SC transition. We thus restrict ourselves to the case of the -type, spin-triplet, orbital-singlet, pairing, defined by the following relations between SC amplitudes

(3) | |||

(4) | |||

(5) |

where is the rotation by the angle around the axis perpendicular to the lattice plane, and going through the site . The even-parity property follows from the condition (5) for and translational symmetry. It should be noted that an alternative approach of purely real extended -wave pairing in TBG has been presented very recently, where the Eliashberg formalism has been used within which the pairing is induced due to the many-body spin- and charge-fluctuations. That paper represents a complementary weak-correlation perspective.

To investigate this scenario, instead of employing the original model (Unconventional topological superconductivity and phase diagram for a two-orbital model of twisted bilayer graphene) with the general kinetic exchange (Unconventional topological superconductivity and phase diagram for a two-orbital model of twisted bilayer graphene), we resort to a simpler effective Hamiltonian favoring spin-triplet pairing, defined as

(6) |

where is the total-spin operator on lattice site . The effective pairing coupling has been introduced. The Hamiltonian (Unconventional topological superconductivity and phase diagram for a two-orbital model of twisted bilayer graphene) reproduces correctly attractive interaction in the spin-triplet channel and thus is applicable as long as solely paramagnetic and superconducting state of symmetry defined by Eqs. (3)-(5) are considered. The effective model is illustrated in Fig. 1.

A brief methodological remark is in place. Namely, in general, the exchange part in Eq. (Unconventional topological superconductivity and phase diagram for a two-orbital model of twisted bilayer graphene) in the metallic state has the formChao et al. (1977); Spałek and Chao (1980)

(7) |

where , with and . Such a form of interaction reduces to the ferromagnetic exchange in the Kugel-Khomskii limitKugel and Khomskii (1982) () and to the antiferromagnetic one in the Mott insulating () limit. For and we can expect the change of the pairing symmetry to that with spin singlet and orbital triplet, but then the whole formalism can be formed in a similar manner.Zegrodnik and Spałek (2014) A detailed analysis of the symmetry change requires separate analysis as also the frustration effects may be an important factor in the analysis. The change of symmetry can be also obscured by the phase separation (see below).

Solution and phase diagram—We employ the statistically consistent Gutzwiller approximation (SGA) that has proven to be effective for various classes of correlated electron systems, including the high- cuprates Jędrak and Spałek (2011) and spin-triplet ferromagnetic SC. Zegrodnik et al. (2014); Spałek and Zegrodnik (2013); Kądzielawa-Major et al. (2017) At zero temperature, SGA reduces to optimization of the ground state energy within the class of trial wave functions in the formal limit of infinite spatial dimensions. Here denotes a Slater determinant (describing uncorrelated electrons) and is a product of local correlators that modify the local many particle electronic configurations by means of variational coefficients . This allows to include the effect of strong correlations on top of the renormalized quasiparticle picture.

Since the width of the narrow bands arising in the magic-angle graphene ,Cao et al. (2018a, b) the role of the local correlations is expected to be crucial. To consider this scenario in detail, we fix the parameters as , , and . The calculations are performed at small nonzero temperature for numerical purposes. This value maps onto absolute temperature scale of less than .

The calculated phase diagram for the model (Unconventional topological superconductivity and phase diagram for a two-orbital model of twisted bilayer graphene) as a function of electron concentration per superlattice site is shown in Fig. 2. In panel (a) we display the correlated SC amplitude component which is one of the principal results of the present contribution. Around the half-filling (), we obtain two asymmetric SC domes, with stronger SC correlations on the lower-concentration side of the phase diagram. Such an asymmetry is expected as triangular lattice is not bipartite and the electron-hole symmetry cannot be expected. We point out a small hump is SC amplitude emerging near the integer filling , where the correlations are again enhanced. This feature is fairly weak for the present choice of parameters and is obscured by the larger dome closer to . The obtained two-dome structure exhibits a remarkable agreement with the recent experimental data for the magic-angle bilayer graphene.Cao et al. (2018a) We emphasize that the agreement is semi-quantitative as when the dimensionless SC amplitude is scaled by the characteristic energy , the maximal gap parameter matches well the measured critical temperatures. The shaded area in the phase diagram marks the phase-separation region, where the SC state coexists with an increasing fraction of the Mott insulating phase as the half-filling is approached. We point out that, even though the maximum of SC correlations for lies within the phase separation regime, the sample-averaged SC amplitude still exhibits a two-dome structure with a proper particle-hole asymmetry. Panel (b) details doping-evolution of the hopping, measured by the nearest-neighbor correlation function , where is the renormalization factor. For kinematic reasons, these correlations drop to zero for an empty and completely filled system ( and , respectively). Due to strong correlations, this happens also close to the half-filling, where the Mott transition is reached. As is shown in panel (c), the probability of double occupancies normalized to its uncorrelated (Hartree-Fock) value is reduced accordingly. To complete the discussion of the electronic correlations, in Fig. 2(d) we plot the ratio of renormalized to uncorrelated kinetic energy which roughly describes the mass enhancement . As expected, the kinetic energy is suppressed as the Mott transition is approached.

Next, we discuss the phase separation occurring in our model. Figure 2(e) shows the doping-dependence of the chemical potential. As the half-filling is approached from either above or below, the chemical potential eventually levels off as a function of electron concentration, to become extremely steep in the close vicinity of . This behavior is reminiscent of that observed for the one-band Hubbard model,Tandon et al. (1999); Kotliar et al. (2002); Eckstein et al. (2007) where phase separation occurs as well. The value of the chemical potential in this regime is determined by Maxwell construction, which is illustrated in the insets. The squares are computational data points and red color marks unstable solutions. At present, we are unable to numerically approach part of the curve sufficiently close to to observe the full S-shaped function . The jump in for is a signature of the opening of the Hubbard gap.Hubbard (1964)

For completeness, we plot in Fig. 3 the principal characteristics for , i.e., in the weak-correlation limit. There is no sign of Mott insulating behavior, so indeed appreciable correlations are required to reproduce experimental data. With the increasing a second dome appears and becomes predominant at large (cf. Fig. 2).

Topological properties—It is established that pairing symmetry might render the system a topological SC.Chern (2016) For the model (Unconventional topological superconductivity and phase diagram for a two-orbital model of twisted bilayer graphene) we have explicitly computed the Chern number by the efficient method of Brillouin-zone triangulation Fukui et al. (2005) with the result , depending on the direction of phase winding of the order parameter. In this situation, a distinct set of topologically-protected edge states is expected for finite-size sample. To investigate the latter we have considered the lattice slab of dimensions sites with open- and periodic-boundary conditions along the shorter and longer ends, respectively. This results in the -band, one-dimensional system. The parameters were set to the same values as in previous section and the electron concentration was fixed at to stay away from the phase-separation regime that would hinder the analysis (cf. Fig. 2).

In Fig. 4(a) we plot the band structure of the effective quasiparticles, calculated within the SGA approach. Most of the bands are gapped due the superconductivity, but gapless modes crossing the Fermi level appear as well. To elucidate the nature of those states, we calculate the zero-temperature spectral functions

(8) |

where the index enumerates the sites along the shorter edge of the system ( and correspond to the two opposite sides of the sample). The operator creates a spin- electron in the -th orbital at the site . The one-dimensional wave vector results from periodic boundary conditions in the longer direction. Moreover, denote eigenvalues of the renormalized one-particle Hamiltonian emerging within the SGA calculations. In Fig. 4(b) and (c) we plot the spectral function for the , spin-up electrons at two opposite ends of the sample, and , respectively. It is apparent that the edges of the system contribute substantially to the states crossing the Fermi energy, with the opposite signs of the Fermi-velocities. The nature of these modes is finally settled in panel (d), where we display the bulk contribution to the spectral weight . The latter exhibits no intensity close to Fermi energy.

Finally, we address the evolution of the electronic correlations as one moves from the edge to the sample bulk region. In Fig. 5 we plot the probability of site double occupancy normalized to its Hartree-Fock value as a function of the site index . In the central part of the system, remains practically constant. At the edges, where the topologically protected states are located, is substantially reduced. This implies that, even though the edge states are robust, it is energetically beneficial to suppress their weight by strongly correlating orbitals near the boundary of the system.

Summary—We have provided a quantitative analysis of the two-orbital model of bilayer graphene with extremely narrow bands. The onset of the Mott insulating state at the half-filling requires the presence of relatively strong correlations, here . The absence of the electron-hole symmetry for the triangular lattice assumed here leads to the asymmetric SC domes on the electron and hole sides, with more prominent character on the low band-filling side. This feature differs from that in the case of two-dimensional model of high-temperature SC. Zegrodnik and Spałek (2017) The topological states appear naturally in the gapful state of the superconductor. Further studies would require the inclusion of electron-concentration dependence (and sign reversal) of the effective intersite exchange integral, which should lead to spin-singlet pairing with essentially the same type of behavior.

Acknowledgments—The discussions with Adam Rycerz are gratefully acknowledged. This work was supported financially by the Grant No. DEC-2012/04/A/ST3/00342 from National Science Centre (NCN) of Poland.

## References

- Pfleiderer (2009) C. Pfleiderer, “Superconducting phases of -electron compounds,” Rev. Mod. Phys. 81, 1551 (2009).
- Hott et al. (2015) R. Hott, R. Kleiner, T. Wolff, and G. Zwicknagl, “Handbook of applied superconductivity,” (Wiley-VCH, Berlin, 2015) Chap. Review on Superconducting Materials, pp. 1–50.
- Lewenstein et al. (2012) M. Lewenstein, A. Sanpera, and V. Ahufinger, Ultracold Atoms in Optical Lattices. Simulating Quantum Many-Body Systems (Oxford University Press, 2012).
- Gebgard (1997) F. Gebgard, The Mott Metal-Insualtor Transition, Springer Tracts in Modern Physics (Springer, Berlin, 1997).
- Ogata and Fukuyama (2008) M. Ogata and H. Fukuyama, “The - model for the oxide high- superconductors,” Rep. Prog. Phys. 71, 036501 (2008).
- Spałek (2006) J. Spałek, “Magnetic properties of almost localized fermions revisited: spin dependent masses and quantum critical behavior,” Phys. Stat. Sol. (b) 243, 78 (2006).
- Carr (2011) L. D. Carr, ed., Understanding Quantum Phase Transitions (CRC Press, Taylor & Francis Group, Boca Raton, 2011).
- Cao et al. (2018a) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, “Unconventional superconductivity in magic-angle graphene superlattices,” Nature 556, 43 (2018a).
- Cao et al. (2018b) Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, Ray C. Ashoori, and P. Jarillo-Herrero, “Correlated insulator behaviour at half-filling in magic-angle graphene superlattices,” Nature 556, 80 (2018b).
- Chen et al. (2018) G. Chen, L. Jiang, S. Wu, B. Lv, H. Li, K. Watanabe, T. Taniguchi, Z. Shi, Y. Zhang, and F. Wang, “Gate-Tunable Mott Insulator in Trilayer Graphene-Boron Nitride Moiré Superlattice,” (2018), arXiv:1803.01985 .
- Yuan and Fu (2018) N. F. Q. Yuan and L. Fu, “A Model for Metal-Insulator Transition in Graphene Superlattices and Beyond,” (2018), arXiv:1803.09699 .
- Zegrodnik and Spałek (2017) M. Zegrodnik and J. Spałek, “Universal properties of high-temperature superconductors from real-space pairing: Role of correlated hopping and intersite Coulomb interaction within the -- model,” Phys. Rev. B 96, 054511 (2017).
- Jędrak and Spałek (2011) J. Jędrak and J. Spałek, “Renormalized mean-field - model of high- superconductivity: Comparison to experiment,” Phys. Rev. B 83, 104512 (2011).
- Zegrodnik et al. (2014) M. Zegrodnik, J. Bünemann, and J. Spałek, “Even-parity spin-triplet pairing by purely repulsive interactions for orbitally degenerate correlated fermions,” New J. Phys. 16, 033001 (2014).
- Spałek and Zegrodnik (2013) J. Spałek and M. Zegrodnik, “Spin-triplet paired state induced by Hunds rule coupling and correlations: a fully statistically consistent Gutzwiller approach,” J. Phys.: Condens. Matter 25, 435601 (2013).
- Xu and Balents (2018) C. Xu and L. Balents, “Topological Superconductivity in Twisted Multilayer Graphene,” (2018), arXiv:1803.08057 .
- Spałek et al. (1985) J. Spałek, D. K. Ray, and M. Acquarone, “A hybridized basis for simple band structures,” Sol. State Commun. 56, 909 – 915 (1985).
- Spałek and Chao (1980) J. Spałek and K. A. Chao, “Kinetic exchange interaction in a doubly degenerate narrow band and its application to and ,” J. Phys. C 13, 5241 (1980).
- Chao et al. (1977) K. A. Chao, J. Spałek, and A. M. Oleś, “The kinetic exchange interaction in doubly degenerate narrow bands,” Phys. Stat. Sol. (b) 84, 747 (1977).
- Zegrodnik and Spałek (2014) M. Zegrodnik and J. Spałek, “Spontaneous appearance of nonzero-momentum Cooper pairing: Possible application to the iron-pnictides,” Phys. Rev. B 90, 174507 (2014).
- Chen et al. (2013) K. S. Chen, Z. Y. Meng, U. Yu, S. Yang, M. Jarrell, and J. Moreno, ‘‘Unconventional superconductivity on the triangular lattice Hubbard model,” Phys. Rev. B 88, 041103 (2013).
- Kugel and Khomskii (1982) K. I. Kugel and D. I. Khomskii, “The Jahn-Teller effect and magnetism: transition metal compounds,” Phys. Usp. 25, 231–256 (1982).
- Kądzielawa-Major et al. (2017) E. Kądzielawa-Major, M. Fidrysiak, P. Kubiczek, and J. Spałek, ‘‘Spin-triplet paired phases inside ferromagnet induced by Hund’s rule coupling and electronic correlations: Application to ,” (2017), arXiv:1712.08028 .
- Tandon et al. (1999) A. Tandon, Z. Wang, and G. Kotliar, “Compressibility of the Two-Dimensional Infinite- Hubbard Model,” Phys. Rev. Lett. 83, 2046 (1999).
- Kotliar et al. (2002) G. Kotliar, S. Murthy, and M. J. Rozenberg, “Compressibility Divergence and the Finite Temperature Mott Transition,” Phys. Rev. Lett. 89, 046401 (2002).
- Eckstein et al. (2007) M. Eckstein, M. Kollar, M. Potthoff, and D. Vollhardt, “Phase separation in the particle-hole asymmetric Hubbard model,” Phys. Rev. B 75, 125103 (2007).
- Hubbard (1964) J. Hubbard, “Electron correlations in narrow energy bands III. An improved solution,” Proc. Roy. Soc. London A 281, 401 (1964).
- Chern (2016) T. Chern, ‘‘ and wave topological superconductors and new mechanisms for bulk boundary correspondence,” AIP Advances 6, 085211 (2016).
- Fukui et al. (2005) T. Fukui, Y. Hatsugai, and H. Suzuki, “Chern Numbers in Discretized Brillouin Zone: Efficient Method of Computing (Spin) Hall Conductances,” J. Phys. Soc. Japan 74, 1674 (2005).