Unconventional topological superconductivity and phase diagram for an effective two-orbital model as applied to twisted bilayer graphene

We consider the superconducting and Mott-insulating states for twisted bilayer graphene, modeled as a two-narrow-band system of electrons with appreciable intra-atomic Coulomb interactions. The interaction induces kinetic exchange which leads to real space, either triplet- or singlet-spin pairing, in direct analogy to heavy fermions and high-temperature superconductors. By employing the statistically consistent Gutzwiller method, we construct explicitly the phase diagram as a function of electron concentration for the spin-triplet $d_{x^2-y^2}+i{d}_{xy}$ paired case, as well as determine the topological edge states. The model reproduces principal features observed experimentally in a semiquantitative manner. The essential role of electronic correlations in driving both the Mott-insulating and superconducting transitions is emphasized. The transformation of the spin-triplet state into its spin-singlet analog is also analyzed, as well as the appearance of the phase-separated superconducting+Mott-insulating state close to the half-filling.

Figure 1

Illustration of the effective two-band model of twisted bilayer graphene. The green arrows mark intraorbital repulsion $U$, interorbital repulsion $U^{\prime }$, intersite exchange interaction $J$, as well as the single-particle hopping integral $t$.

Figure 2

Phase diagram for $t<0,U=U^{\prime }=18\left|t\right|,J=-\left|t\right|$, and temperature $T={10}^{-4}\left|t\right|/k_B$, obtained for $512\times 512$ lattice. (a) Doping dependence of the dimensionless superconducting gap amplitude component ${\langle c_{i\uparrow }^{\left(1\right)}{c}_{j\uparrow }^{\left(2\right)}\rangle }_G$. The shaded area marks the phase-separation region, where the superconducting state appears separated spatially from the Mott-insulating phase emerging near the half-filling. (b) Hopping probability ${\langle c_{i\sigma }^{\left(l\right)†}{c}_{j\sigma }^{\left(l\right)}\rangle }_G$ which represents the electron itineracy. (c) Probability $d^2$ of the double-site occupancy, normalized to its Hartree-Fock value $d_{\mathrm{HF}}^2$. (d) The ratio of kinetic energies calculated within the SGA and Hartree-Fock approximations. (e) Chemical potential $\mu$ as a function of electron density. Note that $\mu$ is constant in the phase-separation region and its value is determined by Maxwell construction. The closeups of the phase-separation regions below and above the half-filling are displayed in the insets. The squares are computational data points and the red solid lines mark unstable spatially homogeneous solutions.

Figure 3

Same as in Fig. 2 except for smaller value of $U=5\left|t\right|$. The superconducting gap appears only on the overdoped side [panel (a)], and both the double occupancy (b) and the chemical potential (c) evolve in a continuous manner.

Figure 4

Band structure and spectral properties of a $40\times 256$ lattice slab with periodic and open boundary conditions along the longer and shorter edges, respectively. The model parameters are $U=U^{\prime }=18\left|t\right|,J=-\left|t\right|$, and temperature $T={10}^{-4}\left|t\right|/k_B\left(t<0\right)$. (a) Band structure presented in the one-dimensional Brillouin zone. The bulk gap arises due to $d+id$ SC, whereas the levels crossing the Fermi energy originate from the topologically protected edge states. (b), (c) Spectral functions calculated for the orbitals on two opposite shorter ends of the sample. Note that they contribute substantially to the levels crossing the Fermi energy. (d) The same as in (b) and (c), but for the bulk states at the system center. The varying intensity reflects the difference in the spectral density values.

Figure 5

Probability of the double occupancy of the lattice site $d^2$, normalized to its Hartree-Fock value $d_{\mathrm{HF}}^2$ across the transverse cross section of the $40\times 256$ lattice slab (the sites are enumerated from 1 to 40). The model parameters are $U=U^{\prime }=18\left|t\right|,J=-\left|t\right|$, temperature $T={10}^{-4}\left|t\right|/k_B$, and $t<0$.

Figure 6

Exemplary processes involving triplet states, contributing to the second-order kinetic exchange close to the quarter-filling (a) and half-filling (b). The left panels show the initial configurations of two involved sites, $i$ and $i$, whereas the right panels represent the situation after hopping takes place. In the quarter-filling case, the Hund's rule exchange reduces the energy cost of the hopping process favoring triplet configurations. On the contrary, close to the half-filling, breaking of the initial spin triplets by hopping increases the energy cost of such a process, disfavoring the triplet formation. Green arrows illustrate action of the hopping term.

Figure 7

Exemplary diagram resulting from the Wicks decomposition of $\langle {\mathrm{\Psi }}_0\left|c_{i\sigma }^{†}{\widehat{n}}_{i\overline{\sigma }}^{\mathrm{HF}}{c}_{j\sigma }{\widehat{n}}_{j\overline{\sigma }}^{\mathrm{HF}}{\widehat{d}}_{l_1}^{\mathrm{HF}}{\widehat{d}}_{l_2}^{\mathrm{HF}}\right|{\mathrm{\Psi }}_0\rangle$, which is one of the four second-order expansion terms from (C2), when ${\widehat{o}}_i=c_{i\sigma }^{†}$ and ${\widehat{o}}_j^{\prime }=c_{j\sigma }$ (i.e., when calculating the hopping contribution in the correlated state $\langle {\mathrm{\Psi }}_G\left|c_{i\sigma }^{†}{c}_{j\sigma }\right|{\mathrm{\Psi }}_G\rangle )$. The black lines connecting the vertices correspond to $P_{mn}$ expectation values. During the summation procedure in real space, the so-called internal vertices $l_1$ and $l_2$ are attached to all possible lattice sites within the region determined by the real-space cutoff $R_{\max }(|{\mathbf{R}}_l-{\mathbf{R}}_i|\le R_{\max }$ and $|{\mathbf{R}}_l-{\mathbf{R}}_j|\le R_{\max })$.

Figure 8

Ground-state energy (a), double occupancies (b), and chemical potential (c), all as a function of band filling for the case of two-band Hubbard model with $U=U^{\prime }=11.5\left|t\right|$ [cf. Eq. (1)] and $J=0$, calculated within the SGA (red squares and black triangles) and DE-GWF methods (blue dots and green diamonds, respectively). Only the paramagnetic phase has been included (no SC phase considered here). The DE-GWF calculation has been carried out up to the third diagrammatic order.

Figure 9

The same as in Fig. 8, but for $U=U^{\prime }=15.5\left|t\right|$.