# Topological Superconductivity in Twisted Multilayer Graphene

###### Abstract

We study a minimal Hubbard model for electronically driven superconductivity in a correlated flat mini-band resulting from the superlattice modulation of a twisted graphene multilayer. The valley degree of freedom drastically modifies the nature of the preferred pairing states, favoring spin triplet order with a valley singlet structure. We identify two candidates in this class, which are both topological superconductors. These states support half-vortices carrying half the usual superconducting flux quantum , and have topologically protected gapless edge states.

Recent experimentsCao et al. (2018a, b); Chen et al. (2018) demonstrate remarkable correlation phenomena in twisted multi-layer graphene with small twist angles, for which the resulting moiré pattern induces an effective triangular superlattice with a unit cell much larger than the microscopic one. The superlattice generally induces mini-bands with a reduced superlattice Brillouin zone. It was theoretically predicted that flat mini-bands should exist in such systems, an effect especially pronounced near “magic angles” in bilayer systems Bistritzer and MacDonald (2011); Suárez Morell et al. (2010); Fang and Kaxiras (2016); Trambly de Laissardière et al. (2012). When the mini-band at the Fermi energy is much narrower than the effective Coulomb interaction energy per electron, then correlation effects may be expected. Experiments on bilayersCao et al. (2018a) and trilayersChen et al. (2018) find evidence for a correlated Mott insulating state when such a mini-band contains an integer number of electrons per superlattice unit cell. Furthermore, gate tuning the charge density away from the half-filling bilayer moiré Mott insulator with 2 electrons per unit cell led to superconductivity with strong coupling characteristicsCao et al. (2018b). Many features are strikingly similar to those of the cuprate high-T materials, for which superconductivity also occurs in close proximity to a Mott insulator. This raises the intriguing possibility of graphene moiré superlattices serving as a new platform for unconventional superconductivity with unprecedented in-situ tunability. The goal of the current work is to understand the nature of the observed superconducting phases. We argue that even in the simplest situation, the valley degree of freedom of graphene leads to dramatic modifications to the superconductivity: the preferred states are topological superconductors with a valley singlet structure.

Our results are based on the minimal description of a correlated flat band in terms of a Hubbard model, with a single “site” per unit cell. This is valid when the superlattice period is large, and when the inter-band mixing may be neglected. For such a flat band, the (weak) tunneling between nearest-neighbor unit cells dominates the kinetic energy. The large period suppresses interactions beyond nearest-neighbor sites. Furthermore, each unit cell effectively hosts two degenerate orbital wave functions for electrons, which correspond to the two original valleys at the Brillouin zone corners, since the large unit cell moiré modulation cannot mix these states due to their large momentum space separation. Our starting point is therefore a two-orbital Hubbard model on the triangular lattice, with in total four flavors of single-electron states on each site, including both the spin and orbital degrees of freedom:

(1) |

Eq. 1 has an SU(4) symmetry which corresponds to the rotation between the four flavors of electron states.

This symmetry is justified as follows. For the hopping term, SU(2) spin-rotation invariance requires the hopping to be spin-independent. Mixing between different orbital states is prohibited by the large valley separation in momentum space. The reality and equality of the hopping amplitudes for the two different orbitals follows, at least for twisted bilayer graphene (Fig. 1), by careful consideration of rotation, reflection (which exchanges the valleys) and reflection . Thus the SU(4) symmetry of the hopping term in Eq. 1 should be an excellent approximation. The SU(4) symmetry of the term follows from its dependence only on the total charge of a site, which physically represents the capacitive energy of a superlattice unit cell due to “medium range” Coulomb interactions, i.e. on scales large compared to the microscopic lattice spacing but small compared to the screening length. Corrections to this SU(4) symmetry arising from short-range interactions do exist and will be considered later, but are weaker than the dominant SU(4) part by a factor proportional to , where is the superlattice spacing and is the microscopic lattice spacing.

When the number of electrons per site of Eq. (1) is or , and when is sufficiently large, the system becomes a Mott insulator. An effective SU(4) Heisenberg model for the Mott insulator can be derived using the standard perturbation theory based on the Hubbard model Eq. 1:

(2) |

where , and with are fifteen Hermitian matrices that form the fundamental representation of the SU(4) Lie-algebra (we choose ). The SU(4) spin model Eq. 2 itself is already an interesting subject to study, and compared with SU(2) spin systems, it is more likely to support exotic spin liquid ground states Wu et al. (2003); Wu (2005); Xu and Wu (2008); Hermele and Gurarie (2011); Savary (2015); Hermele et al. (2009, 2005); Mishra et al. (2002); Gorshkov et al. (2010); Kaul (2011); Yamada et al. (2017); Natori et al. (2018). But in this work we will focus on the superconductor phase next to the Mott insulator after doping.

The Heisenberg interaction in Eq. 2 can be rewritten in a different form (a Fierz identitySavary et al. (2017)):

(3) |

where we defined the 6 component and 10 component pairing fields symmetric and anti-symmetric, respectively, in . Obviously, the anti-ferromagnetic interactions () appropriate near half-filling favor condensing the operators , which are “even parity” in this sense, and we henceforth neglect the odd parity channel. In fact, transforms as an SO(6) vector (SU(4)SO(6)), when written in an appropriate basis:

(4) |

where , and . The six-component vector can be decomposed into a three component spin-triplet and orbital-singlet pairing vector , and another three component spin-singlet and orbital-triplet pairing vector . With the SU(4) symmetry of Eq. 1 and Eq. 2, these two sets of three-component vectors are exactly degenerate.

Upon doping (for instance hole-doping), one can turn on a kinetic term on Eq. 3, and the large limit becomes a t-J model with a projection that prohibits more than two particles per site. In an intermediate coupling scenario we can simply add to the Hamiltonian to represent the effects of antiferromagnetic fluctuations. Then a standard mean field theory leads to a condensate of , i.e. superconductivity. Due to Fermi statistics, this even parity pairing is antisymmetric in SU(4) flavor space, which is the essence of superexchange that favors antiferromagnetism. Amongst the even parity channels, wave pairing is penalized by the large on-site Hubbard interaction, and we expect wave pairing to be favored. Previous studies for SU(2) superconductors on the triangular lattice found that, to ensure the entire Fermi surface is gapped, pairing is often favored Wang et al. (2004); Zhou and Wang (2008); Chen et al. (2013); Watanabe et al. (2007); Weber et al. (2006); Watanabe et al. (2006, 2004).

Now let us consider the effects of SU(4) symmetry-breaking perturbations to the Hubbard model. The dominent effects arise from interactions, which are analogous to Kanamori terms multi-orbital Hubbard systems. As is usually the case for transition metal ions, we assume that the most important of these is the Hunds coupling

(5) |

where and is the total spin on site . The Hunds coupling is expected to further prefer the three-component spin-triplet and orbital singlet pairing vector over the other three components of the SO(6) vector . To see this, consider two nearest neighbor sites that are both doped with one hole, each site is occupied by one electron, at second order perturbation theory in . Suppose the two electrons form a spin-singlet and orbital triplet state, then the virtual intermediate state contains one doubly occupied site which increases the energy relative to two singly occupied sites by ; while if the two electrons form a spin-triplet state, then the virtual intermediate state has energy , which is lower than the previous case due to the Hunds interaction Eq. 5 (Fig. 2). Thus the Hunds coupling will select the spin-triplet and orbital-singlet components from the SO(6) vector , and the energy splitting is at the order of . The analysis of the electron-doped case leads to the same conclusion. Instead of the two site argument, one may alternatively just consider the modification of the super-exchange interaction of Eq. (3) by the term. This leads to a ferromagnetic contribution purely in the spin sector , which by a similar Fierz identity favors triplet pairing.

It is noteworthy that the valley degree of freedom allows the formation of an even parity (d-wave) spin triplet state, which is impossible due to Fermi statistics for a single orbital model. Here it occurs because the orbital singlet is anti-symmetric. However, in our discussion we defined the parity and angular momentum of the pair with respect to the two-orbital Hubbard model. Microscopically, parity also exchanges the two valleys, so in terms of the large microscopic Brillouin zone, the even parity d-wave state becomes an odd-parity f-wave one. We stick with the former convention for concreteness.

Knowing that the system favors spin-triplet or pairing, most generally we can write the spin triplet Cooper pair matrix in the BdG Hamiltonian as

(6) |

where

(7) |

Here , , are all three-component real vectors under spin SO(3) rotation, and . Other types of spin-triplet superconductors, for example , with real vectors , do not have a uniform maximal gap on the Fermi surface, and are thus less favorable within mean field theory than types A and B.

Type A and B states are degenerate within the standard BCS mean field theory. This is apparent from comparing for example the type A state with and the type B state and . In the former, both spin up and down electrons experience pairing, while in the latter, the pair field for up spin electrons is , and the pair field for down electrons is . The gap magnitudes are everywhere identical in the two cases, and hence they have the same mean field energy. This is the consequence of an artificial symmetry in the mean field formalism: a reflection symmetry on spin down electrons only, which interchanges the two types of pairings. Taking the most general form of the pairing order parameter with complex vectors , , the BCS mean field theory generates a Landau-Ginzburg free energy

(8) |

The last term maintains the degeneracy between type A and type B pairing, but disfavors other types of pairings. In general, effects beyond the BCS treatment will generate additional terms in the Landau-Ginzburg free energy and lift the degeneracy between type A and B. We will not attempt to resolve which state is favored here, but simply discuss the properties of both candidate states.

Consider time-reversal symmetry, which flips spin and exchanges the two valleys, hence . It also induces complex conjugation, so it acts on the order parameter as . Thus type A pairing breaks time-reversal symmetry because under complex conjugation, while type B pairing is time-reversal invariant.

Now consider the topology of the order parameter. Within a single time-reversal sector, the type A state has the ground state manifold . Here corresponds to the configuration of the spin SO(3) vector , corresponds to the configuration of . The full order parameter is invariant under a transformation

(9) |

Due to this, type A pairing supports a half-vortex, analogous to that in the polar state of spin-1 bosons in cold atom systems Ho (1998); Ohmi and Machida (1998); Mukerjee et al. (2006). After tracing along a full circle around the half-vortex core, both and acquire a minus sign (while remains single valued). The half-vortex carries a quantized magnetic flux

(10) |

which is half of the magnetic flux quantum of ordinary superconductors. Moreover, as was discussed in Ref. Mukerjee et al. (2006), in this purely two dimensional superconductor, the Mermin-Wagner theorem dictates that SO(3) vector is disordered at infinitesimal temperature due to thermal fluctuations. Hence the system no longer has long-range or even quasi-long-range order of . Instead, what persists are power-law correlations of a spin-singlet charge- order parameter . The Kosterlitz-Thouless transition out of this algebraic charge- superconducting phase is driven by unbinding of half-vortices, which leads to a universal superconductor phase stiffness jump at the transition Mukerjee et al. (2006).

While type B pairing does not break time-reversal symmetry, it has similar finite temperature behavior. The vectors and are disordered immediately by infinitesimal temperature, and the system effectively becomes an algebraic charge superconductor with a half-vortex that carries magnetic flux.

Both type A and type B superconductors are topological, in the sense that they both have gapless edge states at their boundary. In the type A superconductor, the boundary has eight channels of chiral Majorana fermions, which in the ideal case leads to a thermal Hall conductance

(11) |

The edge states of the type-A superconductor are stable against any disorder and interaction because they are chiral and hence no backscattering can occur. In the type A superconductor, because the spin symmetry is spontaneously broken down to U(1), one spin component is still conserved: for , this is . In this case, it is convenient to introduce a new basis of fermion, for orbital (valley) , define ; for orbital , define , . Then the entire mean field Bogoliubov-de Gennes Hamiltonian for quasiparticles reads

(12) |

In this basis, spin-up and spin-down fermions , both have Hall conductivity , which is visible in Eq. (12) because the pair field acts in the orbital space (second index of ) as a vector in the plane which winds twice around the origin in momentum space. Hence the eight channels of chiral Majorana fermion edge states can be reorganized into two channels of chiral edge states each for and . Thus the system also has a “spin quantum Hall” conductance : namely, if we couple the system to a “spin gauge field” , and spin-up, spin-down electrons carry gauge charge under the spin gauge field , then after integrating out all the fermions, the system generates a level-4 Chern-Simons term for the background spin gauge field:.

In the type B superconductor, the boundary has four channels of counter propagating non-chiral Majorana fermions, and there is no thermal Hall effect. The stability of the edge states of type-B superconductor deserves a bit more discussion. Let us again take , and , then this superconductor can be simply viewed as spin-up electrons and spin-down electrons forming and topological superconductors separately, and its edge state Hamiltonian reads

(13) |

The order of and fully breaks SO(3) spin symmetry, while the symmetry in Eq. 9 (a product of rotation in the spin and charge sectors) is preserved. The symmetry acts on the quasiparticles of the superconductor as a fermion parity for the right-moving fermion only: , , which also prohibits any mixing between left and right moving modes. Without interactions, the classification of this topological superconductor is obviously . Even including interactions that preserve this symmetry, the edge state in Eq. 13 with four channels of nonchiral Majorana fermions is still topologically stable, namely it cannot be gapped out without breaking the symmetry Fidkowski and Kitaev (2010, 2011); Ryu and Zhang (2012); Qi (2013); Yao and Ryu (2013); Gu and Levin (2014).

In this work we considered electronically driven superconductivity in graphene moiré superlattices within a minimal single band triangular lattice Hubbard model with spin+orbital degeneracy. We found that the valley degree of freedom of graphene has qualitative effects on the superconductivity compared to single-orbital Hubbard systems, favoring topological paired spin-triplet states. With SU(2) spin-rotation symmetry, these states support exotic charge pairing and half-vortices at non-zero temperature. Thus graphene may become not only a venue for strong correlation physics, but also topological superconductivity.

Subsequent to the posting of the first version of this preprint, a number of papers appeared emphasizing the physics related to Dirac band crossings between a pair of flat mini-bands appearing in some models of the moiré superlatticePo et al. (2018); Yuan and Fu (2018). In this situation the two low energy bands are intertwined and there is no obvious separation between them. Consequently, in the strong coupling limit, the minimal description is two band honeycomb lattice model, for which both the kinetic energy and interactions are more complicated than in Eq. (1). In the intermediate correlation regime, states near the Fermi energy dominate the superconductivity, and it is not clear that either the additional band far from the Fermi energy or any symmetry protected Dirac point plays a major role. The major difference between Ref.Po et al. (2018) and our own work is that the former invokes () ferromagnetism, while we rely upon the effective anti-ferromagnetic interaction, Eq. (2). Naïve arguments beginning from the former perspectivePo et al. (2018) appear to favor more conventional singlet superconductivity, which need not be topological. However, ferromagnetism in Hubbard models is notoriously hard to find, and anti-ferromagnetism may be more robust even when the honeycomb description is more appropriate. The compelling simplicity of the triangular framework suggests that graphene moiré heterostructures which realize the single band triangular regime are favorable for realizing topological physics. This constitutes a design goal which is realizable theoretically and experimentally.

Further studies should address these states quantitatively, the possibility of quantum spin liquid physics in the Mott states, and the effects of perturbations to the minimal Hubbard descriptions such as disorder, magnetic fields, and more. The low energy scale of these graphene superlattices allows vastly larger tuning of doping and magnetic field axes in comparison to conventional correlated transition metal compounds, and their pure two-dimensionality makes probing strictly Zeeman effects also possible. Our results may serve as guidance for such future studies.

###### Acknowledgements.

LB thanks Andrea Young for explanations of flat band physics in graphene moiré superlattices, and Lucile Savary for discussions of group theory, pairing and Fierz identities in multi-orbital systems. The authors are supported by the NSF materials theory program through grant DMR1506119 (LB), and the David and Lucile Packard Foundation and NSF Grant No. DMR-1151208 (CX).## References

- Cao et al. (2018a) Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, et al., Nature p. doi:10.1038/nature26154 (2018a).
- Cao et al. (2018b) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature p. doi:10.1038/nature26160 (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, arXiv preprint arXiv:1803.01985 (2018).
- Bistritzer and MacDonald (2011) R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011), ISSN 0027-8424, eprint http://www.pnas.org/content/108/30/12233.full.pdf, URL http://www.pnas.org/content/108/30/12233.
- Suárez Morell et al. (2010) E. Suárez Morell, J. D. Correa, P. Vargas, M. Pacheco, and Z. Barticevic, Phys. Rev. B 82, 121407 (2010), URL https://link.aps.org/doi/10.1103/PhysRevB.82.121407.
- Fang and Kaxiras (2016) S. Fang and E. Kaxiras, Phys. Rev. B 93, 235153 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.93.235153.
- Trambly de Laissardière et al. (2012) G. Trambly de Laissardière, D. Mayou, and L. Magaud, Phys. Rev. B 86, 125413 (2012), URL https://link.aps.org/doi/10.1103/PhysRevB.86.125413.
- Wu et al. (2003) C. Wu, J.-p. Hu, and S.-c. Zhang, Phys. Rev. Lett. 91, 186402 (2003), URL https://link.aps.org/doi/10.1103/PhysRevLett.91.186402.
- Wu (2005) C. Wu, Phys. Rev. Lett. 95, 266404 (2005), URL https://link.aps.org/doi/10.1103/PhysRevLett.95.266404.
- Xu and Wu (2008) C. Xu and C. Wu, Phys. Rev. B 77, 134449 (2008), URL https://link.aps.org/doi/10.1103/PhysRevB.77.134449.
- Hermele and Gurarie (2011) M. Hermele and V. Gurarie, Phys. Rev. B 84, 174441 (2011), URL https://link.aps.org/doi/10.1103/PhysRevB.84.174441.
- Savary (2015) L. Savary, arXiv preprint arXiv:1511.01505 (2015).
- Hermele et al. (2009) M. Hermele, V. Gurarie, and A. M. Rey, Phys. Rev. Lett. 103, 135301 (2009), URL https://link.aps.org/doi/10.1103/PhysRevLett.103.135301.
- Hermele et al. (2005) M. Hermele, T. Senthil, and M. P. A. Fisher, Phys. Rev. B 72, 104404 (2005), URL https://link.aps.org/doi/10.1103/PhysRevB.72.104404.
- Mishra et al. (2002) A. Mishra, M. Ma, and F.-C. Zhang, Phys. Rev. B 65, 214411 (2002), URL https://link.aps.org/doi/10.1103/PhysRevB.65.214411.
- Gorshkov et al. (2010) A. V. Gorshkov, M. Hermele, V. Gurarie, C. Xu, P. S. Julienne, J. Ye, P. Zoller, E. Demler, M. D. Lukin, and A. M. Rey, Nature Physics 289, 6 (2010).
- Kaul (2011) R. K. Kaul, Phys. Rev. B 84, 054407 (2011), URL https://link.aps.org/doi/10.1103/PhysRevB.84.054407.
- Yamada et al. (2017) M. G. Yamada, M. Oshikawa, and G. Jackeli, arXiv:1709.05252 (2017).
- Natori et al. (2018) W. M. H. Natori, E. C. Andrade, and R. G. Pereira, arXiv:1802.00044 (2018).
- Savary et al. (2017) L. Savary, J. Ruhman, J. W. Venderbos, L. Fu, and P. A. Lee, Physical Review B 96, 214514 (2017).
- Wang et al. (2004) Q.-H. Wang, D.-H. Lee, and P. A. Lee, Phys. Rev. B 69, 092504 (2004), URL https://link.aps.org/doi/10.1103/PhysRevB.69.092504.
- Zhou and Wang (2008) S. Zhou and Z. Wang, Phys. Rev. Lett. 100, 217002 (2008), URL https://link.aps.org/doi/10.1103/PhysRevLett.100.217002.
- Chen et al. (2013) K. S. Chen, Z. Y. Meng, U. Yu, S. Yang, M. Jarrell, and J. Moreno, Phys. Rev. B 88, 041103 (2013), URL https://link.aps.org/doi/10.1103/PhysRevB.88.041103.
- Watanabe et al. (2007) T. Watanabe, H. Yokoyama, Y. Tanaka, and J. ichiro Inoue, Physica C: Superconductivity and its Applications 463-465, 152 (2007), ISSN 0921-4534, proceedings of the 19th International Symposium on Superconductivity (ISS 2006), URL http://www.sciencedirect.com/science/article/pii/S09214534070%10751.
- Weber et al. (2006) C. Weber, A. Läuchli, F. Mila, and T. Giamarchi, Phys. Rev. B 73, 014519 (2006), URL https://link.aps.org/doi/10.1103/PhysRevB.73.014519.
- Watanabe et al. (2006) T. Watanabe, H. Yokoyama, Y. Tanaka, J. ichiro Inoue, and M. Ogata, Journal of Physics and Chemistry of Solids 67, 95 (2006), ISSN 0022-3697, spectroscopies in Novel Superconductors 2004, URL http://www.sciencedirect.com/science/article/pii/S00223697050%04129.
- Watanabe et al. (2004) T. Watanabe, H. Yokoyama, Y. Tanaka, J. ichiro Inoue, and M. Ogata, Journal of the Physical Society of Japan 73, 3404 (2004), eprint https://doi.org/10.1143/JPSJ.73.3404, URL https://doi.org/10.1143/JPSJ.73.3404.
- Ho (1998) T.-L. Ho, Phys. Rev. Lett. 81, 742 (1998), URL https://link.aps.org/doi/10.1103/PhysRevLett.81.742.
- Ohmi and Machida (1998) T. Ohmi and K. Machida, Journal of the Physical Society of Japan 67, 1822 (1998), eprint https://doi.org/10.1143/JPSJ.67.1822, URL https://doi.org/10.1143/JPSJ.67.1822.
- Mukerjee et al. (2006) S. Mukerjee, C. Xu, and J. E. Moore, Phys. Rev. Lett. 97, 120406 (2006), URL https://link.aps.org/doi/10.1103/PhysRevLett.97.120406.
- Fidkowski and Kitaev (2010) L. Fidkowski and A. Kitaev, Phys. Rev. B 81, 134509 (2010).
- Fidkowski and Kitaev (2011) L. Fidkowski and A. Kitaev, Phys. Rev. B 83, 075103 (2011).
- Ryu and Zhang (2012) S. Ryu and S.-C. Zhang, Phys. Rev. B 85, 245132 (2012).
- Qi (2013) X.-L. Qi, New J. Phys. 15, 065002 (2013).
- Yao and Ryu (2013) H. Yao and S. Ryu, Phys. Rev. B 88, 064507 (2013).
- Gu and Levin (2014) Z.-C. Gu and M. Levin, Phys. Rev. B 89, 201113 (2014), URL https://link.aps.org/doi/10.1103/PhysRevB.89.201113.
- Po et al. (2018) H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, arXiv preprint arXiv:1803.09742 (2018).
- Yuan and Fu (2018) N. F. Yuan and L. Fu, arXiv preprint arXiv:1803.09699 (2018).