Competing phases of interacting electrons on triangular lattices in moir\'e heterostructures

TL;DR

This paper investigates the competing electronic phases in moiré heterostructures on triangular lattices, revealing tendencies towards density waves, Chern insulators, and topological superconductivity through a comprehensive renormalization group analysis.

Contribution

It introduces an unbiased functional renormalization group study of an extended Hubbard model for moiré heterostructures, identifying novel phases and their dependence on filling and interactions.

Findings

Hund-like couplings induce density wave instabilities.

SU(4) exchange promotes a Chern insulator phase.

Topological $d\u00b1 id$ and $f$-wave superconductivity are found.

Abstract

We study the quantum many-body instabilities of interacting electrons with SU(2)$\times$SU(2) symmetry in spin and orbital degrees of freedom on the triangular lattice near van-Hove filling. Our work is motivated by effective models for the flat bands in hexagonal moir\'e heterostructures like twisted bilayer boron nitride and trilayer graphene-boron nitride systems. We consider an extended Hubbard model including onsite Hubbard and Hund's couplings, as well as nearest-neighbor exchange interactions and analyze the different ordering tendencies with the help of an unbiased functional renormalization group approach. We find three classes of instabilities controlled by the filling and bare interactions. For a nested Fermi surface at van-Hove filling, Hund-like couplings induce a weak instability towards spin or orbital density wave phases. An SU(4) exchange interaction moves the system towards a Chern insulator, which is robust with respect to perturbations from Hund-like interactions or deviations from perfect nesting. Further, in an extended range of fillings and interactions, we find topological $d\pm id$ and (spin-singlet)-(orbital-singlet) $f$-wave superconductivity.

Equations

$$H_0=-t(\alpha )\sum _{⟨ij⟩}\sum _{\sigma =\uparrow ,\downarrow }\sum _{o=1,2}{c}_{i\sigma o}^{†}{c}_{j\sigma o}+h.c..$$

$$H_I=\frac{U}{2}\sum _{i,\nu ,{\nu }^{\prime }}{n}_{i\nu }{n}_{i{\nu }^{\prime }}+J\sum _{⟨ij⟩}\sum _{a=1}^{15}{\hat{T}}_i^a{\hat{T}}_j^a,$$

$$H_I^{\prime }=-V_h\sum _i{\mathbf{S}}_i^2-K\sum _i{\mathbf{L}}_i^2,$$

$${\tilde{H}}_I=$$

$$\displaystyle+W(k_{1},k_{2},k_{3})c^{\dagger}_{k_{3}\sigma o}c^{\dagger}_{k_{4}\sigma^{\prime}o^{\prime}}c_{k_{2}\sigma^{\prime}o}c_{k_{1}\sigma o^{\prime}}\big{]}\,,$$

$$H_{\mathrm{QAH}}=\bar{g}\sum _a{b}_a^{†}{b}_a,$$

$$H_{\mathrm{SDW}}^a$$

$$H_{\mathrm{SC}}=-\hat{g}{\mathbf{\Delta }}_d^{†}{\mathbf{\Delta }}_d,$$

$${\partial }_{\Lambda }V(q_1,q_2,q_3)$$

$${\partial }_{\Lambda }W(q_1,q_2,q_3)$$

$$T_{pp}^V(q_1,q_2,q_3)$$

$$T_{pp}^W(q_1,q_2,q_3)$$

$$T_{ph,d}^V(q_1,q_2,q_3)$$

$$+V(q_1,k,q_3)V(q_2,q_{h,d},k)+W(k,q_1,q_3)W(q_{ph,d},q_2,k)+W(q_1,k,q_3)W(q_2,q_{ph,d},k)$$

$$-2V(q_1,k,q_3)W(q_{ph,d},q_2,k)+2V(q_1,k,q_3)W(q_2,q_{ph,d},k)-2W(q_1,k,q_3)V(q_{ph,d},q_2,k)$$

$$\displaystyle+2W(k,q_{1},q_{3})V(q_{ph,d},q_{2},k)\Big{]}$$

$$T_{ph,d}^W(q_1,q_2,q_3)$$

$$\displaystyle+W(q_{1},k,q_{3})V(q_{2},q_{ph,d},k)\Big{]}$$

$$T_{ph,cr}^V(q_1,q_2,q_3)$$

$$T_{ph,cr}^W(q_1,q_2,q_3)$$

$$\displaystyle+W(k,q_{2},q_{3})V(q_{1},q_{ph,cr},k)\Big{]}\,,$$

$$V^{\mathrm{SC}}$$

$$+\frac{{\lambda }^W}{4}{d}^W(k{)}^{\ast }{d}^W(k^{\prime })c_{ka\alpha }^{†}{c}_{-kb\beta }^{†}{c}_{-k^{\prime }c\gamma }{c}_{k^{\prime }d\delta }\left[-(G_0^{†}{)}_{ab}(G_0{)}_{cd}+({\mathbf{G}}^{†}{)}_{ab}(\mathbf{G}{)}_{cd}\right]\otimes \left[({\Gamma }_0^{†}{)}_{\alpha \beta }({\Gamma }_0{)}_{\gamma \delta }+({\mathbf{\Gamma }}^{†}{)}_{\alpha \beta }(\mathbf{\Gamma }{)}_{\gamma \delta }\right].$$

$$H_{\mathrm{ODW}}^a=g_{\mathrm{ODW}}{\mathbf{L}}_{M_a}\cdot {\mathbf{L}}_{-M_a}$$

$$H_{\mathrm{OSDW}}^a=g_{\mathrm{SODW}}{\mathbf{J}}_{M_a}\cdot {\mathbf{J}}_{-M_a}$$

Full text

Correlated insulating and superconducting behavior have recently been discovered in twisted bilayer graphene (TBG) Cao et al. (2018a, b); Yankowitz et al. (2018) triggering ample excitement due to their potential to shed new light on the problem of unconventional superconductivity. More generally, experiments with moiré superlattices of two-dimensional van-der-Waals heterostructures – further including, for example, hexagonal boron nitride (hBN) layers Chen et al. (2018) – are established as an experimental platform for studies of correlated electron physics. These systems allow for a high degree of control, e.g., in the regulation of the twist angle, a low level of disorder and gate-tunable effective bandwidths or filling factors. Indeed, signatures of tunable insulating and superconducting states have been reported in trilayer graphene/hBN heterostructures Chen et al. (2018, 2019). Further, it has been suggested that a wider class of “magic-angle” systems can be realized in present cold-atom setups Fu et al. (2018).

The appearance of strong correlations in these systems is usually ascribed to the emergence of low-lying flat bands and an increased density of states amplifying the impact of electronic interactionsTrambly de Laissardière et al. (2010); Bistritzer and MacDonald (2011); Lopes dos Santos et al. (2012); Li et al. (2010); Brihuega et al. (2012); Kim et al. (2016). However, in conjunction with the experimental findings, the nature of correlated states in moiré superlattices has yet to be identified and described by appropriate models and methods. One way, recently pursued, is the construction of effective multi-orbital tight-binding models for the nearly flat bands based on Wannier states Koshino et al. (2018); Kang and Vafek (2018); Po et al. (2018). While this in itself is a non-trivial task, some universal aspects are shared by all models, including the focus on the emergent superlattice, the presence of several orbitals inherited from the valleys of the original bands and sizable further-neighbor interactions. Hence, a basic understanding of the correlated behavior may be developed by analyzing phenomenological models that capture the qualitative features of moiré flatbands Xu and Balents (2018); Yuan and Fu (2018); Dodaro et al. (2018); Kennes et al. (2018); Tang et al. (2018); Isobe et al. (2018); Lin and Nandkishore (2019); Rademaker and Mellado (2018); Huang et al. (2018); Fidrysiak et al. (2018); Guo et al. (2018); Roy and Juricic (2018); Liu et al. (2018); Da Liao et al. (2019). In addition, this approach allows one to study multi-orbital effects in hexagonal systems in general. This opens the possibility to new types of interactions and correlations. For example, spin-singlet(triplet) pairing is not bound to even(odd)-parity because the anti-symmetry can be guaranteed through the orbital degrees of freedom.

Regarding the role of interactions, no consensus has been reached on whether the experimentally observed insulating states in TBG or trilayer graphene/hBN are due to Mott physics in the strong coupling limit or due to interaction-driven symmetry-breaking, e.g., based on Fermi-surface nesting close to van-Hove singularities Li et al. (2010); Brihuega et al. (2012); Kim et al. (2016); Isobe et al. (2018); Lin and Nandkishore (2019). For the lattice structure of the effective model, triangular and honeycomb geometries have been discussed. While for TBG most recent suggestions favor a honeycomb superlattice Koshino et al. (2018); Kang and Vafek (2018); Po et al. (2018), moiré structures with hBN are argued to support a triangular superlattice, e.g., trilayer graphene/hBN Chen et al. (2018); Po et al. (2018) and twisted bilayer hBN Xian et al. (2018). In these systems, the flat bands are suggested to be well isolated Chittari et al. (2019); Chen et al. (2018); Xian et al. (2018) so that application of an effective model for the flat bands seems appropriate.

In this work, we consider a phenomenological model of electrons on a triangular lattice with two spin and two orbital states as well as onsite and nearest-neighbor interactions. We follow earlier arguments that the combined spin and orbital degrees of freedom form an approximate SU(4) flavor symmetry Xu and Balents (2018); Dodaro et al. (2018); Yuan and Fu (2018); Kennes et al. (2018) and perturb it by Hund-like interactions down to SU(2)$\times$SU(2). We assume that the experimental observations correspond to fillings close to a van-Hove singularity and a (near-)nested Fermi surface and determine potential ordered states induced by moderate interactions. In this situation, electronic correlations in several channels are important, which is why we approach the problem with the help of the fermion functional renormalization group (fRG) method Metzner et al. (2012); Platt et al. (2013).

Several other works have analyzed multi-orbital models on hexagonal lattices close to van-Hove filling with a nested Fermi surface and repulsive interactions. We extend different aspects of these works: in contrast to random phase approximation (RPA) approaches Liu et al. (2018), we take into account the essential coupling between particle-particle and particle-hole channels. While such competition between the different channels is included in parquet RG studies Isobe et al. (2018); Lin and Nandkishore (2019), our calculations are not restricted to a small number of patches at the van-Hove singularity, but resolve the momentum dependence of the entire Fermi surface. This allows us to draw a direct connection to microscopic models and upon inclusion of model parameters to provide quantitative estimates on doping ranges and energy scales. Further, the fRG equations are well defined rooting in an exact flow equation Wetterich (1993) and include not only the leading $\ln^{2}$-diverging channels, but also sub-leading $\ln$-divergent contributions. We extend previous fRG studies Kennes et al. (2018); Tang et al. (2018); Xian et al. (2018) by inclusion of other interaction types and a more comprehensive exploration of the induced ordering tendencies. The exploitation of SU(2)$\times$SU(2) symmetry facilitates high Fermi-surface resolution, which is important near van-Hove singularities.

We find that the interplay between (approximate) nesting, a large density of states and interaction terms results in three classes of strongly-growing correlations: (1) spin/orbital density wave states, (2) $d\pm id$ and $f$-wave superconductivity and (3) a Haldane-like loop-current phase. The density waves are generated right around perfect nesting and involve the different SU(2)$\times$SU(2) (pseudo-)spins. It depends on the SU(4)-symmetry-breaking interactions if a density wave in the separate SU(2) spin or orbital sector, or in the combined SU(2)$\times$SU(2) sector develops. Due to the larger number of fermion flavors and available nesting vectors in the hexagonal geometry, however, the density wave instabilities are fragile. Instead, superconducting and loop-current phases are more robust, see also Ref. Lin and Nandkishore, 2019. Which pairing instability emerges depends on the fermiology: for a closed Fermi surface we find topological $d\pm id$ superconductivityXu and Balents (2018); Lin and Nandkishore (2019), for an open Fermi surface pairing correlations in either $d$- or $f$-channel grow strong. The loop-current state also develops around the perfectly nested situation as soon as small nearest-neighbor SU(4) exchange interactions are present. This quantum-anomalous-Hall state breaks time-reversal symmetry and has a non-zero Chern number together with a fully gapped spectrumVenderbos (2016); Lin and Nandkishore (2019).

Model and parameters. We follow early phenomenological approaches Xu and Balents (2018); Dodaro et al. (2018); Lin and Nandkishore (2019) and model the superlattice of twisted hexagonal multi-layered systems by a two-dimensional triangular lattice populated by electrons with SU(2) spin and orbital degrees of freedom. In particular for twisted hBN, this seems to be an appropriate starting point, where the dispersion close to the Fermi level is well-matched by a two-orbital tight-binding model with nearest-neighbor hopping on the superlattice with hopping amplitude $t(\alpha)$ depending on twist angle $\alpha$ Xian et al. (2018),

[TABLE]

Here, the $c_{i\sigma o}^{(\dagger)}$ are the fermion annihilation (creation) operators on site $i$ with spin projection $\sigma$ and orbital (or band) index $o$. We add a chemical potential $\mu\sum_{i,\nu}n_{i\nu}$ with density $n_{i\nu}=c^{\dagger}_{i\nu}c_{i\nu}$. Combining the indices $\sigma$ and $o$ as $\nu=(\sigma,o)$ into four flavors of fermion states makes an SU(4) symmetry of Eq. (1) apparent. In TBG and trilayer graphene/hBN, the index $o$ corresponds to the two valleys of the untwisted lattice systemKoshino et al. (2018); Kang and Vafek (2018); Po et al. (2018). But the model can also be considered in the more general context of multiorbital effects. We expect spin-SU(2) symmetry to be a good approximation because of the small spin-orbit coupling of the light atoms composing graphene- or boron-nitrite-based systems. Approximate orbital degeneracy is assumed because in twisted graphene heterostructures, mixing of different valleys is suppressed by their large separation in momentum space Koshino et al. (2018). In twisted hBN, ab initio calculations even show an exact degeneracyXian et al. (2018).

To study correlations, we start with SU(4)-invariant Hubbard repulsion $U$ and nearest-neighbor exchange $J$

[TABLE]

with $\hat{T}^{a}_{i}=c_{i\nu}^{\dagger}T^{a}_{\nu\nu^{\prime}}c_{i\nu^{\prime}}$. The $T^{a}_{\nu\nu^{\prime}}$ are 4$\times$4 matrices forming the fundamental representation of the SU(4) Lie algebra. Our motivation to consider the exchange $J$ is twofold. First, it was shown that nearest-neighbor interactions can be sizable despite the large moiré lattice spacing Isobe et al. (2018). Second, the exchange term is introduced in the strong coupling limit of the onsite term. Being limited to weak and intermediate interactions by our method, we can use $J$ as a hint for the effect of stronger couplings.

We account for corrections to the approximate SU(4) symmetry by adding a Hund’s coupling $V_{h}$ and its equivalent in orbital space $K$, which can be induced from integrating out higher-energy phonons in an anti-adiabatic limit due to the small band widths of the flat bands Dodaro et al. (2018),

[TABLE]

where $\bm{S}_{i}=\frac{1}{2}c^{\dagger}_{i\sigma o}\bm{\sigma}_{\sigma\sigma^{\prime}}c_{i\sigma^{\prime}o}$ and $\bm{L}_{i}=\frac{1}{2}c^{\dagger}_{i\sigma o}\bm{\tau}_{oo^{\prime}}c_{i\sigma o^{\prime}}$ with the Pauli matrices $\sigma_{a},\tau_{a},a\in\{1,2,3\}$. Summation over repeated spin and orbital indices is implied here and in the following 111The couplings in Eqs. (2), (3) relate to the common spin-SU(2) symmetric Hubbard and Hund terms $H_{I}=1/2\sum_{\sigma\sigma^{\prime}}[\sum_{o}U_{\mathrm{SU(2)}}n_{io\sigma}n_{io\sigma^{\prime}}+U^{\prime}_{\mathrm{SU(2)}}\sum_{o\neq o^{\prime}}n_{io\sigma}n_{io^{\prime}\sigma^{\prime}}+J_{\mathrm{SU(2)}}\sum_{o\neq o^{\prime}}c^{\dagger}_{io\sigma}c^{\dagger}_{io^{\prime}\sigma^{\prime}}c_{io\sigma^{\prime}}c_{io^{\prime}\sigma^{\prime}}]$ via $U_{\mathrm{SU(2)}}=U+V/2+K/2$, $U^{\prime}_{\mathrm{SU(2)}}=U_{\mathrm{SU(2)}}$ and $J_{\mathrm{SU(2)}}=V-K$. Pair-hopping $J_{\mathrm{SU(2)}}^{\prime}$ and $U_{\mathrm{SU(2)}}\neq U^{\prime}_{\mathrm{SU(2)}}$ would break the SU(2)$\times$SU(2) symmetry. The term $\propto V$ can be absorbed by a shift into the terms $\propto U,K$ Dodaro et al. (2018): $U\rightarrow U-V$, $K\rightarrow K+V$.. Alternatively, we consider nearest-neighbor spin or orbital SU(2) exchange couplings $H_{I}^{\prime\prime}=\sum_{\langle ij\rangle}(J_{s}\bm{S}_{i}\bm{S}_{j}+K_{n}\bm{L}_{i}\bm{L}_{j})\,$ to break the SU(4) symmetry and amplify the effect of nearest-neighbor exchange. Upon inclusion of $H_{I}^{\prime},H_{I}^{\prime\prime}$, a SU(2)$\times$SU(2) invariance remains, which originates from the SU(2) spin-rotational invariance and the two-orbital structure of the Hamiltonian. Below, we assume that $V_{h},K,J_{s}$ and $K_{n}$ are small compared to $U,J$, but we will comment on larger corrections to SU(4) in the discussion.

Method. We employ a functional renormalization group (fRG) approach for the one-particle-irreducible interaction vertices of the fermionic many-body system in the Fermi-surface patching scheme Metzner et al. (2012); Platt et al. (2013). Here, we specifically use that any SU(2)$\times$SU(2)-invariant interaction can be decomposed into two parts

[TABLE]

where the first interaction vertex $\propto V$ is SU(4) invariant and the second one $\propto W$ breaks it down to SU(2)$\times$SU(2). The $k_{i},\,i\in\{1,2,3,4\}$ are wavevectors in the first Brillouin zone (BZ) and $k_{4}$ is fixed by wavevector conservation. The fRG approach introduces an infrared cutoff $\Lambda$ and determines the renormalization of the system with respect to $\Lambda$. We derive fRG flow equations for the scale-dependent vertices $V$ and $W$ in the SU(2)$\times$SU(2)-invariant case, see App. A. The fRG flow is initialized at the bandwidth of the system and integrated out toward the Fermi level upon lowering $\Lambda$. This scheme amounts to an infinite-order summation of coupled particle-particle and particle-hole channels of second order in the effective interactions. It enables an unbiased investigation of the competing correlations through the analysis of the components of $V(k_{1},k_{2};k_{3})$ and $W(k_{1},k_{2};k_{3})$, which signal Fermi-liquid instabilities by flowing to large values at a critical fRG scale $\Lambda_{c}$. We use $\Lambda_{c}$ as an estimate for transition temperatures $\Lambda_{c}=T^{*}$. In a RPA-like summation, $T^{*}$ would equal the mean-field transition temperature, whereas it can be altered in our case due to the mutual feedback of the different channels. For the numerical computation, we discretize the wavevector dependence of $V$ and $W$. To this end, we project all wavevectors onto the Fermi surface and resolve the angular dependence by dividing the BZ into $N$ patches as shown in Fig. 1, i.e. $V(k_{1},k_{2},k_{3})\rightarrow V(\varphi_{1},\varphi_{2},\varphi_{3})$ with $|k_{i}|=k_{F}(\varphi_{i})$ being the corresponding Fermi vector. In the following, we discuss our findings from varying $\mu$ and $J$, whilst keeping fixed $U,V_{h}$ and $K$ (or $J_{s}$ and $K_{n}$). Phase diagrams suggested by the diverging correlations are shown in Figs. 1, 2 and 4.

*Instabilities of the nested Fermi surface. * In our model, a chemical potential of $\mu=2t$ (three electrons per site) leads to a perfectly nested Fermi surface with three inequivalent nesting vectors $M_{a}$, $a\in\{1,2,3\}$ (see Fig. 1). As a result, we find both real and imaginary density wave instabilities near van-Hove filling depending on the initial interactions. We assume that the onsite interaction $U$ is always dominating.

The SU(4) exchange coupling $J$ induces an ordering tendency in the imaginary charge density wave channel, which corresponds to a quantum anomalous Hall (QAH) state. It robustly appears in the tentative fRG phase diagram for varying initial couplings and persists when doping slightly away from van-Hove filling, $\mu\neq 2t$, cf. Fig. 1. In particular, it is stable upon inclusion of small SU(4)-breaking terms. The characteristic vertex structure of the QAH instability exhibits a special sign modulation, cf. Fig. 2. The divergence only occurs in vertex $V$ for scattering processes connected by the different $M_{a}$ vectors. At van-Hove filling, the vertex can be fitted by $V(k_{1},k_{2},k_{3})=\bar{g}\,d_{a}(k_{1})d_{a}(k_{2})\delta(k_{3}-k_{1}-M_{a})$ with three $d$-wave functions for each of the $M$ points $d_{1}=2\sin(k_{x}/2)\sin(\sqrt{3}k_{y}/2),d_{2,3}=\mp\cos(k_{x})\pm\cos((k_{x}\pm\sqrt{3}k_{y})/2)$. This leads to the effective Hamiltonian

[TABLE]

where $b_{a}=\sum_{k}d_{a}(k)c_{k+M_{a}\sigma o}^{\dagger}c_{k\sigma o}$ and $\bar{g}>0$. Away from $\mu=2t$, the QAH divergencies are more limited to the vicinity of the $M$ points. $H_{\mathrm{QAH}}$ suggests a purely imaginary mean-field configuration with $\langle b_{a}\rangle=-\langle b_{a}^{\dagger}\rangle$, leading to loop currents in real space Venderbos (2016); Lin and Nandkishore (2019). We find simultaneous ordering tendencies of all three components $b_{a}$. This is expected to be energetically favored as it leads to a fully gapped state, which breaks time-reversal symmetry and has a non-zero Chern invariantVenderbos (2016). We show the corresponding flux configuration in Fig. 2 (b).

For spin or orbital Hund couplings $V_{h},K\neq 0$ larger than the exchange $J$, strong correlations of real SU(2)$\times$SU(2) spin/orbital-density waves (SODW) are induced. This is signaled by sharp horizontal and vertical lines in the vertices as illustrated in Fig. 2 (c). It depends on the SU(4)-breaking terms (3), which spin and orbital correlations are chosen from the allowed SU(2)$\times$SU(2) possibilities: $\bm{S}\cdot\bm{S}$, $\bm{L}\cdot\bm{L}$ or $(c^{\dagger}\bm{\sigma}\otimes\bm{\tau}c)\cdot(c^{\dagger}\bm{\sigma}\otimes\bm{\tau}c)$. For the example in Fig. 2 (b), we choose $U=3t,V_{h}=0.4t$ and $K=0$, which favors $\bm{S}\cdot\bm{S}$. We show other examples in the appendix. The divergent lines in Fig. 2 again correspond to scattering processes that involve points on the Fermi surface connected by one of the three nesting vectors $M_{a}$, cf. Fig. 1. We read off the strength of the vertices as $V(k_{1},k_{2},k_{3})=g\,\delta(k_{3}-k_{1}-M_{a})$ and $W(k_{1},k_{2},k_{3})=2g\,\delta(k_{3}-k_{2}-M_{a})$ with $g>0$. This tensor structure corresponds to the effective Hamiltonian

[TABLE]

with $\bm{S}_{M_{a}}=\sum_{k}c^{\dagger}_{k+M_{a}o\sigma}\bm{\sigma}_{\sigma\sigma^{\prime}}c_{ko\sigma^{\prime}}$. The instability involves all three nesting vectors $M_{a}$ depending on which of them connects the respective part of the Fermi surface. We interpret this as the itinerant triple-$M$ state of Ref. Nandkishore et al., 2012a for the triangular lattice which is conjectured to occur before an insulating chiral SDW sets in at lower temperatures Li (2012); Wang et al. (2012); Martin and Batista (2008). We sketch the real-space spin pattern in the inset of Fig. 2 (c). In contrast to the QAH, the SDW is relatively weak and quickly suppressed by tuning $\mu$ away from perfect nesting or increasing $J$. We attribute this to the increased available phase space due to the three inequivalent nesting vectors. Moreover, the higher flavor number favors the QAH, cf. Ref. Lin and Nandkishore, 2019.

Unconventional superconductivity. When we tune the chemical potential further away from perfect nesting, we find that the system is susceptible to superconductivity (SC) in the $d$-wave channel. In addition, we see strong pairing correlations at $\mu=2t$ for small onsite $U$ or large exchange $J$. Interestingly, the orbital degrees of freedom open the possibility to $d$-wave pairing with (spin-triplet)-(orbital-singlet) or (spin-singlet)-(orbital-triplet) symmetry. In our calculation it depends on the size of Hund’s and orbital-Hund’s couplings which of them is selected. We find that for $V_{h}>K$, superconductivity is mediated by fluctuations of the effective spin Hamiltonian Eq. (6). Mean-field decoupling of Eq. (6) in the pairing channel with even parity leads to attraction in the (spin-singlet)-(orbital-triplet) channel. This is different to the strong-coupling scenario, where the Hund’s coupling favors (spin-triplet)-(orbital-singlet) pairing. The reason is that in this case, superconductivity is mediated by purely ferromagnetic spin-fluctuations, which is qualitatively different from the mediation in terms of the weak-coupling spin-density wave as found here. Note, however, that one still can get (spin-triplet)-(orbital-singlet) $d$-wave SC from weak coupling, e.g., for $K>V_{h}>0$, an orbital density wave is induced for small $J$ mediating (spin-triplet)-(orbital-singlet) $d$-wave pairing.

A snapshot of the SC vertex is shown in Fig. 3 for the example of $\mu=1.9t$, $U=3t$, $J=0.4t$ and $V_{h}=0.1t$. It exhibits a dominant diagonal feature with $d$-wave form factor $d(k)$, i.e. $V(k_{1},k_{2},k_{3})=-\hat{g}\,d(k_{1})d(k_{3})\delta(k_{1}+k_{2})$ where $\hat{g}>0$. As we explain in App. B, the relative sign between vertices $V$ and $W$ gives information about the spin and orbital pairing configuration. In the above example, the divergent structure in $W(\varphi,\varphi^{\prime},\varphi_{0})$ (not shown) is weaker and has the same sign as compared to $V(\varphi,\varphi^{\prime},\varphi_{0})$. Therefore, the leading SC instability corresponds to (spin-singlet)-(orbital-triplet) pairing as expected. We extract the effective Hamiltonian

[TABLE]

with $\bm{\Delta}_{d}=\sum_{k}d(k)c_{-k\sigma o}(i\bm{\tau}\tau^{y})_{oo^{\prime}}(i\sigma^{y})_{\sigma\sigma^{\prime}}c_{k\sigma^{\prime}o^{\prime}}$. From our calculation, the form factor $d(k)$ is predicted to be a superposition of the two $d$-wave functions $d_{xy}=-2\sqrt{3}\sin k_{x}/2\cos\sqrt{2}k_{y}/2$ and $d_{x^{2}-y^{2}}=2\cos k_{x}-2\cos k_{x}/2\cos\sqrt{3}k_{y}/2$ (see Fig. 3). Eventually, we expect the $d_{x^{2}-y^{2}}\pm id_{xy}$ states of Ref. Xu and Balents, 2018 – or rather their spin-singlet, orbital-triplet equivalent – to be favored because this maximizes the pairing gap on the Fermi surface Black-Schaffer and Honerkamp (2014) (see also Refs. Honerkamp, 2003; Nandkishore et al., 2012b). As outlined in Ref. Xu and Balents, 2018, this leaves two types of $d\pm id$ superconductors degenerate at the mean-field level, which are both topological with protected edge states. Note that nematic superconductivity has also been discussed as an alternativeVenderbos and Fernandes (2018).

Regarding the interplay with the QAH, we observe two different scenarios. In the case of a closed Fermi surface, $\mu<2t$, there appears to be a competition between both phases. This is suggested by the behavior shown in Fig. 1 which exhibits a downturn of the superconducting instability temperature just before the steep increase very close to $\mu=2t$ when the QAH tendencies grow strong. In contrast, for $\mu\geq 2t$, the diverging structure of the vertex suggest a crossover from a QAH- to SC-dominated regime when increasing $\mu$. In the intermediate regime, correlations in both channels grow large and one has to go beyond our approach to decide whether or not there will be a phase of coexistence.

We also find the possibility of $f$-wave SC with (spin-singlet)-(orbital-singlet) configuration. This pairing instability occurs for an open free Fermi surface, $\mu>2t$, and dominant $J_{s}>0$ or dominant $K_{n}>0$, see Fig. 4. The unconventional combination of an odd form factor with spin-singlet is due to mediation by spin fluctuations: mean-field decoupling of Eq. (6) in the odd parity channels leads to attraction in the (spin-singlet)-(orbital-singlet) channel. Our vertex data is well fitted by $f_{x(x^{2}-3y^{2})}=\sin k_{x}-2\sin k_{x}/2\cos\sqrt{3}k_{y}/2$ (see Fig. 5) with nodes along the $\Gamma-M_{a}$ directions. Thus, they do not coincide with the Fermi lines which are centered around the $K,K^{\prime}$ points. Since the gap function follows the form factor, this suggests a nodeless $f$-wave SC state.

Conclusion. In this work, we have studied the quantum many-body instabilities of a multi-orbital model for interacting electrons on the triangular lattice. We have incorporated various universal aspects which were suggested to be relevant for the correlated behavior in moiré heterostructures, such as twisted bilayer boron nitride and trilayer graphene-boron nitride systems. We used the $N$-patch functional renormalization group approach for a vertex structure with SU(2)$\times$SU(2) symmetry. This formulation of the fRG is a suitable method to investigate the interplay between (approximate) nesting and a large density of states in an unbiased way because it takes into account all correlation channels and their mutual couplings. We have identified three classes of strongly-growing correlations: spin/orbital density wave states, two types of unconventional superconductivity and a Haldane-like loop-current phase with non-zero Chern number.

While the density waves are fragile and found only very close to perfect nesting, the superconducting and loop-current phases are robust, and occupy an extended range in parameter space. The type of pairing that emerges depends on the Fermi surface topology. We provided further evidence for the previously discussed topological $d\pm id$ superconductivity scenario Xu and Balents (2018); Lin and Nandkishore (2019). On the other hand, for an open Fermi surface the pairing correlations can be of $d$- or $f$-wave form in our computation, depending on the ratio and size of SU(4)-breaking Hund’s or spin/orbital exchange couplings. For larger corrections to the SU(4) invariance we find other types of instabilities. For example, slightly away from van Hove filling for $V_{h}>t$ and the other interactions set to zero, an (orbital-singlet)-(spin-triplet) SC instability with $s$-wave form factors emerges. For negative $J_{s},K_{n}$, the flow is towards a ferromagnetic spin or ferro-orbital order.

As a perspective for future studies we note that the fRG approach has the advantage that all modes of the effective model are taken into account. Hence, in principle, the resulting instability scales $T^{*}$ can serve as estimates for actual (short-range) ordering scales in the experimental system if the model parameters are known. At present, there are only rough estimates for the interaction parameters in TBGKoshino et al. (2018) and twisted bilayer boron nitrideXian et al. (2018) that point to potentially very large onsite repulsions $\sim 30\,\mathrm{meV}$ compared to effective hoppings $t\sim 0.33\,\mathrm{meV}$Koshino et al. (2018). The effect of these large onsite terms might, however, be reduced by the competition with nonlocal interactions Schüler et al. (2013), which are also expected to be sizeable. A study of this interplay is beyond the scope of our paper. Here we only state that the mentioned $t$-value, our critical scales for $U=3t$ end up in the sub-Kelvin range, $T^{*}\sim 0.03t\approx 0.1$K. Note, however, that these scales depend exponentially on the chosen couplings such that energy scales of a few Kelvin can be reached already for moderate interaction strengths.

As the discussion about the role of interactions is still ongoing, we think that the weak-to-intermediate coupling perspective given here, together with the use of the universal aspects of moiré heterostructures can provide essential insights into the general ordering tendencies in such multi-orbital systems. In particular, our prediction of the appearance of robust topological/chiral interaction-induced phases is accessible to experimental examination.

Acknowledgments.

We thank Alexei Tsvelik for discussions. M.M.S. was supported by the DFG through the Collaborative Research Center SFB1238, TP C04. L.C. acknowledges funding from the Alexander-von-Humboldt foundation. Work at BNL is supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Contract No. DE-SC0012704. C.H. acknowledges funding via the DFG research training group 1995 Quantum Many-Body Methods in Condensed Matter Systems.

Appendix A Fermion functional RG setup

We study many-body instabilities of the model defined in Eqs. (1), (2) and (3) by means of the fRG method Wetterich (1993). To that end, we employ its specific formulation for correlated fermion systems Honerkamp et al. (2001); Salmhofer and Honerkamp (2001), see Refs. Metzner et al., 2012; Platt et al., 2013 for reviews. The fRG-$N$-patch scheme starts with a fermionic action corresponding to our model Hamiltonian $S[\bar{\psi},\psi]=-(\bar{\psi},G_{0}^{-1}\psi)+\mathcal{V}[\bar{\psi},\psi]\,.$ The quadratic first term with the free propagator $G_{0}(\omega_{n},k)=1/(i\omega_{n}-\epsilon(k))$ includes Matsubara frequencies $\omega_{n}$ and wavevectors $k$. The dispersion $\epsilon(k)$ is measured relative to the chemical potential. In the present case, the free propagator is diagonal with respect to spin and the orbital quantum numbers, $\sigma$ and $o$, and we have suppressed the according indices. The general interaction contribution $\mathcal{V}[\bar{\psi},\psi]$ is quartic in the fermion fields. Its specific form can be directly inferred from the Hamiltonians in Eqs. (2) and (3).

The bare propagator is then regularized by an infrared cutoff with energy scale $\Lambda$: $G_{0}(\omega_{n},k)\to G_{0}^{\Lambda}(\omega_{n},k)=\theta_{\varepsilon}^{\Lambda}(\epsilon(k))/(i\omega_{n}-\epsilon(k))\,,$ where $\theta_{\varepsilon}^{\Lambda}$ is an approximate step function with smoothening scale $\varepsilon$ cutting off fluctuations with energies $|\epsilon(k)|\lesssim\Lambda$. The regularized propagator $G_{0}^{\Lambda}$ is then used to set up the functional integral for the scale-dependent effective action $\Gamma^{\Lambda}$, which generates the one-particle irreducible vertex functions $\Gamma^{(2n)\Lambda}$.

The fRG flow is generated upon variation of $\Lambda$ producing a hierarchy of differential equations for the vertex functions $\Gamma^{(2n)\Lambda}$. Integration of the flow towards the infrared $\Lambda\rightarrow 0$ yields the full effective action $\Gamma$. We employ the standard truncation for analyzing many-body instabilities in two-dimensional fermion systems, where the fRG flow of all higher $n$-point functions with $n\geq 6$ and self-energy feedback are neglected. In our SU(2)$\times$SU(2)-invariant system, this corresponds to following the scale dependence of the two generally frequency- and momentum-dependent effective interaction vertices $V$ and $W$ as defined in Eq. (4). Their flow consists of three contributions, i.e. the particle-particle, direct and crossed particle-hole channels:

[TABLE]

The expression for the different channels are

[TABLE]

where we defined $q_{pp}=q_{1}+q_{2}-k$, $q_{ph,d}=q_{1}-q_{3}+k$, $q_{ph,cr}=q_{2}-q_{3}+k$ and $L_{\pm}(q,k)=S^{\Lambda}(k)G_{0}^{\Lambda}(q\pm k)+G_{0}^{\Lambda}(k)S^{\Lambda}(q\pm k)$ with the single-scale propagator $S^{\Lambda}=\partial_{\Lambda}G_{0}^{\Lambda}$. In Eqs. (8)-(15), the arguments are combined frequency and momentum vectors $q=(\omega,\bm{q})$. The most singular part of the vertices is expected to come from zero Matsubara frequency. Thus, for the numerical evaluation of instabilities, we do not resolve the frequency dependence of the two-particle vertices and only consider the zero-frequency limit. The wavevector dependence is taken into account via a patching of the Fermi surface which resolves the angular direction. This truncation has been shown to successfully describe the Fermi surface instabilities of numerous systems. In particular, it goes beyond the random phase approximation because the coupling of the different channels is taken into account. We discretize the Fermi surface by patch points, each representing one of $N$ patches covering the BZ, as shown in Fig. 1. We use up to $N=96$ to check the convergence of our results, which is particularly important at van-Hove filling. For $\mu>2t$, the Fermi surface is open and we choose patches that are centered around the $K,K^{\prime}$ points. For details on the patching, we refer to Ref. Hesselmann et al., 2018 which depicts both geometries used here.

We initialize our numerics with an RG scale equal to the bandwidth. An instability towards an ordered state is signaled by a divergence in the two-particle vertices $V,W$ during the flow towards the infrared and we stop the RG evolution when this occurs. In practice, our stopping condition is that one vertex component exceeds $40t$. When $V,W$ remain finite, we stop the flow when $\Lambda<10^{-9}t$.

Appendix B Extraction of effective Hamiltonians

The snapshots of the vertices in Figs. 2 (b)-(d) show that the fRG predicts specific momentum dependencies for the many-body instabilities, which allows us to characterize the occuring ordering tendencies in some detail. Therefore, we systematically extract Scherer et al. (2014) the effective Hamiltonians, cf. Eqs. (5), (6), (7).

In particular, for the SC instability, the procedure is as follows. The diverging part of the vertex suggests $V(k_{1},k_{2},k_{3})\rightarrow V(k,-k,k^{\prime})$. We compute the eigensystem of this effective $N\times N$ matrix as defined on the patch points. The eigenvector with the largest absolute value is expected to provide the order parameter with the highest transition temperature in a meanfield approach and we therefore examine its momentum dependence in more detail. In the case of the $d$-wave SC instability, we find two eigenvectors with degenerate largest eigenvalue. Introducing the momentum-space representations of the $d$-wave form factors on the triangular lattice, $d_{x^{2}-y^{2}}=2\cos(k_{x})-2\cos(k_{x}/2)\cos(\sqrt{3}k_{y}/2)$ and $d_{xy}=-2\sqrt{3}\sin(k_{x}/2)\sin(\sqrt{3}k_{y})/2)$, we can fit the extracted momentum dependencies very well by $d_{+}(k)=-\sin(\phi)d_{x^{2}-y^{2}}+\cos(\phi)d_{xy}$ and $d_{-}(k)=\cos(\phi)d_{x^{2}-y^{2}}+\sin(\phi)d_{xy}$ with $\phi=5\pi/4$ as shown in the Fig. 3. We also show the $f$-wave form factor extracted from the vertices in Fig 4 in Fig. 5.

We now only keep these leading channels, i.e. $V(k,-k,k^{\prime})\rightarrow\lambda^{V}d(k)^{*}d(k^{\prime})$ with corresponding form factor superposition $d(k)$. To connect to the standard notation, we use the Fierz identity $\delta_{ad}\delta_{bc}=(\bm{\sigma}_{ab}\bm{\sigma}_{cd}+\delta_{ab}\delta_{cd})/2=[(\Gamma_{0}^{\dagger})_{ab}(\Gamma_{0})_{cd}+(\bm{\Gamma}^{\dagger})_{ab}(\bm{\Gamma}^{\dagger})_{cd}]/2$ with $\Gamma_{0}=i\sigma_{2}$ and $\Gamma_{i}=\sigma_{i}i\sigma_{2}$. Accordingly, we define the matrices $G_{0}=i\tau_{2}$ and $\bm{G}=\bm{\tau}i\tau_{2}$ in orbital space. We rewrite the interaction

[TABLE]

We used Greek symbols for spin and Latin symbols for orbital space. The minus sign in the expression from the $W$ vertex is due to a transpose that one has to take because of the exchanged orbital index. Now, we see if $d^{i}(k),\ i\in\{V,W\}$ is even, only $G_{0}^{\dagger}G_{0}\otimes\bm{\Gamma}^{\dagger}\bm{\Gamma}$ or $\bm{G}^{\dagger}\bm{G}\otimes\Gamma_{0}^{\dagger}\Gamma_{0}$ are allowed by the antisymmetry of the pairing function. This correspond to (orbital-singlet)-(spin-triplet) and (orbital-triplet)-(spin-singlet) SC, respectively. Due to the minus sign in the expression, we find for $\lambda^{V}=\lambda^{W}$ (orbital-triplet)-(spin-singlet) and for $\lambda^{V}=-\lambda^{W}$ (orbital-singlet)-(spin-triplet).

For the QAH instability we proceed accordingly and find that the form factor $d_{a}$ is determined by the momentum transfer $M_{a}$ as given in the main text. We plot $d_{3}(k)$ as an example, see top panel in Fig. 5.

Appendix C Other instabilities

As we mentioned in the main text, we find an orbital-singlet, spin-triplet $s$-wave instability for small onsite repulsion and dominant Hund’s coupling. The reason is that there is an attractive onsite interaction for this combination. We show the corresponding vertices in Fig. 6 (a).

We also note that besides the spin-density wave discussed in the main text, the system can flow towards an orbital density wave or a combined spin-orbital density wave at $\mu=2t$. It depends on the Hund-like $V,K$ and SU(2) exchange $J_{s},K_{n}$ couplings which density-wave instability of the allowed SU(2)$\times$SU(2) combinations is chosen. We show examples for the vertex of an orbital density wave in Fig. 6 (b) with initial couplings $U=3t$, $K=0.4t$. The leading vertex structure corresponds to

[TABLE]

with $\bm{L}_{M_{a}}=\sum_{k}c^{\dagger}_{k+M_{a}o\sigma}\bm{\tau}_{oo^{\prime}}c_{ko^{\prime}\sigma}$ and $g_{\mathrm{ODW}}>0$. For $U=3t$ and $V=K=0$, but $J_{s},K_{n}\neq 0$, the density-wave instability corresponds to a combined spin-orbital density wave. That is, we find $V(k_{1},k_{2},k_{3})=g_{\mathrm{SODW}}[4\delta(k_{3}-k_{2}-M_{a})-\delta(k_{3}-k_{1}-M_{a})]$ and $W(k_{1},k_{2},k_{3})=2g_{\mathrm{SODW}}[-\delta(k_{3}-k_{2}-M_{a})+\delta(k_{3}-k_{1}-M_{a})]$ with $g_{\mathrm{SODW}}>0$ and the effective Hamiltonian becomes

[TABLE]

with $\bm{J}_{M_{a}}=\sum_{k}c^{\dagger}_{k+M_{a}\sigma o}\bm{\sigma}_{\sigma\sigma^{\prime}}\bm{\tau}_{oo^{\prime}}c_{k\sigma^{\prime}o^{\prime}}$.