# Nematic, chiral and topological superconductivity in transition metal
dichalcogenides

Constantin Schrade  Center for Quantum Devices, Niels Bohr Institute,
University of Copenhagen, 2100 Copenhagen, Denmark  Liang Fu  Department of
Physics, Massachusetts Institute of Technology, 77 Massachusetts Ave.,
Cambridge, MA 02139

###### Abstract

We introduce and study a realistic model for superconductivity in twisted
bilayer WSe2, where electron pairing arises from spin-valley fluctuations in
the weak-coupling regime. Our model comprises both the full continuum model
moiré bandstructure and a short-ranged repulsive interaction that accounts for
the Coulomb interaction projected onto the localized Wannier orbitals. By
calculating the spin-valley susceptibility, we identify a Fermi surface
nesting feature near half-filling of the top-most moiré band, which induces
significantly enhanced spin-valley fluctuations. We then analyze the dominant
Kohn-Luttinger pairing instabilities due to these spin-valley fluctuations and
show that the leading instability can induce nematic, chiral and topological
superconductivity. As our findings are asymptotically exact for small
interaction strengths, they provide a viable starting point for future studies
of superconductivity in twisted transition metal dichalcogenide bilayers.

###### pacs:

68.65.Cd; 68.65.Ac; 71.10.Fd, 74.20.-z

Two-dimensional transition metal dichalcogenides (TMDs) have in recent years
emerged as a promising material platform for realizing a plethora of new
electronic phases bib:Manzeli2017 . Prominent examples of these electronic
phases in monolayer TMDs include the quantum spin Hall effect bib:Qian2014 ;
bib:Tang2017 ; bib:SWu2018 ; bib:Shi2019 , arising due to the strongly spin-
orbit coupled band structure, and Ising superconductivity bib:Lu2015 ;
bib:Xi2016 ; bib:Sohn2018 ; bib:Sergio2018 , which is realized as a result of
the substantial effective electron mass that enhances the importance of
interaction effects.

Beyond the TMD monolayers, another research frontier that is currently
evolving at rapid pace are moiré lattices realized with multilayer TMDs
bib:LeRoy2020 ; bib:Regan2020 ; bib:Yang2020 ; bib:Jin2021 ; bib:Ghiotto2021 ;
bib:Wang2020 ; bib:Li2021 ; bib:Li2021_2 ; bib:Millis2021 . In these systems,
a lattice mismatch or a rotational misalignment generates an effective
superlattice that can be utilized for simulating strongly-correlated electron
states bib:Wu2018 ; bib:Wu2019 ; bib:Schrade2019 ; bib:Pan2020 ; bib:Pan2020_2
; bib:Zhang2020 ; bib:Zhang2020_2 ; bib:Devakul2021 ; bib:Zhang2021 ;
bib:Hsu2021 . For example, in WSe2/WS2 heterobilayers, experiments have
identified Mott insulators bib:Regan2020 ; bib:Yang2020 , Wigner crystals
bib:Regan2020 , and stripe-ordered states bib:Jin2021 . Moreover, in twisted
WSe2 (tWSe2) homobilayers, experiments have found evidence for quantum
criticality bib:Ghiotto2021  and correlated insulating states bib:Wang2020 .
Notably, the discovered correlated insulators in tWSe2 appeared at half-
filling for broad twist angle range, $\theta\sim 4^{\circ}$-$5.1^{\circ}$, for
which the bandwidth is likely comparable to the interaction strength. This
weak- to intermediate coupling scenario is further supported by a resistivity
enhancement upon tuning a van-Hove singularity to half-filling via an external
displacement field. Both features highlight the importance of bandstructure
effects for understanding correlated electron states in tWSe2.

In addition to the discovery of correlated insulators, experiments on tWSe2
also found a zero-resistance state when doping away from half-filling at
$\theta\sim 5.1^{\circ}$ bib:Wang2020 . The appearance of this zero-resistance
state indicates that, similar to graphene-based superlattices bib:Cao20182 ;
bib:Yankowitz2019 ; bib:Chen20192 ; bib:Liu2020 ; bib:Park2021 ; bib:Hao2021 ;
bib:Cao2021 , tWSe2 can also host superconductivity. However, it is so far
unclear what properties of the bandstructure can effect the candidate
superconducting state in tWSe2. Moreover, it is also an open question on what
constitutes the pairing mechanism and pairing symmetry.

Figure 1: (Color online) (a) WSe2 layers at twist angle $\theta$ forming a
triangular moiré lattice. (b) Bandstructure along high-symmetry lines for
$(\theta,V_{z},V,\psi,w)=(5.1^{\circ},0\,\text{meV},9\,\text{meV},128^{\circ},18\,\text{meV})$.
Inset shows the moiré Brillouin zone at the $+K$ valleys. (c) Same as (b) but
for $V_{z}=42\,\text{meV}$. The spin-up (blue) and spin-down (red) bands at
$\pm K$ valleys are split by the finite displacement field. (d) DOS with a
maximum near $n_{\text{max}}\approx 0.9$.

Here, we address these questions by introducing a model for superconductivity
in tWSe2, where electron pairing arises from spin-valley fluctuations. Our
model comprises both the moiré bandstructure of tWSe2 and a short-ranged
repulsive interaction, obtained by projecting the Coulomb interaction onto the
localized moiré Wannier orbitals. Near half-filling of the top-most moiré
band, we demonstrate that the Fermi surfaces of the $\pm K$ valleys exhibit a
strong nesting feature. This nesting maps the occupied states of one valley
onto the other valley’s unoccupied states, leading to strong spin-valley
fluctuations. We show that, interestingly, electron pairing mediated by such
enhanced spin-valley fluctuations can give rise to nematic, chiral and
topological superconductivity. Our findings are asymptotically exact in the
weak coupling limit bib:Raghu2010  and provide a starting point for studies of
superconductivity in TMD systems.

Model. We consider two layers of WSe2 that are initially aligned in the AA-
stacking configuration and, subsequently, rotated by a small twist angle
$\theta$, see Fig. 1(a). In this situation, the $\pm K$ valleys of the layers
have comparable energy so that inter-layer tunneling becomes effective, which
leads to a layer-hybridized moiré bandstructure. For the discussion of this
moiré bandstructure, we focus separately on the spin-up (spin-down) valence
bands of the $+K$ ($-K$) valleys. This approach is justified by the
substantial valley-dependent spin splitting in the valence bands and the
effective decoupling of the valleys due to their large momentum space
separation. The Hamiltonian for the continuum model moiré bands of spin-up
electrons at the $+K$ valleys then reads bib:Wu2019 ,

$H_{\uparrow}=\int\mathrm{d}{{\boldsymbol{r}}}\,\Psi_{\uparrow}^{{\dagger}}({{\boldsymbol{r}}})\begin{pmatrix}h^{t}_{{{\boldsymbol{k}}}\uparrow}({{\boldsymbol{r}}})&\Delta_{T}({{\boldsymbol{r}}})\\\
\Delta^{{\dagger}}_{T}({{\boldsymbol{r}}})&h^{b}_{{{\boldsymbol{k}}}\uparrow}({{\boldsymbol{r}}})\end{pmatrix}\Psi_{\uparrow}({{\boldsymbol{r}}}),$
(1)

where
$\Psi_{\uparrow}({{\boldsymbol{r}}})=(c^{t}_{\uparrow}({{\boldsymbol{r}}}),c^{b}_{\uparrow}({{\boldsymbol{r}}}))^{T}$
is an electron spinor with components corresponding to the top/bottom $(t/b)$
layer. The Hamiltonian for spin-down electrons at the $-K$ valleys is obtained
from $H_{\uparrow}$ via a time-reversal operation.

First, we discuss the diagonal components of $H_{\uparrow}$, which represent
the individual layers. They are given by,

$h^{t/b}_{{{\boldsymbol{k}}}\uparrow}({{\boldsymbol{r}}})=-\frac{\hbar^{2}({{\boldsymbol{k}}}-{\boldsymbol{\kappa}}_{\pm})^{2}}{2m^{*}}\pm\frac{V_{z}}{2}+\sum_{j=1,3,5}2V\cos({\boldsymbol{b}}_{j}\cdot{{\boldsymbol{r}}}\pm\psi).$
(2)

Here, the first term describes the kinetic energy with effective mass
$m^{*}=0.43m_{0}$ ($m_{0}$ is the bare electron mass) and momentum shifts
${\boldsymbol{\kappa}}_{\pm}=\left[4\pi|\theta|/(3a_{0})\right](-\sqrt{3}/2,\mp
1/2)$ ($a_{0}=3.317\,\text{\AA}$ is the monolayer lattice constant), that
account for for the layer rotation in momentum space. The second term
corresponds to a layer potential difference $V_{z}$ due to an out-of-plane
displacement field. The third term models the moiré potential with amplitude
$V$, phase offset $\psi$, and reciprocal lattice vectors
${\boldsymbol{b}}_{j}=C^{j-1}_{3}(4\pi|\theta|/\sqrt{3}a_{0},0)$ ($C_{3}$ is a
$2\pi/3$ rotation).

Besides the individual layer terms, $H_{\uparrow}$ also includes off-diagonal
terms for inter-layer coupling. They read,

$\Delta_{\text{T}}({{\boldsymbol{r}}})=w(1+e^{-i{\boldsymbol{b}}_{2}\cdot{{\boldsymbol{r}}}}+e^{-i{\boldsymbol{b}}_{3}\cdot{{\boldsymbol{r}}}}),$
(3)

where $w$ corresponds to the inter-layer coupling strength.

Having defined the continuum model Hamiltonian, we now proceed by analyzing
its spatial symmetries. We find that the continuum model exhibits a three-fold
rotation symmetry given by
$U^{{\dagger}}_{C_{3}}({{\boldsymbol{r}}})H_{C_{3}({{\boldsymbol{k}}})\uparrow}(C_{3}({{\boldsymbol{r}}}))U_{C_{3}}({{\boldsymbol{r}}})=H_{{{\boldsymbol{k}}}\uparrow}({{\boldsymbol{r}}})$
with the Hamiltonian density
$H_{{{\boldsymbol{k}}}\uparrow}({{\boldsymbol{r}}})$,
$U_{C_{3}}({{\boldsymbol{r}}})=e^{i(1+\tau_{z}/2){\boldsymbol{b}}_{2}\cdot{{\boldsymbol{r}}}_{2}}e^{i{\boldsymbol{b}}_{3}\cdot{{\boldsymbol{r}}}}$,
and a layer-space Pauli matrix $\tau_{z}$. Besides the rotation symmetry,
which is also a microscopic symmetry of the moiré lattice, the continuum model
exhibits an emergent mirror symmetry along $y=0$ given by
$H_{M_{x}({{\boldsymbol{k}}})\uparrow}(M_{x}({{\boldsymbol{r}}}))=H_{{{\boldsymbol{k}}}\downarrow}({{\boldsymbol{r}}})$.
Interestingly, this mirror symmetry flips the two valleys.

As a final step, we plot the moiré bands along high-symmetry lines of the
Brillouin zone for
$(\theta,V,\psi,w)=(5.1^{\circ},9\,\text{meV},128^{\circ},18\,\text{meV})$
bib:Devakul2021  with $V_{z}=0\,\text{meV}$ and $V_{z}=42\,\text{meV}$, see
Fig. 1(b) and (c). We find that finite displacement field induces an
approximate saddle point in the top-most moiré bands
$\xi_{\uparrow/\downarrow}({{\boldsymbol{k}}})$ at
${\boldsymbol{\kappa}}_{\pm}$. At these points, the Fermi velocity is greatly
reduced, which leads to an enhancement of the density of states (DOS) for a
filling of $n\approx 0.9$ electrons per moiré unit cell, see Fig. 1(d). The
emergence of an enhanced DOS near $n\approx 1$ is in accordance with
experiments bib:Wang2020  and signals that correlation effects due to spin-
valley fluctuations may be amplified. In our further analysis, we will fix
above choice for $(\theta,V,\psi,w)$ that we obtained from density functional
calculations bib:Devakul2021  and the field value $V_{z}=42\,\text{meV}$.

Spin-valley fluctuations. We now discuss the ordering instabilities that arise
from the spin-valley fluctuations. To identify the spin-valley fluctuations
with dominant modulation vectors, we will analyze the maxima of the spin-
valley susceptibility,

$\chi_{\alpha\beta}({\boldsymbol{q}},\omega)=-\frac{1}{\beta
N}\sum_{{{\boldsymbol{k}}}}\text{Tr}\left[\sigma_{\alpha}G({{\boldsymbol{k}}},i\omega)\sigma_{\beta}G({{\boldsymbol{k}}}+{\boldsymbol{q}},i\omega)\right].$
(4)

In this definition,
$\sigma_{\alpha,\beta}\in\\{\sigma_{z},\sigma_{\pm}=(\sigma_{x}\pm
i\sigma_{y})/2\\}$ are spin/valley-space Pauli matrices, ‘Tr’ denotes the
trace over the spin/valley indices, $\beta$ is the inverse temperature, and
$N$ is the number of superlattice unit cells. Furthermore, the non-interacting
Green’s function is given by,

$G_{s_{1}s_{2}}({{\boldsymbol{k}}},i\omega)=\delta_{s_{1}s_{2}}[i\omega-\xi_{s_{1}}({{\boldsymbol{k}}})+\mu]^{-1}.$
(5)

Here, $\xi_{s}({{\boldsymbol{k}}})$ are the dispersions for the top-most moiré
band of the $\pm K$ valleys with
$\xi_{\downarrow}({{\boldsymbol{k}}})=\xi_{\uparrow}(-{{\boldsymbol{k}}})$ due
to time-reversal symmetry, $\mu$ is the chemical potential, and $\omega$
corresponds to a fermionic Matsubara frequency.

Figure 2: (Color online) (a) Dispersion $\xi_{\uparrow}({\boldsymbol{q}})$ of
the top-most band at the $+K$ valleys. The $+K$ ($-K$) valley Fermi surfaces
are shown with solid (dashed) black lines. (b) Transversal spin-valley
susceptibility $\chi_{+-}({\boldsymbol{q}})$ at $1/\beta=100\,$meV. (c)
Longitudinal spin-valley susceptibility $\chi_{zz}({\boldsymbol{q}})$ at
$1/\beta=100\,$meV.

Having defined the spin-valley susceptibility, we will now investigate its
static longitudinal and transversal components,
$\chi_{zz}({\boldsymbol{q}})=\chi_{zz}({\boldsymbol{q}},0)$ and
$\chi_{+-}({\boldsymbol{q}})=\chi_{+-}({\boldsymbol{q}},0)$. In our model,
these two components of the spin-valley susceptibility will show different
dependencies on the modulation vector ${\boldsymbol{q}}$ as the external
displacement fields lifts the $SU(2)$ valley rotation symmetry bib:Pan2020 .

To explain these differences, we focus on the filling of $n\approx 1$ electron
per moiré unit cell because it is the experimentally most relevant case. In
this situation, the Fermi surfaces have the shape of triangles centered around
the $\gamma$ point of the Brillouin zone as shown in Fig. 2(a). Due to the
triangle-like shapes, an approximate Fermi surface nesting feature arises,
mapping the occupied states of one valley onto the other valley’s unoccupied
states via a momentum shift
${{\boldsymbol{Q}}}_{j}=C^{j-1}_{3}{\boldsymbol{\kappa}}_{+}$. As shown in
Fig. 2(b), this approximate nesting enhances $\chi_{+-}({\boldsymbol{q}})$
near ${\boldsymbol{q}}\approx{{\boldsymbol{Q}}}_{j}$, indicating an
instability towards an inter-valley excitonic insulator that pairs electrons
and holes of opposite valleys bib:Bi2021 . However, due to the broken $SU(2)$
valley rotation symmetry, the same features are not present in the
longitudinal valley susceptibility. Instead, as shown in Fig. 2(c),
$\chi_{zz}({\boldsymbol{q}})$ is enhanced near
${\boldsymbol{q}}\approx{\boldsymbol{0}}$, suggesting that the excitonic
instability competes with a spin/valley-polarized state bib:Bultinck2020 .

We now compare these results with the situation away from $n\approx 1$. For
fillings $n\lesssim 0.9$, the Fermi surfaces shrink to smaller triangles,
which modifies the dominant modulation vectors. For example, at $n\approx
0.5$, $\chi_{+-}({\boldsymbol{q}})$ is enhanced near
${\boldsymbol{q}}\approx{{\boldsymbol{Q}}}^{\prime}_{j}$ with
${{\boldsymbol{Q}}}^{\prime}_{j}=C^{j-1}_{3}{\boldsymbol{\kappa}}_{-}$, while
$\chi_{zz}({\boldsymbol{q}})$ shows maxima near the Brillouin zone boundaries.
For $n>1$, the Fermi surfaces are center at the Brillouin zone corners. At
$n\approx 1.5$, this change in Fermi surface topology enhances
$\chi_{+-}({\boldsymbol{q}})$ near ${\boldsymbol{q}}\approx{\boldsymbol{0}}$,
while $\chi_{zz}({\boldsymbol{q}})$ is again enhanced near the Brillouin zone
edges.

Superconductivity. We will now demonstrate that spin-valley fluctuations can
mediate superconductivity by a Kohn-Luttinger mechanism that realizes Cooper
pairing from a nominally repulsive interaction term bib:Kohn1965 ;
bib:Raghu2010 .

To start, we introduce the Hubbard Hamiltonian,

$H=\sum_{{{\boldsymbol{k}}}s}\xi_{s}({{\boldsymbol{k}}})c^{{\dagger}}_{{{\boldsymbol{k}}}s}c_{{{\boldsymbol{k}}}s}+\frac{U}{N}\sum_{\begin{subarray}{c}{{\boldsymbol{k}}},{\boldsymbol{p}},{\boldsymbol{q}}\\\
s,s^{\prime}\end{subarray}}c^{{\dagger}}_{{\boldsymbol{q}}s}c^{{\dagger}}_{{\boldsymbol{p}}+{{\boldsymbol{k}}}-{\boldsymbol{q}}s^{\prime}}c_{{\boldsymbol{p}}s^{\prime}}c_{{{\boldsymbol{k}}}s},$
(6)

where the quadratic part describes a tight-binding model of electrons
$c_{{{\boldsymbol{k}}}s}$ with momentum ${{\boldsymbol{k}}}$ and spin $s$ on a
triangular lattice bib:Pan2020 . The dispersion $\xi_{s}({{\boldsymbol{k}}})$
is given by the top-most moiré band of the continuum model. In particular, as
we omit multi-band effects, there is an emergent mirror symmetry $M_{y}$ along
$x=0$ with $\xi_{s}(M_{y}{{\boldsymbol{k}}})=\xi_{s}({{\boldsymbol{k}}})$. The
interaction term in the Hamiltonian is obtained by projecting the Coulomb
interaction onto the Wannier orbitals of the continuum model’s top-most moiré
band. Since the Wannier orbitals are well-localized, the effective interaction
is short-ranged and repulsive, $U>0$.

Next, we compute the Cooper channel interaction vertex, which describes
scattering of electrons with spin polarization $s$ and $s^{\prime}$ from the
Fermi surface momenta
$(\hat{{\boldsymbol{k}}}^{\prime},-\hat{{\boldsymbol{k}}}^{\prime})$ to
$(\hat{{\boldsymbol{k}}},-\hat{{\boldsymbol{k}}})$. In general, the
interaction vertex comprises a contribution from opposite-spin and equal-spin
scattering. However, the equal-spin scattering leads to pairing with a
‘${{\boldsymbol{d}}}$-vector’ pointed orthogonal to the spin polarization
bib:Sigrist2005 . Such an equal-spin pairing exhibits a suppressed critical
temperature bib:Frigeri2004 ; bib:Frigeri2004_2  and, hence, our focus will be
on opposite-spin scattering. To second order in $U$, we find that the
interaction vertex bib:supplemental ,

$\Gamma(\hat{{\boldsymbol{k}}},\hat{{\boldsymbol{k}}}^{\prime})=U+U^{2}\chi_{+-}(\hat{{\boldsymbol{k}}}+\hat{{\boldsymbol{k}}}^{\prime}).$
(7)

Notably, this result is different for TMD monolayers, where a parabolic
dispersion induces a momentum dependence only at order $U^{3}$ bib:Hsu2014 .

Next, we perform a mean-field decoupling of the effective pairing interaction
and compute the superconducting order parameter from the self-consistency
equation bib:supplemental ,

$\int_{\text{FS}_{\uparrow}}\frac{d\hat{{\boldsymbol{k}}}^{\prime}}{v_{\uparrow}(\hat{{\boldsymbol{k}}}^{\prime})}\Gamma(\hat{{\boldsymbol{k}}},\hat{{\boldsymbol{k}}}^{\prime})\Delta(\hat{{\boldsymbol{k}}}^{\prime})=\lambda\Delta(\hat{{\boldsymbol{k}}}),$
(8)

Here,
$v_{\uparrow}({{\boldsymbol{k}}})=||\partial\xi_{\uparrow}({{\boldsymbol{k}}})/\partial{{\boldsymbol{k}}}||$
is the Fermi velocity and $\Delta(\hat{{\boldsymbol{k}}})$ is the
superconducting order parameter along the $+K$ valley Fermi surface,
$\xi_{\uparrow}({{\boldsymbol{k}}})=\mu$. The order parameters at the $-K$
Fermi surface, $\xi_{\downarrow}({{\boldsymbol{k}}})=\mu$, is
$\Delta^{\prime}(\hat{{\boldsymbol{k}}})=-\Delta(-\hat{{\boldsymbol{k}}})$.
Three remarks about Eq. (8) are in order:

Figure 3: (Color online) (a) Pairing strength given by the most negative
eigenvalue $\lambda^{(0)}$ of the the self-consistency equation as a function
of filling $n$ with $t=\hbar^{2}\theta^{2}/2m^{*}a^{2}_{0}$. For fillings
$1.4\lesssim n\lesssim 1.45$, the most negative eigenvalue is in the $1d$
$A_{2}$-representation of the $C_{3v}$ point group. At all other fillings, the
most negative eigenvalue is in the $2d$ $E$-representation. The Chern numbers
in the marked regions are
$(|\mathcal{C}_{1}|,|\mathcal{C}_{2}|,|\mathcal{C}_{3}|,|\mathcal{C}_{4}|)=(4,2,4,4)$.
(b) Critical temperature normalized by the bandwidth $W$ for $(U/t)^{2}=0.5$,
showing a pronounced enhancement near $n\approx n_{\text{max}}=0.9$.

(1) The negative eigenvalues $\lambda<0$ correspond to sectors with an
attractive interaction in the Cooper channel and grow in magnitude under
renormalization. In particular, the most negative eigenvalue $\lambda^{(0)}$
is most relevant in the renormalization group sense, triggering the leading
superconducting instability at $T_{c}\propto\exp(-1/|\lambda^{(0)}|)$
bib:Raghu2010 .

(2) The eigenvectors $\Delta({{\boldsymbol{k}}})$ correspond to form factors
of the superconducting order parameters and transform as irreducible
representations $\\{A_{1},A_{2},E\\}$ of the point group $C_{3v}$ that is
generated by $C_{3}$ and $M_{y}$. In particular, since $C_{3v}$ has no
inversion center, singlet and triplet pairings mix in tWSe2 bib:Hsu2014 ;
bib:Yuan2014 .

(3) In the weak coupling limit, the term $\propto U$ in
$\Gamma(\hat{{\boldsymbol{k}}},\hat{{\boldsymbol{k}}}^{\prime})$ is
parametrically larger than the term $\propto U^{2}$. For an attractive
interaction, the eigenvectors thus need to annihilate the term $\propto U$
exactly. Such an annihilation implies that $s$-wave pairing,
$\Delta^{(s)}({{\boldsymbol{k}}})\propto 1$, is suppressed. Moreover, since
$s$-wave pairing mixes with $f$-wave pairing (and with higher harmonics of
$A_{1}$), also these pairings will be unfavorable bib:Hsu2014 . Hence, the
dominant superconducting instability will arise from the $A_{2}$ or $E$
representation.

To determine the leading eigenvalue $\lambda^{(0)}$ and the representation of
the associated superconducting order parameter, we solve Eq. (8) numerically
for a broad range of fillings $n\in[0.4,1.6]$. As shown in Fig. 3 and 4, we
find for $1.4\lesssim n\lesssim 1.45$ that $\lambda^{(0)}$ is in the $A_{2}$
representation with a single-component order parameter
$\Delta({{\boldsymbol{k}}})\propto\Delta_{1}({{\boldsymbol{k}}})$. In
contrast, for all other $n$, we find that $\lambda^{(0)}$ is in the $E$
representation with a two-component order parameter,
$\Delta({{\boldsymbol{k}}})=\eta_{1}\Delta_{1}({{\boldsymbol{k}}})+\eta_{2}\Delta_{2}({{\boldsymbol{k}}})$.
Notably, we also observe that the magnitude of $\lambda^{(0)}$, which we
associate with the pairing strength, is greatly enhanced for $n\approx 0.9$.
This enhanced pairing strength leads to a significant enhancement in the
critical temperature $T_{c}$ and suggest a pronounced tendency towards
superconductivity at $n=0.9$.

To understand the origin of the $E$ representation superconductivity, we note
that the DOS of $\xi_{+}({{\boldsymbol{k}}})$ is maximal near the saddle
points ${{\boldsymbol{Q}}}^{\prime}_{j}=C^{j-1}_{3}{\boldsymbol{\kappa}}_{-}$
for $n\approx 1$. The amplitude of the scattering processes between these
points is
$\propto\chi\equiv\chi_{+-}({{\boldsymbol{Q}}}^{\prime}_{i}+{{\boldsymbol{Q}}}^{\prime}_{j})$,
which are exactly the points where $\chi_{+-}({\boldsymbol{q}})$ is maximal.
We thus construct a three-patch model for the order parameters
$\boldsymbol{\Delta}=[\Delta({{\boldsymbol{Q}}}^{\prime}_{1}),\Delta({{\boldsymbol{Q}}}^{\prime}_{2}),\Delta({{\boldsymbol{Q}}}^{\prime}_{3})]^{T}$
that (analogous to Eq. (8)) is given by the eigenvalue equation
$\chi(\boldsymbol{P}-\boldsymbol{1})\boldsymbol{\Delta}=\lambda\boldsymbol{\Delta}$
with $(\boldsymbol{P})_{ij}=1$ and $(\boldsymbol{1})_{ij}=\delta_{ij}$. The
most negative eigenvalue $\lambda^{(0)}=-\chi$ is indeed found to be doubly
degenerate with eigenvectors $(-1,1,0)^{T}$ and $(1,1,-2)^{T}$ in the $E$
representation. In contrast, the positive eigenvalue $\lambda^{(1)}=2\chi$ is
in the $A$ representation with eigenvector $(1,1,1)^{T}$.

Finally, we want to identify the superconducting state
$\boldsymbol{\eta}=(\eta_{1},\eta_{2})$ realized in the $E$ representation of
our model. For that purpose, we write down the ($C_{3v}$-symmetric) Ginzburg-
Landau free energy for our system,

$F_{\text{GL}}=a(T-T_{c})\,\boldsymbol{\eta}\cdot\boldsymbol{\eta}^{*}+b_{1}(\boldsymbol{\eta}\cdot\boldsymbol{\eta}^{*})^{2}+b_{2}|\boldsymbol{\eta}\cdot\boldsymbol{\eta}|^{2}+...,$
(9)

where $a>0$ induces the superconducting phase if $T<T_{c}$ and $b_{1}>0$
ensures thermodynamic stability. The sign of $b_{2}$ fixes the form of the
superconducting state: If $b_{2}<0$, $F_{\text{GL}}$ is minimized by the
nematic state with $\boldsymbol{\eta}=(\cos\varphi,\sin\varphi)$ and
$\varphi\in[0,2\pi)$. The latter spontaneously breaks the three-fold rotation
symmetry of the lattice but preserves time-reversal symmetry. If $b_{2}>0$,
$F_{\text{GL}}$ is minimized by the chiral state with
$\boldsymbol{\eta}=(1,\pm i)$. The chiral state breaks time-reversal symmetry
spontaneously and is characterized by a finite Chern number $\mathcal{C}$
counting the number of chiral Majorana edge modes, see Fig. 3. If inversion
symmetry is strongly broken
$|\xi_{\uparrow}({{\boldsymbol{k}}})-\xi_{\downarrow}({{\boldsymbol{k}}})|\gg|\Delta({{\boldsymbol{k}}})|$,
the chiral state should be energetically favorable bib:Scheurer2017 . However,
if inversion symmetry breaking is weak or the sample is strained, the nematic
and chiral state could be competitive and the two states could be
distinguished by anisotropies in the upper critical field bib:Venderbos2016 ;
bib:Cao2021_2 .

Figure 4: (Color online) Linear independent gap functions of the leading
superconducting state in the $E$ representation.

Conclusion. We have studied a realistic model for superconductivity in tWSe2
based on spin-valley fluctuations. We have shown that the enhancement of spin-
valley fluctuations at half-filling can provide a mechanism for Cooper pairing
and explain the experimental observation of a candidate superconducting state
bib:Wang2020 . Since our results are asymptotically exact for small
interaction strength, we hope that our work will provide a starting point for
superconductivity studies in twisted TMDs.

## References

  * (1) S. Manzeli, D. Ovchinnikov, D. Pasquier, O. V. Yazyev, and A. Kis, Nat. Rev. Mat. 2, 17033 (2017).
  * (2) X. Qian, J. Liu, L. Fu, and J. Li, Science 346, 1344 (2014).
  * (3) S. Tang, C. Zhang, D. Wong, Z. Pedramrazi, H.-Z. Tsai, C. Jia, B. Moritz, M. Claassen, H. Ryu, S. Kahn, J. Jiang, H. Yan, M. Hashimoto, D. Lu, R. G. Moore, C.-C. Hwang, C. Hwang, Z. Hussain, Y. Chen, M. M. Ugeda, Z. Liu, X. Xie, T. P. Devereaux, M. F. Crommie, S.-K. Mo, and Z.-X. Shen, Nat. Phys. 13, 683 (2017).
  * (4) S. Wu, V. Fatemi, Q. D. Gibson, K. Watanabe, T. Taniguchi, R. J. Cava, and P. Jarillo-Herrero, Science 359, 76 (2018).
  * (5) Y. Shi, J. Kahn, B. Niu, Z. Fei, B. Sun, X. Cai, B. A. Francisco, D. Wu, Z.-X. Shen, X. Xu, D. H. Cobden, and Y.-T. Cui, Sci. Adv. 5, eaat8799 (2019).
  * (6) J. M. Lu, O. Zeliuk, I. Leermakers, N. F. Q. Yuan, U. Zeitler, K. T. Law, J. T. Ye, Science 350, 1353 (2015).
  * (7) X. Xi, Z. Wang, W. Zhao, J.-H. Park, K. T. Law, H.Berger, L. Forró, J. Shan, and K. F. Mak, Nat. Phys. 12, 139 (2016).
  * (8) E. Sohn, X. Xi, W.-Y. He, S. Jiang, Z. Wang, K. Kang, J.-H. Park, H. Berger, L. Forró, K. T. Law, J. Shan, and K. F. Mak, Nat. Mater. 17, 504 (2018).
  * (9) C. Sergio, M. R. Sinko, D. P. Gopalan, N. Sivadas, K. L. Seyler, K. Watanabe, T. Taniguchi, A.W. Tsen, X. Xu, D. Xiao, and B. Hunt, Nat. Commun. 9, 1427 (2018).
  * (10) E. C. Regan, D. Wang, C. Jin, M. I. B. Utama, B. Gao, X. Wei, S. Zhao, W. Zhao, K. Yumigeta, M. Blei, J. Carlstroem, K. Watanabe, T. Taniguchi, S. Tongay, M. Crommie, A. Zettl, and F. Wang, Nature 579, 359 (2020).
  * (11) Y. Tang, L. Li, T. Li, Y. Xu, S. Liu, K. Barmak, K. Watanabe, T. Taniguchi, A. H. MacDonald, J. Shan, and K. F. Mak, Nature 579, 353 (2020).
  * (12) C. Jin, Z. Tao, T. Li, Y. Xu, Y. Tang, J. Zhu, S. Liu, K. Watanabe, T. Taniguchi, J. C. Hone, L. Fu, J. Shan, and K. F. Mak, Nat. Mater. 20, 940 (2021)
  * (13) A. Ghiotto, E.-M. Shih, G. S. S. G. Pereira, D. A. Rhodes, B. Kim, J. Zang, A. J. Millis, K. Watanabe, T. Taniguchi, J. C. Hone, L. Wang, C. R. Dean, and Abhay N. Pasupathy, arXiv:2103.09796 (2021).
  * (14) L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M. Kennes, Y. Bai, B. Kim, K. Watanabe, T. Taniguchi, X. Zhu, J. Hone, A. Rubio, A. N. Pasupathy, and C. R. Dean, Nature Materials 19, 861 (2020).
  * (15) Z. Zhang, Y. Wang, K. Watanabe, T. Taniguchi, K. Ueno, E. Tutuc, and B. J. LeRoy, Nat. Phys. 16, 1093 (2020).
  * (16) T. Li, S. Jiang, L. Li, Y. Zhang, K. Kang, J. Zhu, K. Watanabe, T. Taniguchi, D. Chowdhury, L. Fu, J. Shan, and K. F. Mak, arXiv:2103.09779 (2021).
  * (17) T. Li, S. Jiang, B. Shen, Y. Zhang, L. Li, T. Devakul, K. Watanabe, T. Taniguchi, L. Fu, J. Shan, and K. F. Mak, arXiv:2107.01796 (2021).
  * (18) J.Zang, J. Wang, J. Cano, and A. J. Millis, Phys. Rev. B 104, 075150 (2021).
  * (19) F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Phys. Rev. Lett. 121, 026402 (2018).
  * (20) F. Wu, T. Lovorn, E. Tutuc, I. Martin, and A. H. MacDonald, Phys. Rev. Lett. 122, 086402 (2019).
  * (21) C. Schrade and L. Fu, Phys. Rev. B 100, 035413 (2019).
  * (22) H. Pan, F. Wu, and S. Das Sarma, Phys. Rev. Research 2, 033087 (2020).
  * (23) H. Pan, F. Wu, and S. Das Sarma, Phys. Rev. B 102, 201104 (2020).
  * (24) Y. Zhang, N. F. Q. Yuan, and L. Fu, Phys. Rev. B 102, 201115(R) (2020).
  * (25) Y. Zhang, T. Liu, and L. Fu, Phys. Rev. B 103, 155142 (2021).
  * (26) T. Devakul, V. Crépel, Y. Zhang, and L. Fu, arXiv:2106.11954 (2021).
  * (27) Y. Zhang, T. Devakul, and L. Fu, arXiv:2107.02167 (2021).
  * (28) Y.-T. Hsu, F. Wu, S. Das Sarma, arXiv:2105.02415 (2021).
  * (29) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018).
  * (30) M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Science 363, 1059 (2019).
  * (31) G. Chen, A. L. Sharpe, P. Gallagher, I. T. Rosen, E. J. Fox, L. Jiang, B. Lyu, H. Li, K. Watanabe, T. Taniguchi, J. Jung, Z. Shi, D. Goldhaber-Gordon, Y. Zhang, and F. Wang, Nature 572, 215 (2019).
  * (32) X. Liu, Z. Hao, E. Khalaf, J. Y. Lee, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, Nature 583, 221-225 (2020).
  * (33) J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 590, 249 (2021).
  * (34) Z. Hao, A. Zimmerman, P. Ledwith, E. Khalaf, D. H. Najafabadi, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, Science 371, 1133 (2021).
  * (35) Y. Cao, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, arXiv:2103.12083 (2021).
  * (36) S. Raghu, S. A. Kivelson, and D. J. Scalapino, Phys. Rev. B 81, 224505 (2010).
  * (37) Z. Bi and L. Fu, Nat. Commun. 12, 642 (2021).
  * (38) N. Bultinck, S. Chatterjee, and M. P. Zaletel, Phys. Rev. Lett. 124, 166601 (2020).
  * (39) M. Sigrist, AIP Conference Proceedings 789, 165 (2005).
  * (40) P. A. Frigeri, D. F. Agterberg, A. Koga, and M. Sigrist, Phys. Rev. Lett. 92, 097001 (2004).
  * (41) P. A. Frigeri, D.F. Agterberg, M. Sigrist, New J. Phys. 6, 115 (2004).
  * (42) W. Kohn and J. M. Luttinger, Phys. Rev. Lett. 15, 524 (1965).
  * (43) In the Supplemental Material, we provide details on the derivation of the interaction vertex, the self-consistency equation, and the Ginzburg-Landau free energy. We also discuss an effective triangular-lattice Hubbard model.
  * (44) Y.-T. Hsu, A. Vaezi, M. H. Fischer, E.-A. Kim, Nat. Commun. 8, 14985 (2017).
  * (45) N. F. Q. Yuan, K. F. Mak, K. T. Law, Phys. Rev. Lett. 113, 097001 (2014).
  * (46) W. Cho, R. Thomale, S. Raghu, S. A. Kivelson, Phys. Rev. B 88, 064505 (2013).
  * (47) M. S. Scheurer, D. F. Agterberg, and J. Schmalian, npj Quant. Mater. 2, 9 (2017).
  * (48) J. W. F. Venderbos, V. Kozii, and L. Fu, Phys. Rev. B 94, 094522 (2016).
  * (49) Y. Cao, D. Rodan-Legrain, J. M. Park, F. N. Yuan, K. Watanabe, T. Taniguchi, R. M. Fernandes, L. Fu, and P. Jarillo-Herrero, Science 372, 264 (2021).

Supplemental Material to ‘Nematic, chiral and topological superconductivity in
transition metal dichalcogenides’

Constantin Schrade1 and Liang Fu2

1 Center for Quantum Devices, Niels Bohr Institute, University of Copenhagen,
2100 Copenhagen, Denmark
2 Department of Physics, Massachusetts Institute of Technology, 77
Massachusetts Ave., Cambridge, MA 02139

In the Supplemental Material, we provide details on the derivation of the
interaction vertex, the self-consistency equation, and the Ginzburg-Landau
free energy. We also discuss an additional effective triangular-lattice
Hubbard model.

##  I Effective pairing interaction

In this first section of the Supplemental Material (SM), we provide details on
the derivation of the effective pairing interaction from a Kohn-Luttinger
renormalization of the bare repulsive interaction.
\- As a first step in this derivation, we write the partition function of the
Hamiltonian $H$ (as given in Eq. (6) of the main text) in terms of a coherent
state Feynman path integral over Grassman variables,

$\begin{split}Z&=\int\mathcal{D}[c^{*},c]\,e^{-S_{0}-S_{I}},\\\
S_{0}&=\sum_{s}\sum_{{{\boldsymbol{k}}},\omega_{n}}c^{*}_{{{\boldsymbol{k}}},s}(i\omega_{n})\left[-i\omega_{n}+\xi_{s}({{\boldsymbol{k}}})\right]c_{{{\boldsymbol{k}}},s}(i\omega_{n}),\\\
S_{I}&=\frac{U}{\beta
N}\sum_{s,s^{\prime}}\sum_{\begin{subarray}{c}{{\boldsymbol{k}}},{{\boldsymbol{k}}}^{\prime},{\boldsymbol{q}}\\\
\omega_{n},\omega^{\prime}_{n},\omega_{\nu}\end{subarray}}c^{*}_{{{\boldsymbol{k}}}+{\boldsymbol{q}},s}(i\omega_{n}+i\omega_{\nu})c^{*}_{{{\boldsymbol{k}}}^{\prime}-{\boldsymbol{q}},s^{\prime}}(i\omega^{\prime}_{n}-i\omega_{\nu})c_{{{\boldsymbol{k}}}^{\prime},s^{\prime}}(i\omega^{\prime}_{n})c_{{{\boldsymbol{k}}},s}(i\omega_{n}).\end{split}$
(10)

Next, we introduce a low-energy effective action for modes in the vicinity of
the Fermi surfaces bib:Shankar1994 ; bib:Vafek2014 ,

$S^{\text{eff}}_{I}=\langle S_{I}\rangle-\frac{1}{2}\left(\langle
S^{2}_{I}\rangle-\langle S_{I}\rangle^{2}\right),$ (11)

where $\langle...\rangle$ denotes the averaging over ‘high’ energy modes with
action $S_{0}$. We note that the first term in the expression for
$S^{\text{eff}}_{I}$ corresponds to the bare repulsive interaction and is
$\propto U$. In contrast, the second term renormalizes the bare repulsive
interaction and is $\propto U^{2}$.

To evaluate the renormalization of the bare repulsive interaction, we write it
in the more compact form,

$S^{\text{eff}}_{I}-\langle
S_{I}\rangle=-\frac{1}{2}\left[\sum_{1,2,3,4}\sum_{1^{\prime},2^{\prime},3^{\prime},4^{\prime}}\langle
c^{*}_{1}c^{*}_{2}c_{3}c_{4}c^{*}_{1^{\prime}}c^{*}_{2^{\prime}}c_{3^{\prime}}c_{4^{\prime}}\rangle-\langle
S_{I}\rangle^{2}\right]=-\frac{1}{2}\sum_{1,2,3,4}\sum_{1^{\prime},2^{\prime},3^{\prime},4^{\prime}}\langle
c^{*}_{1}c^{*}_{2}c_{3}c_{4}c^{*}_{1^{\prime}}c^{*}_{2^{\prime}}c_{3^{\prime}}c_{4^{\prime}}\rangle_{\text{con}}.$
(12)

Here, we have introduced the short-hand notation $\sum_{1,2,3,4}\equiv(U/\beta
N)\sum_{s,s^{\prime}}\sum_{{{\boldsymbol{k}}},{{\boldsymbol{k}}}^{\prime},{\boldsymbol{q}}}\sum_{\omega_{n},\omega^{\prime}_{n},\omega_{\nu}}$.
Moreover, in the second equality, we have noted that the Wick contractions
which correspond to disconnected Feynman diagrams in $\langle
S^{2}_{I}\rangle$ are cancelled by $\langle S_{I}\rangle^{2}$. In particular,
the average $\langle...\rangle_{\text{con}}$ only involves Wick contractions
that correspond to connected Feynman diagrams.

Within the Kohn-Luttinger approach, the corrections to the bare repulsive
interactions are given by the Wick contractions in the particle-hole channel.
These Wick contractions are classified into three types,

$\begin{split}\text{Type
1:}\qquad&\langle\wick{\c{1}c^{*}_{1}c^{*}_{2}c_{3}\c{2}c_{4}\c{2}c^{*}_{1^{\prime}}c^{*}_{2^{\prime}}\c{1}c_{3^{\prime}}c_{4^{\prime}}}\rangle,\quad\langle\wick{\c{1}c^{*}_{1}c^{*}_{2}c_{3}\c{2}c_{4}c^{*}_{1^{\prime}}\c{2}c^{*}_{2^{\prime}}c_{3^{\prime}}\c{1}c_{4^{\prime}}}\rangle,\quad\langle\wick{c^{*}_{1}\c{1}c^{*}_{2}\c{2}c_{3}c_{4}\c{2}c^{*}_{1^{\prime}}c^{*}_{2^{\prime}}\c{1}c_{3^{\prime}}c_{4^{\prime}}}\rangle,\quad\langle\wick{c^{*}_{1}\c{1}c^{*}_{2}\c{2}c_{3}c_{4}c^{*}_{1^{\prime}}\c{2}c^{*}_{2^{\prime}}c_{3^{\prime}}\c{1}c_{4^{\prime}}}\rangle\\\
&\langle\wick{\c{1}c^{*}_{1}c^{*}_{2}\c{2}c_{3}c_{4}\c{2}c^{*}_{1^{\prime}}c^{*}_{2^{\prime}}c_{3^{\prime}}\c{1}c_{4^{\prime}}}\rangle,\quad\langle\wick{c^{*}_{1}\c{1}c^{*}_{2}c_{3}\c{2}c_{4}\c{2}c^{*}_{1^{\prime}}c^{*}_{2^{\prime}}c_{3^{\prime}}\c{1}c_{4^{\prime}}}\rangle,\quad\langle\wick{\c{1}c^{*}_{1}c^{*}_{2}\c{2}c_{3}c_{4}c^{*}_{1^{\prime}}\c{2}c^{*}_{2^{\prime}}\c{1}c_{3^{\prime}}c_{4^{\prime}}}\rangle,\quad\langle\wick{c^{*}_{1}\c{1}c^{*}_{2}c_{3}\c{2}c_{4}c^{*}_{1^{\prime}}\c{2}c^{*}_{2^{\prime}}\c{1}c_{3^{\prime}}c_{4^{\prime}}}\rangle\\\
\\\ \text{Type
2:}\qquad&\langle\wick{c^{*}_{1}\c{1}c^{*}_{2}\c{2}c_{3}c_{4}\c{2}c^{*}_{1^{\prime}}c^{*}_{2^{\prime}}c_{3^{\prime}}\c{1}c_{4^{\prime}}}\rangle,\quad\langle\wick{c^{*}_{1}\c{1}c^{*}_{2}\c{2}c_{3}c_{4}c^{*}_{1^{\prime}}\c{2}c^{*}_{2^{\prime}}\c{1}c_{3^{\prime}}c_{4^{\prime}}}\rangle,\quad\langle\wick{\c{1}c^{*}_{1}c^{*}_{2}c_{3}\c{2}c_{4}c^{*}_{1^{\prime}}\c{2}c^{*}_{2^{\prime}}\c{1}c_{3^{\prime}}c_{4^{\prime}}}\rangle,\quad\langle\wick{\c{1}c^{*}_{1}c^{*}_{2}c_{3}\c{2}c_{4}\c{2}c^{*}_{1^{\prime}}c^{*}_{2^{\prime}}c_{3^{\prime}}\c{1}c_{4^{\prime}}}\rangle\\\
\\\ \text{Type
3:}\qquad&\langle\wick{c^{*}_{1}\c{1}c^{*}_{2}c_{3}\c{2}c_{4}\c{2}c^{*}_{1^{\prime}}c^{*}_{2^{\prime}}\c{1}c_{3^{\prime}}c_{4^{\prime}}}\rangle,\quad\langle\wick{\c{1}c^{*}_{1}c^{*}_{2}\c{2}c_{3}c_{4}c^{*}_{1^{\prime}}\c{2}c^{*}_{2^{\prime}}c_{3^{\prime}}\c{1}c_{4^{\prime}}}\rangle,\quad\langle\wick{\c{1}c^{*}_{1}c^{*}_{2}\c{2}c_{3}c_{4}\c{2}c^{*}_{1^{\prime}}c^{*}_{2^{\prime}}\c{1}c_{3^{\prime}}c_{4^{\prime}}}\rangle,\quad\langle\wick{c^{*}_{1}\c{1}c^{*}_{2}c_{3}\c{2}c_{4}c^{*}_{1^{\prime}}\c{2}c^{*}_{2^{\prime}}c_{3^{\prime}}\c{1}c_{4^{\prime}}}\rangle\end{split}$
(13)

For each of these three types of Wick contractions, we show example Feynman
diagrams in Fig. 5. We now proceed with the evaluation of the different types
of Wick contractions:

###  I.1 Wick contractions of type 1 and 2

We begin by evaluation the Wick contractions of type 1 and, for that purpose,
consider the following example,

$\displaystyle\ \ \
-\frac{1}{2}\sum_{1,2,3,4}\sum_{1^{\prime},2^{\prime},3^{\prime},4^{\prime}}\langle\wick{\c{1}c^{*}_{1}c^{*}_{2}c_{3}\c{2}c_{4}\c{2}c^{*}_{1^{\prime}}c^{*}_{2^{\prime}}\c{1}c_{3^{\prime}}c_{4^{\prime}}}\rangle$
(14) $\displaystyle=-\frac{1}{2(\beta
N)^{2}}\sum_{\begin{subarray}{c}{{\boldsymbol{k}}},{{\boldsymbol{k}}}^{\prime},{\boldsymbol{q}}\\\
\omega_{n},\omega^{\prime}_{n},\omega_{\nu}\end{subarray}}G_{\uparrow}({\boldsymbol{q}},i\omega_{\nu})G_{\uparrow}({\boldsymbol{q}}+{{\boldsymbol{k}}}-{{\boldsymbol{k}}}^{\prime},i\omega_{\nu}+i\omega_{n}+i\omega^{\prime}_{n})\langle
c^{*}_{{{\boldsymbol{k}}},\uparrow}(i\omega_{n})c^{*}_{-{{\boldsymbol{k}}},\downarrow}(-i\omega_{n})c_{-{{\boldsymbol{k}}}^{\prime},\downarrow}(-i\omega^{\prime}_{n})c_{{{\boldsymbol{k}}}^{\prime},\uparrow}(i\omega^{\prime}_{n})\rangle$
$\displaystyle\quad-\frac{1}{2(\beta
N)^{2}}\sum_{\begin{subarray}{c}{{\boldsymbol{k}}},{{\boldsymbol{k}}}^{\prime},{\boldsymbol{q}}\\\
\omega_{n},\omega^{\prime}_{n},\omega_{\nu}\end{subarray}}G_{\downarrow}({\boldsymbol{q}},i\omega_{\nu})G_{\downarrow}({\boldsymbol{q}}+{{\boldsymbol{k}}}-{{\boldsymbol{k}}}^{\prime},i\omega_{\nu}+i\omega_{n}+i\omega^{\prime}_{n})\langle
c^{*}_{{{\boldsymbol{k}}},\uparrow}(i\omega_{n})c^{*}_{-{{\boldsymbol{k}}},\downarrow}(-i\omega_{n})c_{-{{\boldsymbol{k}}}^{\prime},\downarrow}(-i\omega^{\prime}_{n})c_{{{\boldsymbol{k}}}^{\prime},\uparrow}(i\omega^{\prime}_{n})\rangle$

In the second line of the example, we have defined the Green’s function
$G_{s}({{\boldsymbol{k}}},i\omega_{n})=\langle
c^{{\dagger}}_{{{\boldsymbol{k}}},s}(i\omega_{n})c_{{{\boldsymbol{k}}},s}(i\omega_{n})\rangle$
and only kept the terms that lead to a pairing interaction. The remaining Wick
contractions of type 1 give the same result and, therefore, we do not display
them here.

Next, we want to compare this finding with a particular example for a Wick
contraction of type 2,

$\displaystyle-\frac{1}{2}\sum_{1,2,3,4}\sum_{1^{\prime},2^{\prime},3^{\prime},4^{\prime}}\langle\wick{c^{*}_{1}\c{1}c^{*}_{2}\c{2}c_{3}c_{4}\c{2}c^{*}_{1^{\prime}}c^{*}_{2^{\prime}}c_{3^{\prime}}\c{1}c_{4^{\prime}}}\rangle$
(15) $\displaystyle\equiv\frac{1}{(\beta
N)^{2}}\sum_{\begin{subarray}{c}{{\boldsymbol{k}}},{{\boldsymbol{k}}}^{\prime},{\boldsymbol{q}}\\\
\omega_{n},\omega^{\prime}_{n},\omega_{\nu}\end{subarray}}G_{\uparrow}({\boldsymbol{q}},i\omega_{\nu})G_{\uparrow}({\boldsymbol{q}}+{{\boldsymbol{k}}}-{{\boldsymbol{k}}}^{\prime},i\omega_{\nu}+i\omega_{n}+i\omega^{\prime}_{n})\langle
c^{*}_{{{\boldsymbol{k}}},\uparrow}(i\omega_{n})c^{*}_{-{{\boldsymbol{k}}},\downarrow}(-i\omega_{n})c_{-{{\boldsymbol{k}}}^{\prime},\downarrow}(-i\omega^{\prime}_{n})c_{{{\boldsymbol{k}}}^{\prime},\uparrow}(i\omega^{\prime}_{n})\rangle$
$\displaystyle+\frac{1}{(\beta
N)^{2}}\sum_{\begin{subarray}{c}{{\boldsymbol{k}}},{{\boldsymbol{k}}}^{\prime},{\boldsymbol{q}}\\\
\omega_{n},\omega^{\prime}_{n},\omega_{\nu}\end{subarray}}G_{\downarrow}({\boldsymbol{q}},i\omega_{\nu})G_{\downarrow}({\boldsymbol{q}}+{{\boldsymbol{k}}}-{{\boldsymbol{k}}}^{\prime},i\omega_{\nu}+i\omega_{n}+i\omega^{\prime}_{n})\langle
c^{*}_{{{\boldsymbol{k}}},\uparrow}(i\omega_{n})c^{*}_{-{{\boldsymbol{k}}},\downarrow}(-i\omega_{n})c_{-{{\boldsymbol{k}}}^{\prime},\downarrow}(-i\omega^{\prime}_{n})c_{{{\boldsymbol{k}}}^{\prime},\uparrow}(i\omega^{\prime}_{n})\rangle$

Here, we observe that this finding differs from the previous one by a minus
sign and a prefactor of 2. Since the remaining Wick contractions of type 2 all
give the same result as this example, we find that the type 1 and type 2
contributions cancel each other and, therefore, do not contribute to the
effective action.

Figure 5: (Color online) (a) Feynman diagrams corresponding to Wick
contractions of type 1 and type 2 within the Kohn-Luttinger renormalization of
the bare repulsive interaction. The closed fermion loop in the Feynman diagram
of type 2 produces an additional minus sign and, therefore, the Feynman
diagrams of type 1 and 2 cancel each other. (b) Feynman diagram corresponding
to Wick contractions of type 3, which yield a finite contribution to the
effective pairing interaction.

###  I.2 Wick contractions of type 3

For illustrating the evaluation of a Wick contraction of type 3, we again pick
a specific example,

$\displaystyle\ \ \
-\frac{1}{2}\sum_{1,2,3,4}\sum_{1^{\prime},2^{\prime},3^{\prime},4^{\prime}}\langle\wick{c^{*}_{1}\c{1}c^{*}_{2}c_{3}\c{2}c_{4}\c{2}c^{*}_{1^{\prime}}c^{*}_{2^{\prime}}\c{1}c_{3^{\prime}}c_{4^{\prime}}}\rangle$
(16) $\displaystyle\equiv-\frac{1}{(\beta
N)^{2}}\sum_{\begin{subarray}{c}{{\boldsymbol{k}}},{{\boldsymbol{k}}}^{\prime},{\boldsymbol{q}}\\\
\omega_{n},\omega^{\prime}_{n},\omega_{\nu}\end{subarray}}G_{\downarrow}({\boldsymbol{q}},i\omega_{\nu})G_{\uparrow}({\boldsymbol{q}}+{{\boldsymbol{k}}}+{{\boldsymbol{k}}}^{\prime},i\omega_{\nu}+i\omega_{n}+i\omega^{\prime}_{n})\langle
c^{*}_{{{\boldsymbol{k}}},\uparrow}(i\omega_{n})c^{*}_{-{{\boldsymbol{k}}},\downarrow}(-i\omega_{n})c_{-{{\boldsymbol{k}}}^{\prime},\downarrow}(-i\omega^{\prime}_{n})c_{{{\boldsymbol{k}}}^{\prime},\uparrow}(i\omega^{\prime}_{n})\rangle\rangle.$

Similar to the previous subsection, we find that all type 3 Wick contractions
give the same contribution as this example. Hence, we find that the amplitude
of the effective pairing interaction reads,

$\begin{split}\Gamma({{\boldsymbol{k}}},{{\boldsymbol{k}}}^{\prime})=U+U^{2}\chi_{+-}({{\boldsymbol{k}}}+{{\boldsymbol{k}}}^{\prime})\quad\text{with}\quad\chi_{+-}({\boldsymbol{q}})&\equiv\chi_{+-}({\boldsymbol{q}},i\omega=0)\\\
&=-\frac{1}{\beta
N}\sum_{{{\boldsymbol{k}}}}G_{\downarrow}({\boldsymbol{q}},i\omega=0)G_{\uparrow}({\boldsymbol{q}}+{{\boldsymbol{k}}},i\omega=0)\\\
&=-\frac{1}{N}\sum_{{{\boldsymbol{k}}}}\frac{f(\xi_{\downarrow}({{\boldsymbol{k}}})-\mu)-f(\xi_{\uparrow}({\boldsymbol{q}}+{{\boldsymbol{k}}})-\mu)}{\xi_{\downarrow}({{\boldsymbol{k}}})-\xi_{\uparrow}({\boldsymbol{q}}+{{\boldsymbol{k}}})}.\end{split}$
(17)

Here, we have set the Matsubara frequencies at the external legs of the
Feynman diagrams to zero and omitted numerical prefactors in $\chi_{+-}$ as
they will not be relevant for the weak coupling study of superconductivity.
Furthermore, we have introduced the Fermi distribution $f$ and carried out the
summation over the Matsubara frequencies by using
$G_{s}({{\boldsymbol{k}}},i\omega_{n})=1/(i\omega_{n}-\xi_{s}({{\boldsymbol{k}}})+\mu)$.
This concludes our derivation of the effective pairing interaction.

##  II Mean-field theory and self-consistency equation

In this second section of the SM, we provide details on the derivation of the
pairing Hamiltonian within mean-field theory and determine the self-
consistency equation for the superconducting order parameter.

In a first step, we consider the effective pairing interaction in the
Hamiltonian form,

$H_{I}^{\text{eff}}=\frac{1}{N}\sum_{{{\boldsymbol{k}}},{{\boldsymbol{k}}}^{\prime}}\Gamma({{\boldsymbol{k}}},{{\boldsymbol{k}}}^{\prime})\,c^{{\dagger}}_{{{\boldsymbol{k}}},\uparrow}c^{{\dagger}}_{-{{\boldsymbol{k}}},\downarrow}c_{-{{\boldsymbol{k}}}^{\prime},\downarrow}c_{{{\boldsymbol{k}}}^{\prime},\uparrow}\quad\text{with}\quad\Gamma({{\boldsymbol{k}}},{{\boldsymbol{k}}}^{\prime})=U+U^{2}\chi_{+-}({{\boldsymbol{k}}}+{{\boldsymbol{k}}}^{\prime}),$
(18)

and perform a mean-field decoupling,

$\begin{split}H_{I}^{\text{eff}}\rightarrow
H_{\Delta}&=\frac{1}{2N}\sum_{{{\boldsymbol{k}}},{{\boldsymbol{k}}}^{\prime}}\Gamma({{\boldsymbol{k}}},{{\boldsymbol{k}}}^{\prime})\,\left[\langle
c^{{\dagger}}_{{{\boldsymbol{k}}},\uparrow}c^{{\dagger}}_{-{{\boldsymbol{k}}},\downarrow}\rangle
c_{-{{\boldsymbol{k}}}^{\prime},\downarrow}c_{{{\boldsymbol{k}}}^{\prime},\uparrow}+c^{{\dagger}}_{{{\boldsymbol{k}}},\uparrow}c^{{\dagger}}_{-{{\boldsymbol{k}}},\downarrow}\langle
c_{-{{\boldsymbol{k}}}^{\prime},\downarrow}c_{{{\boldsymbol{k}}}^{\prime},\uparrow}\rangle\right]\\\
&=-\sum_{{{\boldsymbol{k}}}}\Delta({{\boldsymbol{k}}})\,c_{-{{\boldsymbol{k}}},\downarrow}c_{{{\boldsymbol{k}}},\uparrow}-\sum_{{{\boldsymbol{k}}}}\Delta({{\boldsymbol{k}}})^{*}\,c^{{\dagger}}_{{{\boldsymbol{k}}},\uparrow}c^{{\dagger}}_{-{{\boldsymbol{k}}},\downarrow},\end{split}$
(19)

where we have defined the superconducting order parameter as,

$\begin{split}\Delta({{\boldsymbol{k}}})&=-\frac{1}{2N}\sum_{{{\boldsymbol{k}}}^{\prime}}\Gamma({{\boldsymbol{k}}},{{\boldsymbol{k}}}^{\prime})\,\langle
c^{{\dagger}}_{{{\boldsymbol{k}}}^{\prime},\uparrow}c^{{\dagger}}_{-{{\boldsymbol{k}}}^{\prime},\downarrow}\rangle.\end{split}$
(20)

In a second step, we write the pairing Hamiltonian $H_{\Delta}$ and the normal
state Hamiltonian $H_{0}$ as a Bogoliubov-de Gennes Hamiltonian matrix,

$\begin{split}H&=H_{0}+H_{\Delta}\\\
&=\sum_{{{\boldsymbol{k}}}}\xi_{\uparrow}({{\boldsymbol{k}}})\,c^{{\dagger}}_{{{\boldsymbol{k}}},\uparrow}c_{{{\boldsymbol{k}}},\uparrow}+\xi_{\downarrow}({{\boldsymbol{k}}})\,c^{{\dagger}}_{{{\boldsymbol{k}}},\downarrow}c_{{{\boldsymbol{k}}},\downarrow}-\Delta({{\boldsymbol{k}}})\,c_{-{{\boldsymbol{k}}},\downarrow}c_{{{\boldsymbol{k}}},\uparrow}-\Delta({{\boldsymbol{k}}})^{*}\,c^{{\dagger}}_{{{\boldsymbol{k}}},\uparrow}c^{{\dagger}}_{-{{\boldsymbol{k}}},\downarrow}\\\
&=\sum_{{{\boldsymbol{k}}}}\xi_{\uparrow}({{\boldsymbol{k}}})\left[c^{{\dagger}}_{{{\boldsymbol{k}}},\uparrow}c_{{{\boldsymbol{k}}},\uparrow}+c^{{\dagger}}_{-{{\boldsymbol{k}}},\downarrow}c_{-{{\boldsymbol{k}}},\downarrow}\right]-\Delta({{\boldsymbol{k}}})\,c_{-{{\boldsymbol{k}}},\downarrow}c_{{{\boldsymbol{k}}},\uparrow}-\Delta({{\boldsymbol{k}}})^{*}\,c^{{\dagger}}_{{{\boldsymbol{k}}},\uparrow}c^{{\dagger}}_{-{{\boldsymbol{k}}},\downarrow}\\\
&=\sum_{{{\boldsymbol{k}}}}\begin{pmatrix}c^{{\dagger}}_{{{\boldsymbol{k}}},\uparrow}&c_{-{{\boldsymbol{k}}},\downarrow}\end{pmatrix}\begin{pmatrix}\xi_{\uparrow}({{\boldsymbol{k}}})&-\Delta({{\boldsymbol{k}}})^{*}\\\
-\Delta({{\boldsymbol{k}}})&-\xi_{\uparrow}({{\boldsymbol{k}}})\end{pmatrix}\begin{pmatrix}c_{{{\boldsymbol{k}}},\uparrow}\\\
c^{{\dagger}}_{-{{\boldsymbol{k}}},\downarrow}\end{pmatrix}+E_{0}\end{split}$
(21)

where we have used that
$\xi_{\downarrow}(-{{\boldsymbol{k}}})=\xi_{\uparrow}({{\boldsymbol{k}}})$ and
defined the energy offset $E_{0}$. We now introduce the Bogoliubov
transformations,

$\begin{split}\begin{pmatrix}c_{{{\boldsymbol{k}}},\uparrow}\\\
c^{{\dagger}}_{-{{\boldsymbol{k}}},\downarrow}\end{pmatrix}&=\begin{pmatrix}u^{*}_{{{\boldsymbol{k}}}}&v_{{{\boldsymbol{k}}}}\\\
-v_{{{\boldsymbol{k}}}}&u_{{{\boldsymbol{k}}}}\end{pmatrix}\begin{pmatrix}\gamma_{{{\boldsymbol{k}}},1}\\\
\gamma^{{\dagger}}_{-{{\boldsymbol{k}}},2}\end{pmatrix}\ \ \text{with}\ \
u_{{{\boldsymbol{k}}}}=\frac{\Delta({{\boldsymbol{k}}})}{|\Delta({{\boldsymbol{k}}})|}\sqrt{\frac{1}{2}\left(1+\frac{\xi_{\uparrow}({{\boldsymbol{k}}})}{E({{\boldsymbol{k}}})}\right)},\quad
v_{{{\boldsymbol{k}}}}=\sqrt{\frac{1}{2}\left(1-\frac{\xi_{\uparrow}({{\boldsymbol{k}}})}{E({{\boldsymbol{k}}})}\right)},\end{split}$
(22)

and apply it to the Bogoliubov-de Gennes Hamiltonian matrix. This yields the
diagonalized Hamiltonian,

$\begin{split}H&=\sum_{{{\boldsymbol{k}}}}\begin{pmatrix}\gamma^{{\dagger}}_{{{\boldsymbol{k}}},1}&\gamma_{-{{\boldsymbol{k}}},2}\end{pmatrix}\begin{pmatrix}E({{\boldsymbol{k}}})&0\\\
0&-E({{\boldsymbol{k}}})\end{pmatrix}\begin{pmatrix}\gamma_{{{\boldsymbol{k}}},1}\\\
\gamma^{{\dagger}}_{-{{\boldsymbol{k}}},2},\end{pmatrix}=\sum_{{{\boldsymbol{k}}}}E({{\boldsymbol{k}}})\left[\gamma^{{\dagger},1}_{{{\boldsymbol{k}}}}\gamma_{{{\boldsymbol{k}}},1}+\gamma^{{\dagger}}_{-{{\boldsymbol{k}}},2}\gamma_{-{{\boldsymbol{k}}},2}\right]+E^{\prime}_{0},\end{split}$
(23)

where we have defined the quasiparticle energy
$E({{\boldsymbol{k}}})=\sqrt{\xi_{\uparrow}({{\boldsymbol{k}}})^{2}+|\Delta({{\boldsymbol{k}}})|^{2}}$
and another energy offset $E^{\prime}_{0}$.

In a third step, we evaluate the expectation value,

$\begin{split}\langle
c^{{\dagger}}_{{{\boldsymbol{k}}},\uparrow}c^{{\dagger}}_{-{{\boldsymbol{k}}},\downarrow}\rangle&=\frac{1}{2}\left[-u_{{{\boldsymbol{k}}}}v_{{{\boldsymbol{k}}}}\langle\gamma^{{\dagger}}_{{{\boldsymbol{k}}},1}\gamma_{{{\boldsymbol{k}}},1}\rangle+v_{{{\boldsymbol{k}}}}u_{{{\boldsymbol{k}}}}\langle\gamma_{-{{\boldsymbol{k}}},2}\gamma^{{\dagger}}_{-{{\boldsymbol{k}}},2}\rangle\right]\\\
&=\frac{1}{2}u_{{{\boldsymbol{k}}}}v_{{{\boldsymbol{k}}}}\left[1-\langle\gamma^{{\dagger}}_{{{\boldsymbol{k}}},1}\gamma_{{{\boldsymbol{k}}},1}\rangle-\langle\gamma^{{\dagger}}_{-{{\boldsymbol{k}}},2}\gamma_{-{{\boldsymbol{k}}},2}\rangle\right]\\\
&=\frac{1}{2}u_{{{\boldsymbol{k}}}}v_{{{\boldsymbol{k}}}}\left[1-2f(E({{\boldsymbol{k}}}))\right]\\\
&=\frac{1}{2}u_{{{\boldsymbol{k}}}}v_{{{\boldsymbol{k}}}}\tanh\left(\frac{E({{\boldsymbol{k}}})}{2k_{B}T}\right)\\\
&=\frac{1}{4}\frac{\Delta({{\boldsymbol{k}}})}{E({{\boldsymbol{k}}})}\tanh\left(\frac{E({{\boldsymbol{k}}})}{2k_{B}T}\right).\end{split}$
(24)

Here, $f$ denotes the Fermi-Dirac distribution. We now reinsert the expression
of the expectation value in the definition of the superconducting order
parameter. Under the assumption that the superconducting order parameter is
sufficiently small, we find the self-consistency equation,

$\begin{split}\Delta({{\boldsymbol{k}}})&=-\frac{1}{8N}\sum_{{{\boldsymbol{k}}}^{\prime}}\Gamma({{\boldsymbol{k}}},{{\boldsymbol{k}}}^{\prime})\,\frac{\Delta({{\boldsymbol{k}}}^{\prime})}{\xi_{\uparrow}({{\boldsymbol{k}}}^{\prime})}\tanh\left(\frac{\xi_{\uparrow}({{\boldsymbol{k}}}^{\prime})}{2k_{B}T}\right).\end{split}$
(25)

In a final step, we rewrite this self-consistency equation as an integral
eigenvalue equation,

$\int_{\text{FS}_{\uparrow}}\frac{\mathrm{d}\hat{{\boldsymbol{k}}}^{\prime}}{v_{\uparrow}(\hat{{\boldsymbol{k}}}^{\prime})}\Gamma(\hat{{\boldsymbol{k}}},\hat{{\boldsymbol{k}}}^{\prime})\Delta(\hat{{\boldsymbol{k}}}^{\prime})=\lambda\Delta(\hat{{\boldsymbol{k}}}),$
(26)

which determines the superconducting order parameter on the $+K$ Fermi
surface, $\xi_{\uparrow}({{\boldsymbol{k}}})-\mu=0$, with
$v_{\uparrow}({{\boldsymbol{k}}})=||\partial\xi_{\uparrow}({{\boldsymbol{k}}})/\partial{{\boldsymbol{k}}}||$
denoting the Fermi velocity. We note that the superconducting order parameters
at the $+K$ and $-K$ Fermi surfaces are given by $\Delta({{\boldsymbol{k}}})$
and $\Delta^{\prime}({{\boldsymbol{k}}})=-\Delta(-{{\boldsymbol{k}}})$,
respectively. This concludes our derivation of the self-consistency equation
of the superconducting order parameter.

##  III Numerical evaluation of the self-consistency equation

In this third section of the SM, we give details on the numerical evaluation
of the self-consistency equation.

To start, we write the line integral in the self-consistency equation as an
area integral over a small neighborhood of the $+K$ Fermi surface given by the
level set $\xi_{\uparrow}({{\boldsymbol{k}}})-\mu=\pm\Omega$,

$\int_{\text{FS}_{\uparrow}}\frac{\mathrm{d}\hat{{\boldsymbol{k}}}^{\prime}}{v_{\uparrow}(\hat{{\boldsymbol{k}}}^{\prime})}\Gamma(\hat{{\boldsymbol{k}}},\hat{{\boldsymbol{k}}}^{\prime})\Delta(\hat{{\boldsymbol{k}}}^{\prime})\approx\frac{1}{2\Omega}\int_{|\xi_{\uparrow}({{\boldsymbol{k}}}^{\prime})-\mu|<\Omega}\mathrm{d}{{\boldsymbol{k}}}^{\prime}\,\Gamma(\hat{{\boldsymbol{k}}},\hat{{\boldsymbol{k}}}^{\prime})\Delta(\hat{{\boldsymbol{k}}}^{\prime})\approx\frac{1}{2\Omega}\sum_{\begin{subarray}{c}{{\boldsymbol{k}}}^{\prime}\\\
|\xi_{\uparrow}({{\boldsymbol{k}}}^{\prime})-\mu|<\Omega\end{subarray}}\delta{{\boldsymbol{k}}}^{\prime}\,\Gamma(\hat{{\boldsymbol{k}}},\hat{{\boldsymbol{k}}}^{\prime})\Delta(\hat{{\boldsymbol{k}}}^{\prime}).$
(27)

Here, we have noted that the width of the area defined by the level set is
roughly
$2\Omega|\nabla_{{{\boldsymbol{k}}}}\xi_{\uparrow}({{\boldsymbol{k}}})|^{-1}=2\Omega\,v_{\uparrow}(\hat{{\boldsymbol{k}}}^{\prime})^{-1}$
and divided the integrand by this factor. Subsequently, we have written the
area integral in terms of a Riemann sum.

There are now two alternative approaches for numerically solving the
eigenvalue equation:

A first approach is to fix a set of vectors $\\{{{\boldsymbol{k}}}_{i}\\}$ in
momentum space within the window $|\xi_{+}({{\boldsymbol{k}}})-\mu|<\Omega$
and write the Riemann sum as a matrix $\boldsymbol{A}$ with entries
$[\boldsymbol{A}]_{ij}=\Gamma({{\boldsymbol{k}}}_{i},{{\boldsymbol{k}}}^{\prime}_{j})\,\delta{{\boldsymbol{k}}}^{\prime}_{j}/2\Omega$.
Diagonalization of this matrix yields the desired eigenvalues $\lambda$ and
eigenvectors $\Delta({{\boldsymbol{k}}})$. However, it is important to note
that within this approach the resulting eigenvectors are not automatically
classified in terms of the irreducible representations of C3v.

A second approach is to fix an irreducible representation $D$ of C3v and
expand the superconducting order parameter in terms of basis functions that
are orthogonal with respect to the Fermi surface inner product,

$\Delta^{D}({{\boldsymbol{k}}})=\sum_{i}c^{D}_{i}\psi^{D}_{i}({{\boldsymbol{k}}})\quad\text{with}\quad\left(\frac{1}{2\Omega}\right)^{2}\sum_{\begin{subarray}{c}{{\boldsymbol{k}}},{{\boldsymbol{k}}}^{\prime}\\\
|\xi_{\uparrow}({{\boldsymbol{k}}})-\mu|<\Omega\\\
|\xi_{\uparrow}({{\boldsymbol{k}}}^{\prime})-\mu|<\Omega\end{subarray}}\psi_{i}^{D}({{\boldsymbol{k}}})^{*}\psi_{j}^{D}({{\boldsymbol{k}}}^{\prime})\,\delta{{\boldsymbol{k}}}\,\delta{{\boldsymbol{k}}}^{\prime}=\delta_{ij}.$
(28)

This expansion allows us to write the self-consistency equation in the matrix
form,

$\boldsymbol{A}^{D}\boldsymbol{c}^{D}=\lambda^{D}\boldsymbol{c}^{D}\quad\text{with}\quad[\boldsymbol{A}^{D}]_{ij}=\left(\frac{1}{2\Omega}\right)^{2}\sum_{\begin{subarray}{c}{{\boldsymbol{k}}},{{\boldsymbol{k}}}^{\prime}\\\
|\xi_{\uparrow}({{\boldsymbol{k}}})-\mu|<\Omega\\\
|\xi_{\uparrow}({{\boldsymbol{k}}}^{\prime})-\mu|<\Omega\end{subarray}}\psi_{i}^{D}({{\boldsymbol{k}}})^{*}\Gamma(\hat{{\boldsymbol{k}}},\hat{{\boldsymbol{k}}}^{\prime})\psi_{j}^{D}({{\boldsymbol{k}}}^{\prime})\,\delta{{\boldsymbol{k}}}\,\delta{{\boldsymbol{k}}}^{\prime},$
(29)

so that we can solve for $\boldsymbol{c}^{D}=(c_{1},c_{2},...)^{T}$ to obtain
the eigenvalues $\lambda^{D}$ and the eigenvectors
$\Delta^{D}({{\boldsymbol{k}}})$.

Finally, we outline a two-step protocol bib:Platt2013  to generate the
orthogonal basis functions for a $d$-dimensional representation $D$: As a
first step, we fix a set of linear-independent vectors
$\\{{\boldsymbol{\delta}}_{n,i}\\}_{i=1,..,d}$ that connect a site on the
superlattice with an $n$-th nearest neighbor site. This set of vectors allows
us to define the lattice harmonics,

$\tilde{\psi}^{D}_{n,i}({{\boldsymbol{k}}})=\sum_{G}\text{char}_{D}([G])\,e^{i{{\boldsymbol{k}}}\cdot
G{\boldsymbol{\delta}}_{n,i}},$ (30)

where the summation is over all the group elements $G$ of $C_{3v}$ and
$\text{char}_{D}([G])$ denotes the character of the conjugacy class $[G]$ in
the representation $D$. In a second step, we use a Gram-Schmidt procedure to
orthogonalize the lattice harmonics, leading to a desired set of orthogonal
basis functions $\psi_{i}^{D}({{\boldsymbol{k}}})$. This concludes our
discussion on the numerical evaluation of the self-consistency equation.

##  IV Ginzburg-Landau free energy

In this fourth section of the Supplemental Material, we present details on the
analysis of the phenomenological Ginzburg-Landau free energy, which allows to
us to identify the candidate pairing states that is realized in the
$E$-representation of our model for superconductivity in tWSe2.

In a first step, we introduce a general linear combination of the
superconducting order parameter,

$\Delta({{\boldsymbol{k}}})=\eta_{1}\Delta_{1}({{\boldsymbol{k}}})+\eta_{2}\Delta_{2}({{\boldsymbol{k}}}),$
(31)

where $\\{\Delta_{1}({{\boldsymbol{k}}}),\Delta_{2}({{\boldsymbol{k}}})\\}$
are two linear-independent basis functions of the $E$-representation and
$\\{\eta_{1},\eta_{2}\\}$ are complex expansion coefficients.

In a second step, we expand the Ginzburg-Landau free energy in terms of
$\boldsymbol{\eta}=(\eta_{1},\eta_{2})$,

$F_{\text{GL}}=a\,\boldsymbol{\eta}\cdot\boldsymbol{\eta}^{*}+b_{1}(\boldsymbol{\eta}\cdot\boldsymbol{\eta}^{*})^{2}+b_{2}|\boldsymbol{\eta}\cdot\boldsymbol{\eta}|^{2}.$
(32)

Here, we have required that the Ginzburg-Landau free energy is (1) a real
function, (2) invariant under all symmetry operations of $C_{3v}$, (3) an even
function in $\eta_{1}$ and $\eta_{2}$. Moreover, we have denoted the Ginzburg-
Landau expansion coefficients by $\\{a,b_{1},b_{2}\\}$. We point out that
$a>0$ induces the superconducting phase if $T<T_{c}$. Also, we note that
$|\boldsymbol{\eta}^{*}\times\boldsymbol{\eta}|$ is not an independent quartic
term because
$|\boldsymbol{\eta}\cdot\boldsymbol{\eta}|^{2}=(\boldsymbol{\eta}\cdot\boldsymbol{\eta}^{*})^{2}-|\boldsymbol{\eta}^{*}\times\boldsymbol{\eta}|$.

In a third step, we perform a variation of the Ginzburg-Landau free energy
with respect to $\boldsymbol{\eta}^{*}$. This yields the requirement,

$\begin{split}\frac{\delta
F_{\text{GL}}}{\delta\boldsymbol{\eta}^{*}}=a\,\boldsymbol{\eta}+2b_{1}\left(\boldsymbol{\eta}\cdot\boldsymbol{\eta}^{*}\right)\boldsymbol{\eta}+2b_{2}\left(\boldsymbol{\eta}\cdot\boldsymbol{\eta}\right)\boldsymbol{\eta}^{*}\stackrel{{\scriptstyle!}}{{=}}0.\end{split}$
(33)

To ensure thermodynamic stability, we further require that
$b_{1}\stackrel{{\scriptstyle!}}{{>}}0$.

In a fourth step, we note that if $b_{2}>0$, the superconducting order
parameter needs to satisfy $\boldsymbol{\eta}\cdot\boldsymbol{\eta}=0$ for
minimizing the Ginzburg-Landau free energy. From this constraint, we see that
the superconducting order parameter needs to take on the chiral form,

$\Delta({{\boldsymbol{k}}})\propto\Delta_{1}({{\boldsymbol{k}}})\pm
i\Delta_{2}({{\boldsymbol{k}}}).$ (34)

In a final step, we note that if $b_{2}<0$, the superconducting order
parameter needs to maximize $\boldsymbol{\eta}\cdot\boldsymbol{\eta}$ for
minimizing the Ginzburg-Landau free energy. This condition implies that the
superconducting order parameter needs to take on the nematic form,

$\Delta({{\boldsymbol{k}}})\propto\cos\varphi\,\Delta_{1}({{\boldsymbol{k}}})+\sin\varphi\,\Delta_{2}({{\boldsymbol{k}}}),$
(35)

where $\varphi$ is an arbitrary phase. This concludes our derivation of the
form of the superconducting order parameter.

##  V Chern number

In this sixth section of the Supplemental Material, we present the details on
the computation of the Chern number for the Bogoliubov-de Gennes Hamiltonian
of superconducting tWSe2.

As a first step, we write the Bogoliubov-de Gennes Hamiltonian in the form,

$H=\frac{1}{2}\sum_{{{\boldsymbol{k}}}}\begin{pmatrix}c^{{\dagger}}_{{{\boldsymbol{k}}},\uparrow}&c^{{\dagger}}_{{{\boldsymbol{k}}},\downarrow}&c_{-{{\boldsymbol{k}}},\uparrow}&c_{-{{\boldsymbol{k}}},\downarrow}\end{pmatrix}\begin{pmatrix}\xi_{\uparrow}({{\boldsymbol{k}}})&0&0&-\Delta({{\boldsymbol{k}}})^{*}\\\
0&\xi_{\downarrow}({{\boldsymbol{k}}})&\Delta(-{{\boldsymbol{k}}})^{*}&0\\\
0&\Delta(-{{\boldsymbol{k}}})&-\xi_{\downarrow}({{\boldsymbol{k}}})&0\\\
-\Delta({{\boldsymbol{k}}})&0&0&-\xi_{\uparrow}({{\boldsymbol{k}}})\end{pmatrix}\begin{pmatrix}c_{{{\boldsymbol{k}}},\uparrow}\\\
c_{{{\boldsymbol{k}}},\downarrow}\\\
c^{{\dagger}}_{-{{\boldsymbol{k}}},\uparrow}\\\
c^{{\dagger}}_{-{{\boldsymbol{k}}},\downarrow}\end{pmatrix},$ (36)

and assume a superconducting order parameter of the chiral form,
$\Delta({{\boldsymbol{k}}})=\Delta_{1}({{\boldsymbol{k}}})+i\Delta_{2}({{\boldsymbol{k}}})$.

As a second step, we define the Chern number as bib:Yuan2014 ,

$\mathcal{C}=\frac{1}{2\pi}\sum_{E_{s}({{\boldsymbol{k}}})<0}\int\mathrm{d}{{\boldsymbol{k}}}\,\\{\partial_{k_{x}}a^{s}_{y}({{\boldsymbol{k}}})-\partial_{k_{y}}a^{s}_{x}({{\boldsymbol{k}}})\\}.$
(37)

In this expression, the filled Bogoliubov bands are given by,

$E_{\uparrow}({{\boldsymbol{k}}})=-\sqrt{\xi_{\uparrow}({{\boldsymbol{k}}})^{2}+|\Delta({{\boldsymbol{k}}})|^{2}}\quad\text{with}\quad
E_{\downarrow}({{\boldsymbol{k}}})=-\sqrt{\xi_{\downarrow}({{\boldsymbol{k}}})^{2}+|\Delta(-{{\boldsymbol{k}}})|^{2}}.$
(38)

The corresponding eigenstates $|\uparrow,{{\boldsymbol{k}}}\rangle$ and
$|\downarrow,{{\boldsymbol{k}}}\rangle$ of the Hamiltonian density then
define,

$a^{s}_{i}({{\boldsymbol{k}}})=-i\langle
s,{{\boldsymbol{k}}}|\partial_{k_{i}}|s,{{\boldsymbol{k}}}\rangle.$ (39)

Finally, we note that for the numerical evaluation of the Chern number, we use
the method described in bib:Fukui2005 . This concludes our discussion on the
Chern numbers for the Bogoliubov bandstructure of superconducting tWSe2.

Figure 6: (Color online) (a) Triangular lattice with directed bonds. Electrons
of the $+K$ valley, pick up a tunneling phase $e^{i\pi/6}$ when hopping along
the bond direction. Electrons of the $-K$ valley pick up the opposite phase.
(b) Dispersion $\tilde{\xi}_{\uparrow}({{\boldsymbol{k}}})/2\tilde{t}$ along
high-symmetry lines of the moiré Brillouin zone for the effective
‘$\pi/6$’-model. (c) Density of states of the band structure of the effective
‘$\pi/6$’-model.

Figure 7: (Color online) (a) Dispersion
$\tilde{\xi}_{\uparrow}({\boldsymbol{q}})$ of the top-most band at the $+K$
valleys. The $+K$ ($-K$) valley Fermi surfaces are shown with solid (dashed)
black lines. (b) Transversal spin-valley susceptibility
$\tilde{\chi}_{+-}({\boldsymbol{q}})$ at $1/\beta=100\,$meV. (c) Longitudinal
spin-valley susceptibility $\tilde{\chi}_{zz}({\boldsymbol{q}})$ at
$1/\beta=100\,$meV.

##  VI Effective ‘$\pi/6$’-model

In this seventh section of the Supplemental Material, we introduce and discuss
an effective model that captures some of the features of the full continuum
model of tWSe2 for the parameter choices given in the main text.

###  VI.1 Hamiltonian

As a starting point, we introduce the Hamiltonian of the effective model,

$\tilde{H}=\sum_{{{\boldsymbol{k}}},s}[\tilde{\xi}_{s}({{\boldsymbol{k}}})-\mu]c^{{\dagger}}_{{{\boldsymbol{k}}},s}c_{{{\boldsymbol{k}}},s}+\frac{U}{N}\sum_{\begin{subarray}{c}{{\boldsymbol{k}}},{\boldsymbol{p}},{\boldsymbol{q}}\\\
s,s^{\prime}\end{subarray}}c^{{\dagger}}_{{\boldsymbol{q}},s}c^{{\dagger}}_{{\boldsymbol{p}}+{{\boldsymbol{k}}}-{\boldsymbol{q}},s^{\prime}}c_{{\boldsymbol{p}},s^{\prime}}c_{{{\boldsymbol{k}}},s}$
(40)

In this expression, we have defined the dispersion relations, which describe
(as shown in Fig. 6(a)) electron hopping on a triangular lattice,

$\tilde{\xi}_{s}({{\boldsymbol{k}}})=-2\tilde{t}\sum_{i=1,2,3}\cos\left(\frac{\pi}{6}\right)\cos({{\boldsymbol{k}}}\cdot{\boldsymbol{\delta}}_{i})+s\sin\left(\frac{\pi}{6}\right)\sin({{\boldsymbol{k}}}\cdot{\boldsymbol{\delta}}_{i}),$
(41)

Here, the lattice vectors are given by
${\boldsymbol{\delta}}_{1}=a_{0}/|\theta|(0,1)^{T}$,
${\boldsymbol{\delta}}_{2}=a_{0}/|\theta|(\sqrt{3}/2,-1/2)^{T}$, and
${\boldsymbol{\delta}}_{3}=a_{0}/|\theta|(-\sqrt{3}/2,-1/2)^{T}$. Moreover,
the tunneling amplitude is denoted by $\tilde{t}$. As shown in Fig. 6(b), the
dispersion relation $\tilde{\xi}_{s}({{\boldsymbol{k}}})$ is qualitatively
similar to the dispersion relation $\xi_{s}({{\boldsymbol{k}}})$ obtained from
the full continuum model. In particular, the dispersion relation shows a
saddle points at $\kappa_{\pm}$ that lead to a vHs at half-filling, as
depicted in Fig. 6(c).

###  VI.2 Comparing the effective model to the continuum model

Next, want to compare the effective model to the continuum model of . For that
purpose

we have computed the valley susceptibility from the expression,

$\tilde{\chi}_{\alpha\beta}({\boldsymbol{q}})=-\frac{1}{\beta
N}\sum_{{{\boldsymbol{k}}},i\omega}\text{Tr}\left[\sigma_{\alpha}\tilde{G}({{\boldsymbol{k}}},i\omega)\sigma_{\beta}\tilde{G}({{\boldsymbol{k}}}+{\boldsymbol{q}},i\omega)\right]\quad\text{with}\quad\tilde{G}_{s_{1}s_{2}}({{\boldsymbol{k}}},i\nu)=\delta_{s_{1}s_{2}}[i\nu-\tilde{\xi}_{s_{1}}({{\boldsymbol{k}}})+\mu]^{-1}.$
(42)

The results of the computation are shown in Fig. 7. Most notably, we have
found that the transversal susceptibility
$\tilde{\chi}_{+-}({\boldsymbol{q}})$ has a qualitatively similar structure as
in the full continuum model. In particular, at half-filling, the transversal
susceptibility shows the same characteristics maxima at the wavevectors
${{\boldsymbol{Q}}}_{j}=R_{2\pi(j-1)/3}{\boldsymbol{\kappa}}_{+}$.

Finally, we compare the possible superconducting orders arising from inter-
valley fluctuations in the ‘$\pi/6$’-model to the results of the full
continuum model. As shown in Fig. 8 and 9, the most negative eigenvalue
$\tilde{\lambda}^{(0)}$ of the Cooper channel interaction vertex
$\tilde{\Gamma}(\hat{{\boldsymbol{k}}},\hat{{\boldsymbol{k}}}^{\prime})=U+U^{2}\tilde{\chi}_{+-}(\hat{{\boldsymbol{k}}}+\hat{{\boldsymbol{k}}}^{\prime})$
is in the $E$-representation of the $C_{3v}$ point group at half-filling,
which agrees with the continuum model findings. However, we find that the
leading superconducting order in the ‘$\pi/6$’-model and the full continuum
model differ at quarter filling. Specifically, we find that the leading
superconducting instability at quarter-filling is not in the $A_{2}$
representation as in the full continuum model, but rather in the
$E$-representation.

Figure 8: (Color online) Pairing strengths of the effective ‘$\pi/6$’-model
given by the most negative eigenvalue $|\tilde{\lambda}^{(0)}|$ of the the
self-consistency equation (within a fixed representation) as a function of the
filling $n$. For fillings $0.27\lesssim n\lesssim 0.29$ and $0.71\lesssim
n\lesssim 0.73$, we find that the most negative eigenvalue is in the one-
dimensional $A_{2}$ representation of the $C_{3v}$ point group. At all other
fillings, we find that the most negative eigenvalue is in the two-dimensional
$E$ representation of $C_{3v}$.

Figure 9: (Color online) Gap functions of the leading superconducting states
of the effective ‘$\pi/6$’-model (normalized by their maximum value). At
$n\approx 1/4,\,1/2\,3/4$, the most negative eigenvalue is two-fold degenerate
and the gap basis functions $\tilde{\Delta}_{1,2}({{\boldsymbol{k}}})$
transform in the $2d$ $E$-representation of the $C_{3v}$ point group.

## References

  * (1) R. Shankar, Reviews of Modern Physics 66, 129 (1994).
  * (2) L. Wang and O. Vafek, Physica C 497, 6 (2014).
  * (3) C. Platt, W. Hanke, and R. Thomale, Advances in Physics 62, 453 (2013).
  * (4) N. F. Q. Yuan, K. F. Mak, and K. T. Law, Phys. Rev. Lett. 113, 097001 (2014).
  * (5) T. Fukui, Y. Hatsugai, and H. Suzuki, J. Phys. Soc. Jpn. 74, 1674 (2005).


# Modularity benefits reinforcement learning agents with competing homeostatic
drives

Zack Dulberg
Princeton University
Princeton, NJ 08544
<EMAIL_ADDRESS>
&Rachit Dubey
Princeton University
Princeton, NJ 08544
<EMAIL_ADDRESS>
&Isabel M. Berwian
Princeton University
Princeton, NJ 08544
<EMAIL_ADDRESS>
&Jonathan D. Cohen
Princeton University
Princeton, NJ 08544
<EMAIL_ADDRESS>

###### Abstract

The problem of balancing conflicting needs is fundamental to intelligence.
Standard reinforcement learning algorithms maximize a scalar reward, which
requires combining different objective-specific rewards into a single number.
Alternatively, different objectives could also be combined at the level of
_action_ value, such that specialist modules responsible for different
objectives submit different action suggestions to a decision process, each
based on rewards that are independent of one another. In this work, we explore
the potential benefits of this alternative strategy. We investigate a
biologically relevant multi-objective problem, the continual homeostasis of a
set of variables, and compare a monolithic deep Q-network to a modular network
with a dedicated Q-learner for each variable. We find that the modular agent:
a) requires minimal exogenously determined exploration; b) has improved sample
efficiency; and c) is more robust to out-of-domain perturbation.

Keywords:

modular reinforcement learning, homeostasis, conflict, multi-objective
decision-making, exploration

#### Acknowledgements

This project / publication was made possible through the support of a grant
from the John Templeton Foundation.

##  1 Introduction

Humans (and other animals) must satisfy a large set of distinct and possibly
conflicting objectives. For example, we must find food, water, shelter,
socialize, maintain our temperature, reproduce, etc.. Artificial agents that
need to function autonomously in natural environments may face similar
problems (e.g., balancing the need to accomplish a specified goal with the
need to recharge, etc.). Finding a way to balance disparate needs is thus an
important challenge for intelligent agents, and can often be a source of
psychological conflict in humans.

In standard reinforcement learning (RL), monolithic agents act in order to
maximize a single future discounted reward [1]. The standard way to generalize
this to problems with multiple objectives is scalarization. For example, in
homeostatically-regulated reinforcement learning (HRRL), an agent is modelled
as having separable homeostatic drives, and is rewarded based on its ability
to maintain all its ”homeostats” at their set points. This is done by
_combining_ deviations from _all_ set-points into a single reward which, when
maximized, minimizes homeostatic deviations overall [2].

This approach faces several challenges typical to RL. First, to avoid settling
on a sub-optimal policy, an agent must trade-off exploitation of knowledge
about its primary objective with some form of exogenous exploration (typically
by acting randomly or according to an exploration-specific bonus). Second,
sample inefficiency follows from the “curse of dimensionality”: as
environmental complexity increases, an agent must learn how exponentially more
states relate to its objective. Third, RL agents tend to over-fit their
environment, performing poorly out-of-domain (i.e. they are not robust to
distribution shifts). In the broader context of balancing multiple objectives,
reward scalarization might be undesirable if the relative importance of
different objectives is unknown or variable [3]. Finally, is not clear that
the brain itself uses a common currency to navigate such trade-offs [4].

What is the alternative? Given that objectives may conflict with each other
due to environmental constraints, and that agents only have one body with
which to act, conflict must be resolved at some point between affordance and
action. The monolithic solution resolves conflict at the level of reward (i.e.
close to affordance). We suggest resolution could occur later; a set of
modules with separate reward functions could submit action values to a
decision process that selects a final action [5, 6]. This approach has the
potential to address the three aforementioned challenges. Exploration might
emerge naturally as a property of the system rather than having to be imposed
or regulated as a separate factor, as specialist modules are “dragged along”
by other modules when those have the “upper hand” on action. Modules might
also have smaller sub-sets of relevant features to learn about, improving
sample efficiency, and be less sensitive to distribution shifts in irrelevant
features, improving robustness.

Here, we report simulations that provide evidence for benefits of such a
modular approach with respect to exploration, sample efficiency, and
robustness using deep RL in the context of homeostatic objectives. We
construct a simple but flexible environment of homeostatic tasks and construct
a deep RL implementation of the HRRL reward function. We then use this
framework to quantify differences between monolithic and modular deep
Q-agents, finding that modular agents seem to explore well on their own,
achieve homeostasis faster, and better maintain it after distribution shift.
Together, these results highlight the potential learning benefits of modular
RL, while at the same time offering a framework through which psychological
conflict and resolution might be better understood.

##  2 Methods

###  2.1 Environment

To study conflicting needs, we constructed a toy grid-world environment
containing multiple different resources. Specifically, each location $(x,y)$
in the environment contained a vector of resources of length $N$ (i.e., there
were $N$ overlaid resource maps). The spatial distribution of each individual
resource was specified by a normalized 2D Gaussian with mean $\mu_{x},\mu_{y}$
and co-variance matrix $\Sigma$ (see also Figure 1(a)).

The agent received as perceptual input a 3x3 egocentric slice of the $N$
resource maps (i.e. it could see all the resource levels at each position in
its local vicinity). In addition to the resource landscape, the agent also
perceived a vector $H_{t}=(h_{1,t},h_{2,t},...,h_{N,t})$ consisting of $N$
internal variables with each representing the agent’s homeostatic need with
respect to the resource. We refer to these variables as ”internal stats” or
just ”stats” (such as osmostat, glucostat, etc.) which we assume are
independent ($h_{i}$ is only affected by acquisition of resource $i$) and have
some desired set-point $h_{i}^{*}$ (see Figure 1(b)). Set-points were fixed at
$H^{*}=(h_{1}^{*},h_{2}^{*},...,h_{N}^{*})$ and did not change over the course
of training.

The agent could move in each of four cardinal directions and, with each step,
the individual stats, $h_{i}$, increased by the amount of resource $i$ at the
agent’s next location. Additionally, each internal stat decayed at a constant
rate to represent the natural depletion of internal resources over time (note:
resources in the environment themselves did not deplete). Thus, if the agent
discovered a location with a high level of resource for a single depleted
stat, staying at that location would optimize that stat toward its set-point,
however others would progressively deplete. Agents were initialized in the
center of the grid, with internal stats below their set-points, and were
trained for $30,000$ steps for a single episode (i.e. agents had to learn in
real time as internal stats depleted). For all experiments, the internal stats
started at the same level and shared the same set-points. Environmental
parameters are summarized in Table 1.

(a) Gridworld environment

(b) Overview of monolithic DQN and modular GmQ agents

Figure 1: Environment and model schematics

###  2.2 Models

##### Monolithic agent

We created a monolithic agent based on the deep Q network (DQN) [7]. The
agent’s perceptual input was a concatenation of all local resource levels
along with all internal stat levels at each time step (we used neural networks
as function approximators since stats were continuous variables). It’s output
was 4 action logits subsequently used for $\epsilon$-greedy action selection.
We used the HRRL reward function [2] which defined reward at each time-step
$r_{t}$ as drive reduction, where drive $D$ was a convex function of set-point
deviations; see equation (1).

$\vspace{-0.2cm}r_{t}=D(H_{t})-D(H_{t+1})\quad\text{where}\quad
D(H_{t})=\sqrt[m]{\sum_{i=1}^{N}|h_{i}^{*}-h_{i,t}|^{n}}$ (1)

##### Modular agent

We created a modular agent based on greatest-mass Q-learning (GmQ) [5], which
consisted of a separate DQN for each of the 4 resources/stats. Here, each
module had the same input as the monolithic model (i.e. the full egocentric
view and all 4 stat levels), but received a separate reward $r_{i,t}$ derived
from only a single stat. The reward function for the $i$th module was
therefore defined as in equation (2), where drive $D$ depended on the $i$th
resource only.

$\vspace{-0.2cm}r_{i,t}=D(h_{i,t})-D(h_{i,t+1})\quad\text{where}\quad
D(h_{i,t})=\sqrt[m]{|h_{i}^{*}-h_{i,t}|^{n}}$ (2)

To select a single action from the suggestions of the multiple modules, we
used a simple additive heuristic. We first summed Q-values for each action
across modules, and then performed standard $\epsilon$-greedy action selection
on the result. More specifically, if $Q_{i}(a)$ was the Q-value of action $a$
suggested by module $i$, greedy actions were selected as
$\underset{a}{\arg\max}\sum\limits_{i}Q_{i}(a)$.

##### Common features

Schematics for both models are shown in Figure 1(b). All Q-networks were
multi-layered perceptrons (MLP) with rectified linear nonlinearities trained
using a standard temporal difference loss function with experience replay and
target networks [7]. The Adam optimizer was used to perform one gradient
update on each step in the environment. For both models, $\epsilon$ was
annealed linearly from its initial to final value at the beginning of training
at a rate that was experimentally manipulated as described below.
Hyperparameters are summarized in Table 1.

Table 1: Parameter settings for environment and models Model | DQN | GmQ | Environment |
---|---|---|---|---
Trainable parameters | 1.09e6 | 1.09e6 | # of resources $N$ | 4
MLP hidden layer units | 1024 | 500 | HRRL exponents $(n,m)$ | $(4,2)$
Learning rate | 1e-3 | 1e-3 | Stat set-points $H^{*}$ | $(5,5,5,5)$
Discount factor $\gamma$ | 0.5 | 0.5 | Initial stat levels $H_{t=0}$ | $(0.5,0.5,0.5,0.5)$
Memory buffer capacity | 30k | 30k | Resource locations $\mu_{x},\mu_{y}$ | {0,10},{0,10}
Target network update frequency | 200 | 200 | Resource covariance $\Sigma$ | $\big{(}\begin{smallmatrix}1&0\\\ 0&1\end{smallmatrix}\big{)}$
Batch size | 512 | 512 | Stat depletion per step | 0.004
Initial $\epsilon$ | 1 | 1 |  |
Final $\epsilon$ | 0.01 | 0.01 |  |

##  3 Results

###  3.1 Optimizing monolithic DQN for homeostasis

We first characterized whether a standard DQN could reliably perform the task
of homeostasis in our environment. Figure 2(a) summarizes the mean internal
stat levels of 10 models averaged over all 4 stats and over the final 1k steps
of training for different desired set-points. The model reliably achieved each
set-point by the end of training, indicated by points tracking the identity
line. The slight under-shooting (i.e. points slightly below the diagonal line)
may reflect a bias from stats being initialized far below their set-points.

We then fixed the set-points $h_{i}^{*}=5$ for all stats, and optimized the
performance of the DQN baseline by performing a search over performance-
relevant hyper-parameters, such as the discount factor $\gamma$. To quantify
performance, we calculated the average homeostatic deviation per step $\Delta$
after the exploration annealing phase (using $t_{1}=15k$ and $t_{2}=30k$) as
in equation (3). Lower $\Delta$ indicates better homeostatic performance.

$\vspace{-0.2cm}\Delta=\frac{\sum\limits_{t=t_{1}}^{t_{2}}\sum\limits_{i}|h_{i}^{*}-h_{i,t}|}{t_{2}-t_{1}}$
(3)

Baseline DQN performance over a range of discount factors $\gamma$ is shown in
Figure 2(b). We selected the best performing setting of $\gamma=0.5$ and
matched this and other parameters (see Table 1) between DQN and GmQ to compare
them in the following two head-to-head experiments.

(a) Homeostasis is achieved for a range of set-points ($n=10$)

(b) Performance over a range of discount factors ( $n=50$)

Figure 2: DQN baseline learnability and performance; Boxplots display inter-
quartile range and outliers for $n$ models

###  3.2 Modularity provides an exploration benefit

To investigate the impact of modularity on the need for exploration, we
systematically varied the number of steps used to anneal $\epsilon$ in
$\epsilon$-greedy exploration from its initial to final value. We varied the
$\epsilon$ annealing time for both models from 1 (i.e. minimal exploration of
$\epsilon=0.01$ only) to 10k (i.e. annealing from $\epsilon=1$ to
$\epsilon=0.01$ over 10k steps).

Figure 3(a) shows the results of varying the amount of exploration annealing
for both models. While DQN gains incremental performance benefits from
increasing periods of initial exploration, GmQ displays a striking
indifference to the exploration period; with only 1 step of annealing, it
performs as well or better than the best DQN models. In other words, DQN
requires careful tuning of an appropriate exploration annealing period, but
GmQ achieves good performance with effectively no exogenously specified
exploration.

###  3.3 Modularity provides robustness in the face of perturbation

Finally, we tested how robust each model was to a perturbation out-of-domain
that occurred halfway through training, i.e. at time-step 15k. At that point,
a single internal stat variable (i.e. $h_{4}$) was clamped to a value of 20 (a
value previously unseen by the network), and did not change (thus contributing
no drive reduction and therefore 0 reward). We tested how well homeostasis was
maintained for remaining stats after this perturbation.

Figure 3(b) shows the time-course of the four stats from the beginning of
training, through a perturbation at time-step 15k, using $5000$
$\epsilon$-annealing steps. First, it can be seen that GmQ achieves stable
homeostasis first, whereas DQN over-shoots set-points initially and takes
longer to stabilize. Second, when stat 4 is clamped, DQN displays a
significant disturbance to homeostasis of remaining stats, without clear
recovery, whereas GmQ is robust in the face of the perturbation.

(a) Varying exploration period ($n=50$)

(b) Stat time-courses ($n=50$) perturbed at $t=15k$ by clamping $h_{4}$. Green
line shows set-point. Shading reflects s.d. accross models.

Figure 3: Experiments comparing DQN and GmQ with respect to (a) exploration
and (b) perturbation

##  4 Discussion

We have shown that in a grid-world task with competing homeostatic drives, a
simple modular agent based on greatest-mass Q-learning (GmQ) requires less
hand-coded exploration, learns faster, and is more robust to environmental
perturbations compared to a traditional monolithic deep Q-network (DQN). Our
findings in the context of competing drives complement work showing mixture of
expert systems display improved sample efficiency and generalization [8]. We
also believe the exploration benefits we observed are novel. The problem of
exploration in RL is fundamental, and existing solutions make use of noise,
explicit exploratory drives/bonuses, and/or other forms of auto-annealing [9,
10]. We suggest an additional class of strategies, namely, exploration as an
added benefit of having multiple independent drives, since exploitation from
the perspective of one module is exploration from the perspective of another.
We hypothesize that the ability of modules to suggest conflicting actions may
provide modular agents with an implicit source of exploration.

##### Future Work

Our toy environment has highlighted some initial benefits of modularity
(exploration, sample efficiency and robustness), and we predict that these
advantages will be amplified in more complex environments, or as the number of
drives/objectives increases, due to the curse of dimensionality (more states
to explore, learn about, or perturb). We aim to test our agents in rich 3D
environments with homeostatic objectives. Next, modular drives immediately
pose the problem of coordination. While we simply summed Q-values, more
complex arbitrators (such as an additional RL agent that dynamically re-
weights individual drives) might better exploit the benefits of modularity in
the context of multiple objectives. Finally, humans experience psychological
conflict, with various resolution mechanisms long described by psychodynamic
theories [11]. Modular RL, with its implicit conflicts and resolutions, could,
for the first time, offer a formal, computationally-explicit, and normative
explanatory framework that could undergird and/or replace elements of
psychodynamic theory for understanding the mechanisms responsible for conflict
and resolution in the human brain.

## References

  * [1] Richard S Sutton and Andrew G Barto  “Reinforcement learning: An introduction”  MIT press, 2018
  * [2] Mehdi Keramati and Boris Gutkin  “Homeostatic reinforcement learning for integrating reward collection and physiological stability”  In _Elife_ 3 eLife Sciences Publications Limited, 2014, pp. e04811
  * [3] Diederik M Roijers, Peter Vamplew, Shimon Whiteson and Richard Dazeley  “A survey of multi-objective sequential decision-making”  In _Journal of Artificial Intelligence Research_ 48, 2013, pp. 67–113
  * [4] Benjamin Y Hayden and Yael Niv  “The case against economic values in the orbitofrontal cortex (or anywhere else in the brain).”  In _Behavioral Neuroscience_ 135.2 American Psychological Association, 2021, pp. 192
  * [5] Stuart J Russell and Andrew Zimdars  “Q-decomposition for reinforcement learning agents”  In _Proceedings of the 20th International Conference on Machine Learning (ICML-03)_ , 2003, pp. 656–663
  * [6] Harm Van Seijen et al.  “Hybrid reward architecture for reinforcement learning”  In _Advances in Neural Information Processing Systems_ 30, 2017
  * [7] Volodymyr Mnih et al.  “Playing atari with deep reinforcement learning”  In _arXiv preprint arXiv:1312.5602_ , 2013
  * [8] Robert A Jacobs, Michael I Jordan, Steven J Nowlan and Geoffrey E Hinton  “Adaptive mixtures of local experts”  In _Neural computation_ 3.1 MIT Press, 1991, pp. 79–87
  * [9] Tianpei Yang et al.  “Exploration in deep reinforcement learning: a comprehensive survey”  In _arXiv:2109.06668_ , 2021
  * [10] Samuel M McClure, Mark S Gilzenrat and Jonathan D Cohen  “An exploration-exploitation model based on norepinepherine and dopamine activity”  In _Advances in neural information processing systems_ 18, 2005
  * [11] Anna Freud  “The ego and the mechanisms of defence”  Routledge, 2018