Universal mechanism of Ising superconductivity in twisted bilayer, trilayer and quadrilayer graphene

TL;DR

This paper reveals a universal mechanism for superconductivity in twisted multilayer graphene, driven by valley symmetry breaking and Fermi line nesting, which explains the emergence and fading of superconductivity with doping.

Contribution

It demonstrates that valley symmetry breaking and Fermi line nesting are fundamental to superconductivity across twisted bilayer, trilayer, and quadrilayer graphene, providing a unified theoretical framework.

Findings

Superconductivity arises from valley symmetry breaking in moiré systems.

Fermi line nesting induces a strong Kohn-Luttinger pairing instability.

Superconductivity diminishes as the Fermi line becomes more isotropic at higher doping.

Abstract

We show that the superconducivity in twisted graphene multilayers originates from a common feature, which is the strong valley symmetry breaking characteristic of these moir\'e systems at the magic angle. This leads to a breakdown of the rotational symmetry of the flat moir\'e bands down to $C_3$, and to ground states in which the time-reversal symmetry is broken for a given spin projection. However, this symmetry can be recovered upon exchange of spin-up and spin-down electrons, as we illustrate by means of a self-consistent microscopic Hartree-Fock resolution where the states for the two spin projections acquire opposite sign of the valley polarization. There is then a spin-valley locking by which the Fermi lines for the two spin projections are different and related by inversion symmetry. This effect represents a large renormalization of the bare spin-orbit coupling of the graphene multilayers, lending protection to the superconductivity against in-plane magnetic fields. In the twisted bilayer as well as in trilayer and quadrilayer graphene, the pairing glue is shown to be given by the nesting between parallel segments of the Fermi lines which arise from the breakdown of symmetry down to $C_3$. This leads to a strong Kohn-Luttinger pairing instability, which is relevant until the Fermi line recovers gradually a more isotropic shape towards the bottom of the second valence band, explaining why the superconductivity fades away beyond three-hole doping of the moir\'e unit cell.

Tables

Table 1. Table 1: Eigenvalues of the Cooper-pair vertex with largest magnitude and dominant harmonics grouped according to the irreducible representations of the approximate C 3 ​ v subscript 𝐶 3 𝑣 C_{3v} symmetry, for the Fermi line of TBG shown in Fig. 1 .

Eigenvalue $\lambda$	harmonics	Irr. Rep.
4.25	$1$
2.14	$\{\cos \left(\varphi \right),\sin \left(\varphi \right)\}$	E
2.09
$-0.43$	$\{\cos \left(7\varphi \right),\sin \left(7\varphi \right),$	E
$-0.40$	$\cos \left(8\varphi \right),\sin \left(8\varphi \right)\}$
$-0.36$	$\cos \left(6\varphi \right)$	$A_1$
0.33	$\cos \left(3\varphi \right)$	$A_1$
0.26	$\{\cos \left(2\varphi \right),\sin \left(2\varphi \right),$	E
0.21	$\cos \left(5\varphi \right),\sin \left(5\varphi \right)\}$

Table 2. Table 2: Eigenvalues of the Cooper-pair vertex with largest magnitude and dominant harmonics grouped according to the irreducible representations of the approximate C 3 ​ v subscript 𝐶 3 𝑣 C_{3v} symmetry, for the Fermi line of TQG shown in Fig. 4 .

Eigenvalue $\lambda$	harmonics	Irr. Rep.
3.54	$1$
1.40	$\{\cos \left(\varphi \right),\sin \left(\varphi \right),$	E
1.39	$\cos \left(2\varphi \right),\sin \left(2\varphi \right)\}$
0.37	$\cos \left(3\varphi \right)$	$A_1$
$-0.31$	$\{\cos \left(5\varphi \right),\sin \left(5\varphi \right),$	E
$-0.30$	$\cos \left(7\varphi \right),\sin \left(7\varphi \right)\}$
$-0.30$	$\cos \left(6\varphi \right)$	$A_1$
0.27	$\sin \left(6\varphi \right)$	$A_2$

Equations

$$H_{\mathrm{int}}=\frac{1}{2}\sum _{i,j,\sigma ,{\sigma }^{\prime }}{a}_{i\sigma }^{†}{a}_{i\sigma }{v}_{\sigma {\sigma }^{\prime }}({\mathbf{r}}_i-{\mathbf{r}}_j)a_{j{\sigma }^{\prime }}^{†}{a}_{j{\sigma }^{\prime }},$$

$$G^{-1}=G_0^{-1}-\Sigma$$

$${\left(G\right)}_{i\sigma ,j\sigma }=-\sum _a\frac{1}{{\epsilon }_{a\sigma }}{\varphi }_{a\sigma }({\mathbf{r}}_i){\varphi }_{a\sigma }({\mathbf{r}}_j{)}^{\ast }.$$

$${\left(\Sigma \right)}_{i\sigma ,j\sigma }=$$

$$-v_{\sigma \sigma }({\mathbf{r}}_i-{\mathbf{r}}_j)\text{sideset}{}^{\prime }{\sum }_a{\varphi }_{a\sigma }({\mathbf{r}}_i){\varphi }_{a\sigma }({\mathbf{r}}_j{)}^{\ast },$$

$$h_{ij}^{(\sigma )}=\text{sideset}{}^{\prime }{\sum }_a{\varphi }_{a\sigma }({\mathbf{r}}_i){\varphi }_{a\sigma }({\mathbf{r}}_j{)}^{\ast }.$$

$$P_{\pm }^{(\sigma )}=\mathrm{Im}\left(\sum _{i\in A}{\left(h_{i_1{i}_2}^{(\sigma )}{h}_{i_2{i}_3}^{(\sigma )}{h}_{i_3{i}_1}^{(\sigma )}\right)}^{\frac{1}{3}}\pm \sum _{i\in B}{\left(h_{i_1{i}_2}^{(\sigma )}{h}_{i_2{i}_3}^{(\sigma )}{h}_{i_3{i}_1}^{(\sigma )}\right)}^{\frac{1}{3}}\right)$$

$$C^{(\sigma )}=$$

$$V(\varphi ,{\varphi }^{\prime })=V_0(\varphi ,{\varphi }^{\prime })-$$

$$\frac{1}{(2\pi {)}^2}{\int }^{{\Lambda }_0}\frac{d\epsilon }{\epsilon }{\int }_0^{2\pi }d{\varphi }^{\prime \prime }\frac{\partial k_{\perp }}{\partial \epsilon }\frac{\partial k_{\parallel }}{\partial {\varphi }^{\prime \prime }}{V}_0(\varphi ,{\varphi }^{\prime \prime })V({\varphi }^{\prime \prime },{\varphi }^{\prime })$$

$$\epsilon \frac{\partial \hat{V}(\varphi ,{\varphi }^{\prime })}{\partial \epsilon }=\frac{1}{2\pi }{\int }_0^{2\pi }d{\varphi }^{\prime \prime }\hat{V}(\varphi ,{\varphi }^{\prime \prime })\hat{V}({\varphi }^{\prime \prime },{\varphi }^{\prime })$$

$$V_0(\varphi ,{\varphi }^{\prime })=v_{\mathbf{k}-{\mathbf{k}}^{\prime }}+\frac{v_{\mathbf{Q}}^2{\tilde{\chi }}_{\mathbf{k}+{\mathbf{k}}^{\prime }}}{1-v_{\mathbf{Q}}{\tilde{\chi }}_{\mathbf{k}+{\mathbf{k}}^{\prime }}},$$

$${\epsilon }_c={\Lambda }_0{e}^{-1/|\lambda |}$$

Full text

Introduction.— The origin of superconductivity (SC) in twisted multilayer graphene samplesCao et al. (2018a, b); Yankowitz et al. (2019); Kerelsky et al. (2019) is still highly debated and poses one of the major theoretical challenges of condensed matter physics.Codecido et al. (2019); Shen et al. (2020); Lu et al. (2019); Chen et al. (2019); Xu and Balents (2018); Volovik (2018); Yuan and Fu (2018); Po et al. (2019); Roy and Jurivcić (2019); Guo et al. (2018); Dodaro et al. (2018); Liu et al. (2018); Slagle and Kim (2019); Peltonen et al. (2018); Kennes et al. (2018); Koshino et al. (2018); Kang and Vafek (2018); Isobe et al. (2018); You and Vishwanath ; Wu et al. (2018); Zhang et al. (2019); Ochi et al. (2018); Thomson et al. (2018); Carr et al. (2018); Guinea and Walet (2018); Zou et al. (2018); Laksono et al. (2018); González and Stauber (2019); Kang and Vafek (2019); González and Stauber (2020); Carr et al. (2020); Lopez-Bezanilla and Lado (2020); Cao et al. (2021); Fischer et al. (2022); Liu et al. (2022); Christos et al. ; Scammell et al. Whereas initially also electron-phonon pairing was intensively discussed, it has now become an experimental evidence that the electron-electron interaction is responsible for the unconventional low-temperature behavior of twisted bilayer graphene (TBG) at the magic angle.

In this paper, we establish a common ground to understand the interaction effects in TBG, twisted trilayer graphene (TTG), and twisted quadrilayer graphene (TQG). We elucidate that the SC originates in all of them from a common feature, which is the strong valley symmetry breaking (VSB) characteristic of these moiré systems at the magic angle. Such a strong symmetry breaking pattern has profound consequences and, in TBG as well as for TTG and TQG, it leads to the breakdown of the $C_{6}$ symmetry of the quasi-flat moiré bands down to $C_{3}$, which is the natural symmetry group operating in a single valley.

The breakdown of inversion symmetry gives rise to a strong distortion of the quasi-flat valence bands (VBs) in the Brillouin zone, with a strong modulation of the Coulomb interaction along the anisotropic Fermi line which reflects in the development of some negative (attractive) couplings in the decomposition over the different harmonics. This type of Kohn-Luttinger (KL) pairing instability is relevant until the Fermi line recovers a more isotropic shape towards the bottom of the second VB, which explains the reduction of the critical temperature of SC beyond three-hole doping per moiré unit cell.

The strong VSB can be confirmed from the experimental observations in TBG and TTG. It offers an explanation for the reset of the Hall density around filling fraction $\nu=-2$, which otherwise would be an ordinary filling level in the middle of the lowest-energy VBs. At $\nu=-2$, the strong VSB splits the two valleys, placing the Fermi level at the vertices of the Dirac cones in the energetically lower valley. The Dirac nodes are destabilized by a sufficiently strong Coulomb interaction, opening a gap through a mechanism of dynamical symmetry breaking. This subsequent pattern turns out to be time-reversal symmetry breaking (TRSB) driving into an insulator phase in the case of TTG, thus explaining the experimental observation of a Chern number $C=-2$ in the trilayer at $\nu=-2$. However, the subsequent metallic or insulating character (after TRSB) is less robust in TBG and it may fluctuate depending on the twist angle and the screening environment, making more variable the experimental observations.

The lack of inversion symmetry in the regime of VSB implies also that singlet Cooper pairs can only be formed if electronic states with opposite spin projection are allowed to have opposite sign of the valley polarization, so that the two spin projections are attached to opposite valleys. This type of solution for the two spin sectors arises naturally in the self-consistent resolution of the interacting system in the Hartree-Fock approximation, as we show below.

VSB leads then to spin-valley locking, where the Fermi lines for the two spin projections are different and related by inversion symmetry, as shown in Fig. 1. This leads to a particular realization of the so-called Ising SC, as the spin-valley locking represents a large renormalization of the small spin-orbit coupling of the graphene multilayers, lending protection to the SC against in-plane magnetic fields. This may explain the large in-plane fields required to destroy the SC and the Pauli limit violation observed in TTG, and which would be absent in TBG due to the large orbital Zeeman coupling characteristic of the bilayer system.

Symmetry breaking in Hartree-Fock approximation.— We carry out a real space description of the twisted multilayers, taking into account all the carbon atoms in the supercell of the moiré systems. In this respect, we are going to consider some representatives from the sequence of commensurate superlattices with twist angle $\theta_{i}=\arccos((3i^{2}+3i+1/2)/(3i^{2}+3i+1))$. In the case of TBG we will focus on the system with $i=30$, with a twist angle $\theta=1.08^{\circ}$ which matches well the magic-angle condition. We will recall the results obtained for TTG in Ref. Gonzalez and Stauber, , which considered the configuration with alternating twist angle between layers with $\theta=1.61^{\circ}$, very close to the magic-angle condition for this symmetric TTG. In the case of TQG, we will assemble two twisted bilayers with $i=18$ and alternating twist angle $\theta=1.79^{\circ}$, which place the system very close to the flat-band condition.

The starting point in the real space approach is the tight-binding approximation for the noninteracting system, in which the Hamiltonian $H_{0}$ is taken with a typical parametrization of the exponential decay of hopping matrix elements, as given in Ref. Moon and Koshino, 2013. To this, we add the interacting Hamiltonian $H_{\rm int}$ accounting for the Coulomb interaction, which is expressed in terms of creation (annihilation) operators $a_{i\sigma}^{\dagger}$ ($a_{i\sigma}$) for electrons at each carbon site $i$ with spin $\sigma$ as

[TABLE]

For ${\bm{r}}_{i}\neq{\bm{r}}_{j}$, we take $v_{\sigma\sigma^{\prime}}({\bm{r}}_{i}-{\bm{r}}_{j})=v({\bm{r}}_{i}-{\bm{r}}_{j})$, $v$ being the extended Coulomb potential. This is supposed to be screened by metallic gates above and below the graphene compound, which introduces an exponential decay with screening length $\xi$. Moreover, the strength of the potential is also reduced by a dielectric constant $\epsilon$, accounting mainly for internal screening from electron-hole excitations. For ${\bm{r}}_{i}={\bm{r}}_{j}$, we have the Hubbard interaction $v_{\sigma\sigma^{\prime}}=U\delta_{\sigma,-\sigma^{\prime}}$. The precise value of this coupling is not very important, as long as it is nonvanishing, but it plays an essential role to constrain the relative orientation of the spin projections in the two valleys of the multilayers.

We adopt a Hartree-Fock (HF) approach to the many-body problem, which implies the self-consistent resolution of the Dyson equation involving the interacting (noninteracting) electron propagator $G$ ($G_{0}$) and the self-energy $\Sigma$

[TABLE]

The HF approximation proceeds by assuming that the static limit of the propagator $G$ admits a representation similar to that of $G_{0}$, but with a set of eigenvalues $\varepsilon_{a\sigma}$ and eigenvectors $\phi_{a\sigma}({\bm{r}}_{i})$ modified by the interaction:

[TABLE]

The resolution of (2) is then feasible since the self-energy is given in terms of the eigenvectors by

[TABLE]

where the prime denotes that the sum is only over the occupied levels Fetter and Walecka (1971).

The HF approach is well-suited to our discussion, since the Fock potential is the essential contribution triggering the different symmetry-breaking order parameters. These can be written from the matrix elements

[TABLE]

Thus, the order parameters for TRSB are given by

[TABLE]

where the sums run over the loops made of three nearest neighbors $i_{1},i_{2}$ and $i_{3}$ of each atom $i$ in graphene sublattices $A$ and $B$. A nonvanishing $P_{+}$ corresponds to the usual Chern insulator phase with Haldane mass, while $P_{-}\neq 0$ signals instead an uneven rigid shift in the energies of the two valleys of the moiré system. Moreover, the other pattern likely to be found corresponds to chiral symmetry breaking, measured by the order parameter

[TABLE]

The resolution of Eq. (2) allows us to map the different symmetry breaking order parameters as the strength of the Coulomb interaction evolves with $\epsilon$, from weak to strong coupling. The different phases obtained at filling fraction $\nu=-2$ for TBG and TQG at their respective magic angles are shown in Fig. 2. The phase diagram for TTG with $\theta=1.61^{\circ}$ at the same hole doping can be found in Ref. Gonzalez and Stauber, .

It is remarkable that the picture of symmetry breaking patterns is qualitatively the same in the three multilayers. The critical point for the opening of a gap at the Dirac nodes is shifted towards weaker coupling in Fig. 2(a), but this is due to the fact that TBG at $\theta=1.08^{\circ}$ starts from a more critical condition, with a width of the noninteracting quasi-flat bands smaller than in the other two cases. On the other hand, we observe that another transition takes place to the right of the phase diagrams. The order parameter of VSB vanishes at a certain interaction strength, which is the signal that VSB is giving way to a different phase characterized by intervalley coherence in the strong coupling limit. However, in the absence of valley polarization, this new pattern cannot account for the jumps observed in the Hall density of TBG and TTG at $\nu=-2$, so we will not consider such a strong coupling phase in the subsequent discussion.

The common structure of the phases in the three multilayers is indeed ubiquitous in the moiré systems, as it is also found in TBG at larger twist angles where the magic-angle condition is tuned by applying hydrostatic pressure. These systems with smaller supercells afford an exact HF treatment, including all the valence bands. In the real systems with 5000-10000 atoms in the moiré, we only include up to 128 bands around the neutrality point in the self-consistent resolution, due to technical reasons. Nevertheless, the exact resolutions of the smaller systems show a similar structure of the phases as in Figs. 2(a)-(b) which lends further support to the results obtained here.

We remark that the gap plotted in Figs. 2(a)-(b) corresponds to the separation between valence and conduction bands at the $K$ point in the lower valley, in the presence of valley polarization. We consider that TBG, as well as TTG and TQG, must be placed in the phase with VSB, since the splitting of the conduction and the valence band in the lower valley is the feature that explains the jump in the Hall density observed experimentally at $\nu=-2$. In the phase with VSB, TBG and TQG show indeed at $\nu=-2$ a splitting of the bands at the Fermi level, as can be seen in Fig. 3. It has been shown in Ref. Gonzalez and Stauber, that a similar effect takes place in the case of TTG, and that the splitting of the valence and the conduction band reproduces the reset of the Hall density found experimentally near the integer filling.

Let us finally comment on the behavior of the band-structure of TBG close to the $\Gamma$ point. As can be seen in Fig. 3(a), the eigenenergies are shifted to lower energies in contrary to the expectation that the Coulomb energy should push them upward. The reason is that the Hartree contribution is not fully taken into account as only 128 bands are considered. Nevertheless, it is the exchange energy that is responsible for the symmetry breaking which is treated adequately.

Ising superconductivity.— The strong VSB leads to ground states in which the time-reversal symmetry is broken for a given spin projection, but where this symmetry is recovered upon exchange of the two spin projections. This is indeed realized in the HF approximation as this approach allows for solutions where the states for the two spin sectors have opposite sign of the VSB order parameter. The bands are then symmetric under the transformation ${\bm{k}}\rightarrow-{\bm{k}}$ and concurrent exchange of the spin projections, as shown in Figs. 1 and 4. We observe in the figures that the two spin sectors have different Fermi lines, related by inversion symmetry in momentum space.

The main consequence of the spin-valley locking is that the electrons with opposite momenta of a Cooper pair are forced to live on different Fermi lines attached to opposite valleys. In this regard, this effect represents a strong reinforcement of the small spin-orbit coupling in the graphene layers. A similar renormalization was considered in Ref. Kane and Mele, 2005, taking into account diagrammatic contributions which coincide with those encoded in the Fock part of our approximation. This implies that the SC in the twisted multilayers must be of Ising type, as the spin-valley locking favors a preferred orientation of the spins of a Cooper pair, normal to the graphene planes.

What remains to be discussed is the pairing instability in the system, which demands the analysis of the Cooper pair vertex $V$ for electrons with zero total momentum. This vertex can be parametrized in terms of the angles $\phi$ and $\phi^{\prime}$ of the respective momenta of the spin-up incoming and outgoing electrons on each contour line of energy $\varepsilon$. The instabilities of the vertex can be obtained by solving the equation encoding the iteration of the scattering of Cooper pairs:

[TABLE]

where $k_{\parallel},k_{\perp}$ are the longitudinal and transverse components of the momentum for each energy contour line while $V_{0}(\phi,\phi^{\prime})$ is the bare vertex at a high-energy energy cutoff $\Lambda_{0}$.

Eq. (8) can be simplified by differentiating with respect to the cutoff, which leads to

[TABLE]

with $\widehat{V}(\phi,\phi^{\prime})=F(\phi)F(\phi^{\prime})V(\phi,\phi^{\prime})$ and $F(\phi)=\sqrt{(\partial k_{\perp}/\partial\varepsilon)(\partial k_{\parallel}/\partial\phi)/2\pi}$. Eq. (9) implies that the vertex is a function of the variable $\varepsilon/\Lambda_{0}$. If the initial condition $V_{0}(\phi,\phi^{\prime})$ has a negative eigenvalue for any of its projections (harmonics) over the Fermi line, the solutions of (9) will display then a divergence in the magnitude of the negative coupling as $\varepsilon\rightarrow 0$, which is the signature of the pairing instability.

The crucial question is to start with a sensible representation of the initial vertex $V_{0}(\phi,\phi^{\prime})$ at the high-energy cutoff. For this we rely on the iteration of diagrams in the particle-hole channel shown in Fig. 5, which are not included in the RPA diagrams already accounted for by the internal screening of the Coulomb potential $v$. We have therefore

[TABLE]

where ${\bm{k}},{\bm{k}}^{\prime}$ are the respective momenta for angles $\phi,\phi^{\prime}$ and $\widetilde{\chi}_{{\bm{q}}}$ is particle-hole susceptibility appearing in Fig. 5.

The final step is to project the vertex onto the harmonics $\cos(n\phi),\sin(n\phi)$ which build up the different contributions to $\widehat{V}(\phi,\phi^{\prime})$ at the high-energy cutoff. We have carried out this operation along the Fermi lines of TBG shown in Fig. 1, at filling fraction $\nu=-2.4$. The eigenvalues for the different harmonics can be grouped according to the irreducible representations of the approximate symmetry group $C_{3v}$. The results can be seen in Table 1.

We find that, among the eigenvalues in Table 1, there is a pair of negative couplings with relatively large magnitude $|\lambda|\approx 0.4$. We recall that, for a negative coupling $\lambda$ in a given representation, the solution of Eq. (9) develops a divergence in that channel at the energy scale

[TABLE]

In TBG, the magnitude of $\Lambda_{0}$ is constrained by the reduced bandwidth of the second VB in Fig. 1. We can assign to $\Lambda_{0}$ the value of half the bandwidth, so that $\Lambda_{0}\sim 1$ meV. In this case, the small magnitude of the cutoff is compensated by the relatively large value of $|\lambda|$, which leads to a critical temperature $T_{c}\sim 1$ K.

The same analysis of Fourier decomposition can be done in the case of TQG, with the results shown in Table 2 for the Fermi line in Fig. 4 at $\nu=-2.8$. The value of the dominant negative coupling is somewhat smaller than in TBG, $|\lambda|\approx 0.3$, but the estimate of the critical temperature gives the same order of magnitude $T_{c}\sim 1$ K. In both cases, we observe that the main limitation in $T_{c}$ may come from the reduced space for the scattering of the Cooper pairs, which is inherent to the small bandwidth of the second VB.

Conclusion.— In both TBG and TQG, the large values of the dominant attractive couplings are a consequence of the approximate nesting between parallel segments in the Fermi lines shown in Figs. 1 and 4. In this respect, the Kohn-Luttinger mechanism operates in rather similar way as in TTG, where there is approximate nesting between the sides of the triangular patches of the Fermi line (as shown in Ref. Gonzalez and Stauber, ). In the three cases, the nesting features disappear as the hole doping increases and the Fermi lines recover a more isotropic shape, explaining why the Kohn-Luttinger instability fades away before reaching the bottom of the band.

We have seen that the Fermi lines for the two spin projections are different and related by inversion symmetry in the twisted multilayers. This leads to a particular realization of the so-called Ising SC, as the spin-valley locking represents a large renormalization of the small spin-orbit coupling of the graphene multilayers, lending protection to the SC against in-plane magnetic fields. This may explain the large in-plane fields required to destroy the SC and the Pauli limit violation observed in TTG, and which would be absent in TBG due to the large orbital Zeeman coupling characteristic of the bilayer system.

Acknowledgements. This work was supported by the project PID2020-113164GB-I00 financed by MCIN/ AEI/10.13039/501100011033. The access to computational resources of CESGA (Centro de Supercomputación de Galicia) is also gratefully acknowledged.