8\frac{\sin^{2}\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(1-\frac{\sin^{2}\delta_{f}}{\pi\widetilde{\Delta}_{0}}(-\widetilde{U}_{1})\right)\
.$ (129)

###  B.4 Ward identities and exact susceptibilities

We consider a perturbation

$\Delta S=\sum_{1}f^{\dagger}(1)V(1)f(1)$ (130)

where $V(1)=V_{\alpha_{1}\eta_{1}s_{1}}$ is a diagonal matrix and is
independent of frequency. Then the first order correction to the Green’s
function can be written as

$\displaystyle\delta\mathcal{G}(1)=$ $\displaystyle-\langle
f(1)f^{\dagger}(1)(-\Delta S)e^{-S_{I}}\rangle_{0C}$ $\displaystyle=$
$\displaystyle\sum_{2}\wick{\c{1}f(1)\c{2}f^{\dagger}(1)\c{1}f^{\dagger}(2)V(2)\c{2}f(2)}-4\times\frac{T}{4}\sum_{23456}\wick{\c{1}f(1)\c{2}f^{\dagger}(1)\c{1}f^{\dagger}(3)\c{3}f^{\dagger}(4)\c{4}f(5)\c{2}f(6)\c{4}f^{\dagger}(2)\c{3}f(2)}\cdot\Gamma^{0}(3,4;5,6)\cdot
V(2)+\cdots$ $\displaystyle=$ $\displaystyle[\mathcal{G}^{0}(1)]^{2}\cdot
V(1)+T\sum_{2}[\mathcal{G}^{0}(1)]^{2}\cdot\Gamma^{0}(1,2;2,1)\cdot[\mathcal{G}^{0}(2)]^{2}\cdot
V(2)+\cdots$ (131)

The subscript $0C$ in the first row represents the connected diagrams, and the
factor $4$ in the second term of the second row comes from four equivalent
contractions. Summing over all the terms to infinite order, we obtain the
exact result

$\delta\mathcal{G}(1)=[\mathcal{G}(1)]^{2}\cdot
V(1)+T\sum_{2}[\mathcal{G}(1)]^{2}\cdot\Gamma(1,2;2,1)\cdot[\mathcal{G}(2)]^{2}\cdot
V(2)\ ,$ (132)

which can be represented by the skeleton diagram

(133)

Here bare Green’s functions and vertex function are replaced by the full
Green’s functions and vertex function respectively. Writing
$\delta\mathcal{G}(1)=-[\mathcal{G}(1)]^{2}\cdot(-V(1)-\delta\Sigma(1))$, we
obtain the perturbed self-energy

$\delta\Sigma_{I_{1}}(\omega_{1})=\sum_{I_{2}}\int\frac{d\omega_{2}}{2\pi}\cdot\Gamma_{I_{1},I_{2};I_{2},I_{1}}(\omega_{1},\omega_{2};\omega_{2},\omega_{1})\cdot[\mathcal{G}(\omega_{2})]^{2}\cdot
V_{I_{2}}\ .$ (134)

Here $I_{i}=(\alpha_{i},\eta_{i},s_{i})$ is a composite index. We consider the
perturbation term $V_{I}=-h^{O}\cdot O_{I}$ with $O$ being a diagonal matrix,
then there is

$\sum_{I_{1}}\frac{\partial\Sigma_{I_{1}}(\omega_{1})}{\partial
h^{O}}O_{I_{1}}=-\sum_{I_{1}I_{2}}\int\frac{d\omega_{2}}{2\pi}\cdot
O_{I_{1}}O_{I_{2}}\cdot\Gamma_{I_{1},I_{2};I_{2},I_{1}}(\omega_{1},\omega_{2};\omega_{2},\omega_{1})\cdot[\mathcal{G}(\omega_{2})]^{2}\
.$ (135)

We then consider a gauge transformation
$f_{I}^{\prime}(\omega)=f_{I}(\omega+\nu\cdot O_{I})$. On the one hand, there
must be

$\displaystyle\mathcal{G}_{I_{1}}^{\prime}(\omega_{1})-\mathcal{G}(\omega_{1})=\mathcal{G}(\omega_{1}+\nu\cdot
O_{I_{1}})-\mathcal{G}(\omega_{1})=-[\mathcal{G}(\omega_{1})]^{2}\cdot\frac{\partial}{\partial\omega_{1}}\left([\mathcal{G}^{0}(\omega_{1})]^{-1}-\Sigma(\omega_{1})\right)\cdot\nu\cdot
O_{I_{1}}+\mathcal{O}(\nu^{2})$ $\displaystyle=$
$\displaystyle-[\mathcal{G}(\omega_{1})]^{2}\cdot\frac{\partial[\mathcal{G}^{0}(\omega_{1})]^{-1}}{\partial\omega_{1}}\cdot\nu\cdot
O_{I_{1}}+[\mathcal{G}(\omega_{1})]^{2}\frac{\partial\Sigma(\omega)}{\partial
i\omega}\cdot i\nu\cdot O_{I_{1}}+\mathcal{O}(\nu^{2})$ (136)

Here
$\mathcal{G}^{0}(\omega_{1})=(i\omega_{1}+i\Delta_{0}\cdot\mathrm{sgn}(\omega_{1})-\epsilon_{f})^{-1}$
is the free propagator. On the other hand, one can rewrite the action
$S[f^{\dagger},f]$ in terms of $S[f^{\prime\dagger},f^{\prime}]+\Delta
S[f^{\prime\dagger},f^{\prime}]$, with the perturbation

$\Delta S=\sum_{1}f^{\prime\dagger}(1)(-i\nu\cdot O(1)-i\nu\cdot
2\Delta_{0}\delta(\omega_{1})O(1))f^{\prime}(1)=\sum_{1}f^{\prime\dagger}(1)\left(-\nu\cdot
O(1)\partial_{\omega_{1}}[\mathcal{G}^{0}(1)]^{-1}\right)f^{\prime}(1)\ .$
(137)

There is no other perturbation term in $\Delta S$ as long as $O$ is a
generator of the continuous symmetry such that it leaves all the instantaneous
terms ($\omega$-independent terms) invariant. Substituting $V(1)=-\nu\cdot
O(1)\partial_{\omega_{1}}[\mathcal{G}^{0}(1)]^{-1}$ into Eq. 132 and comparing
it to Sec. B.4, one obtains

$\displaystyle\frac{\partial\Sigma(\omega_{1})}{\partial i\omega_{1}}\cdot
i\nu\cdot O_{I_{1}}=$
$\displaystyle-\sum_{I_{2}}\int\frac{d\omega_{2}}{2\pi}\Gamma_{I_{1},I_{2};I_{2},I_{1}}(\omega_{1},\omega_{2};\omega_{2},\omega_{1})\cdot[\mathcal{G}(\omega_{2})]^{2}\cdot\nu\cdot
O_{I_{2}}\cdot\partial_{\omega_{2}}[\mathcal{G}^{0}(\omega_{2})]^{-1}$ (138)

We do not explicitly write
$\partial_{\omega_{2}}[\mathcal{G}^{0}(\omega_{2})]^{-1}$ as
$i+2i\Delta_{0}\delta(\omega_{2})$ because it will lead to the ill-defined
term $[\mathcal{G}(\omega_{2})]^{2}\delta(\omega_{2})$, where
$\mathcal{G}(\omega_{2})$ is discontinuous at $\omega_{2}=0$. This problem
comes from abuse of the $\delta$-function. Here we avoid this problem by
rewriting
$[\mathcal{G}(\omega_{2})]^{2}\partial_{\omega_{2}}[\mathcal{G}^{0}(\omega_{2})]^{-1}$
as

$\displaystyle[\mathcal{G}(\omega_{2})]^{2}\partial_{\omega_{2}}[\mathcal{G}^{0}(\omega_{2})]^{-1}=[\mathcal{G}(\omega_{2})]^{2}\left(\partial_{\omega_{2}}[\mathcal{G}(\omega_{2})]^{-1}+\partial_{\omega_{2}}\Sigma(\omega_{2})\right)=-\partial_{\omega_{2}}\mathcal{G}(\omega_{2})+[\mathcal{G}(\omega_{2})]^{2}\partial_{\omega_{2}}\Sigma(\omega_{2})$
(139)

The derivative $\partial_{\omega_{2}}\mathcal{G}(\omega_{2})$ consists of a
continuous part at $\omega_{2}\neq 0$ and a $\delta$-function peak at
$\omega_{2}=0$ that is responsible for the discontinuity of
$\mathcal{G}(\omega_{2})$, i.e.,

$\partial_{\omega_{2}}\mathcal{G}(\omega_{2})=-[\mathcal{G}(\omega_{2})]^{2}(i-\partial_{\omega_{2}}\Sigma(\omega_{2}))+(\mathcal{G}(0^{+})-\mathcal{G}(0^{-}))\delta(\omega_{2})\
.$ (140)

This expression is correct in the sense that integrating it gives the correct
antiderivative. Then there is

$[\mathcal{G}(\omega_{2})]^{2}\partial_{\omega_{2}}[\mathcal{G}^{0}(\omega_{2})]^{-1}=-(\mathcal{G}(0^{+})-\mathcal{G}(0^{-}))\delta(\omega_{2})+i[\mathcal{G}(\omega_{2})]^{2}\
.$ (141)

Substituting this expression into Eq. 138 yields

$\displaystyle\frac{\partial\Sigma(\omega_{1})}{\partial i\omega_{1}}\cdot
O_{I_{1}}=$ $\displaystyle\frac{\mathcal{G}(0^{+})-\mathcal{G}(0^{-})}{2\pi
i}\sum_{I_{2}}\Gamma_{I_{1},I_{2};I_{2},I_{1}}(\omega_{1},0;0,\omega_{1})\cdot
O_{I_{2}}-\sum_{I_{2}}\int\frac{d\omega_{2}}{2\pi}\Gamma_{I_{1},I_{2};I_{2},I_{1}}(\omega_{1},\omega_{2};\omega_{2},\omega_{1})\cdot[\mathcal{G}(\omega_{2})]^{2}\cdot
O_{I_{2}}$ (142)

According to the Friedel sum rule (Eq. 86), there is

$\frac{\mathcal{G}(0^{+})-\mathcal{G}(0^{-})}{2\pi
i}=-A(0)=-\frac{\sin^{2}(\delta_{f})}{\pi\Delta_{0}}$ (143)

and hence

$\frac{\partial\Sigma(\omega_{1})}{\partial i\omega_{1}}\cdot
O_{I_{1}}=-\frac{\sin^{2}(\delta_{f})}{\pi\Delta_{0}}\sum_{I_{2}}\Gamma_{I_{1},I_{2};I_{2},I_{1}}(\omega_{1},0;0,\omega_{1})\cdot
O_{I_{2}}-\sum_{I_{2}}\int\frac{d\omega_{2}}{2\pi}\Gamma_{I_{1},I_{2};I_{2},I_{1}}(\omega_{1},\omega_{2};\omega_{2},\omega_{1})\cdot[\mathcal{G}(\omega_{2})]^{2}\cdot
O_{I_{2}}\ .$ (144)

Comparing Eq. 144 to Eq. 135, we have

$\frac{\partial\Sigma(\omega_{1})}{\partial
i\omega_{1}}\sum_{I_{1}}O_{I_{1}}^{2}=\sum_{I_{1}}\frac{\partial\Sigma_{I_{1}}(\omega_{1})}{\partial
h^{O}}O_{I_{1}}-\frac{\sin^{2}(\delta_{f})}{\pi\Delta_{0}}\sum_{I_{1}I_{2}}\Gamma_{I_{1},I_{2};I_{2},I_{1}}(\omega_{1},0;0,\omega_{1})\cdot
O_{I_{1}}O_{I_{2}}\ .$ (145)

Now we consider the response $\overline{O}=\sum_{I}O_{I}\langle
f^{\dagger}_{I}f_{I}\rangle$

$\displaystyle\overline{O}=-\sum_{I}O_{I}\langle
f_{I}(-0^{+})f_{I}^{\dagger}(0)\rangle=T\sum_{I}O_{I}\sum_{\omega}e^{i\omega
0^{+}}\mathcal{G}_{I}(\omega)$ (146)

to the external field $-h^{O}O_{I}$.

$\displaystyle\chi^{O}=\frac{\partial\overline{O}}{\partial
h^{O}}=-T\sum_{I}O_{I}\sum_{\omega}e^{i\omega
0^{+}}[\mathcal{G}(\omega)]^{2}\frac{\partial}{\partial
h^{O}}[h^{O}O_{I}-\Sigma_{I}(\omega)]=T\sum_{I}O_{I}\sum_{\omega}e^{i\omega
0^{+}}[\mathcal{G}(\omega)]^{2}[-O_{I}+\frac{\partial\Sigma_{I}(\omega)}{\partial
h^{O}}]\ .$ (147)

Making use of the Ward identity Eq. 145, we can replace the $h$-derivative by
the $\omega$-derivative

$\displaystyle\chi^{O}=$ $\displaystyle
T\sum_{I}O_{I}^{2}\sum_{\omega}e^{i\omega
0^{+}}[\mathcal{G}(\omega)]^{2}\left[-1+\frac{\partial\Sigma(\omega)}{\partial
i\omega}\right]+T\frac{\sin^{2}(\delta_{f})}{\pi\Delta_{0}}\sum_{II_{2}}\sum_{\omega}e^{i\omega
0^{+}}[\mathcal{G}(\omega)]^{2}\cdot\Gamma_{I,I_{2};I_{2},I}(\omega,0;0,\omega)O_{I}O_{I_{2}}$
(148)

$[\mathcal{G}_{I}(\omega)]^{2}\left[-1+\frac{\partial\Sigma(\omega)}{\partial
i\omega}\right]$ is the continuous part of the Green’s function derivative
(Eq. 140), and it can be rewritten as

$[\mathcal{G}_{I}(\omega)]^{2}\left[-1+\frac{\partial\Sigma(\omega)}{\partial
i\omega}\right]=\frac{\partial\mathcal{G}(\omega)}{\partial
i\omega}-\frac{\mathcal{G}(0^{+})-\mathcal{G}(0^{-})}{i}\delta(\omega)=\frac{\partial\mathcal{G}(\omega)}{\partial
i\omega}+2\pi\frac{\sin^{2}\delta_{f}}{\pi\Delta_{0}}\delta(\omega)\ .$ (149)

Thus,

$\chi^{O}=T\sum_{I}O_{I}^{2}\sum_{\omega}e^{i\omega
0^{+}}\frac{\partial\mathcal{G}(\omega)}{\partial
i\omega}+\frac{\sin^{2}\delta_{f}}{\pi\Delta_{0}}\sum_{I}O_{I}^{2}+T\frac{\sin^{2}(\delta_{f})}{\pi\Delta_{0}}\sum_{II_{2}}\sum_{\omega}e^{i\omega
0^{+}}[\mathcal{G}(\omega)]^{2}\cdot\Gamma_{I,I_{2};I_{2},I}(\omega,0;0,\omega)O_{I}O_{I_{2}}\
,$ (150)

where the first term vanishes. Making use of the identity Eq. 135, we rewrite
the above equation as

$\displaystyle\chi^{O}=$
$\displaystyle\frac{\sin^{2}\delta_{f}}{\pi\Delta_{0}}\sum_{I}O_{I}^{2}-\frac{\sin^{2}(\delta_{f})}{\pi\Delta_{0}}\cdot\sum_{I}\frac{\partial\Sigma_{I}(0)}{\partial
h^{O}}O_{I}=\frac{\sin^{2}(\delta_{f})}{\pi\Delta_{0}}\sum_{I}\left(O_{I}^{2}-\frac{\partial\Sigma_{I}(0)}{\partial
h^{O}}O_{I}\right)\ .$ (151)

Now we have expressed the exact susceptibility in terms of the derivatives of
zero-frequency self-energy. To relate the susceptibility to the vertex
function, we replace the $\partial_{h^{O}}\Sigma_{I}(0)$ by the
$\omega$-derivative via Eq. 145

$\displaystyle\chi^{O}=\frac{\sin^{2}(\delta_{f})}{\pi\Delta_{0}}\sum_{I}O_{I}^{2}\left(1-\frac{\Sigma(\omega)}{\partial
i\omega}\right)_{\omega=0}-\left(\frac{\sin^{2}(\delta_{f})}{\pi\Delta_{0}}\right)^{2}\sum_{I_{1}I_{2}}\Gamma_{I_{1}I_{2};I_{2}I_{1}}(00;00)\cdot
O_{I_{1}}O_{I_{2}}\ .$ (152)

Notice that $\left(1-\frac{\Sigma(\omega)}{\partial
i\omega}\right)_{\omega=0}$ equals inverse of the quasi-particle weight
$z^{-1}$, there is

$\chi^{O}=\frac{\sin^{2}(\delta_{f})}{z\cdot\pi\Delta_{0}}\left(\sum_{I}O_{I}^{2}-\frac{\sin^{2}(\delta_{f})}{z\cdot\pi\Delta_{0}}\sum_{I_{1}I_{2}}z^{2}\Gamma_{I_{1}I_{2};I_{2}I_{1}}(00;00)\cdot
O_{I_{1}}O_{I_{2}}\right)\ .$ (153)

We now have expressed the exact susceptibility in terms of quasi-particle
weight $z$, hybridization $\Delta_{0}$, filling $\delta_{f}=\pi(\nu_{f}+4)/8$,
and full vertex function at zero-frequency. It is worth emphasizing that this
identity applies only when $O$ is a generator of the continuous symmetry.
Applying it to the $O$ operators discussed in the last subsection, we find the
five susceptibilities in Eqs. (120) to (124) are exact.

###  B.5 Pairing susceptibilities and irreducible vertex in Cooper channel

Following the same procedure in Sec. B.3, we can calculate the pairing
susceptibilities $\chi_{\rm p}$.

One can diagonalize the pairing susceptibilities by enumerating distinct two-
electron representations, which have been discussed in detail in Sec. A.5.
There are five non-degenerate two-electron levels (Eqs. 61, 62, 63, 64 and
65). To proceed, we choose the following five operators to represent the five
channels

$\sum_{I_{1}I_{2}}f_{I_{1}}^{\dagger}O^{0xy}_{I_{1}I_{2}}f_{I_{2}}^{\dagger}\qquad\text{(intra-
orbital inter-valley singlet)}\ ,$ (154)
$\sum_{I_{1}I_{2}}f_{I_{1}}^{\dagger}O^{x0y}_{I_{1}I_{2}}f_{I_{2}}^{\dagger}\qquad\text{(inter-
orbital intra-valley singlet)}\ ,$ (155)
$\sum_{I_{1}I_{2}}f_{I_{1}}^{\dagger}O^{00y}_{I_{1}I_{2}}f_{I_{2}}^{\dagger}\qquad\text{(intra-
orbital intra-valley singlet)}\ ,$ (156)
$\sum_{I_{1}I_{2}}f_{I_{1}}^{\dagger}O^{xxy}_{I_{1}I_{2}}f_{I_{2}}^{\dagger}\qquad\text{(inter-
orbital inter-valley singlet)}\ ,$ (157)
$\sum_{I_{1}I_{2}}f_{I_{1}}^{\dagger}O^{xy0}_{I_{1}I_{2}}f_{I_{2}}^{\dagger}\qquad\text{(inter-
orbital inter-valley triplet)}\ ,$ (158)

where $O^{\mu\nu\rho}=\sigma^{\mu}\tau^{\nu}\varsigma^{\rho}$. Each of them is
chosen as a linear combination of a few degenerate two-electron states in Eqs.
61, 62, 63, 64 and 65.

To apply Eq. 95 to pairing susceptibility, we should define
$\hat{A}=\sum_{I_{1}I_{2}}f_{I_{2}}O^{*}_{I_{1},I_{2}}f_{I_{1}}$,
$\hat{B}=\sum_{I_{1}I_{2}}f_{I_{1}}^{\dagger}O_{I_{1},I_{2}}f^{\dagger}_{I_{2}}$,
for $O$ matrices given in the above equations. Notice that $O$ is always anti-
symmetric. Fouriering Eq. (95),

$\chi^{O}_{\rm
p}(\nu)=T\sum_{1234}\delta_{\omega_{1}+\omega_{2},\nu}\delta_{\omega_{3}+\omega_{4},\nu}\left\langle
e^{-S_{I}}f(2)O^{*}(1,2)f(1)f^{\dagger}(3)O(3,4)f^{\dagger}(4)\right\rangle_{0C}\
.$ (159)

At the zeroth order of $S_{I}$,

$\displaystyle\chi^{O(0)}_{\rm p}(\nu)$
$\displaystyle=T\sum_{1234}\delta_{\omega_{1}+\omega_{2},\nu}\delta_{\omega_{3}+\omega_{4},\nu}\left(\wick{\c{1}f(2)O^{*}(1,2)\c{2}f(1)\c{2}f^{\dagger}(3)O(3,4)\c{1}f^{\dagger}(4)}-\wick{\c{1}f(2)O^{*}(1,2)\c{2}f(1)\c{1}f^{\dagger}(3)O(3,4)\c{2}f^{\dagger}(4)}\right)$
$\displaystyle=\sum_{12}\left[O^{*}(1,2)O(1,2)-O^{*}(1,2)O(2,1)\right]T\sum_{\omega_{1}}\mathcal{G}^{0}(\omega_{1})\mathcal{G}^{0}(-\omega_{1}+\nu)$
$\displaystyle\overset{\nu\to 0^{+}}{=}\sum_{12}2\cdot
O^{*}(1,2)O(1,2)\cdot\frac{\tan\delta_{f}}{\pi\Delta_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)$
(160)

where we used

$\displaystyle
T\sum_{\omega}\frac{1}{i\omega-\epsilon_{f}+i\Delta(\omega)}\frac{1}{-i\omega-\epsilon_{f}+i\Delta(-\omega)}$
$\displaystyle=$ $\displaystyle
2T\sum_{\omega>0}\frac{1}{(\omega+\Delta_{0})^{2}+(\epsilon_{f})^{2}}=\frac{1}{\pi}\int_{0}^{\infty}\mathrm{d}\omega\frac{1}{(\omega+\Delta_{0})^{2}+(\epsilon_{f})^{2}}$
$\displaystyle=$
$\displaystyle\frac{1}{\pi\epsilon_{f}}\left(\frac{\pi}{2}-\arctan\frac{\Delta_{0}}{\epsilon_{f}}\right)=\frac{\tan\delta_{f}}{\pi\Delta_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)$

and $O(1,2)=-O(2,1)$. Thus the zeroth order susceptibility can be written as

$\chi^{O(0)}_{\rm
p}=2\cdot\mathrm{tr}\left[O^{\dagger}O\right]\cdot\frac{\tan\delta_{f}}{\pi\Delta_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\
.$ (162)

At the first order of $S_{I}$,

$\displaystyle\chi^{O(1)}_{\rm p}(\nu)$ $\displaystyle=$
$\displaystyle-T^{2}\sum_{1234}\delta_{\omega_{1}+\omega_{2},\nu}\delta_{\omega_{3}+\omega_{4},\nu}O^{*}(1,2)O(3,4)\sum_{5678}\frac{1}{4}\delta_{\omega_{5}+\omega_{6},\omega_{7}+\omega_{8}}\Gamma^{0}(5,6;7,8)$
$\displaystyle\times\left\langle
f(2)f(1)f^{\dagger}(3)f^{\dagger}(4)f^{\dagger}(5)f^{\dagger}(6)f(7)f(8)\right\rangle_{0C}$
$\displaystyle=$
$\displaystyle-T^{2}\sum_{1234}\delta_{\omega_{1}+\omega_{2},\nu}\delta_{\omega_{3}+\omega_{4},\nu}O^{*}(1,2)O(3,4)\sum_{5678}\frac{1}{4}\delta_{\omega_{5}+\omega_{6},\omega_{7}+\omega_{8}}\Gamma^{0}(5,6;7,8)$
$\displaystyle\times
4\times\wick{\c{1}f(2)\c{2}f(1)\c{3}f^{\dagger}(3)\c{4}f^{\dagger}(4)\c{1}f^{\dagger}(5)\c{2}f^{\dagger}(6)\c{3}f(7)\c{4}f(8)}$

where the $\times 4$ factor is contributed by equivalent contractions. The
contraction gives

$\displaystyle(-1)^{2}\delta_{25}\delta_{16}\delta_{37}\delta_{48}\mathcal{G}^{0}(\omega_{1})\mathcal{G}^{0}(\omega_{2})\mathcal{G}^{0}(\omega_{3})\mathcal{G}^{0}(\omega_{4})\
,$ (164)

thereby,

$\displaystyle\chi^{O(1)}_{\rm p}(0^{+})$ $\displaystyle=$
$\displaystyle-T^{2}\sum_{1234}O^{*}(1,2)O(3,4)\Gamma^{0}(2,1;3,4)\sum_{\omega_{1},\omega_{3}}\mathcal{G}^{0}(\omega_{1})\mathcal{G}^{0}(-\omega_{1})\mathcal{G}^{0}(\omega_{3})\mathcal{G}^{0}(-\omega_{3})$
$\displaystyle=$
$\displaystyle-\left[\sum_{1234}O^{*}(1,2)\Gamma^{0}(2,1;3,4)O(3,4)\right]\left(\frac{\tan\delta_{f}}{\pi\Delta_{0}}(\frac{\pi}{2}-\delta_{f})\right)^{2}\
.$

Replacing ${\Delta}_{0}$, $\Gamma^{0}$ by the effective parameters of quasi-
particles $\widetilde{\Delta}_{0}$, $\widetilde{\Gamma}^{0}$, we obtain the
quasi-particle contributed pairing susceptibilities

$\widetilde{\chi}^{O}_{\rm p}\approx
2\cdot\mathrm{tr}[O^{\dagger}O]\cdot\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\left[1-\kappa_{\rm
p}[O,\widetilde{\Gamma}^{0}]\cdot\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\right]$
(166)

where

$\kappa_{\rm
p}[O,\widetilde{\Gamma}^{0}]=\frac{1}{2\cdot\mathrm{tr}[O^{\dagger}O]}\sum_{I_{1,2,3,4}}O^{*}_{I_{2},I_{1}}\widetilde{\Gamma}^{0}_{I_{1},I_{2};I_{3},I_{4}}(0,0;0,0)O_{I_{3},I_{4}}\
.$ (167)

One should be aware that $\widetilde{\Gamma}^{0}$ equals the exact full vertex
function $\widetilde{\Gamma}$ of quasi-particles at zero frequency, but the
frequency dependence of $\widetilde{\Gamma}$ has been omitted in the above
equation. We have also omitted the frequency dependence the quasi-particle
self-energy $\widetilde{\Sigma}(\omega)$ (beyond quasi-particle weight
correction), e.g., $\omega^{2}$ terms in $\widetilde{\Sigma}(\omega)$.
Therefore, Eq. 166 is an approximate result.

We find $\kappa_{\rm p}[O,\widetilde{\Gamma}^{0}]$ is nothing but the energy
(in terms of bare $\widetilde{\Gamma}^{0}$) of the two-electron states defined
by $O$. Readers may directly verify this. Here we prove this statement. $O$
defines a two-electron state
$|\Phi\rangle=\sum_{IJ}O_{IJ}f_{I}^{\dagger}f_{J}^{\dagger}|0\rangle$. Its
norm is given by
$\langle\Phi|\Phi\rangle=\sum_{IJ}O_{IJ}^{*}O_{IJ}-O_{IJ}^{*}O_{JI}=2\cdot\mathrm{tr}[O^{\dagger}O]$.
Then $\kappa_{\rm p}[O,\widetilde{\Gamma}^{0}]$ is by definition
$\langle\Phi|H_{\rm eff}|\Phi\rangle/\langle\Phi|\Phi\rangle$ with $H_{\rm
eff}$ being the effective Hamiltonian defined by $\widetilde{\Gamma}^{0}$. For
the two-electron states given in Eqs. 61, 62, 63, 64 and 65, where $O$ is
chosen as $O^{0xy}$, $O^{x0y}$, $O^{00y}$, $O^{xxy}$, $O^{xy0}$, respectively,
the corresponding pairing susceptibilities are

$\widetilde{\chi}_{\rm p}^{0xy}\approx
16\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\left[1-\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\widetilde{U}_{2}\right]\
,$ (168) $\widetilde{\chi}_{\rm p}^{x0y}\approx
16\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\left[1-\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\widetilde{U}_{3}\right]\
,$ (169) $\widetilde{\chi}_{\rm p}^{00y}\approx
16\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\left[1-\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\widetilde{U}_{4}\right]\
,$ (170) $\widetilde{\chi}_{\rm p}^{xxy}\approx
16\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\left[1-\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\left(\widetilde{U}_{1}-\widetilde{\mathcal{J}}\right)\right]\
,$ (171) $\widetilde{\chi}_{\rm p}^{xy0}\approx
16\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\left[1-\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\left(\widetilde{U}_{1}+\widetilde{\mathcal{J}}\right)\right]\
.$ (172)

We consider to decompose $\widetilde{\chi}_{\rm p}$ into a summation of ladder
diagrams

(173)

where $\Gamma^{\rm p}$ is the irreducible vertex in the pairing channel.
$\Gamma^{\rm p}$ cannot be divided into two by cutting two left-going Green’s
function lines. Due to the $\rm
U(1)\\!\times\\!U(1)\\!\times\\!U(2)\\!\times\\!U(2)$ symmetry, $\Gamma^{\rm
p}$ must have the same form as $\Gamma$, and can be parameterized by
$\Gamma^{\rm p}_{U1,2,3,4}$, ${\Gamma}^{\rm p}_{\mathcal{J}}$. We define the
effective interaction in Cooper channel as $U_{1,2,3,4}^{\rm
p}=z^{2}\Gamma^{\rm p}_{U1,2,3,4}$, $\mathcal{J}^{\rm p}=z^{2}\Gamma^{\rm
p}_{\mathcal{J}}$. Neglecting the frequency dependencies of $\Gamma^{\rm p}$
and the renormalized self-energy, the pairing susceptibilities can be written
in terms of effective interactions as

$\widetilde{\chi}_{\rm p}^{0xy}\approx
16\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\left[1+\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)U_{2}^{\rm
p}\right]^{-1}\ ,$ (174) $\widetilde{\chi}_{\rm p}^{x0y}\approx
16\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\left[1+\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)U_{3}^{\rm
p}\right]^{-1}\ ,$ (175) $\widetilde{\chi}_{\rm p}^{00y}\approx
16\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\left[1+\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)U_{4}^{\rm
p}\right]^{-1}\ ,$ (176) $\widetilde{\chi}_{\rm p}^{xxy}\approx
16\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\left[1+\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\left(U_{1}^{\rm
p}-\mathcal{J}^{\rm p}\right)\right]^{-1}\ ,$ (177) $\widetilde{\chi}_{\rm
p}^{xy0}\approx
16\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\left[1+\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\left(U_{1}^{\rm
p}+\mathcal{J}^{\rm p}\right)\right]^{-1}\ .$ (178)

Equaling these equations to Eqs. (168) to (172), we obtain the effective
interactions (in the pairing channel) in terms of the renormalized
interactions

$U^{\rm
p}_{2}=\frac{\widetilde{U}_{2}}{1-\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\widetilde{U}_{2}}\
,\qquad U^{\rm
p}_{3}=\frac{\widetilde{U}_{3}}{1-\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\widetilde{U}_{3}}\
,\qquad U^{\rm
p}_{4}=\frac{\widetilde{U}_{4}}{1-\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)\widetilde{U}_{4}}\
,$ (179) $U^{\rm p}_{1}-{\mathcal{J}}^{\rm
p}=\frac{\widetilde{U}_{1}-\widetilde{\mathcal{J}}}{1-\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)(\widetilde{U}_{1}-\widetilde{\mathcal{J}})}\
,\qquad U^{\rm p}_{1}+{\mathcal{J}}^{\rm
p}=\frac{\widetilde{U}_{1}+\widetilde{\mathcal{J}}}{1-\frac{\tan\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(\frac{\pi}{2}-\delta_{f}\right)(\widetilde{U}_{1}+\widetilde{\mathcal{J}})}\
.$ (180)

###  B.6 Asymptotic behavior of exact vertex function

####  B.6.1 The $U\gg\mathcal{J}\gg T_{\rm K}$ limit

##### The $\nu_{f}=\pm 2$ states.

We focus on the doped the correlated insulators at the total fillings
$\nu=\pm(2+\delta\nu)$, where $0\\!<\\!\delta\nu\\!<\\!1$. According to the
(approximate) particle-hole symmetry of the problem, we only consider the
$\nu=-2$ case. According to the calculations in Ref. Zhou  _et al._ (2023); Hu
_et al._ (2023b); Datta  _et al._ (2023); Rai  _et al._ (2023), the strong
repulsion interaction will fix $\nu_{f}$ around $-2$ for $-3<\nu<-2$. The two-
electron states are already discussed in Sec. A.5. In the $T_{\rm
K}\sim\widetilde{\Delta}_{0}\ll\mathcal{J}$ limit, only the ground states (Eq.
64) of $H_{I}$

$\frac{f_{\alpha\eta\uparrow}^{\dagger}f_{\overline{\alpha}\overline{\eta}\downarrow}^{\dagger}-f_{\alpha\eta\downarrow}^{\dagger}f_{\overline{\alpha}\overline{\eta}\uparrow}^{\dagger}}{\sqrt{2}}|0\rangle,\qquad(\alpha\eta=1+,1-),\qquad
E=U-2\mathcal{J}\ .$ (181)

will participate in the Kondo screening Nozières and Blandin (1980). (See
section 6.4 of Ref. Nozières and Blandin (1980) for the discussion about
splittings in atomic levels.) They transform as $d_{x^{2}-y^{2}}$ and $d_{xy}$
orbitals under the symmetry operations in $D_{6}$ and form the $E_{2}$
representation. We have shown that Eqs. (120) to (124) are the exact
susceptibilities of bare particles. The charge
($\sigma^{0}\tau^{0}\varsigma^{0}$), the valley-charge
($\sigma^{0}\tau^{z}\varsigma^{0}$), the orbital-charge
($\sigma^{z}\tau^{0}\varsigma^{0}$), and the total spin
($\sigma^{0}\tau^{0}\varsigma^{z}$) take fixed values 2, 0, 0, 0 in the ground
state manifold, respectively. Therefore, these degrees of freedom are frozen
at the Kondo energy scale $\widetilde{\Delta}_{0}$, and the corresponding
susceptibilities are not contributed by quasi-particles. For example, in the
single-orbital Anderson impurity model at half-filling, the charge
susceptibility given by the Bethe ansatz is Kawakami and Okiji (1982)

$\chi^{c}=\frac{1}{\pi\cdot U}\int_{-\infty}^{\infty}dy\cdot
e^{-\pi\frac{U}{2\Delta_{0}}y^{2}}\frac{(U/2\Delta_{0})^{\frac{3}{2}}}{(y+\frac{1}{2})^{2}(U/2\Delta_{0})^{2}+\frac{1}{4}}\
.$ (182)

One can see that this $\chi^{c}$ is much smaller than the quasi-particle
density of states at Fermi level ($\widetilde{\Delta}_{0}^{-1}$). In fact,
$\chi^{c}\cdot\widetilde{\Delta}_{0}\ll 1$ is a universal behavior of Fermi
liquid in the Kondo regime. For Eq. 120 to reproduce this correct behavior of
$\chi^{c}$, the renormalized interactions in the bracket must cancel the
$\widetilde{\Delta}_{0}^{-1}$ divergence, implying the constraint

$2\widetilde{U}_{1}+2\widetilde{U}_{2}+2\widetilde{U}_{3}+\widetilde{U}_{4}+\widetilde{\mathcal{J}}=\frac{\pi\widetilde{\Delta}_{0}}{\sin^{2}\delta_{f}}\
.$ (183)

The same argument also applies to valley, orbital, and spin degrees of freedom
because they are also frozen at the Kondo energy scale and are not contributed
by quasi-particles. The constraints
$\chi^{c,v,o,s}\cdot\widetilde{\Delta}_{0}\ll 1$ imply

$\widetilde{U}_{1}=-\frac{\pi\widetilde{\Delta}_{0}}{\sin^{2}\delta_{f}},\qquad\widetilde{U}_{2}=\widetilde{U}_{3}=\frac{\pi\widetilde{\Delta}_{0}}{\sin^{2}\delta_{f}}-\frac{1}{2}\widetilde{\mathcal{J}},\qquad\widetilde{U}_{4}=-\frac{\pi\widetilde{\Delta}_{0}}{\sin^{2}\delta_{f}}+\widetilde{\mathcal{J}}\
.$ (184)

The two-electron energies of a generic $H_{I}$ are given in Sec. A.5.
Replacing the bare interaction parameters with the renormalized interaction
parameters, we can obtain eigenvalues of $\widetilde{\Gamma}^{0}$ as

$\widetilde{U}_{2}=\frac{\pi\widetilde{\Delta}_{0}}{\sin^{2}\delta_{f}}-\frac{1}{2}\widetilde{\mathcal{J}}\;\text{(intra-
orbital inter-
valley)},\qquad\widetilde{U}_{3}=\frac{\pi\widetilde{\Delta}_{0}}{\sin^{2}\delta_{f}}-\frac{1}{2}\widetilde{\mathcal{J}}\;\text{(inter-
valley intra-orbital)}\ ,$ (185)
$\widetilde{U}_{4}=-\frac{\pi\widetilde{\Delta}_{0}}{\sin^{2}\delta_{f}}+\widetilde{\mathcal{J}}\;\text{(intra-
orbital intra-valley singlet)}\ ,$ (186)
$\widetilde{U}_{1}-\widetilde{\mathcal{J}}=-\frac{\pi\widetilde{\Delta}_{0}}{\sin^{2}\delta_{f}}-\widetilde{\mathcal{J}}\;\text{(inter-
valley inter-orbital
singlet)},\qquad\widetilde{U}_{1}+\widetilde{\mathcal{J}}=-\frac{\pi\widetilde{\Delta}_{0}}{\sin^{2}\delta_{f}}+\widetilde{\mathcal{J}}\;\text{(inter-
valley inter-orbital triplet)}\ .$ (187)

Here $\delta_{f}=\frac{\pi}{4}$. One of the last two must be negative.
Therefore, we prove the statement that the renormalized interaction has at
least one negative channel.

Unlike the charge, valley, orbital, and spin, the angular momentum
($\sigma^{z}\tau^{0}\varsigma^{z}$) is not quenched in the ground state of
$H_{I}$. Thus, the quasi-particle contributed part of $\chi^{a}$ should
diverge at the order $\widetilde{\Delta}_{0}^{-1}$ in the
$\widetilde{\Delta}_{0}\to 0$ limit. $\chi^{a}$ expressed in terms of
$\widetilde{\mathcal{J}}$ is

$\chi^{a}=8\frac{\sin^{2}\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(8-4\frac{\sin^{2}\delta_{f}}{\pi\widetilde{\Delta}_{0}}\widetilde{\mathcal{J}}\right)\
.$ (188)

Requiring $\chi^{a}$ non-negative leads to the condition
$\widetilde{\mathcal{J}}\leq
2\frac{\pi\widetilde{\Delta}_{0}}{\sin^{2}\delta_{f}}$.

We can also extract useful conditions from the other quasi-particle
susceptibilities in Eqs. (125) to (129) and (168) to (172). Even though they
are not related to the exact susceptibilities of the bare particles, they
should also be positive for the Fermi liquid theory to be valid. Then we have
the inequalities

$\widetilde{\chi}^{0zz}\approx
8\frac{\sin^{2}\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(0+\frac{\sin^{2}\delta_{f}}{\pi\widetilde{\Delta}_{0}}2\widetilde{\mathcal{J}}\right)\gtrsim
0,\qquad\widetilde{\chi}^{x00}=\widetilde{\chi}^{xz0}=\widetilde{\chi}^{0x0}=\widetilde{\chi}^{zx0}\approx
8\frac{\sin^{2}\delta_{f}}{\pi\widetilde{\Delta}_{0}}\cdot\left(2-\frac{\sin^{2}\delta_{f}}{\pi\widetilde{\Delta}_{0}}\frac{1}{2}\widetilde{\mathcal{J}}\right)\gtrsim
0\ ,$ (189) $\widetilde{\chi}^{xx0}\approx
8\frac{\sin^{2}\delta_{f}}{\pi\widetilde{\Delta}_{0}}\left(0+\frac{\sin^{2}\delta_{f}}{\pi\widetilde{\Delta}_{0}}2\widetilde{\mathcal{J}}\right)\gtrsim
0,\qquad\widetilde{\chi}^{xxz}=8\frac{\sin^{2}\delta_{f}}{\pi\widetilde{\Delta}_{0}}\cdot
0\gtrsim 0\ ,$ (190) $\widetilde{\chi}_{\rm p}^{0xy}=\widetilde{\chi}_{\rm
p}^{x0y}\approx\frac{4}{\widetilde{\Delta}_{0}}\left(1-\frac{1}{4\widetilde{\Delta}_{0}}\left(2\pi\widetilde{\Delta}_{0}-\frac{1}{2}\widetilde{\mathcal{J}}\right)\right)\gtrsim
0,\qquad\widetilde{\chi}_{\rm
p}^{00y}\approx\frac{4}{\widetilde{\Delta}_{0}}\left(1-\frac{1}{4\widetilde{\Delta}_{0}}\left(-2\pi\widetilde{\Delta}_{0}+\widetilde{\mathcal{J}}\right)\right)\gtrsim
0,$ (191) $\widetilde{\chi}_{\rm
p}^{xxy}\approx\frac{4}{\widetilde{\Delta}_{0}}\left(1-\frac{1}{4\widetilde{\Delta}_{0}}\left(-2\pi\widetilde{\Delta}_{0}-\widetilde{\mathcal{J}}\right)\right)\gtrsim
0,\qquad\widetilde{\chi}_{\rm
p}^{xy0}\approx\frac{4}{\widetilde{\Delta}_{0}}\left(1-\frac{1}{4\widetilde{\Delta}_{0}}\left(-2\pi\widetilde{\Delta}_{0}+\widetilde{\mathcal{J}}\right)\right)\gtrsim
0.$ (192)

As the susceptibilities are approximate, we use “$\gtrsim$” rather than
“$\geq$” in these inequalities. This requirement leads to

$4.6\widetilde{\Delta}_{0}\approx
8\left(\frac{\pi}{2}-1\right)\widetilde{\Delta}_{0}\lesssim\widetilde{\mathcal{J}}\lesssim\left(1+\frac{\pi}{2}\right)\widetilde{\Delta}_{0}\approx
10.3\widetilde{\Delta}_{0}\ .$ (193)

One can directly verify that the inter-orbital inter-valley singlet pairing
fluctuation ($\chi^{xxy}_{\rm p}$) is strongest among all the pairings. The
universality hypothesis states that $T_{\rm K}$, or $\widetilde{\Delta}_{0}$,
is the only energy scale in the Kondo regime, implying that the ratio
$\widetilde{\mathcal{J}}/\widetilde{\Delta}_{0}$ is a universal constant. Our
analyses above show that this ratio is in the range from 4.6 to 10.3.

For later convenience, here we also express the effective interactions in the
pairing channel as functions of $\widetilde{\mathcal{J}}$ (Eqs. 179 and 180)

$U^{\rm p}_{2}=U^{\rm
p}_{3}\approx\frac{2\pi\widetilde{\Delta}_{0}-\frac{1}{2}\widetilde{\mathcal{J}}}{1-\frac{1}{4\widetilde{\Delta}_{0}}(2\pi\widetilde{\Delta}_{0}-\frac{1}{2}\widetilde{\mathcal{J}})},\qquad
U^{\rm
p}_{4}\approx\frac{-2\pi\widetilde{\Delta}_{0}+\widetilde{\mathcal{J}}}{1-\frac{1}{4\widetilde{\Delta}_{0}}(-2\pi\widetilde{\Delta}_{0}+\widetilde{\mathcal{J}})}\
,$ (194) $U^{\rm p}_{1}-\mathcal{J}^{\rm
p}\approx\frac{-2\pi\widetilde{\Delta}_{0}-\widetilde{\mathcal{J}}}{1-\frac{1}{4\widetilde{\Delta}_{0}}(-2\pi\widetilde{\Delta}_{0}-\widetilde{\mathcal{J}})},\qquad
U^{\rm p}_{1}+\mathcal{J}^{\rm
p}\approx\frac{-2\pi\widetilde{\Delta}_{0}+\widetilde{\mathcal{J}}}{1-\frac{1}{4\widetilde{\Delta}_{0}}(-2\pi\widetilde{\Delta}_{0}+\widetilde{\mathcal{J}})}\
.$ (195)

##### The $\nu_{f}=0$ states.

The flattened interaction (Eq. 56) is designed for the $\nu_{f}=\pm 2$ states
in such a way that it reproduces the correct two-electron ground states, as
explained at the beginning of Sec. A.5. Thus, it may not apply to the
$\nu_{f}=0$ states. However, here we still discuss its physics at $\nu_{f}=0$
for theoretical interests. The four-particle ground state is unique. It has
$N_{\alpha\eta}=1$, $\mathbf{S}_{\alpha\eta}^{2}=\frac{3}{4}$ for all
$\alpha\eta$ flavors, and
$(\mathbf{S}_{1+}+\mathbf{S}_{2-})^{2}=(\mathbf{S}_{1-}+\mathbf{S}_{2+})^{2}=0$.
Its energy can be read from Sec. A.5 as $E=6U-2\mathcal{J}$. Apart from
$\chi^{c,s,v,o}\cdot\widetilde{\Delta}_{0}\ll 1$, there is also
$\chi^{a}\cdot\widetilde{\Delta}_{0}\ll 1$ in the
$\widetilde{\Delta}_{0}\ll\mathcal{J}$ limit, which thus fixes
$\widetilde{J}=2\frac{\pi\widetilde{\Delta}}{\sin^{2}\delta_{f}}=2\pi\widetilde{\Delta}_{0}$,
where $\delta_{f}=\frac{\pi}{2}$ at $\nu_{f}=0$. Substituting the renormalized
interactions and $\delta_{f}\to\frac{\pi}{2}$ into Eqs. 168, 169, 170, 171 and
172, we have $\tan\delta_{f}\cdot(\frac{\pi}{2}-\delta_{f})=1$
$\widetilde{U}_{2}=\widetilde{U}_{3}=0$,
$\widetilde{U}_{4}=\pi\widetilde{\Delta}_{0}$,
$\widetilde{U}_{1}-\widetilde{\mathcal{J}}=-3\pi\widetilde{\Delta}_{0}$,
$\widetilde{U}_{1}+\widetilde{\mathcal{J}}=\pi\widetilde{\Delta}_{0}$ and
hence

$\widetilde{\chi}_{\rm p}^{0xy}=\widetilde{\chi}_{\rm
p}^{x0y}\approx\frac{16}{\pi\widetilde{\Delta}_{0}}\left(1+0\right),\qquad\widetilde{\chi}_{\rm
p}^{00y}\approx\frac{16}{\pi\widetilde{\Delta}_{0}}\left(1-1\right),$ (196)
$\widetilde{\chi}_{\rm
p}^{xxy}\approx\frac{16}{\pi\widetilde{\Delta}_{0}}\left(1+3\right),\qquad\widetilde{\chi}_{\rm
p}^{xy0}\approx\frac{16}{\pi\widetilde{\Delta}_{0}}\left(1-1\right)\ .$ (197)

One can see that the inter-orbital inter-valley singlet fluctuation is the
only one that is promoted by interaction. However, as $\widetilde{\Delta}_{0}$
is negligible at $\nu_{f}=0$ Zhou  _et al._ (2023); Hu  _et al._ (2023a), the
energy scale of renormalized attractive interactions is also negligible.

####  B.6.2 The $U\gg T_{\rm K}\gg\mathcal{J}$ limit

We only focus the states at $\nu_{f}=-2$ in this subsubsection. On the one
hand, only the charge degree of freedom is frozen in the $U\gg T_{\rm
K}\gg\widetilde{\Delta}_{0}\sim T$ limit, implying (Eq. 183)

$2\widetilde{U}_{1}+2\widetilde{U}_{2}+2\widetilde{U}_{3}+\widetilde{U}_{4}+\widetilde{\mathcal{J}}=\frac{\pi\widetilde{\Delta}_{0}}{\sin^{2}\delta_{f}}\
.$ (198)

$\widetilde{\mathcal{J}}$ is barely renormalized in the Fermi liquid phase
because it is much smaller than the Kondo energy scale. Thus, there must be
$\widetilde{\mathcal{J}}\approx\mathcal{J}$. On the other hand, the model has
an U(8) symmetry in the $\widetilde{\mathcal{J}}\to 0$ limit. This approximate
U(8) symmetry implies
$\widetilde{U}_{1}\approx\widetilde{U}_{2}\approx\widetilde{U}_{3}\approx\widetilde{U}_{4}$.
Therefore, we have

$\widetilde{U}_{1}\approx\widetilde{U}_{2}\approx\widetilde{U}_{3}\approx\widetilde{U}_{4}\approx\frac{2\pi}{7}\widetilde{\Delta}_{0},\qquad\widetilde{\mathcal{J}}=\mathcal{J}\
.$ (199)

The ground state of the renormalized interaction is still the $d$-wave $E_{2}$
states (Eq. 64). It has the energy

$\widetilde{U}_{1}-\widetilde{\mathcal{J}}=\frac{2\pi}{7}\widetilde{\Delta}_{0}-\mathcal{J}\
.$ (200)

This result implies that the pairing potential should start becoming negative
when $\mathcal{J}\sim T_{\rm K}$. Combined with the discussions in the last
subsubsection, this shows that the pairing potential is negative when $T_{\rm
K}\ll\mathcal{J}$ and remains negative to the regime $T_{\rm K}>\mathcal{J}$
as $T_{\rm K}$ increases, as shown in Fig. 2(e).

###  B.7 The counter term

In this section we demonstrate how the counter term (Eq. 91) is determined
perturbatively in terms of $\widetilde{\Gamma}^{0}$.

The first order $\widetilde{\Gamma}$ (Fig. 4(a)) is

$\widetilde{\Gamma}^{(1)}(00;00)=\widetilde{\Gamma}^{0}(00;00)+\Lambda^{(1)}\
,$ (201)

where $\Lambda^{(1)}$ is the first-order counter term.
$\widetilde{\Gamma}^{(1)}(00;00)$ equals $\widetilde{\Gamma}^{0}(00;00)$ by
definition (Eq. 94). Thus, $\Lambda^{(1)}=0$.

The first order self-energy is given by the diagrams in Fig. 4(b), where the
first term is contributed by $\widetilde{\Gamma}^{0}$ and the second term is
contributed by the first-order counter terms $\lambda_{0,1}^{(1)}$. To derive
this self-energy, we consider its correction to the Green’s function

$\displaystyle\delta\widetilde{\mathcal{G}}^{(1)}(1)=$
$\displaystyle-\left\langle
f(1)\left(-\sum_{2}f^{\dagger}(2)(-\lambda_{0}^{(1)}\omega_{2}+\lambda_{1}^{(1)})f(2)-\frac{T}{4}\sum_{2345}\delta_{\omega_{2}+\omega_{3},\omega_{4}+\omega_{5}}\widetilde{\Gamma}^{0}(23;45)f^{\dagger}(2)f^{\dagger}(3)f(4)f(5)\right)f^{\dagger}(1)\right\rangle_{0C}$
$\displaystyle=$
$\displaystyle[\widetilde{\mathcal{G}}^{0}(1)]^{2}\times\left(-\lambda_{0}^{(1)}\omega_{1}+\lambda_{1}^{(1)}+4\frac{T}{4}\sum_{2}\widetilde{\Gamma}^{0}(12;21)\widetilde{\mathcal{G}}^{0}(2)\right)$
(202)

where the factor 4 in the second term of the second row comes from 4
equivalent contractions, $\widetilde{\mathcal{G}}^{0}$ is the free propagator
of the quasi-particles, and the subscript $0C$ means connected diagrams. The
self-energy can be read as

$\widetilde{\Sigma}^{(1)}(1)=-i\lambda_{0}^{(1)}\omega_{1}+\lambda_{1}^{(1)}+T\sum_{2}\widetilde{\Gamma}^{0}(12;21)\widetilde{\mathcal{G}}^{0}(2)=-i\lambda_{0}^{(1)}\omega_{1}+\lambda_{1}^{(1)}+\sum_{I_{2}}\widetilde{\Gamma}^{0}_{I_{1}I_{2};I_{2}I_{1}}\int\frac{d\omega_{2}}{2\pi}\widetilde{\mathcal{G}}^{0}(\omega_{2})e^{i0^{+}\omega_{2}}\
,$ (203)

where the factor $e^{i0^{+}\omega_{2}}$ factor is due to the time-ordering in
the Green’s function $-\langle f(-0^{+})f^{\dagger}(0)\rangle$ such that the
creation operator is before the annihilation operator. The interaction factor
is (Eq. 109)

$\displaystyle\sum_{I_{2}}\widetilde{\Gamma}^{0}_{I_{1}I_{2};I_{2}I_{1}}=\frac{1}{8}\sum_{I_{1}I_{2}}\widetilde{\Gamma}^{0}_{I_{1}I_{2};I_{2}I_{1}}=$
$\displaystyle
2\widetilde{U}_{1}+2\widetilde{U}_{2}+2\widetilde{U}_{3}+\widetilde{U}_{4}+\widetilde{\mathcal{J}}\
.$ (204)

The frequency integral is

$\int_{-\infty}^{\infty}\frac{d\omega_{2}}{2\pi}e^{i0^{+}\omega_{2}}\frac{1}{i\omega_{2}+i\widetilde{\Delta}_{0}\cdot\mathrm{sgn}(\omega_{2})-\widetilde{\epsilon}_{f}}=\int_{0}^{\infty}\frac{d\omega_{2}}{2\pi}e^{i0^{+}\omega_{2}}\frac{1}{i\omega_{2}+i\widetilde{\Delta}_{0}-\widetilde{\epsilon}_{f}}+c.c.$
(205)

Define $z=i\omega_{2}$, the function
$\frac{1}{z+i\widetilde{\Delta}_{0}-\widetilde{\epsilon}_{f}}$ is analytical
in the half plane $\mathrm{Im}[z]>0$. We can continuously deform the integral
line from $(0,i\infty)$ to $(0,-\infty)$, then we have

$\int_{0}^{-\infty}\frac{dz}{2\pi
i}e^{0^{+}z}\frac{1}{z+i\widetilde{\Delta}_{0}-\widetilde{\epsilon}_{f}}+c.c.=\frac{1}{\pi}\int_{-\infty}^{0}dz\cdot
e^{0^{+}z}\frac{\widetilde{\Delta}_{0}}{(z-\widetilde{\epsilon}_{f})^{2}+\widetilde{\Delta}_{0}^{2}}=\pi\delta_{f}\
.$ (206)

Therefore, there is

$\widetilde{\Sigma}^{(1)}(\omega)=-i\lambda_{0}^{(1)}\omega+\lambda_{1}^{(1)}+\left(2\widetilde{U}_{1}+2\widetilde{U}_{2}+2\widetilde{U}_{3}+\widetilde{U}_{4}+\widetilde{\mathcal{J}}\right)\pi\delta_{f}$
(207)

The condition Eq. 94 implies

$\lambda_{0}^{(1)}=0,\qquad\lambda_{1}^{(1)}=-\left(2\widetilde{U}_{1}+2\widetilde{U}_{2}+2\widetilde{U}_{3}+\widetilde{U}_{4}+\widetilde{\mathcal{J}}\right)\pi\delta_{f}\
.$ (208)

This also implies $\widetilde{\Sigma}^{(1)}(\omega)=0$.

Figure 4: Diagrams to determine the counter term. (a), (b), (c), (d) are
diagrams for the first order vertex, first order self-energy, second order
vertex, and second order self-energy, respectively.

The second order vertex is contributed by diagrams in Fig. 4(c). The last two
diagrams cancel each other because $\widetilde{\Sigma}^{(1)}(\omega)=0$. Then
$\Lambda^{(2)}$ should cancel the first two diagrams at zero frequency. We do
not attempt to calculate $\Lambda^{(2)}$ in this work.

We also do not attempt to calculate the full second order self-energy
$\widetilde{\Sigma}^{(2)}$. Instead, we only show how to calculate its linear
frequency-dependence, i.e., $\partial_{\omega}\widetilde{\Sigma}^{(2)}$. The
condition $\partial_{\omega}\widetilde{\Sigma}^{(2)}=0$ (Eq. 94) will
determine the counter term $\lambda_{0}^{(2)}$. $\widetilde{\Sigma}^{(2)}$ is
contributed by the diagrams shown in Fig. 4(d). The last two diagrams cancel
each other because $\widetilde{\Sigma}^{(1)}(\omega)=0$ according to the
second paragraph above. The third diagram is contributed by $\Lambda^{(2)}$,
which by definition is $\omega$-independent (Eq. 91) and hence does not
contribute to $\partial_{\omega}\widetilde{\Sigma}^{(2)}$. Thus, we only need
to consider the first two diagrams, dubbed as $\Sigma^{(2,1)}$ and
$\Sigma^{(2,2)}$, respectively. The first diagram gives the correction

$\displaystyle\delta\widetilde{\mathcal{G}}^{(2,1)}(1)=$
$\displaystyle-\left\langle
f(1)\left(\frac{1}{2}\frac{T^{2}}{4^{2}}\sum_{\begin{subarray}{c}2345\\\
5^{\prime}4^{\prime}3^{\prime}2^{\prime}\end{subarray}}\delta_{\omega_{2}+\omega_{3},\omega_{4}+\omega_{5}}\delta_{\omega_{5}^{\prime}+\omega_{4}^{\prime},\omega_{3}^{\prime}+\omega_{2}^{\prime}}\widetilde{\Gamma}^{0}(23;45)\widetilde{\Gamma}^{0}(5^{\prime}4^{\prime};3^{\prime}2^{\prime})f^{\dagger}(2)f^{\dagger}(3)f(4)f(5)f^{\dagger}(5^{\prime})f^{\dagger}(4^{\prime})f(3^{\prime})f(2^{\prime})\right)f^{\dagger}(1)\right\rangle_{0C}$
(209)

Applying Wick’s theorem, there is

$\displaystyle\delta\widetilde{\mathcal{G}}^{(2,1)}(1)$
$\displaystyle=-2^{4}\frac{1}{2}\frac{T^{2}}{4^{2}}\sum_{\begin{subarray}{c}2345\\\
5^{\prime}4^{\prime}3^{\prime}2^{\prime}\end{subarray}}\delta_{\omega_{2}+\omega_{3},\omega_{4}+\omega_{5}}\delta_{\omega_{5}^{\prime}+\omega_{4}^{\prime},\omega_{3}^{\prime}+\omega_{2}^{\prime}}\widetilde{\Gamma}^{0}(23;45)\widetilde{\Gamma}^{0}(5^{\prime}4^{\prime};3^{\prime}2^{\prime})$
$\displaystyle\times\wick{\c{1}f(1)\c{1}f^{\dagger}(2)\c{3}f^{\dagger}(3)\c{2}f(4)\c{1}f(5)\c{1}f^{\dagger}(5^{\prime})\c{2}f^{\dagger}(4^{\prime})\c{3}f(3^{\prime})\c{1}f(2^{\prime})\c{1}f^{\dagger}(1)}\
.$ (210)

The factor $2^{4}$ comes from equivalent contractions. It equals

$\displaystyle\delta\widetilde{\mathcal{G}}^{(2,1)}(1)$
$\displaystyle=-[\mathcal{G}^{0}(1)]^{2}\cdot\frac{1}{2}\cdot
T^{2}\sum_{345}\delta_{\omega_{1}+\omega_{3},\omega_{4}+\omega_{5}}\widetilde{\Gamma}^{0}(13;45)\widetilde{\Gamma}^{0}(54;31)\widetilde{\mathcal{G}}^{0}(3)\widetilde{\mathcal{G}}^{0}(4)\widetilde{\mathcal{G}}^{0}(5)$
(211)

The self-energy can be extracted as

$\displaystyle\widetilde{\Sigma}^{(2,1)}(\omega_{1})=-\frac{1}{2}T^{2}\sum_{I_{3,4,5}}|\widetilde{\Gamma}^{0}_{I_{1}I_{3};I_{4}I_{5}}|^{2}\sum_{\omega_{3,4}}\widetilde{\mathcal{G}}^{0}(\omega_{3})\widetilde{\mathcal{G}}^{0}(\omega_{4})\widetilde{\mathcal{G}}^{0}(\omega_{1}+\omega_{3}-\omega_{4})\
,$ (212)

where the interaction factor can be calculated as

$u^{2}=\frac{1}{2}\sum_{I_{3,4,5}}|\widetilde{\Gamma}^{0}_{I_{1}I_{2};I_{3}I_{4}}|^{2}=\frac{1}{16}\sum_{I_{1,2,3,4}}|\widetilde{\Gamma}^{0}_{I_{1}I_{3};I_{4}I_{5}}|^{2}=\frac{1}{16}\mathrm{tr}[\widetilde{\Gamma}^{0}\cdot\widetilde{\Gamma}^{0\dagger}]$
(213)

Here we regard $I_{1,2}$ and $I_{3,4}$ as the left and right matrix indices
respectively.
$\mathrm{tr}[\widetilde{\Gamma}^{0}\cdot\widetilde{\Gamma}^{0\dagger}]$ equals
the sum of squares of the eigenvalues. $\widetilde{\Gamma}^{0}$ is a
64-dimensional matrix. Its eigenspace with nonzero eigenvalues consists of the
56 anti-symmetric states. Each of the 28 states discussed in Sec. A.5 appears
twice in the eigenspace because states related by permutation have been
regarded physically identical. Suppose $\Phi_{1}$ is a two-electron state with
the energy $E$. Define $\Phi_{2}$ as the permutation of $\Phi_{1}$, it
represents the same state. The physical energy is given by the generalized
eigenvalue problem defined by the Hamiltonian matrix
$H_{nm}=\langle\Phi_{n}|\hat{H}|\Phi_{m}\rangle$ and overlap matrix
$S_{nm}=\langle\Phi_{n}|\Phi_{m}\rangle$. There must be
$H=E(\sigma_{0}-\sigma_{x})$ and $S=\sigma_{0}-\sigma_{x}$ such that the
generalized eigenvalue is $E$. However, the Hamiltonian matrix $H$ itself has
two eigenvalues: $2E$ and 0. Therefore, each of 28 two-electron state
contribute to an eigenvalue $2E$ and an eigenvalue 0. Thereby,

$\displaystyle u^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{16}\sum_{E}|2E|^{2}=\frac{1}{4}\left(8\widetilde{U}_{2}^{2}+8\widetilde{U}_{3}^{2}+4\widetilde{U}_{4}^{2}+2(\widetilde{U}_{1}-\widetilde{\mathcal{J}})^{2}+6(\widetilde{U}_{1}+\mathcal{J})^{2}\right)$
(214) $\displaystyle=$ $\displaystyle
2\widetilde{U}_{1}^{2}+2\widetilde{U}_{2}^{2}+2\widetilde{U}_{3}^{2}+\widetilde{U}_{4}^{2}+2\widetilde{U}_{1}\widetilde{\mathcal{J}}+2\widetilde{\mathcal{J}}^{2}\
.$

Now we have

$\displaystyle\widetilde{\Sigma}^{(2,1)}(\omega_{1})=$
$\displaystyle-u^{2}\int\frac{d\omega_{3}d\omega_{4}}{(2\pi)^{2}}\cdot\frac{1}{i\omega_{3}+i\widetilde{\Delta}_{0}\mathrm{sgn}(\omega_{3})-\widetilde{\epsilon}_{f}}\frac{1}{i\omega_{4}+i\widetilde{\Delta}_{0}\mathrm{sgn}(\omega_{4})-\widetilde{\epsilon}_{f}}\frac{1}{i(\omega_{1}+\omega_{3}-\omega_{4})+i\widetilde{\Delta}_{0}\mathrm{sgn}(\omega_{1}+\omega_{3}-\omega_{4})-\widetilde{\epsilon}_{f}}$
$\displaystyle=$ $\displaystyle-u^{2}\int\frac{d\nu}{2\pi}\cdot
F(\nu)\cdot\frac{1}{i(\omega_{1}+\nu)+i\widetilde{\Delta}_{0}\mathrm{sgn}(\omega_{1}+\nu)-\widetilde{\epsilon}_{f}}$
(215)

where $\nu=\omega_{3}-\omega_{4}$ and

$F(\nu)=\int\frac{d\omega}{2\pi}\cdot\frac{1}{i\omega+i\widetilde{\Delta}_{0}\mathrm{sgn}(\omega)-\widetilde{\epsilon}_{f}}\cdot\frac{1}{i(\omega-\nu)+i\widetilde{\Delta}_{0}\mathrm{sgn}(\omega-\nu)-\widetilde{\epsilon}_{f}}\
.$ (216)

By definition $F(\nu)=F(-\nu)$. We hence have

$\displaystyle F(\nu)=$
$\displaystyle\int_{-\infty}^{0}\frac{d\omega}{2\pi}\frac{1}{i\omega-i\widetilde{\Delta}_{0}-\widetilde{\epsilon}_{f}}\frac{1}{i(\omega-|\nu|)-i\widetilde{\Delta}_{0}-\widetilde{\epsilon}_{f}}+\int_{0}^{|\nu|}\frac{d\omega}{2\pi}\frac{1}{i\omega+i\widetilde{\Delta}_{0}-\widetilde{\epsilon}_{f}}\frac{1}{i(\omega-|\nu|)-i\widetilde{\Delta}_{0}-\widetilde{\epsilon}_{f}}$
$\displaystyle+$
$\displaystyle\int_{|\nu|}^{\infty}\frac{d\omega}{2\pi}\frac{1}{i\omega+i\widetilde{\Delta}_{0}-\widetilde{\epsilon}_{f}}\frac{1}{i(\omega-|\nu|)+i\widetilde{\Delta}_{0}-\widetilde{\epsilon}_{f}}\
,$ (217)

Define $z=i\omega$, the above integral along imaginary axes of $z$ can be
continuously deformed to integral along real axes:

$\displaystyle F(\nu)=$ $\displaystyle\int_{0}^{\infty}\frac{dz}{2\pi
i}\left(\frac{1}{z+i\widetilde{\Delta}_{0}-\widetilde{\epsilon}_{f}}-c.c.\right)\frac{1}{z-i|\nu|-i\widetilde{\Delta}_{0}-\widetilde{\epsilon}_{f}}+\int_{0}^{\infty}\frac{dz}{2\pi
i}\frac{1}{z+i|\nu|+i\widetilde{\Delta}_{0}-\widetilde{\epsilon}_{f}}\left(\frac{1}{z+i\widetilde{\Delta}_{0}-\widetilde{\epsilon}_{f}}-c.c.\right)$
$\displaystyle=$
$\displaystyle\frac{2}{\pi}\int_{0}^{\infty}dz\cdot\frac{\widetilde{\Delta}_{0}}{(z-\widetilde{\epsilon}_{f})^{2}+\widetilde{\Delta}_{0}^{2}}\cdot\frac{z-\widetilde{\epsilon}_{f}}{(z-\widetilde{\epsilon}_{f})^{2}+(|\nu|+\widetilde{\Delta}_{0})^{2}}=\frac{\widetilde{\Delta}_{0}}{\pi(\nu^{2}+2\widetilde{\Delta}_{0}|\nu|)}\ln\left(1+\frac{\nu^{2}+2\widetilde{\Delta}_{0}|\nu|}{\widetilde{\Delta}_{0}^{2}+\widetilde{\epsilon}_{f}^{2}}\right)$
(218)

It is direct to verify that
$F(0^{+})=F(0^{-})=\frac{\widetilde{\Delta}_{0}}{\pi(\widetilde{\Delta}_{0}^{2}+\widetilde{\epsilon}_{f}^{2})}=\frac{\sin^{2}\delta_{f}}{\pi\widetilde{\Delta}_{0}}$
is the density of states at the Fermi level. Due to Eq. 140, the leading order
$\omega_{1}$-dependence of $\widetilde{\Sigma}^{(2,1)}$ can be obtained as

$\displaystyle\partial_{\omega}\widetilde{\Sigma}^{(2,1)}(\omega)=$
$\displaystyle-i\cdot u^{2}\int\frac{d\nu}{2\pi}\cdot
F(\nu)\left([\widetilde{\mathcal{G}}^{0}(\nu)]^{2}+2\pi\cdot\frac{\sin^{2}\delta_{f}}{\pi\widetilde{\Delta}_{0}}\delta(\nu)\right)$
$\displaystyle=$
$\displaystyle-i\cdot\left(\frac{u\cdot\sin^{2}\delta_{f}}{\pi\widetilde{\Delta}_{0}}\right)^{2}-i\cdot
u^{2}\int_{0}^{\infty}\frac{d\nu}{2\pi}\cdot
F(\nu)\frac{2\widetilde{\epsilon}_{f}^{2}-2(\nu+\widetilde{\Delta}_{0})^{2}}{(\widetilde{\epsilon}_{f}^{2}+(\nu+\widetilde{\Delta}_{0})^{2})^{2}}$
(219)

We define the dimensionless function

$\xi^{2}(\delta_{f})=\sin^{4}\delta_{f}+\pi^{2}\int_{0}^{\infty}\frac{dx}{2\pi}\cdot\frac{1}{\pi(x^{2}+2x)}\ln\left(1+\frac{x^{2}+2x}{1+\cot^{2}\delta_{f}}\right)\frac{2\cot^{2}(\delta_{f})-2(1+x)^{2}}{[\cot^{2}(\delta_{f})+(1+x)^{2}]^{2}}$
(220)

then we can write $\partial_{\omega}\widetilde{\Sigma}^{(2,1)}(\omega)$ as

$\partial_{\omega}\widetilde{\Sigma}^{(2,1)}(\omega)=-i\left(\frac{u\cdot\xi(\delta_{f})}{\pi\widetilde{\Delta}_{0}}\right)^{2}$
(221)

We find

$\xi\left(\frac{\pi}{2}\right)=\sqrt{3-\frac{\pi^{2}}{4}}\approx
0.7298,\qquad\xi\left(\frac{\pi}{4}\right)=\sqrt{\frac{4-4\mathcal{C}+\pi-\pi\ln
2}{8}}\approx 0.4031\ ,$ (222)

where $\mathcal{C}\approx 0.915966$ is the Catalan’s constant. The result
$\xi(\pi/2)$ can also be found in Ref. Hewson (1993a). The total frequency
derivative $\widetilde{\Sigma}^{(2)}$ also contains the contribution from the
second diagram of Fig. 4(d), i.e., $-i\lambda_{0}^{(2)}$. Imposing
$\partial_{\omega}\widetilde{\Sigma}^{(2)}=0$ gives

$\lambda^{(2)}_{0}=-\left(\frac{u\cdot\xi(\delta_{f})}{\pi\widetilde{\Delta}_{0}}\right)^{2}\
.$ (223)

##  Appendix C Quasi-particle mean-field theory of the superconductivity

###  C.1 Effective BdG Hamiltonian

The free part of topological heavy fermion model is given by Song and Bernevig
(2022)

${H}_{0}=-\mu\sum_{\mathbf{R}\alpha\eta s}f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta s}+\sum_{\eta
s}\sum_{aa^{\prime}}\sum_{|\mathbf{k}|<\Lambda_{c}}(H^{(c,\eta)}_{aa^{\prime}}(\mathbf{k})-\mu\delta_{aa^{\prime}})c^{\dagger}_{\mathbf{k}a\eta
s}c_{\mathbf{k}a\eta s}+\sum_{\eta s\alpha
a}\sum_{|\mathbf{k}|<\Lambda_{c}}\left(e^{-\frac{|\mathbf{k}|^{2}\lambda^{2}}{2}}H^{(cf,\eta)}_{a\alpha}(\mathbf{k})c^{\dagger}_{\mathbf{k}a\eta
s}f_{\mathbf{k}\alpha\eta s}+h.c.\right)$ (224)

Here $\mu$ is the chemical potential, $c_{\mathbf{k}a\eta s}$ is the fermion
operator for the $c$-electron of the momentum $\mathbf{k}$, orbital $a$
($=1,2,3,4$), valley $\eta$ ($=\pm$), and spin $s$ ($=\uparrow,\downarrow$),
$f_{\mathbf{R}\alpha\eta s}$ is the fermion operator for the $f$-electron of
the site $\mathbf{R}$ orbital $\alpha$ (=1,2), valley $\eta$, and spin.
$\mathbf{R}$’s form the triangular lattice generated by the Bravais lattice
basis $\mathbf{a}_{\rm M1}=\frac{2\pi}{3k_{\theta}}(\sqrt{3},1)$,
$\mathbf{a}_{\rm M2}=\frac{2\pi}{3k_{\theta}}(-\sqrt{3},1)$, where
$k_{\theta}=2|\mathbf{K}|\sin\frac{\theta}{2}$, $\theta=1.08^{\circ}$ is the
first magic angle, $|\mathbf{K}|=1.703\mathring{A}^{-1}$. The triangular
lattice corresponds to the AA stacking regions of MATBG (Fig. 1(a)).
$f_{\mathbf{k}\alpha\eta
s}=\frac{1}{N_{M}}\sum_{\mathbf{R}}e^{-i\mathbf{R}\cdot\mathbf{k}}f_{\mathbf{R}\alpha\eta
s}$ is the $f$-electron operator at the momentum $\mathbf{k}$, with $N_{\rm
M}$ being the number of moiré cells and $\mathbf{k}$ taking values in the
moiré Brillouin zone spanned by $\mathbf{b}_{\rm
M1}=k_{\theta}(\frac{\sqrt{2}}{2},\frac{3}{2})$, $\mathbf{b}_{\rm
M2}=k_{\theta}(-\frac{\sqrt{2}}{2},\frac{3}{2})$. However, $c$-electron is
described by a continuous field that cannot be regularized to a lattice due to
the symmetry anomaly Song  _et al._ (2021). Thus, momentum $\mathbf{k}$ of
$c$-electrons is within a cutoff $\Lambda_{c}$ that will be sent to infinity
in the end Song and Bernevig (2022). The matrices $H^{(c,\eta)}$ and
$H^{(cf,\eta)}$ are defined by

$H^{(c,\eta)}(\mathbf{k})=\begin{pmatrix}0_{2\times 2}&v_{\star}(\eta
k_{x}\sigma^{0}+ik_{y}\sigma^{z})\\\ v_{\star}(\eta
k_{x}\sigma^{0}-ik_{y}\sigma^{z})&M\sigma^{x}\end{pmatrix},\qquad
H^{(cf,\eta)}(\mathbf{k})=\begin{pmatrix}\gamma\sigma^{0}+v_{\star}^{\prime}(\eta
k_{x}\sigma^{x}+k_{y}\sigma^{y})\\\ 0_{2\times 2}\end{pmatrix}\ .$ (225)

The parameter $\lambda$ is the spread of the Wannier functions of
$f$-electrons, and it suppresses the hybridization at large momentum. In this
work we adopt the parameters: $\gamma=-24.75$meV,
$v_{\star}=-4.303\mathrm{eV\cdot\mathring{A}}$,
$v_{\star}^{\prime}=1.623\mathrm{eV\cdot\mathring{A}}$,
$\lambda=1.4131/k_{\theta}$, $M=3.697$meV. The resulting band structure is
shown in the inset of Fig. 1(a).

The free Hamiltonian has time-reversal $T$ and $D_{6}$ point group, generated
by $C_{2z}$-rotation, $C_{3z}$-rotation, and $C_{2x}$-rotation symmetries.
Their representations on $f$-electrons and $c$-electrons are given by (Eq. 8)

$T=\sigma^{0}\tau^{x}\varsigma^{0}K,\qquad
C_{2z}=\sigma^{x}\tau^{x}\varsigma^{0},\qquad
C_{3z}=e^{i\frac{2\pi}{3}\sigma^{z}\tau^{z}\varsigma^{0}},\qquad
C_{2x}=\sigma^{x}\tau^{0}\varsigma^{0}\ ,$ (226)

and

$T=(\sigma^{0}\tau^{x}\varsigma^{0}\oplus\sigma^{0}\tau^{x}\varsigma^{0})K,\qquad
C_{2z}=\sigma^{x}\tau^{x}\varsigma^{0}\oplus\sigma^{x}\tau^{x}\varsigma^{0},\qquad
C_{3z}=e^{i\frac{2\pi}{3}\sigma^{z}\tau^{z}\varsigma^{0}}\oplus\sigma^{0}\tau^{0}\varsigma^{0},\qquad
C_{2x}=\sigma^{x}\tau^{0}\varsigma^{0}\oplus\sigma^{x}\tau^{0}\varsigma^{0}\
,$ (227)

respectively. Detailed derivations of these symmetries are given sections S2A
and S2B of the supplementary material Ref. Song and Bernevig (2022). Readers
may directly verify that Eq. 224 is invariant under these symmetries. In this
work we do not distinguish the single-particle and the second-quantized
representations of symmetry operators. The correspondence between the two
representations are given after Eq. 8.

Eq. 224 with Coulomb interaction is solved within the framework of dynamical
mean-field theory (DMFT) in Refs. Zhou  _et al._ (2023); Datta  _et al._
(2023); Rai  _et al._ (2023). Here we are only interested in the doped
correlated insulator at the filling $\nu=-2-\delta\nu$, where the
$\widetilde{\Delta}_{0}\ll\mathcal{J}$ is assumed. In the mean-field
calculations, we will regard $z$ as an additional input to $\delta\nu$. In the
$\widetilde{\Delta}_{0}\ll\mathcal{J}$ limit, the $f$-filling is frozen to
$\nu_{f}=-2$, indicating that the doping is fully undertaken by $c$-electrons,
i.e., $\nu_{c}=\nu-\nu_{f}=-\delta\nu$.

To describe the quasi-particles at $\nu=-2-\delta\nu$, we express $f$-electron
operators in terms of the quasi-particles operators as
$f_{\mathbf{R}\alpha\eta s}=z^{\frac{1}{2}}\widetilde{f}_{\mathbf{R}\alpha\eta
s}$. We also introduce the DMFT self-consistent on-site energies
$\epsilon_{f}$, $\epsilon_{c1}=\epsilon_{c2}$, $\epsilon_{c3}=\epsilon_{c4}$
for $\widetilde{f}$-electrons, $a=1,2$ $c-$electrons, $a=3,4$ $c$-electrons,
respectively Zhou  _et al._ (2023). Then the effective free Hamiltonian for
the quasi-particle can be written as

$\displaystyle{H}_{0}=$ $\displaystyle\epsilon_{f}\sum_{\mathbf{R}\alpha\eta
s}\widetilde{f}_{\mathbf{R}\alpha\eta
s}^{\dagger}\widetilde{f}_{\mathbf{R}\alpha\eta s}+\sum_{\eta
s}\sum_{aa^{\prime}}\sum_{|\mathbf{k}|<\Lambda_{c}}(H^{(c,\eta)}_{aa^{\prime}}(\mathbf{k})+\delta_{aa^{\prime}}\epsilon_{c,a})c^{\dagger}_{\mathbf{k}a\eta
s}c_{\mathbf{k}a\eta s}$ $\displaystyle+z^{\frac{1}{2}}\sum_{\eta s\alpha
a}\sum_{|\mathbf{k}|<\Lambda_{c}}\left(e^{-\frac{|\mathbf{k}|^{2}\lambda^{2}}{2}}H^{(cf,\eta)}_{a\alpha}(\mathbf{k})c^{\dagger}_{\mathbf{k}a\eta
s}\widetilde{f}_{\mathbf{k}\alpha\eta s}+h.c.\right)\ .$ (228)

One can see that the $cf$ hybridization is suppressed by $z^{\frac{1}{2}}$.
$\epsilon_{f}$, $\epsilon_{c,1}$, $\epsilon_{c,3}$ should be determined self-
consistently in the DMFT plus Hartree-Fock calculations. Here we provide a
quick approach to estimate $\epsilon_{f}$, $\epsilon_{c,1}$, $\epsilon_{c,3}$
for given $z$ and $\nu=-2-\delta\nu$. First, it is found that
$G=\epsilon_{c,3}-\epsilon_{c,1}$ is insensitive to doping and can be
approximated as a fixed quantity (taken as 8meV here) Zhou  _et al._ (2023).
We then search for $\epsilon_{f}$ and $\epsilon_{c,1}$ that give the correct
fillings $\nu_{f}=-2$ and $\nu_{c}=-\delta\nu$.

As detailed in Appendix B.5, we take the irreducible vertex $z^{2}\Gamma^{\rm
p}$ \- multiplied by $z^{2}$ \- in the pairing channel as the effective
interaction in mean-field calculation. (One would over-count Feynman diagrams
if $z^{2}\Gamma$ is used instead Georges  _et al._ (1996).) For the sake of
investigating superconducting orders, it suffices to keep the strongest
channel in $\Gamma^{\rm p}$, i.e., the $d$-wave $E_{2}$ states (Eq. 64). At
$\nu_{f}=-2$, the pairing potential is determined by the renormalized
interaction $\widetilde{\mathcal{J}}$ via (Eq. 195)

$U^{\rm p}_{1}-\mathcal{J}^{\rm
p}\approx\frac{-2\pi\widetilde{\Delta}_{0}-\widetilde{\mathcal{J}}}{1+\frac{\pi}{2}+\widetilde{\mathcal{J}}/(4\widetilde{\Delta}_{0})},\qquad
4.6\widetilde{\Delta}_{0}\lesssim\widetilde{\mathcal{J}}\lesssim
10.3\widetilde{\Delta}_{0}\ .$ (229)

The explicit form of the corresponding interaction can be constructed as
projectors to Eq. 64

$\displaystyle H_{\rm p}=$ $\displaystyle(U_{1}^{\rm p}-\mathcal{J}^{\rm
p})\sum_{\mathbf{R}}\sum_{\eta}\frac{(\widetilde{f}_{\mathbf{R}1\eta\uparrow}^{\dagger}\widetilde{f}_{\mathbf{R}2\overline{\eta}\downarrow}^{\dagger}-\widetilde{f}_{\mathbf{R}1\eta\downarrow}^{\dagger}\widetilde{f}_{\mathbf{R}2\overline{\eta}\uparrow}^{\dagger})(\widetilde{f}_{\mathbf{R}2\overline{\eta}\downarrow}\widetilde{f}_{\mathbf{R}1\eta\uparrow}-\widetilde{f}_{\mathbf{R}2\overline{\eta}\uparrow}\widetilde{f}_{\mathbf{R}1\eta\downarrow})}{2}$
$\displaystyle=$ $\displaystyle\frac{U_{1}^{\rm p}-\mathcal{J}^{\rm
p}}{2}\sum_{\mathbf{R}}\sum_{\eta}(\widetilde{f}_{\mathbf{R}2\overline{\eta}\downarrow}^{\dagger}\widetilde{f}_{\mathbf{R}1\eta\uparrow}^{\dagger}-\widetilde{f}_{\mathbf{R}2\overline{\eta}\uparrow}^{\dagger}\widetilde{f}_{\mathbf{R}1\eta\downarrow}^{\dagger})(\widetilde{f}_{\mathbf{R}1\eta\uparrow}\widetilde{f}_{\mathbf{R}2\overline{\eta}\downarrow}-\widetilde{f}_{\mathbf{R}1\eta\downarrow}\widetilde{f}_{\mathbf{R}2\overline{\eta}\uparrow})\
.$ (230)

We only consider the order parameters

$\Psi_{\eta}=\langle\Phi_{\rm
G}|\widetilde{f}_{\mathbf{R}1\eta\uparrow}\widetilde{f}_{\mathbf{R}2\overline{\eta}\downarrow}|\Phi_{\rm
G}\rangle\ .$ (231)

where $\Phi_{\rm G}$ is the ground state. Due to the spin-SU(2) symmetry
($i\varsigma^{y})$, there must be

$\Psi_{\eta}=-\langle\Phi_{G}|\widetilde{f}_{\mathbf{R}1\eta\downarrow}\widetilde{f}_{\mathbf{R}2\overline{\eta}\uparrow}|\Phi_{G}\rangle\
.$ (232)

The mean-field decomposition of $H_{\rm p}$ reads

$\displaystyle H_{\rm p}\approx$ $\displaystyle\frac{U^{\rm
p}_{1}-\mathcal{J}^{\rm
p}}{2}\sum_{\mathbf{R}\eta}\Big{[}2\Psi_{\eta}(\widetilde{f}^{\dagger}_{\mathbf{R}2\overline{\eta}\downarrow}\widetilde{f}^{\dagger}_{\mathbf{R}1\eta\uparrow}-\widetilde{f}^{\dagger}_{\mathbf{R}2\overline{\eta}\uparrow}\widetilde{f}^{\dagger}_{\mathbf{R}1\eta\downarrow})+2\Psi_{\eta}^{*}(\widetilde{f}_{\mathbf{R}1\eta\uparrow}\widetilde{f}_{\mathbf{R}2\overline{\eta}\downarrow}-\widetilde{f}_{\mathbf{R}1\eta\downarrow}\widetilde{f}_{\mathbf{R}2\overline{\eta}\uparrow})-2|\Psi_{\eta}|^{2}\Big{]}$
(233)

We can see that $(\eta,s)=(+,\uparrow)$ electrons only couples to
$(-,\downarrow)$ electrons, whereas $(-,\uparrow)$ electrons only couples to
$(+,\downarrow)$ electrons. Thus, the BdG Hamiltonian is block diagonal. We
can write the total Hamiltonian as

$\displaystyle H\approx$ $\displaystyle-N_{\rm M}(U^{\rm
p}_{1}-\mathcal{J}^{\rm
p})(|\Psi_{+}|^{2}+|\Psi_{-}|^{2})+\overline{H}_{+}+\overline{H}_{-}$ (234)

$\overline{H}_{+}$ is the BdG Hamiltonian for the $(\eta,s)=(+,\uparrow)$ and
$(-,\downarrow)$ electrons

$\displaystyle\overline{H}_{+}=$
$\displaystyle\sum_{\mathbf{k}}\begin{pmatrix}\widetilde{f}_{\mathbf{k}+\uparrow}^{\dagger}&c_{\mathbf{k}+\uparrow}^{\dagger}&\widetilde{f}_{-\mathbf{k}-\downarrow}&c_{-\mathbf{k}-\downarrow}\end{pmatrix}\begin{pmatrix}\widetilde{\epsilon}_{f}&z^{\frac{1}{2}}\mathcal{H}^{(cf,+)\dagger}(\mathbf{k})&\mathcal{V}&0\\\
z^{\frac{1}{2}}\mathcal{H}^{(cf,+)}(\mathbf{k})&\mathcal{H}^{(c,+)}(\mathbf{k})&0&0\\\
\mathcal{V}^{\dagger}&0&-\widetilde{\epsilon}_{f}&-z^{\frac{1}{2}}\mathcal{H}^{(cf,-)T}(-\mathbf{k})\\\
0&0&-z^{\frac{1}{2}}\mathcal{H}^{(cf,-)*}(-\mathbf{k})&-\mathcal{H}^{(c,-)T}(-\mathbf{k})\end{pmatrix}\begin{pmatrix}\widetilde{f}_{\mathbf{k}+\uparrow}\\\
c_{\mathbf{k}+\uparrow}\\\ \widetilde{f}_{-\mathbf{k}-\downarrow}^{\dagger}\\\
c_{-\mathbf{k}-\downarrow}^{\dagger}\end{pmatrix}$ (235)

where

$\mathcal{V}=-(U^{\rm p}_{1}-\mathcal{J}^{\rm p})\begin{pmatrix}0&\Psi_{+}\\\
\Psi_{-}&0\end{pmatrix}\ ,$ (236)

is a two-by-two matrix. Here $\widetilde{f}_{\mathbf{k}+\uparrow}$ and
$\widetilde{f}_{-\mathbf{k}-\downarrow}^{\dagger}$ are column vectors
consisting of $\widetilde{f}_{\mathbf{k}\alpha+\uparrow}$ and
$\widetilde{f}_{-\mathbf{k}\alpha-\downarrow}^{\dagger}$, respectively.
$c_{\mathbf{k}+\uparrow}$ and $c_{-\mathbf{k}-\downarrow}$ are column vectors
consisting of $c_{\mathbf{k}+\mathbf{G},a+\uparrow}$ and
$c_{-\mathbf{k}-\mathbf{G},a-\downarrow}$, respectively, where $\mathbf{k}$ is
in the moiré Brillouin zone and $\mathbf{G}$, regarded as a matrix index, is
moiré reciprocal lattice. The $\mathcal{H}^{(c,\eta)}$ and
$\mathcal{H}^{(cf,\eta)}$ matrices are defined as

$\mathcal{H}^{(c,\eta)}_{\mathbf{G}a,\mathbf{G}^{\prime}a^{\prime}}(\mathbf{k})=\delta_{\bf
G,G^{\prime}}\left(\epsilon_{c,a}\delta_{aa^{\prime}}+H^{(c,\eta)}_{aa^{\prime}}(\mathbf{k}+\mathbf{G})\right),\qquad\mathcal{H}^{(cf,\eta)}_{\mathbf{G}a,\alpha}(\mathbf{k})=e^{-\frac{|\mathbf{k}+\mathbf{G}|^{2}\lambda^{2}}{2}}H^{(cf,\eta)}_{a\alpha}(\mathbf{k}+\mathbf{G})\
.$ (237)

$\overline{H}_{-}$ is the BdG Hamiltonian for
$(\eta,s)=(+,\downarrow),(-,\uparrow)$ electrons. It can be obtained by
applying the $i\varsigma^{y}$ rotation to $\overline{H}_{+}$, and hence it has
identical spectrum as $\overline{H}_{+}$. The total energy (per moiré unit
cell) can be calculated as

$\frac{E_{\rm tot}}{N_{\rm M}}=-N_{\rm M}(U^{\rm p}_{1}-\mathcal{J}^{\rm
p})(|\Psi_{+}|^{2}+|\Psi_{-}|^{2})+\frac{2}{N_{\rm M}}\sum_{n\in{\rm
occ}}\sum_{\mathbf{k}}\mathcal{E}_{n\mathbf{k}}\ ,$ (238)

where $\mathcal{E}_{n\mathbf{k}}$ is the $n$-th BdG band of
$\overline{H}_{+}$, and “$\sum_{n\in{\rm occ}}$” means summation over the
negative branches.

There are two relevant saddle points of the pairings in the numerical
calculation:

$\mathcal{V}=e^{i\theta_{1}}\mathcal{V}_{1}(\cos\varphi\cdot\sigma^{x}+\sin\varphi\cdot\sigma^{y}),\quad
e^{i\theta_{2}}\mathcal{V}_{2}\frac{\sigma^{x}\pm i\sigma^{y}}{2}\ .$ (239)

Here $\mathcal{V}_{1}$ ($>0$), $\mathcal{V}_{2}$ ($>0$) represent the pairing
mean-field Hamiltonian terms

$e^{i\theta_{1}}\mathcal{V}_{1}\sum_{\mathbf{k}}\left[e^{-i\varphi}\widetilde{f}_{\mathbf{k}1+\uparrow}^{\dagger}\widetilde{f}_{-\mathbf{k}2-\downarrow}^{\dagger}+e^{i\varphi}\widetilde{f}_{\mathbf{k}2+\uparrow}^{\dagger}\widetilde{f}_{-\mathbf{k}1-\downarrow}^{\dagger}\right]-(\uparrow\leftrightarrow\downarrow)\
.$ (240)
$e^{i\theta_{2}}\mathcal{V}_{2}\sum_{\mathbf{k}}\widetilde{f}_{\mathbf{k}1+\uparrow}^{\dagger}\widetilde{f}_{-\mathbf{k}2-\downarrow}^{\dagger}-(\uparrow\leftrightarrow\downarrow)\quad\text{or}\quad
e^{i\theta_{2}}\mathcal{V}_{2}\sum_{\mathbf{k}}\widetilde{f}_{\mathbf{k}2+\uparrow}^{\dagger}\widetilde{f}_{-\mathbf{k}1-\downarrow}^{\dagger}-(\uparrow\leftrightarrow\downarrow)\
,$ (241)

respectively. $\theta_{1,2}$ are the U(1) phases spontaneously chosen in the
symmetry breaking.
$\sum_{\mathbf{k}}\widetilde{f}_{\mathbf{k}1+\uparrow}^{\dagger}\widetilde{f}_{-\mathbf{k}2-\downarrow}^{\dagger}-\widetilde{f}_{\mathbf{k}1+\downarrow}^{\dagger}\widetilde{f}_{-\mathbf{k}2-\uparrow}^{\dagger}$
and
$\sum_{\mathbf{k}}\widetilde{f}_{\mathbf{k}2+\uparrow}^{\dagger}\widetilde{f}_{-\mathbf{k}1-\downarrow}^{\dagger}-\widetilde{f}_{\mathbf{k}2+\downarrow}^{\dagger}\widetilde{f}_{-\mathbf{k}1-\uparrow}^{\dagger}$
transform as $d_{x^{2}-y^{2}}\pm id_{xy}$ orbitals under the symmetry
operations and form the two-dimensional representation $E_{2}$ of $D_{6}$ (Eq.
226). The linear combination in Eq. 240 breaks the $C_{3z}$ symmetry and
respects $C_{2z}$ and time-reversal ($T=\tau^{x}K$) symmetries. The linear
combinations in Eq. 241 break the time-reversal symmetry and preserve
$C_{6z}$. Thus, $\mathcal{V}_{1}$ and $\mathcal{V}_{2}$ correspond to the
nematic $d$-wave pairing and chiral $d$-wave pairing, respectively. As we will
show in next subsubsection, $\mathcal{V}_{1}$ leads to a $p$-wave-like nodal
structure in the BdG spectrum due to the Berry’s phase on Fermi surface.

###  C.2 The $p$-wave-like nodal structure in the nematic $d$-wave phase

We now show that the $\mathcal{V}_{1}$ pairing leads to a $p$-wave-like nodal
structure in the BdG spectrum. We consider the weak pairing limit and project
Eq. 235 onto the Fermi surface basis. As shown in Fig. 3(b), there are two
Fermi surfaces around the $\Gamma_{M}$ point, and both of them are dominated
by the $f$-orbitals. For either Fermi surface, we can approximate the Bloch
band state around the Fermi surface as

$\psi_{\mathbf{k}\eta
s}^{\dagger}=\sum_{\alpha}U^{(\eta)}_{\alpha}(\mathbf{k})\widetilde{f}^{\dagger}_{\mathbf{k}\alpha\eta
s}\ .$ (242)

The time-reversal $T=\tau^{x}K$ and $C_{2z}T=\sigma^{x}K$ symmetries (Eq. 226)
indicate

$T\psi_{\mathbf{k}+s}^{\dagger}T^{-1}=\psi_{-\mathbf{k}-s}^{\dagger}\quad\Rightarrow\quad
U_{\alpha}^{(-)}(-\mathbf{k})=U_{\alpha}^{(+)*}(\mathbf{k})$ (243)
$C_{2z}T\psi_{\mathbf{k}+s}^{\dagger}(C_{2z}T)^{-1}=\psi_{\mathbf{k}+s}^{\dagger}e^{i\phi_{\mathbf{k}}}\quad\Rightarrow\quad
U_{2}^{(+)}(\mathbf{k})=U_{1}^{(+)*}(\mathbf{k})e^{-i\phi_{\mathbf{k}}}\ .$
(244)

Without loss of generality, we can choose the gauge

$U^{(+)}(\mathbf{k})=\frac{1}{\sqrt{2}}\begin{pmatrix}1\\\
e^{-i\phi_{\mathbf{k}}}\end{pmatrix},\qquad
U^{(-)}(-\mathbf{k})=\frac{1}{\sqrt{2}}\begin{pmatrix}1\\\
e^{i\phi_{\mathbf{k}}}\end{pmatrix}\ .$ (245)

According to Ref. Ahn  _et al._ (2019), $\phi_{\mathbf{k}}$ must have an odd
winding number along the Fermi surface, i.e.,

$\mathcal{W}=\frac{1}{2\pi}\oint_{\rm
FS}d\mathbf{k}\cdot\partial_{\mathbf{k}}\phi_{\mathbf{k}}=1\mod 2\ ,$ (246)

if the Fermi surface encloses an odd number of Dirac points. In the symmetric
Fermi liquid, $\phi_{\mathbf{k}}$ winding number is $\pm 1$, as shown in Fig.
3(b). Now we project the $\mathcal{V}_{1}$ pairing onto the Fermi surface

$\displaystyle\sum_{\mathbf{k}}e^{i\theta_{1}}\mathcal{V}_{1}\left(e^{-i\varphi}\widetilde{f}_{\mathbf{k}1+\uparrow}^{\dagger}\widetilde{f}_{-\mathbf{k}2-\downarrow}^{\dagger}+e^{i\varphi}\widetilde{f}_{\mathbf{k}2+\uparrow}^{\dagger}\widetilde{f}_{-\mathbf{k}1-\downarrow}^{\dagger}\right)-(\uparrow\leftrightarrow\downarrow)$
$\displaystyle=$
$\displaystyle\sum_{\mathbf{k}}2e^{i\theta_{1}}\mathcal{V}_{1}\cdot\cos(\phi_{\mathbf{k}}+\varphi)\cdot\psi_{\mathbf{k}+\uparrow}^{\dagger}\psi_{-\mathbf{k}-\downarrow}^{\dagger}-(\uparrow\leftrightarrow\downarrow)\
.$ (247)

As $\phi_{\mathbf{k}}$ winds $\pm 2\pi$ along the Fermi surface, there will be
two nodes.

The nodes are robust against perturbations that preserve the $C_{2z}T$
symmetry. For example, we consider a mixture of $s$-wave pairing and the
nematic $d$-wave pairing. Due to the $D_{6}$ point group symmetry, the
$s$-wave pairing must have the form
$e^{i\theta_{0}}\cdot\mathcal{V}_{0}\sigma^{0}$ ($\mathcal{V}_{0}>0$). The
projected pairing on Fermi surface becomes

$\displaystyle\sum_{\mathbf{k}}e^{i\theta_{0}}\mathcal{V}_{0}\left(\widetilde{f}_{\mathbf{k}1+\uparrow}^{\dagger}\widetilde{f}_{-\mathbf{k}1-\downarrow}^{\dagger}+\widetilde{f}_{\mathbf{k}2+\uparrow}^{\dagger}\widetilde{f}_{-\mathbf{k}2-\downarrow}^{\dagger}\right)+e^{i\theta_{1}}\mathcal{V}_{1}\left(e^{-i\varphi}\widetilde{f}_{\mathbf{k}1+\uparrow}^{\dagger}\widetilde{f}_{-\mathbf{k}2-\downarrow}^{\dagger}+e^{i\varphi}\widetilde{f}_{\mathbf{k}2+\uparrow}^{\dagger}\widetilde{f}_{-\mathbf{k}1-\downarrow}^{\dagger}\right)-(\uparrow\leftrightarrow\downarrow)$
$\displaystyle=$
$\displaystyle\sum_{\mathbf{k}}2\left(e^{i\theta_{0}}\mathcal{V}_{0}+e^{i\theta_{1}}\mathcal{V}_{1}\cdot\cos(\phi_{\mathbf{k}}+\varphi)\right)\psi_{\mathbf{k}+\uparrow}^{\dagger}\psi_{-\mathbf{k}-\downarrow}^{\dagger}-(\uparrow\leftrightarrow\downarrow)\
.$ (248)

The $C_{2z}T$ symmetry requires $\theta_{0}=\theta_{1}$ mod $\pi$ such that
the pairing term gains a uniform phase $e^{-i2\theta_{0}}$ under $C_{2z}T$.
There are two pairing nodes on the Fermi surface as long as
$\mathcal{V}_{0}<\mathcal{V}_{1}$.

The above discussions are based on the weak pairing presumption. If the
pairing energy is much larger than the band separation, then the nodes on
inner and outer Fermi surfaces will merge each other, leading to a gapped
phase.

Refs. Wu  _et al._ (2018); Liu  _et al._ (2023) reported that the $d$-wave
pairing has four nodes on the Fermi surface. The difference is because Refs.
Wu  _et al._ (2018); Liu  _et al._ (2023) used the bare band structure of
MATBG, where $\phi_{\mathbf{k}}$ has a winding number 2.