# Spontaneous spin superconductor state in ABCA-stacked tetralayer graphene

Shuai Li<EMAIL_ADDRESS>School of Physical Science and Technology,
Soochow University, Suzhou 215006, China.  Institute for Advanced Study,
Soochow University, Suzhou 215006, China.  Yuan-Hang Ren  School of Physical
Science and Technology, Soochow University, Suzhou 215006, China.  Ao-Long Li
State Key Laboratory of Surface Physics and Department of Physics, Fudan
University, Shanghai 200433, China  Interdisciplinary Center for Theoretical
Physics and Information Sciences (ICTPIS), Fudan University, Shanghai 200433,
China.  Hua Jiang<EMAIL_ADDRESS>Interdisciplinary Center for
Theoretical Physics and Information Sciences (ICTPIS), Fudan University,
Shanghai 200433, China.

###### Abstract

We theoretically demonstrate a spontaneous spin superconductor (SC) state in
ABCA-stacked tetralayer graphene, under sequential effects of electron-
electron (e-e) and electron-hole (e-h) interactions. First of all, we examine
the ferromagnetic (FM) exchange instability and phase diagram of the system
induced by the long-range e-e interaction. At non- or low-doping levels, the
interaction trends to stabilize a FM phase with the coexisting electron and
hole carriers. Superior to bilayer and trilayer systems, tetralayer graphene
has a larger FM phase region and spin splitting, making it more advantageous
to realize the spin SC state. Subsequently, we prove that the FM phase becomes
unstable when attractive e-h interaction is considered. As a consequence, the
spin SC state can be spontaneously formed at low temperature, where spin-
triplet exciton pairs act as the equivalent of Cooper pairs. We further
develop a consistent BCS-type theory for the spin SC state in ABCA-stacked
graphene. The predicted spin superconducting gap can reach about $7.0$ meV,
with a critical temperature of about 45 K for non-doping system. At last, we
demonstrated a spin-current Josephson effect in the ABCA-stacked graphene spin
SC heterojunction. Our findings enrich the prospective spin SC candidate
materials, illuminating more possibilities for achieving non-dissipative
super-spintronics.

##  I Introduction

Ever since its discovery[1], the development of superconductivity has
witnessed remarkable advancements[1, 2, 3, 4, 5, 6, 7, 8, 9, 10].
Theoretically, the physics of superconductors (SCs) are well understood by the
fundamental BCS theory: Electrons combine into spin-singlet Cooper pairs, then
condense into a charge superfluid at low temperature supporting persistent
charge supercurrents[2]. On the other hand, searching for spin superfluid
possessing spin-polarized supercurrents is a key aspect in spintronics[11, 12,
13]. The combination of superconductivity and spintronics opens a new research
direction: super-spintronics, and offers a promising avenue for the
application of new types of spintronic devices with a substantial decrease in
power consumption[14, 15].

The unusual ferromagnetic (FM) SC is believed to hold a spin superfluid with
equal-spin triplet paring[16, 17, 18, 19, 20]. However, this intrinsic triplet
pairing state is rare and only identified in few materials[18, 19, 20, 21, 22,
23]. Fortunately, the spin-triplet excitonic condensate, dubbed as the spin
SC, constitutes another prospective candidate for the realization of spin
superfluid and non-dissipative spin transport [24, 25, 26, 27, 28, 29, 30, 31,
32, 33, 34, 35, 36]. Although the spin SC is a charge insulator, a spin
supercurrent can flow without dissipation, i.e. zero spin resistance. In
addition, the London-type equation[24, 28] and Ginzburg-Landau-type theory[29]
for spin SCs are well established, unveiling intriguing electric ”Meissner
effect” and the spin-current Josephson effect. Interplayed with topology, the
triplet spin SC might also be exploited as a topological SC bearing potentials
for topological quantum computing without using Majorana fermions[37, 38].

Several promising materials are predicted to host spin-triplet excitonic
condensate, including monolayer and trilayer graphene[24, 25, 26, 27], Kagome
lattice systems[30, 31], perovskites[33, 34], and topological insulators[35,
36]. In these systems, spin splitting and the presence of dual carriers can be
induced externally by magnetic fields and ferromagnets, or intrinsically
through Coulomb interactions. Subsequently, the spin-triplet exciton pairs
could form and condense into a spin superfluid. However, the experimental
confirmation of the spin SC states confronts challenges stemming from inherent
limitations, such as short exciton lifetimes resulted from electron-hole (e-h)
recombination[39, 40], inadequate spin split by weak magnetic exchange
interaction[41], and difficulties in material synthesis [30, 31]. These
constraints complicate the experimental validation of the spin SCs.
Encouragingly, the recently reported ABCA-stacked tetralayer (ABCA) graphene
[42, 43, 44] might serve as a promising platform to overcome these drawbacks.
Compared to bilayer and trilayer graphene, its notable flatter bands and more
quenched kinetic energy can possibly enhance the effects of Coulomb
interaction[45], leading to a larger intrinsic spin split. Additionally, the
generated hole carriers reside above the electronic states, thereby inhibiting
e-h recombination and endowing excitons with long lifetimes[24]. Therefore,
the benefits of ABCA graphene could possibly make it an exceptional platform
for achieving spin superconductivity. Given the experimental accessibility of
high-quality samples, this system merits further theoretical and experimental
explorations of the triplet spin SC state.

In this paper, we theoretically demonstrate that the ABCA graphene can realize
a spontaneous spin SC state by sequentially considering the intrinsic
interactions. At first, we study the FM exchange instability induced by long-
range electron-electron (e-e) interaction using a variational method. In
comparison with the bilayer and trilayer graphene, the phase diagram of
tetralayer system shows an obviously lager FM phase region as a result of its
divergent density of states (DOS) near the Dirac points. For non-doping or
weakly doped systems, the interaction trends to stabilize a FM phase with the
coexistence of electron and hole carriers.

Subsequently, we verify that the FM phase will not be stable with the
consideration of attractive e-h interaction. Electrons and holes spontaneously
bind into spin-triplet exciton pairs and then condense into a spin SC state. A
consistent BCS-type theory of the spin SC state is developed. For non-doping
system, the estimated spin superconducting gap and critical temperature can
reach about $7.0$ meV and $45$ K, respectively. Given the successful
fabrication of ABCA graphene in recent experiments, we expect this spin SC
state to be identified experimentally by finely tuning the system parameters.
Finally, we demonstrate the behavior of spin-current Josephson effect in the
ABCA graphene spin SC heterojunction.

The remainder of this paper is organized as: In Sec. II, we describe the model
Hamiltonian and the variational wave function method to study exchange
instability. In Sec. III, we study the FM instability of the system and derive
its phase diagram. Sec. IV is devoted to the verification of instability of
the FM phase in the presence of e-h interaction. A BCS-type theory is further
established for the spin SC state in ABCA graphene. In Sec. V, we study the
related spin-current Josephson effect. A brief conclusion of our work is given
in Sec. VI. All supplementary materials are relegated into the appendices.

##  II Model and variational method

The lattice structure of ABCA graphene in one primitive cell is illustrated in
Fig. 1(a). We model the kinetic energy of the system by a tight-binding
Hamiltonian

$H_{\mathrm{tb}}=-t\sum_{\langle ij\rangle
l\sigma}c^{\dagger}_{A_{l}i\sigma}c_{B_{l}j\sigma}-\gamma\sum_{il\sigma}c^{\dagger}_{B_{l}i\sigma}c_{A_{l+1}i\sigma}+\mathrm{H.c.},$
(1)

where the operator $c_{\alpha_{l}i\sigma}$ annihilates an electron with spin
$\sigma(\uparrow,\downarrow)$ at the site $i$ of sublattice
$\alpha(\mathrm{A,B})$ on layer $l(1,\cdots,4$), $t\approx 3$ eV and
$\gamma\approx 0.35$ eV are the nearest neighbor (NN) in-plane and out-of-
plane hoppings, respectively[46, 25]. Perform an expansion of
$H_{\mathrm{tb}}$ around $\bm{K}_{+}$ or $\bm{K}_{-}$ point at the corners of
the Brillouin zone (BZ), then the low-energy kinetic Hamiltonian can be
obtained as

$H_{\rm
kin}=\sum_{\bm{p}}\sum_{\eta=+,-}\sum_{\sigma=\uparrow,\downarrow}\Psi^{\dagger}_{\eta\sigma}(\bm{p})K_{\eta}(\bm{p})\Psi_{\eta\sigma}(\bm{p}).$
(2)

Here
$\Psi_{\eta\sigma}(\bm{p})=[\psi_{A_{1}\eta\sigma}(\bm{p}),\psi_{B_{1}\eta\sigma}(\bm{p}),\cdots,\psi_{B_{4}\eta\sigma}(\bm{p})]^{T}$,
where $\psi_{\alpha\eta\sigma}(\bm{p})$ is the operator that annihilates an
electron with spin $\sigma(\uparrow,\downarrow)$ on sublattice
$\alpha(A_{1},\cdots,B_{4})$ of a small momentum $\bm{p}$ measured from valley
$\bm{K}_{\eta}$ ($\eta=\pm$).

$K_{\eta}(\bm{p})=\left(\begin{array}[]{cccc}H^{0}_{\eta}(\bm{p})&t_{\bot}&0&0\\\
t_{\bot}^{\dagger}&H^{0}_{\eta}(\bm{p})&t_{\bot}&0\\\
0&t_{\bot}^{\dagger}&H^{0}_{\eta}(\bm{p})&t_{\bot}\\\
0&0&t_{\bot}^{\dagger}&H^{0}_{\eta}(\bm{p})\\\ \end{array}\right)$ (3)

gives the valley dependent Hamiltonian in momentum space, where
$t_{\bot}=(0~{}0;-\gamma~{}0)$ is the interlayer coupling matrix and
$H^{0}_{\eta}(\bm{p})$ denotes the Hamiltonian of single layer graphene at
long wavelength limit, truncated to the first order in $\bm{p}$ [47]

$H^{0}_{\eta}(\bm{p})=\left(\begin{array}[]{cc}0&\hbar\upsilon_{F}pe^{i\eta\theta(\bm{p})}\\\
\hbar\upsilon_{F}pe^{-i\eta\theta(\bm{p})}&0\\\ \end{array}\right),$ (4)

with $\tan\theta(\bm{p})=\frac{p_{y}}{p_{x}}$ and
$\hbar\upsilon_{F}=\frac{3ta_{0}}{2}$ ($a_{0}\approx 0.142$ nm is the nearest
carbon-carbon distance).

Figure 1: (Color online) (a) Schematic lattice structure of ABCA-stacked
tetralayer graphene. (b) Numerical energy bands (left) along the $p_{x}$ axis
near a valley $\bm{K}_{\eta}$ (set as $p=0$) and the corresponding DOS (right)
calculated by diagonalizing $H_{\mathrm{tb}}$ and $H_{\mathrm{kin}}$,
respectively.

By solving the Schrödinger equation
$K_{\eta}(\bm{p})M_{\eta}(\bm{p})=M_{\eta}(\bm{p})\varepsilon(\bm{p})$, the
kinetic Hamiltonian is diagonalized as

$H_{\rm
kin}=\sum_{\bm{p}\eta\sigma}\Phi^{\dagger}_{\eta\sigma}(\bm{p})\varepsilon(\bm{p})\Phi_{\eta\sigma}(\bm{p})=\sum_{\bm{p}m\eta\sigma}\varepsilon_{m}(\bm{p})\phi^{\dagger}_{\bm{p}m\eta\sigma}\phi_{\bm{p}m\eta\sigma},$
(5)

where
$\Phi_{\eta\sigma}(\bm{p})=M^{\dagger}_{\eta}(\bm{p})\Psi_{\eta\sigma}(\bm{p})=[\cdots,\phi_{\bm{p}m\eta\sigma},\cdots]^{T}$
with $m$ denoting the band index. We numerically plot the energy bands of ABCA
graphene near $\bm{K}_{\eta}$ in Fig. 1(b), where the lowest-energy conduction
and valence bands (labeled 1 and 2) intersect at the Dirac point ($E=0$) with
quartic dispersion
$\varepsilon_{1,2}(\bm{p})\approx\pm\frac{(\hbar\upsilon_{F}p)^{4}}{\gamma^{3}}$
[48].

The low-energy Hamiltonian $H_{\mathrm{kin}}$ is only valid when the momentum
$p$ is smaller than some critical value $k_{c}$ [49, 50]. Thus, we have to
guarantee that the physics of ABCA graphene should mostly depend on the low-
energy region and treat high-energy bands as non-important approximations. The
reason why $H_{\mathrm{kin}}$ is used instead of $H_{\mathrm{tb}}$ is that it
endows a rigid rotation symmetry around valleys $\bm{K}_{\eta}$, thereby
making numerical calculations easier to implement. To explain that
$H_{\mathrm{kin}}$ can depict the low-energy physics as well as
$H_{\mathrm{tb}}$, we numerically calculate their energy bands and DOS, as
shown in Fig. 1(b). In the region of $|E|\leq 0.5\gamma$ and $pa_{0}\leq 0.1$,
the bands (labeled 1 and 2) and DOS given by $H_{\mathrm{kin}}$ coincide well
with that given by $H_{\mathrm{tb}}$, and are accurate in tolerance within 2%.
Since the Fermi energy and momentum invoked in this work are well located in
this region, $H_{\mathrm{kin}}$ is then applicable in calculations.

Below we study the kinetic energy density of a given variational state, which
can be calculated by

$\frac{E_{\rm kin}}{S}=\frac{\langle
H_{\rm{kin}}\rangle}{S}=\int_{D}\frac{d^{2}\bm{p}}{(2\pi)^{2}}~{}\sum_{m\eta\sigma}\varepsilon_{m}(\bm{p})n_{m\eta\sigma}(\bm{p}),$
(6)

where $S$ represents the area of the system,
$n_{m\eta\sigma}(\bm{p})=\langle\phi^{\dagger}_{\bm{p}m\eta\sigma}\phi_{\bm{p}m\eta\sigma}\rangle=0$,
1, $\Theta(Q_{\sigma}-p)$, or $1-\Theta(Q_{\sigma}-p)$ when the band $m$ with
spin $\sigma$ at valley $\bm{K}_{\eta}$ is empty, fully occupied, partially
filled by electrons, or by holes to momentum $Q_{\sigma}$ (pocket size). Here
$Q_{\sigma}$ is the variational parameter, and $\Theta$ is the Heaviside step
function. The integral is performed over a disk area $D$ centered at
$\bm{K}_{\eta}$ with radius $\Lambda=\sqrt{\frac{4\pi}{3\sqrt{3}a^{2}_{0}}}$,
which ensures that the number of states is conserved in the first BZ [46].

Next we consider the exchange interaction between electrons. Because only
weakly doped system is studied, the Coulomb interaction should be slightly
screened and thus long ranged [51, 52]. Therefore, the e-e interaction can be
given by
$V(\bm{r}-\bm{r}^{\prime})=\frac{e^{2}}{\epsilon|\bm{r}-\bm{r}^{\prime}|}$,
where $\epsilon$ is the dielectric constant of the substrate ($\epsilon\approx
1$ for suspended samples). In the framework of second quantization, the e-e
exchange interaction can be written in terms of field operators
$\hat{\varphi}(\bm{r})$ as

$H_{\rm ex}=\frac{1}{2}\int\int
d\bm{r}d\bm{r}^{\prime}~{}\hat{\varphi}^{\dagger}(\bm{r})\hat{\varphi}^{\dagger}(\bm{r}^{\prime})V(\bm{r}-\bm{r}^{\prime})\hat{\varphi}(\bm{r}^{\prime})\hat{\varphi}(\bm{r}).$
(7)

Expanding $\hat{\varphi}(\bm{r})$ in the basis of orbital states depicted by
$\psi_{\alpha a\sigma}(\bm{p})$, $H_{\rm ex}$ can be recast to (see Appendix A
for derivation) [46, 25]

$H_{\rm
ex}=\frac{1}{2S}\sum_{\bm{q}\neq\bm{0}}\sum^{4}_{i,j=1}\hat{\rho}_{i}(\bm{q})V_{ij}(\bm{q})\hat{\rho}_{j}(-\bm{q}),$
(8)

where $V_{ij}(\bm{q})=2\pi e^{2}e^{-qd_{ij}}/\epsilon q=2\pi
g\hbar\upsilon_{F}e^{-qd_{ij}}/q$ gives the interaction between layers $i$ and
$j$, with $d_{ij}$ the vertical distance between the two layers and
$g=e^{2}/(\epsilon\hbar\upsilon_{F})$ defined as the dimensionless interaction
governing the system and
$\hat{\rho}_{i}(\bm{q})=\sum_{\alpha\bm{p}a\sigma}\psi^{\dagger}_{\alpha_{i}a\sigma}(\bm{p}+\bm{q})\psi_{\alpha_{i}a\sigma}(\bm{p})$
is the Fourier component of layer charge density $\hat{\rho}_{i}(\bm{r}_{i})$.

The presence of long-range interaction described by Eq. (8) can lead to the
exchange instability in the noninteracting ground state of ABCA graphene. We
will utilize a variational wave function technique to study the FM
instability, which is important in 2D electron gas with diluted carriers[53,
46, 25, 50, 54]. First, it is convenient to convert $\hat{\rho}_{i}(\bm{q})$
into the diagonalized basis by

$\hat{\rho}_{i}(\bm{q})=\sum_{\bm{p}\eta\sigma}\Phi^{\dagger}_{\eta\sigma}(\bm{p}+\bm{q})~{}\chi^{i}(\bm{p}+\bm{q},\bm{p})~{}\Phi_{\eta\sigma}(\bm{p})$
(9)

where

$\chi^{i}(\bm{p}+\bm{q},\bm{p})=M_{\eta}(\bm{p}+\bm{q})^{\dagger}\mathrm{diag}(\cdots,0,1,1,0,\cdots)~{}M_{\eta}(\bm{p}).$
(10)

Then substituting $\hat{\rho}_{i}(\bm{q})$ in Eq. (8) with Eq. (9) yields the
exchange energy density for a chosen variational state

$\begin{split}&\frac{E_{\mathrm{ex}}}{S}=\frac{\langle
H_{\mathrm{ex}}\rangle}{S}\approx-\frac{1}{2}\int_{D}\frac{d^{2}\bm{p}}{(2\pi)^{2}}\int_{D}\frac{d^{2}\bm{q}}{(2\pi)^{2}}\sum_{i,j}\sum_{m,n}\sum_{\eta,\sigma}\\\
&\chi^{i}_{mn,\eta}(\bm{q},\bm{p})\chi^{j}_{nm,\eta}(\bm{p},\bm{q})V_{ij}(\bm{p}-\bm{q})n_{m\eta\sigma}(\bm{q})n_{n\eta\sigma}(\bm{p}),\end{split}$
(11)

where the scattering between states of the two different valleys has been
dropped since their momenta difference $\bm{p}-\bm{q}$ is very large to render
$V_{ij}(\bm{p}-\bm{q})$ to be extremely small and thus safely neglected. Note
that $E_{\mathrm{kin}}/S$ and $E_{\mathrm{ex}}/S$ are all functions of the
variational parameter $Q_{\sigma}$, as shown by Eqs. (6) and (11).

The total energy density for the ABCA graphene is given by
$E/S=(E_{\mathrm{kin}}+E_{\mathrm{ex}})/S$. To determine the pocket size of
the equilibrium state, we have to minimize $E/S$ with respect to $Q_{\sigma}$
to find $Q_{\sigma,\mathrm{min}}$, for which the energy density is minimized.

Figure 2: (Color online) (a) Half-filling case. Left: noninteracting ground
state with zero pockets. Right: variational state with one electron pocket
$Q_{\uparrow}$ in spin up and one hole pocket $Q_{\downarrow}$ in spin down
channels. (b-c) Low-doping cases. Left: noninteracting ground state with two
electron pockets $Q_{0\uparrow}=Q_{0\downarrow}$. Right in (b): variational
state with two asymmetric electron pockets $Q_{\uparrow}>Q_{\downarrow}>0$;
Right in (c): variational state with one electron pocket $Q_{\uparrow}$ and
one hole pocket $Q_{\downarrow}$. At the two valleys $\bm{K}_{+}$ and
$\bm{K}_{-}$, the pockets are degenerate.

##  III FM exchange instability and phase diagram

As mentioned above, the FM exchange instability plays important roles in
electronic systems with low density of carriers[54]. In this section, we will
study the FM instability in weakly doped ABCA graphene with long-range Coulomb
interaction and acquire its phase diagram at zero temperature, utilizing the
variational method.

First, let us start with a non-doping system with $n=0$ (here $n$ is the
number of carriers in one primitive cell doped away from the Dirac point). In
the noninteracting ground state, the bands of the two spin channels are filled
to the Dirac point with zero pockets $Q_{0\uparrow}=Q_{0\downarrow}=0$. Under
e-e interaction, the FM variational state can be chosen with one electron
pocket in the spin-up branch ($Q_{\uparrow}>0$) and one hole pocket in the
spin-down branch ($Q_{\downarrow}<0$) at each valley, such that the sizes of
the pockets satisfy $Q_{\uparrow}=-Q_{\downarrow}=Q$, because of the
conservation of particle number in the non-doping system. The noninteracting
ground state and the variational trial state are schematically shown in Fig.
2(a). To gain the properties of equilibrium state, we calculate the energy
density variations $\Delta E_{\mathrm{kin}}/S$ and $\Delta E_{\mathrm{ex}}/S$
relative to noninteracting ground state as functions of $Q$, which are potted
in Fig. 3(a). Analytically, the kinetic energy cost by an electron (hole)
pocket is given by the integral $\int_{0}^{E(Q)}\rho(E)EdE$. Since the DOS for
the quartic band is $\rho(E)\sim 1/\sqrt{E}$ and $E\sim Q^{4}$, hence $\Delta
E_{\mathrm{kin}}/S\sim Q^{6}$, as consistent with the numerical results in
Fig. 3(a). On the other hand, $\Delta E_{\mathrm{ex}}/S$ exhibits a
nonmonotonic behavior with $Q$, which first decreases to the minimum and then
goes up. Therefore, the sum of $\Delta E_{\mathrm{kin}}/S$ and $\Delta
E_{\mathrm{ex}}/S$ gives rise to $\Delta E/S$ showing a minimum at some
specific $Q_{\mathrm{min}}$. For the unscreened system ($\epsilon=1$ and
$g=2.1$ inherent of the ABCA graphene), $\Delta E/S$ reaches the minimum
[labeled by the blue dot in Fig. 3(a)] with a spin-split energy
$\hbar\upsilon_{F}Q_{\mathrm{min}}\approx 0.46\gamma$ of the equilibrium
state. Therefore, the e-e interaction stabilizes a spin-polarized FM phase in
the non-doping ABCA graphene where electron and hole pockets form
simultaneously, with a magnetization $m=n_{\uparrow}-n_{\downarrow}\approx
1.10\times 10^{-3}g_{s}\mu_{B}$ in one cell. The spin-split energy and
magnetization significantly exceed that reported in bilayer
($\hbar\upsilon_{F}Q_{\mathrm{min}}\approx 0.05\gamma$ [46]) and trilayer
($\hbar\upsilon_{F}Q_{\mathrm{min}}\approx 0.09\gamma$ [25]) graphene, which
arises from the flatter energy bands ($\sim k^{4}$) and divergent DOS of the
tetralyer graphene system near the Dirac point.

In realistic experiment, the interaction $g$ can be flexibly tuned by external
substrates[55, 56]. Here, we explore the effect of $g$ on the equilibrium
states. As shown in Fig. S1 in Appendix B, $\Delta E_{\mathrm{ex}}/S$ is
lowered while $\Delta E_{\mathrm{kin}}/S$ keeps invariant as $g$ increases.
The pocket size $Q_{\mathrm{min}}$ of the equilibrium state gradually goes up
and saturates for large enough $g$. Since FM magnetization $m\propto
Q^{2}_{\mathrm{min}}$, $m$ displays a similar behavior with $g$ as shown by
the blue line in Fig. 3(d). For $g=3.0$, $m$ can reach a notable value
$1.10\times 10^{-3}g_{s}\mu_{B}$, about two orders larger than that in bilayer
graphene [46]. It is now clear from above study that a FM phase is always
energetically favorable for any small $g$ in the non-doping ABCA graphene
system. This is also confirmed and illustrated by the $n=0$ line in the phase
diagram shown in Fig. 3(c).

The magnetization of equilibrium states should not have a spin preference,
since the two spin channels are degenerate and scattering between them are
forbidden [see Eq. (11)]. We present a quantitative study in Appendix B to
elucidate that the spin-flipped state is degenerate with the FM state in Fig.
2(a) at equilibrium. Note that the selected FM variational state above
preserves $\mathrm{Z}_{2}$ symmetry but with broken $\mathrm{SU}(2)$ symmetry,
hence it is spin polarized[46]. There is another possible state which breaks
both symmetries but without net magnetization: The spin up (down) channels at
the two valleys have reversed pockets and the two spin channels are opposite
as well at the same valley [see Fig. S2(c)]. Because the intervalley
scattering is neglected, this nonmagnetic (NM) state is also degenerate with
the FM state in Fig. 2(a). Thus, it is necessary to incorporate the
intervalley scattering contribution to determine which state is more
energetically favorable. A detailed discussion is given in Appendix B, where
we find that the FM phase has a smaller energy density than the NM phase, and
thus is more stable.

Next, we investigate the FM instability in weakly doped ABCA graphene. Suppose
that the noninteracting ground state is initially doped with electron pockets
$Q_{0\uparrow}=Q_{0\downarrow}=Q_{0}$ for the two spin channels, as shown in
the left panels of Figs. 2(b-c). When considering e-e interaction, FM
instability occurs. There are two types of FM variational states in the doped
case, as displayed in the right panels of 2(b-c): One type with two asymmetric
electron pockets ($Q_{\uparrow}>Q_{\downarrow}>0$), and the other type with
one electron pocket and one hole pocket ($Q_{\uparrow}>0>Q_{\downarrow}$) for
the spin up and spin down channels, respectively. With conserved number of
particles, the pocket sizes are constrained by

$s_{\uparrow}\frac{Q^{2}_{\uparrow}}{4\pi}+s_{\downarrow}\frac{Q^{2}_{\downarrow}}{4\pi}=s_{0\uparrow}\frac{Q^{2}_{0\uparrow}}{4\pi}+s_{0\downarrow}\frac{Q^{2}_{0\downarrow}}{4\pi}=\frac{n}{S_{\mathrm{cell}}}$
(12)

where $s_{(0)\uparrow/\downarrow}$ equals +1 for electron pockets, and -1 for
hole pockets, and $S_{\mathrm{cell}}$ is the area of one primitive cell. From
above, we see that there is only one independent variational parameter that
accounts. Here we choose $Q_{\downarrow}$ without loss of generality. Note
that there is always an electronlike pocket for the FM variational state [see
Figs. 2(b-c)], therefore we generically set $s_{\uparrow}=1$. Also, the spin-
resolved pockets are imbalanced with $Q_{\uparrow}\neq|Q_{\downarrow}|$.

Figure 3: (Color online) (a) The variation of kinetic energy density $\Delta
E_{\mathrm{kin}}/S$, exchange energy density $\Delta E_{\mathrm{ex}}/S$ and
total energy density $\Delta E/S$ vs $\hbar\upsilon_{F}Q/\gamma$ for non-
doping system ($n=0$) with $g=2.1$. Blue dot labels the minimum of $\Delta
E/S$. (b) $\Delta E/S$ vs $x=s_{\downarrow}Q^{4}_{\downarrow}/Q^{4}_{0}$ for
low-doping systems with different $n$ and $g=2.1$. Squares label minima of
$\Delta E/S$. (c) Phase diagrams vs $n$ and $g$ for ABCA, ABC and AB graphene
systems. The three lines are boundaries separating FM and NM phases. (d)
Equilibrium FM magnetization $m$ vs $g$ for half filling ($n=0$) and doped
($n\neq 0$) ABCA graphene.

Below we evaluate the energy density $\Delta E/S$ vs
$x=s_{\downarrow}Q^{4}_{\downarrow}/Q^{4}_{0}$ to study the equilibrium
states. The calculation results for $g=2.1$ with different doping constraints
$n$ are plotted in Fig. 3(b). For slightly doping $n=5.03\times 10^{-4}$ (see
black solid line), the equilibrium state is FM phased with a hole pocket in
the spin down channel (denoted by the black square: $x<0\Rightarrow
Q_{\downarrow}<0$) and also an electron pocket in spin up channel. When doping
is enhanced to $n=5.34\times 10^{-4}$ (see cyan dash-dot line), the system is
transformed to NM phase at equilibrium (denoted by the cyan square:
$x=1\Rightarrow Q_{\downarrow}=Q_{\uparrow}=Q_{0}$). There is a transition
point where the FM phase is degenerate with the NM phase at equilibrium when
$n=5.20\times 10^{-4}$ (see pink dash line and squares). By now, we
numerically find that there exist two types of exchange instabilities, the FM
and NM instabilities, in doped ABCA graphene. To systematically study the
equilibrium states, we have calculated the phase diagram as functions of $g$
and $n$, as given in Fig. 3(c). At small $n$, the system trends to develop a
FM phase even with small $g$. For moderate $n$, the FM phases can be
stabilized by increasing $g$. Whereas it can not stay magnetic anymore when
$n$ exceeds some critical threshold, even $g$ is significantly large. The FM
phases are all characterized with one electron pocket in one spin channel and
one hole pocket in the other, rather than two electronlike pockets [see Figs.
2(b-c)]. We further note that the vanishing magnetization $m$ in the NM state
dramatically changes to finite once $g$ crosses the transition line, as shown
in Fig. 3(d). And $m$ will be suppressed as $n$ increases.

Moreover, we also studied the layer effect on the phase diagram. The results
for AB bilayer and ABC trilayer graphene are displayed in Fig. 3(c) to give
comparison with the ABCA graphene. We see that the FM phase region of the
tetralayer system is much larger than that of the bilayer and trilayer
systems. This is attributed to the peculiar quartic spectrum near the Dirac
point, which gives rise to divergent DOS [see Fig. 1(b)]. FM phase can be more
favorable and easily generated according to the Stoner criteria [54]. Thus,
the ABCA graphene is more flexible to exploit with the FM phases, superior to
bilayer and trilayer graphene systems.

In summary, we have studied the FM exchange instability of the weakly doped
ABCA graphene in the presence of long-range e-e interaction and acquired zero
temperature phase diagram. At non- or low-doping, the interaction prefers to
stabilize spin-polarized FM phases at equilibrium, where an electronlike
pocket and a holelike pocket are formed at the same time for the two spin
channels. The ABCA graphene will be more advantageous than the bilayer and
trilayer systems in the flexibility on manipulating FM phases, due to a larger
FM-phase region of its phase diagram.

##  IV Spontaneous spin superconductor state

###  IV.1 Instability of the FM state

The intrinsic e-e interaction will result in a FM state in ABCA graphene with
the coexistence of electron and hole pockets. Below we will prove that the
naturally occurring attractive interaction between electron and hole carriers,
no matter how weak it is, will lead to the instability of this FM state.
Electrons and holes will spontaneously bound into spin-triplet exciton pairs.

At the low energies of the FM state of ABCA graphene, the noninteracting
Hamiltonian for the spin $\uparrow$ (electron-like) and spin $\downarrow$
(hole-like) subsystems in valley $\bm{K}_{\eta}(\eta=\pm)$ can be
approximately written as

$H^{0}_{\uparrow(\downarrow)}=\left(\begin{array}[]{cc}-M_{e(h)}&\frac{\hbar^{4}\upsilon^{4}_{F}(\hat{k}^{e}_{x}+i\hat{k}^{e}_{y})^{4}}{\gamma^{3}}\\\
\frac{\hbar^{4}\upsilon^{4}_{F}(\hat{k}^{e}_{x}-i\hat{k}^{e}_{y})^{4}}{\gamma^{3}}&-M_{e(h)}\\\
\end{array}\right),$ (13)

where $M_{e}>0$ ($M_{h}<0$) is the chemical potential of the spin $\uparrow$
($\downarrow$) subsystem with respect to the Dirac point, and
$\bm{\hat{k}}^{e}=(\hat{k}^{e}_{x},\hat{k}^{e}_{y})$ is the momentum operator
for the electronic states. We temporarily set $\hbar=1$ hereafter in this
section. Then, the corresponding Schrödinger equations for electrons of the
two subsystems are given by

$\begin{split}\left(\begin{array}[]{cc}-M_{e}&\frac{\upsilon^{4}_{F}(\hat{k}^{e}_{x}+i\hat{k}^{e}_{y})^{4}}{\gamma^{3}}\\\
\frac{\upsilon^{4}_{F}(\hat{k}^{e}_{x}-i\hat{k}^{e}_{y})^{4}}{\gamma^{3}}&-M_{e}\\\
\end{array}\right)\varphi^{0e}_{\bm{k}n\eta\uparrow}(\bm{r}_{e})=&\varepsilon_{\bm{k}n\eta\uparrow}\varphi^{0e}_{\bm{k}n\eta\uparrow}(\bm{r}_{e});\\\
\left(\begin{array}[]{cc}-M_{h}&\frac{\upsilon^{4}_{F}(\hat{k}^{e}_{x}+i\hat{k}^{e}_{y})^{4}}{\gamma^{3}}\\\
\frac{\upsilon^{4}_{F}(\hat{k}^{e}_{x}-i\hat{k}^{e}_{y})^{4}}{\gamma^{3}}&-M_{h}\\\
\end{array}\right)\varphi^{0e}_{\bm{k}n\eta\downarrow}(\bm{r}_{e})=&\varepsilon_{\bm{k}n\eta\downarrow}\varphi^{0e}_{\bm{k}n\eta\downarrow}(\bm{r}_{e}).\end{split}$
(14)

The eigenenergies and eigenstates for the spin up subsystem are solved as

$\begin{split}&\varepsilon_{\bm{k}\bm{\pm}\eta\uparrow}=\pm\frac{\upsilon^{4}_{F}k^{4}}{\gamma^{3}}-M_{e};\\\
&\varphi^{0e}_{\bm{k}\bm{\pm}\eta\uparrow}(\bm{r}_{e})=\frac{e^{i\bm{k}\cdot\bm{r}_{e}}}{\sqrt{2}}\left(\begin{array}[]{cc}1\\\
\pm e^{-i4\theta_{\bm{k}}}\\\ \end{array}\right),\end{split}$ (15)

where $\theta_{\bm{k}}=\arctan\frac{k_{y}}{k_{x}}$. For spin down system,
$\varepsilon_{\bm{k}\bm{\pm}\eta\downarrow}=\bm{\pm}\frac{\upsilon^{4}_{F}k^{4}}{\gamma^{3}}-M_{h}$
and
$\varphi^{0e}_{\bm{k}\bm{\pm}\eta\downarrow}(\bm{r}_{e})=\varphi^{0e}_{\bm{k}\pm\eta\uparrow}(\bm{r}_{e})$.
We show the low-energy bands of the FM state in the top panel of Fig. 4(a).
The spin-up conduction band are partially filled with electron-like carriers,
while the spin-down carriers are hole-like due to the hole pocket in the
valence band. We further take the e-h transformation for the spin down
subsystem. The annihilation of an electron with spin down in valley
$\bm{K}_{\eta}$ is equivalent to the creation of a spin up hole with reversed
momentum and energy in the opposite valley
$\bm{K}_{\bar{\eta}}=\bm{K}_{-\eta}$. By transformation of Eq. (14), the
Schrödinger equation for free holes can be obtained as

$-\left(\begin{array}[]{cc}-M_{h}&\frac{\upsilon^{4}_{F}(-\hat{k}^{h}_{x}+i\hat{k}^{h}_{y})^{4}}{\gamma^{3}}\\\
\frac{\upsilon^{4}_{F}(-\hat{k}^{h}_{x}-i\hat{k}^{h}_{y})^{4}}{\gamma^{3}}&-M_{h}\\\
\end{array}\right)\varphi^{0h}_{\bm{k}n\bar{\eta}\uparrow}(\bm{r}_{h})=\varepsilon_{\bm{k}n\bar{\eta}\uparrow}\varphi^{0h}_{\bm{k}n\bar{\eta}\uparrow}(\bm{r}_{h}).$
(16)

Its solutions are given by

$\begin{split}&\varepsilon_{\bm{k}\bm{\pm}\bar{\eta}\uparrow}=\pm\frac{\upsilon^{4}_{F}k^{4}}{\gamma^{3}}+M_{h},\\\
&\varphi^{0h}_{\bm{k}\bm{\pm}\bar{\eta}\uparrow}(\bm{r}_{h})=\frac{e^{i\bm{k}\cdot\bm{r}_{h}}}{\sqrt{2}}\left(\begin{array}[]{cc}1\\\
\mp e^{i4\theta_{\bm{k}}}\\\ \end{array}\right).\end{split}$ (17)

The bottom panel of Fig. 4(a) shows the energy bands for the holes and
electrons in different valleys. In the FM state of ABCA graphene, the lowest
conduction bands of electrons and holes are filled up to the chemical
potentials $M_{\mathrm{e}}$ and $-M_{\mathrm{h}}$, respectively.

Next, we study the instability of the FM state by taking e-h interaction into
account. Following the theory of Ref. [24], we consider an additional e-h pair
residing out of the FM state composed with free fermions [as illustrated in
the bottom panel of Fig. 4(a)]. The only effect of other free electrons and
holes is to prohibit the e-h pair from occupying states below Fermi energy
$E_{f}=0$. We will subsequently show that with attractive interaction
$-U(\bm{r}_{e},\bm{r}_{h})$, this e-h pair spontaneously forms mixed state
with neutral charge and $2(\hbar/2)$ spin. The Schrödinger equation for the
interacting e-h pair reads

$\left[\left(\begin{array}[]{cc}-M_{e}&\frac{\upsilon^{4}_{F}(\hat{k}^{e}_{x}+i\hat{k}^{e}_{y})^{4}}{\gamma^{3}}\\\
\frac{\upsilon^{4}_{F}(\hat{k}^{e}_{x}-i\hat{k}^{e}_{y})^{4}}{\gamma^{3}}&-M_{e}\\\
\end{array}\right)-\left(\begin{array}[]{cc}-M_{h}&\frac{\upsilon^{4}_{F}(-\hat{k}^{h}_{x}+i\hat{k}^{h}_{y})^{4}}{\gamma^{3}}\\\
\frac{\upsilon^{4}_{F}(-\hat{k}^{h}_{x}-i\hat{k}^{h}_{y})^{4}}{\gamma^{3}}&-M_{h}\\\
\end{array}\right)-U(\bm{r}_{e},\bm{r}_{h})\right]\psi_{\uparrow}(\bm{r}_{e},\bm{r}_{h})=E\psi_{\uparrow}(\bm{r}_{e},\bm{r}_{h}),$
(18)

where $\psi_{\uparrow}(\bm{r}_{e},\bm{r}_{h})$ is a two-particle wave function
and can be expanded by

$\psi_{\uparrow}(\bm{r}_{e},\bm{r}_{h})=\sum_{\bm{k}}g_{\bm{k}}\varphi^{0e}_{\bm{k}\bm{+}\eta\uparrow}(\bm{r}_{e})\varphi^{0h}_{-\bm{k}\bm{+}\bar{\eta}\uparrow}(\bm{r}_{h}),$
(19)

with the coefficient $g_{\bm{k}}=0$ if
$\frac{\upsilon^{4}_{F}|\bm{k}|^{4}}{\gamma^{3}}-max(M_{e},-M_{h})\leq 0$ due
to the Pauli exclusion principle. Note that only free electron and hole states
with opposite momenta are paired since they are energetically favorable, and
the mixed $\psi_{\uparrow}(\bm{r}_{e},\bm{r}_{h})$ has a nonzero spin
$2(\hbar/2)$. Put $\psi_{\uparrow}(\bm{r}_{e},\bm{r}_{h})$ into Eq. (18), we
obtain

$\begin{split}(&2\frac{\upsilon^{4}_{F}|\bm{k}|^{4}}{\gamma^{3}}-M_{e}+M_{h})g_{\bm{k}}-\sum_{\bm{k}^{\prime}}g_{\bm{k}^{\prime}}U_{\bm{k}\bm{k}^{\prime}}=Eg_{\bm{k}}\\\
&\Rightarrow~{}~{}g_{\bm{k}}=\sum_{\bm{k}^{\prime}}\frac{g_{\bm{k}^{\prime}U_{\bm{k}\bm{k}^{\prime}}}}{2\frac{\upsilon^{4}_{F}|\bm{k}|^{4}}{\gamma^{3}}-M_{e}+M_{h}-E}\end{split}$
(20)

where $U_{\bm{k}\bm{k}^{\prime}}=\int\int
d\bm{r}_{e}d\bm{r}_{h}~{}U(\bm{r}_{e},\bm{r}_{h})e^{-i(\bm{k}-\bm{k}^{\prime})(\bm{r}_{e}-\bm{r}_{h})}$.
Because that $U_{\bm{k}\bm{k}^{\prime}}$ decreases slowly with
$|\bm{k}-\bm{k}^{\prime}|$, it is reasonable to assume it as a constant $U$,
independent of $\bm{k}$ and $\bm{k}^{\prime}$. Then the above equation reduces
to

$1=-U\int_{max(M_{e},-M_{h})}^{\varepsilon_{D}}d\varepsilon~{}~{}\frac{\rho(\varepsilon)}{E+M_{e}-M_{h}-2\varepsilon},$
(21)

where $\rho(\varepsilon)\sim 1/\sqrt{\varepsilon}$ is the low-energy DOS of
the bare ABCA graphene without interaction. With positive $U$ and
$M_{e}~{}(-M_{h})>0$, a negative $E$ solution of this equation is always
possible. Therefore, we can prove that the mixed e-h pair is more stable since
it is lower in energy than a free e-h pair (energy
$E=2\frac{\upsilon^{4}_{F}|\bm{k}|^{4}}{\gamma^{3}}-M_{e}+M_{h}>0$). Electron
and hole will spontaneously bind into exciton pair with neutral charge and
spin polarization $2(\hbar/2)$. In conclusion, the FM state in ABCA graphene
will be unstable in the presence of attractive e-h interaction.

Figure 4: (Color online) (a) Top panel: Schematic of low-energy bands of FM
phased ABCA gaphene for spin-up (left) and spin-down (right) electrons in
valley $\bm{K}_{\eta}$. Bottom panel: Low-energy bands of FM ABCA graphene for
spin-up electrons in valley $\bm{K}_{\eta}$ (left) and spin-up holes in
opposite valley $\bm{K}_{\bar{\eta}}$ (right). An additional e-h pair with
attractive interaction $-U$ is labeled. The horizontal dashed lines mark the
Fermi levels $E_{f}=0$. (b) Schematic of low-energy bands of FM phased ABCA
graphene for free fermions with pair potential $\Delta=0$ (dashed lines), and
for excitonic condensate with $\Delta\neq 0$ (solid lines). (c) $\Delta$ vs
$T$ with different $g$ and $n$. (d) $T_{c}$ vs $n$ with different $g$.

###  IV.2 BCS-type theory for the spin superconductor

In FM phased ABCA graphene, we have shown that electrons and holes will
spontaneously bind into more stable spin-polarized exciton pairs. At low
temperature, they could condense into a spin-polarized superconducting state,
namely, a spin SC. Here, we develop a BCS-type theory for the spin SC state in
ABCA graphene.

Let us begin by considering an interacting electron system in FM phased ABCA
graphene at low energies. Its Hamiltonian reads $H=H_{0}+H_{\mathrm{int}}$,
where $H_{0}$ represents the noninteracting Hamiltonian for electrons from the
lowest conduction band (labeled $\bm{+}$) of spin up channel and the highest
valence band (labeled $\bm{-}$) of the spin down channel [see Fig. 4(a)], and
$H_{\mathrm{int}}$ depicts the e-e interaction between the two band branches:

$\begin{split}&H_{0}=\sum_{\eta,\bm{k}}\varepsilon_{\bm{k}\bm{+}\eta\uparrow}\phi^{\dagger}_{\bm{k}\bm{+}\eta\uparrow}\phi_{\bm{k}\bm{+}\eta\uparrow}+\varepsilon_{\bm{k}\bm{-}\eta\downarrow}\phi^{\dagger}_{\bm{k}\bm{-}\eta\downarrow}\phi_{\bm{k}\bm{-}\eta\downarrow},\\\
&H_{\mathrm{int}}=\sum_{\eta,\bm{k},\bm{k}^{\prime},\bm{q}}U^{\eta}_{\bm{k}\bm{k^{\prime}}\bm{q}}\phi^{\dagger}_{\bm{k}-\bm{q}\bm{+}\eta\uparrow}\phi^{\dagger}_{\bm{k}^{\prime}+\bm{q}\bm{-}\eta\downarrow}\phi_{\bm{k}^{\prime}\bm{-}\eta\downarrow}\phi_{\bm{k}\bm{+}\eta\uparrow}.\end{split}$
(22)

Here, the interaction matrix element
$U^{\eta}_{\bm{k}\bm{k^{\prime}}\bm{q}}=\langle\bm{k}-\bm{q}\bm{+},\eta\uparrow|\langle\bm{k}^{\prime}+\bm{q}\bm{-},\eta\downarrow|V(\bm{r}-\bm{r}^{\prime})|\bm{k}^{\prime}\bm{-},\eta\downarrow\rangle|\bm{k}\bm{+},\eta\uparrow\rangle$.
The two band branches
$\varepsilon_{\bm{k}\bm{+}\eta\uparrow}=\varepsilon_{\bm{k}}-M_{e}$ and
$\varepsilon_{\bm{k}\bm{-}\eta\downarrow}=-\varepsilon_{\bm{k}}-M_{h}$ can be
numerically calculated by the diagonalization of Eq. (3), with their chemical
potentials $M_{e}$ and $M_{h}$ determined by the variational method from last
section. Because the spin-down carrier is being hole-like in FM phased system,
the annihilation of an electron by $\phi_{\bm{k}-\eta\downarrow}$ also means
the creation of a hole with opposite momentum and spin. Therefore we can
define the operator for a hole state by
$\alpha^{\dagger}_{-\bm{k}h\bar{\eta}\uparrow}=\phi_{\bm{k}\bm{-}\eta\downarrow}$.
For an electronic state in spin-up band, we define
$\alpha_{\bm{k}e\eta\uparrow}=\phi_{\bm{k}\bm{+}\eta\uparrow}$. Then the
noninteracting Hamiltonian $H_{0}$ can be transformed to

$H_{0}=\sum_{\eta,\bm{k}}\left(\begin{array}[]{cc}\alpha^{\dagger}_{\bm{k}e\eta\uparrow}&\alpha_{-\bm{k}h\bar{\eta}\uparrow}\\\
\end{array}\right)\left(\begin{array}[]{cc}\varepsilon_{\bm{k}e\eta\uparrow}&0\\\
0&-\varepsilon_{-\bm{k}h\bar{\eta}\uparrow}\\\
\end{array}\right)\left(\begin{array}[]{cc}\alpha_{\bm{k}e\eta\uparrow}\\\
\alpha^{\dagger}_{-\bm{k}h\bar{\eta}\uparrow}\\\ \end{array}\right),$ (23)

where we set
$\varepsilon_{\bm{k}e\eta\uparrow}=\varepsilon_{\bm{k}\bm{+}\eta\uparrow}$ and
$-\varepsilon_{-\bm{k}h\bar{\eta}\uparrow}=\varepsilon_{\bm{k}\bm{-}\eta\downarrow}$.
Notably, the nonequal chemical potentials $M_{e}$ and $-M_{h}$ of the electron
band $\varepsilon_{\bm{k}e\eta\uparrow}=\varepsilon_{\bm{k}}-M_{e}$ and hole
band $\varepsilon_{\bm{k}h\bar{\eta}\uparrow}=\varepsilon_{\bm{k}}+M_{h}$
result in a splitted Fermi surface of the electron and hole states, which can
be effectively viewed as a Zeeman split induced by the “magnetic field”.

The interaction part $H_{\mathrm{int}}$ includes terms given by
$\alpha^{\dagger}_{\bm{k}-\bm{q}e\eta\uparrow}\alpha_{-(\bm{k^{\prime}}+\bm{q})h\bar{\eta}\uparrow}\alpha^{\dagger}_{-\bm{k}^{\prime}h\bar{\eta}\uparrow}\alpha_{\bm{k}e\eta\uparrow}$.
In the presence of the splitted Fermi surface, an exotic Fulde-Ferrell-Larkin-
Ovchinnikov superconducting state with nonzero momentum paring was predicted
to appear[57, 58, 59, 60, 61]. Whereas, its formation requires large split of
the spin-resolved chemical potentials, or the BCS state with zero momentum
paring will continue to dominate[58, 59, 60]. In this work, the imbalance of
$M_{e}$ and $-M_{h}$ is small with the chosen parameters. Thus, we can
reasonably reserve terms satisfying $\bm{k}^{\prime}+\bm{q}=\bm{k}$ in
$H_{\mathrm{int}}$, which pair electron and hole with zero momentum
($-\bm{k}h\bar{\eta}\uparrow,\bm{k}e\eta\uparrow$). Under BCS paring,
$H_{\mathrm{int}}$ transforms to the attractive e-h interaction

$H_{\mathrm{int}}=-\sum_{\eta,\bm{k},\bm{k}^{\prime}}U^{\eta}_{\bm{k}\bm{k}^{\prime}}\alpha^{\dagger}_{\bm{k}^{\prime}e\eta\uparrow}\alpha^{\dagger}_{-\bm{k}^{\prime}h\bar{\eta}\uparrow}\alpha_{-\bm{k}h\bar{\eta}\uparrow}\alpha_{\bm{k}e\eta\uparrow},$
(24)

where the interaction matrix element

$U^{\eta}_{\bm{k}\bm{k}^{\prime}}=\frac{1}{N}\sum_{\alpha,\beta}\left[M^{\alpha,\ast}_{\bm{k}^{\prime}\bm{+}\eta\uparrow}M^{\beta,\ast}_{\bm{k}\bm{-}\eta\downarrow}M^{\beta}_{\bm{k}^{\prime}\bm{-}\eta\downarrow}M^{\alpha}_{\bm{k}\bm{+}\eta\uparrow}\right]V^{\alpha\beta}_{\bm{k}\bm{k}^{\prime}},$
(25)

$V^{\alpha\beta}_{\bm{k}\bm{k}^{\prime}}=\sum_{j}V^{\alpha\beta}_{0j}e^{-i(\bm{k}-\bm{k}^{\prime})(\bm{R}_{j}+\bm{\delta}_{\beta}-\bm{\delta}_{\alpha})}$
with the nonlocal real-space Coulomb interaction matrix element
$V^{\alpha\beta}_{0j}=\langle\alpha 0|\langle\beta
j|V(\bm{r}-\bm{r}^{\prime})|\beta j\rangle|\alpha 0\rangle=\langle\alpha
0|\langle\beta j|\frac{g\hbar\nu_{F}}{|\bm{r}-\bm{r}^{\prime}|}|\beta
j\rangle|\alpha 0\rangle$, and $M^{\alpha(\beta)}_{\bm{k}\bm{\pm}\eta\sigma}$
is the bare system eigenstate of spin $\sigma$($\uparrow,\downarrow$) band in
valley $\bm{K}_{\eta}$ with momentum $\bm{k}$, projected on sublattice
$\alpha$ $(\beta)$. The mean-field approximation of $H_{\mathrm{int}}$ is

$H_{\mathrm{int}}\approx-\sum_{\eta,\bm{k}}\Delta^{\eta}_{\bm{k}}\alpha^{\dagger}_{\bm{k}e\eta\uparrow}\alpha^{\dagger}_{-\bm{k}h\bar{\eta}\uparrow}+\Delta^{\eta,\ast}_{\bm{k}}\alpha_{-\bm{k}h\bar{\eta}\uparrow}\alpha_{\bm{k}e\eta\uparrow},$
(26)

where the exciton pair potential is defined as
$\Delta^{\eta}_{\bm{k}}=\sum_{\bm{k}^{\prime}}U^{\eta}_{\bm{k}^{\prime}\bm{k}}\langle\alpha_{-\bm{k}^{\prime}h\bar{\eta}\uparrow}\alpha_{\bm{k}^{\prime}e\eta\uparrow}\rangle$.

Combing Eqs. (23) and (26), we have the total mean-field Hamiltonian $H$ in
the BdG representation

$H=\sum_{\eta,\bm{k}}\left(\begin{array}[]{cc}\alpha^{\dagger}_{\bm{k}e\eta\uparrow}&\alpha_{-\bm{k}h\bar{\eta}\uparrow}\\\
\end{array}\right)\left(\begin{array}[]{cc}\varepsilon_{\bm{k}e\eta\uparrow}&-\Delta^{\eta}_{\bm{k}}\\\
-\Delta^{\eta,\ast}_{\bm{k}}&-\varepsilon_{-\bm{k}h\bar{\eta}\uparrow}\\\
\end{array}\right)\left(\begin{array}[]{cc}\alpha_{\bm{k}e\eta\uparrow}\\\
\alpha^{\dagger}_{-\bm{k}h\bar{\eta}\uparrow}\\\ \end{array}\right).$ (27)

Since the two valleys $\bm{K}_{\pm}$ are degenerate, we only need to deal with
one single valley and omit the index $\eta$ afterwards. When $M_{e}=-M_{h}=M$,
the energy spectrum of $H_{\rm BdG}$ are given by
$E_{\bm{k}\bm{\pm}}=\pm\sqrt{(\varepsilon_{\bm{k}}-M)^{2}+|\Delta_{\bm{k}}|^{2}}$,
which are schematically plotted in Fig. 4(c). Compared to the bands in the
absence of pair potential [see dashed lines in Fig. 4(c)], an energy gap
$\Delta$ is opened. This means that when an electron and a hole attract and
bind into e-h pair, the energy of the system is lowered. Therefore, the ground
state of the FM phased ABCA graphene is a superfluid composed of spin-triplet
exciton pairs($-\bm{k}h\uparrow,\bm{k}e\uparrow$), namely, a spin SC.
Analogous to the normal SC, a spin supercurrent can flow without dissipation
in the spin SC. On the other hand, because of the neutral charge of the
exciton pairs and the opened energy gap, no charge current can flow making the
spin SC a charge insulator[24].

In spin SCs, the spin superconducting gap $\Delta(T)$ can be estimated by
solving the self-consistent equation[60]:

$\Delta_{\bm{k}}=\sum_{\bm{k}^{\prime}}\frac{U_{\bm{k}^{\prime}\bm{k}}\Delta_{\bm{k}^{\prime}}}{2E_{\bm{k}^{\prime}}}[f(E_{\bm{k}^{\prime}\bm{-}})-f(E_{\bm{k}^{\prime}\bm{+}})],$
(28)

where $E_{\bm{k}\bm{\pm}}=-(M_{e}+M_{h})/2\pm E_{\bm{k}}$ with
$E_{\bm{k}}=\sqrt{[\varepsilon_{\bm{k}}-(M_{e}-M_{h})/2]^{2}+|\Delta_{\bm{k}}|^{2}}$,
$f(E)=1/(e^{E/K_{B}T}+1)$ is the Fermi-Dirac distribution function with $T$
the temperature. Due to that $U_{\bm{k}^{\prime}\bm{k}}$ varies slowly with
$\bm{k}^{\prime}-\bm{k}$, let us assume $U_{\bm{k}^{\prime}\bm{k}}\approx
U_{\bm{k}\bm{k}}\theta(k_{D}-|\bm{k}^{\prime}-\bm{k}|)=\frac{U_{0}}{N}\theta(k_{D}-|\bm{k}^{\prime}-\bm{k}|)$,
where $k_{D}$ is the momentum cut-off and

$U_{\bm{k}\bm{k}}=\frac{U_{0}}{N}=\sum_{\alpha,\beta}\left[M^{\alpha,\ast}_{\bm{k}\bm{+}\eta\uparrow}M^{\beta,\ast}_{\bm{k}\bm{-}\eta\downarrow}M^{\beta}_{\bm{k}\bm{-}\eta\downarrow}M^{\alpha}_{\bm{k}\bm{+}\eta\uparrow}\right]\sum_{j}\frac{V^{\alpha\beta}_{0j}}{N}.$
(29)

Note that $V^{\alpha\beta}_{0j}$ depends linearly on the parameter $g$, so
does $U_{0}$. Invoking the values of real-space nonlocal interaction
$V^{\alpha\beta}_{0j}$ given by Ref. [62], we thus can estimate $U_{0}\approx
28.0$ eV in ABCA graphene (referring to $g=2.1$). Under the above assumption,
Eq. (28) is reduced to

$1=U_{0}\int_{0}^{\varepsilon_{D}}d\varepsilon~{}\frac{\rho(\varepsilon)[f(E_{\varepsilon\bm{-}})-f(E_{\varepsilon\bm{+}})]}{\sqrt{(2\varepsilon-
M_{e}+M_{h})^{2}+4\Delta(T)^{2}}}$ (30)

where $\rho(\varepsilon)\sim 1/\sqrt{\varepsilon}$ is the DOS of bare ABCA
graphene for one valley and one spin channel, and the energy cut-off
$\varepsilon_{D}$ is a function of $k_{D}$. Due to that $M_{e}$ and $M_{h}$
are equilibrium properties that are determined by doping $n$ and interaction
$g$, $\Delta(T)$ should also depend on them according to Eq. (30). We
numerically calculate $\Delta$ as a function of $T$ for several $n$ and $g$,
with $\varepsilon_{D}=\gamma$. As shown in Fig. 4(c), $\Delta$ decreases with
$T$ as expected. At $T=0$ K, it is less affected by $n$ with a gap
$\Delta\approx 7.0$ meV for $g=2.1$ (see solid lines), which is about 2 times
larger than that in FM graphene ($\approx 3$ meV [24]) and comparable to that
in trilayer graphene ($\approx 7.8$ meV [25]). At $T\neq 0$ K, $\Delta$ and
thus the spin superconductivity is suppressed by enlarging doping $n$, as the
Fermi surface split and equivalent “magnetic filed” grow larger driving the
system out of the spin SC state. In screened system with $g=1.9$, the spin
split-energy is decreased, thus giving rise to smaller gaps $\Delta$ compared
to that given by $g=2.1$ (see dash lines).

The critical temperature $T_{c}$ can be estimated by bring $\Delta(T_{c})=0$
into Eq. (30). We numerically calculate $T_{c}$ as a function of $n$ for
different $g$, with $\varepsilon_{D}=\gamma$. As shown in Fig. 4(d), $T_{c}$
is lowered by increasing $n$, which is attributed to the effect of the
equivalent “magnetic filed” as explained above. It can as well be suppressed
by decreasing $g$ from 2.1 to 1.9, as illustrated by the solid and dash lines.
For non-doping system with $n=0$ and $g=2.1$, $T_{c}$ reaches about $45$ K.
The spin superconducting gap $\Delta(0)\approx 7.0$ meV and critical
temperature $T_{c}\approx 45$ K is moderate, which should make the spin SC
state in ABCA graphene readily to confirm under the state-of-the-art
experimental conditions.

One more note to stress, the spin SC state in ABCA graphene can not exist if
the system is heavily doped with large $n$, where FM phase disappears [see
Fig. 3(c)].

##  V Spin-current Josephson effect

Josephson effect is a fundamental hallmark of the charge SCs and has widely
applications. As a counterpart of the SC, the spin SC holds a similar spin-
current Josephson effect, which has been demonstrated in the FM graphene [24].
Here, we will study it in a junction composed of FM ABCA graphene spin SC
/quantum dot/FM ABCA graphene spin SC, as shown in Fig. 5(a). The junction
Hamiltonian comprises three part:
$H=\sum_{\beta=L,R}H_{\beta}+H_{\mathrm{dot}}+H_{\mathrm{coup}}$, where the
lead Hamiltonian $H_{\beta}$, the quantum dot $H_{\mathrm{dot}}$ and their
coupling $H_{\mathrm{coup}}$ are, respectively,

$\begin{split}&H_{\beta}=\sum_{\eta,\bm{k}}\left(\begin{array}[]{cc}\alpha^{\dagger}_{\beta\bm{k}\bm{+}\eta\uparrow}&\alpha^{\dagger}_{\beta\bm{k}\bm{-}\eta\downarrow}\\\
\end{array}\right)\left(\begin{array}[]{cc}\varepsilon_{\bm{k}}-M&\Delta
e^{i\phi_{\beta}}\\\ \Delta e^{-i\phi_{\beta}}&-\varepsilon_{\bm{k}}+M\\\
\end{array}\right)\left(\begin{array}[]{cc}\alpha_{\beta\bm{k}\bm{+}\eta\uparrow}\\\
\alpha_{\beta\bm{k}\bm{-}\eta\downarrow}\\\ \end{array}\right),\\\
&H_{\mathrm{dot}}=\sum_{\sigma=\uparrow,\downarrow}(\epsilon_{d}+\sigma
M_{d})d^{\dagger}_{\sigma}d_{\sigma},\\\
&H_{\mathrm{coup}}=\sum_{\beta,\eta,\bm{k}}(t_{\beta}\alpha^{\dagger}_{\beta\bm{k}\bm{+}\eta\uparrow}d_{\uparrow}+t_{\beta}\alpha^{\dagger}_{\beta\bm{k}\bm{-}\eta\downarrow}d_{\downarrow}+\mathrm{H.c.}).\end{split}$
(31)

Here, $H_{\beta}$ is the same as $H_{\mathrm{BdG}}$ in Eq. (27). $\beta=L/R$
denotes the clean spin SC lead ($n=0$) with $M_{e}=-M_{h}=M$, and its order
parameter $\Delta e^{i\phi_{\beta}}$ ($\Delta$ and $\phi_{\beta}$ are SC gap
and phase, respectively) is assumed to be independent of momentum $\bm{k}$.
$\epsilon_{d}$ is the quantum dot level with a spin-split energy $M_{d}$.
$t_{\beta}$ is the coupling of the lead with the quantum dot.

Using the non-equilibrium Green’s function method, the spin-resolved current
flowing from lead $\beta$ to the central region at equilibrium can be
calculated by
$I_{\beta\uparrow(\downarrow)}=\mathrm{Re}\frac{2et_{\beta}}{\hbar}\sum_{\eta,\bm{k}}G^{<}_{d\beta\bm{k}\eta,11(22)}(t,t)$,
where $G^{<}_{d\beta\bm{k}\eta,11(22)}(t,t)$ is the up (down) diagonal element
of the lesser Green’s function matrix defined in the Nambu basis. Applying
Dyson’s equation to expand $G^{<}_{d\beta\bm{k}\eta}(t,t)$, the spin-resolved
currents flowing in the left lead can then be derived as (details can be found
in Appendix C)

$\begin{split}&I_{L\uparrow}\approx\frac{4e}{\hbar}\int\frac{d\varepsilon}{2\pi}~{}\frac{\Gamma^{2}\Delta^{2}}{\varepsilon^{2}-\Delta^{2}}f(\varepsilon)\mathrm{Im}\frac{1}{B^{\ast}(\varepsilon)}\mathrm{sin}(\phi_{R}-\phi_{L});\\\
&I_{L\downarrow}\approx\frac{4e}{\hbar}\int\frac{d\varepsilon}{2\pi}~{}\frac{\Gamma^{2}\Delta^{2}}{\varepsilon^{2}-\Delta^{2}}f(\varepsilon)\mathrm{Im}\frac{1}{B^{\ast}(\varepsilon)}\mathrm{sin}(\phi_{L}-\phi_{R}),\end{split}$
(32)

where $\Gamma=2\pi|t_{\beta}|^{2}\rho_{N}(0)$ is the linewidth function of
$\beta$ lead with $\rho_{N}(0)$ the DOS of the lead in norm state ($\Delta=0$)
at Fermi energy $E_{f}=0$. The expression of $B(\varepsilon)$ is complex and
thus shown in the appendix.

Figure 5: (Color online) (a) Schematic of the spin SC/quantum dot/spin SC
heterojunction. L/R labels the left/right FM-phased ABCA graphene leads and C
is the central quantum dot. Blue arrows denote magnetization direction. (b)
Spin supercurrent $I^{s}_{L}$ versus phase difference $\Delta\phi/2\pi$ for
different $M_{d}$. Here we set $M=0.036\gamma$, $\Delta=7$ meV ($g=2.1$ and
$n=0$), $\Gamma=0.05M$ and $\epsilon_{d}=0$. (c) $I^{s}_{L}$ versus
$\Delta\phi/2\pi$ for different $\Delta$. Here we choose $n=0$ and
$g=1.9,2.1,2.3,2.5$, which gives $M=(0.0356,0.0364,0.0375,0.0381)\gamma$ and
$\Delta=(5.05,6.81,8.80,10.96)$ meV. Other parameters are $\epsilon_{d}=0$,
$M_{d}=0.1M$ and $\Gamma=0.05M$.

Above equations reveal that the tunneling current with spin up is opposite to
that with spin down, rendering a vanishing charge current
$I^{e}_{L}=(I_{L\uparrow}+I_{L\downarrow})=0$. This is in agreement with the
fact that the charge neutral e-h pairs in the spin SC can not transport charge
current. Whereas, there is a net spin supercurrent flowing through the
junction that resembles the Josephson effect in the conventional SCs:

$\begin{split}I^{s}_{L}=&-(I_{L\uparrow}-I_{L\downarrow})\frac{\hbar}{2e}\\\
\approx&\int\frac{d\varepsilon}{2\pi}~{}\frac{4\Gamma^{2}\Delta^{2}}{\varepsilon^{2}-\Delta^{2}}f(\varepsilon)\mathrm{Im}\frac{1}{B^{\ast}(\varepsilon)}\mathrm{sin}(\phi_{L}-\phi_{R}).\end{split}$
(33)

We numerically calculate $I^{s}_{L}$ as functions of phase difference
$\Delta\phi=\phi_{L}-\phi_{R}$ for several $M_{d}$ and $\Delta$, which are
shown in Figs. 5(b-c). The spin supercurrent $I^{s}_{L}$ can exist with non-
zero $\Delta\phi$ and $\Delta$, even in a spin degenerate quantum dot where
$M_{d}=0$. In summary, we have shown that the junction composed of ABCA
graphene spin SCs and quantum dot harbors a net spin supercurrent at
equilibrium, manifesting the spin-current Josephson effect.

##  VI Conclusions

In conclusion, we demonstrate that the ground state of ABCA-stacked tetralayer
graphene can be a spontaneous spin SC state at low temperature. This finding
is quantitatively elucidated by first studying the FM exchange instability and
the equilibrium phase of the system in the presence of long-range e-e
interaction. The phase diagram indicates that at non- or low-doping levels,
the interaction will stabilize a FM phase with the coexistence of electron and
hole pockets in the two spin channels, respectively. However, this FM state
becomes unstable when considering the naturally occurring attractive e-h
interaction. As a result, electrons and holes spontaneously form spin-triplet
excitonic pairs, which further condense into a spin superfluid, termed spin
SC. We developed a consistent BCS-type theory for this spin SC state in ABCA
graphene, estimating the superconducting gap and critical temperature to be
approximately 7.0 meV and 45 K for the non-doping system. Tetralayer graphene,
superior to bilayer and trilayer forms, offers more flexibility in
manipulating the FM phase, potentially enhancing the realization of the spin
SC state. Additionally, a spin-current Josephson effect is demonstrated in the
ABCA graphene spin SC/quantum dot/ABCA graphene spin SC junction.

###### Acknowledgements.

We thank Y.-H. Li, M. Gong and H.-L. Li for fruitful discussions. This work is
financially supported by the National Key R&D Program of China (Grants No.
2019YFA0308403), and the National Natural Science Foundation of China (Grant
No. 12350401). S. Li and Y.-H. Ren contributed equally to this work.

##  Appendix A Derivation of the e-e exchange interaction

Here we provide the derivation of Eq. (8) in the main text starting from the
second quantization formula Eq. (7). First, we expand the field operators
$\varphi(\bm{r})$ in the orbital Bloch wave function basis described by
$\psi_{\alpha\eta\sigma}(\bm{k})$

$\varphi^{\dagger}(\bm{r})=\sum_{\alpha\bm{k}\sigma}\psi^{\dagger}_{\alpha\sigma}(\bm{k})~{}\langle\alpha\bm{k}\sigma|;~{}~{}~{}\varphi(\bm{r})=\sum_{\alpha\bm{k}\sigma}|\alpha\bm{k}\sigma\rangle~{}\psi_{\alpha\sigma}(\bm{k}),$
(S1)

where $\alpha(\mathrm{A}_{1},\cdots,\mathrm{B}_{4})$ is the sublattice index.
The valley index $\eta$ has been omitted for that momentum $\bm{k}$ is summed
over the whole BZ. Bringing above into Eq. (7), we have

$H_{\rm
ex}=\frac{1}{2}\sum_{\alpha_{1}\bm{k}_{1}\sigma_{1}}\cdots\sum_{\alpha_{4}\bm{k}_{4}\sigma_{4}}\langle\alpha_{1}\bm{k}_{1}\sigma_{1},\alpha_{2}\bm{k}_{2}\sigma_{2}|V(\bm{r}-\bm{r}^{\prime})|\alpha_{3}\bm{k}_{3}\sigma_{3},\alpha_{4}\bm{k}_{4}\sigma_{4}\rangle\psi^{\dagger}_{\alpha_{1}\sigma_{1}}(\bm{k}_{1})\psi^{\dagger}_{\alpha_{2}\sigma_{2}}(\bm{k}_{2})\psi_{\alpha_{3}\sigma_{3}}(\bm{k}_{3})\psi_{\alpha_{4}\sigma_{4}}(\bm{k}_{4}),$
(S2)

where the Coulomb interaction
$V(\bm{r}-\bm{r}^{\prime})=\frac{e^{2}}{\epsilon|\bm{r}-\bm{r}^{\prime}|}$.
The interaction matrix elements can be simplified a step further. Suppose that
the sublattices $\alpha_{1}$ and $\alpha_{4}$ belong to layer $j$, and that
$\alpha_{2}$ and $\alpha_{3}$ locate on layer $i$, then

$\begin{split}\langle\alpha_{1}\bm{k}_{1}\sigma_{1},\alpha_{2}\bm{k}_{2}\sigma_{2}|V(\bm{r}_{i}-\bm{r}_{j})|\alpha_{3}\bm{k}_{3}\sigma_{3},\alpha_{4}\bm{k}_{4}\sigma_{4}\rangle&=\frac{\delta_{\sigma_{1},\sigma_{4}}\delta_{\sigma_{2},\sigma_{3}}}{S}\sum_{\bm{q}\neq\bm{0}}V_{ij}(\bm{q})\langle\alpha_{1}\bm{k}_{1},\alpha_{2}\bm{k}_{2}|e^{i\bm{q}\cdot(\bm{r}_{i}-\bm{r}_{j})}|\alpha_{3}\bm{k}_{3},\alpha_{4}\bm{k}_{4}\rangle\\\
&=\frac{\delta_{\sigma_{1},\sigma_{4}}\delta_{\sigma_{2},\sigma_{3}}\delta_{\alpha_{1},\alpha_{4}}\delta_{\alpha_{2},\alpha_{3}}}{S}\sum_{\bm{q}\neq\bm{0}}\sum_{\bm{G,G^{\prime}}}V_{ij}(\bm{q})\delta_{\bm{k}_{3}+\bm{q},\bm{k}_{2}+\bm{G}}\delta_{\bm{k}_{4}-\bm{q},\bm{k}_{1}+\bm{G^{\prime}}},\end{split}$
(S3)

where $S$ denotes the area of the system, $V_{ij}(\bm{q})=\frac{2\pi
e^{2}}{\epsilon q}e^{-qd_{ij}}$ is the Fourier component of
$V(\bm{r}_{i}-\bm{r}_{j})$ with $d_{ij}$ the vertical distance between layers
$i$ and $j$, $\bm{G}$ and $\bm{G}^{\prime}$ are reciprocal lattice vectors and
$\delta$’s are the Kronecker symbols. Here, the $\bm{q}=\bm{0}$ term is
excluded due to the compensation of interaction from the positive ionic
background. To obtain the second line of above equation, we have invoked the
relation that

$\langle\alpha\bm{k}_{1}|e^{i\bm{q}\cdot\bm{r}}|\beta\bm{k}_{2}\rangle=\frac{1}{N}\sum_{i,j}e^{i\bm{k}_{2}\cdot\bm{R}^{\beta}_{j}}e^{-i\bm{k}_{1}\cdot\bm{R}^{\alpha}_{i}}\langle\alpha
i|e^{i\bm{q}\cdot\bm{r}}|\beta
j\rangle=\frac{\delta_{\alpha,\beta}}{N}\sum_{i}e^{i(\bm{k}_{2}+\bm{q}-\bm{k}_{1})\cdot\bm{R}^{\alpha}_{i}}=\sum_{\bm{G}}\delta_{\alpha,\beta}\delta_{\bm{k}_{2}+\bm{q},\bm{k}_{1}+\bm{G}},$
(S4)

where $R^{\alpha}_{i}$ ($R^{\beta}_{j}$) represents the lattice vector of
sublattice $\alpha$ ($\beta$) in the $i$th ($j$th) unit cell. Therefore in
$H_{\mathrm{ex}}$, the contribution from layers $i$ and $j$ can be recast to

$\begin{split}H^{ji}_{\mathrm{ex}}=\frac{1}{2S}\sum_{\alpha_{1}\epsilon
j}\sum_{\alpha_{2}\epsilon
i}\sum_{\bm{k}_{1},\bm{k}_{2}}\sum_{\bm{q}\neq\bm{0}}\sum_{\sigma}V_{ij}(\bm{q})~{}[\psi^{\dagger}_{\alpha_{1}\sigma}(\bm{k}_{1})\psi^{\dagger}_{\alpha_{2}\sigma}(\bm{k}_{2})&\psi_{\alpha_{2}\sigma}(\bm{k}_{2}-\bm{q})\psi_{\alpha_{1}\sigma}(\bm{k}_{1}+\bm{q})+\\\
&\psi^{\dagger}_{\alpha_{1}\sigma}(\bm{k}_{1})\psi^{\dagger}_{\alpha_{2}\bar{\sigma}}(\bm{k}_{2})\psi_{\alpha_{2}\bar{\sigma}}(\bm{k}_{2}-\bm{q})\psi_{\alpha_{1}\sigma}(\bm{k}_{1}+\bm{q})].\end{split}$
(S5)

We can define the Fourier component of the layer-dependent charge density
operator $\hat{\rho}_{i}(\bm{r}_{i})$ as usual:
$\hat{\rho}_{i}(\bm{q})=\sum_{\alpha\bm{p}\sigma}\psi^{\dagger}_{\alpha_{i}\sigma}(\bm{p}+\bm{q})\psi_{\alpha_{i}\sigma}(\bm{p})$.
By summing all the contributions from different combinations of $i$ and $j$,
the exchange interaction can be written in the form that

$H_{\mathrm{ex}}=\sum_{i,j}H^{ji}_{\mathrm{ex}}=\frac{1}{2S}\sum_{\bm{q}\neq\bm{0}}\sum_{i,j}\hat{\rho}_{i}(\bm{q})V_{ij}(\bm{q})\hat{\rho}_{j}(-\bm{q}),$
(S6)

which is the expression shown in Eq. (8) of the main text.

##  Appendix B Equilibrium properties of non-doping ABCA graphene

Here, we give a complementary discussion of the equilibrium properties of non-
doping ABCA graphene in the presence of e-e interaction. First, we investigate
the effects of $g$. As shown in Fig. S1(a), we numerically calculate the
kinetic $\Delta E_{\mathrm{kin}}/S$ and exchange $\Delta E_{\mathrm{ex}}/S$
energy densities of the FM variational states versus $Q$ for $g=1.0,2.1$ and
$3.0$. $\Delta E_{\mathrm{kin}}/S$ is invariant with $g$ for sure because it
is not relevant to interactions. Whereas, $\Delta E_{\mathrm{ex}}/S$ is
firstly lowered with the

Figure S1: (Color online) (a) $\Delta E_{\mathrm{kin}}/S$ and $\Delta
E_{\mathrm{ex}}/S$ of the FM variational state as functions of variational
parameter $Q$ for various $g$. (b) The equilibrium pocket $Q_{\mathrm{min}}$
vs $g$ for half filling ABCA graphene ($n=0$).

increase of $g$, and then goes up by increasing $g$ when $Q$ exceeds some
critical value. The position of the minimum of $\Delta E_{\mathrm{ex}}/S$ is
affected by $g$, thus rendering the equilibrium pocket $Q_{\mathrm{min}}$
dependent on $g$ as well. We show the $Q_{\mathrm{min}}$-$g$ relation in Fig.
S1(b). $Q_{\mathrm{min}}$ first goes up rapidly and then saturates for large
enough $g$. The above study reveals that the FM phase is more stabilized in
the half filling system with increasing $g$.

Figure S2: (Color online) (a) The spin-flipped variational state with respect
to that given in Fig. 2(a). (b) $\Delta E/S$ vs $Q$ for the spin-flipped state
with $\downarrow$ magnetization (dash lines) and for the FM state with
$\uparrow$ magnetization (solid lines). (c) Nonmagnetic variational state with
reversed pockets between the two valleys $\bm{K}_{a}$ and $\bm{K}_{\bar{a}}$,
and between the two spin channels. (d) $\Delta E/S$ vs $Q$ for the NM state
(dash lines) and for the FM stat (solid lines), incorporating intervalley
scattering contribution.

Second, we discuss the magnetization direction of the equilibrium FM phase.
Consider a spin-flipped variational state where the pockets of each spin
channel are opposite to that in Fig. 2(a), as shown in Fig. S2(a). We then
calculate $\Delta E/S$ vs $Q$ for this state with $\downarrow$ magnetization
to compare with the FM state with $\uparrow$ magnetization [see Fig. S2(b)]. A
clear degeneracy in energy is found between the two cases, demonstrating a
degenerate FM state at equilibrium with reversed magnetization. In another
words, FM equilibrium states have no spin preference.

Furthermore, we study the NM variational state where the pockets are reversed
at the two valleys and also opposite in the two spin channels, as displayed in
Fig. S2(c). The FM state preserves $\mathrm{Z}_{2}$ symmetry but with broken
$\mathrm{SU}(2)$ symmetry. Whereas, the NM state breaks both. We then
calculate their energy density $\Delta E/S$ vs $Q$ by taking intervalley
scattering contribution into account, to see which is more energetically
favorable. The exchange energy density contributed by intervalley scattering
is given by

$\frac{E^{\mathrm{inter}}_{\mathrm{ex}}}{S}=-\frac{1}{2}\int_{D}\frac{d^{2}\bm{p}}{(2\pi)^{2}}\int_{D}\frac{d^{2}\bm{q}}{(2\pi)^{2}}\sum_{i,j}\sum_{m,n,\sigma}\sum_{a\neq
a^{\prime}}\chi^{i}_{mn,a^{\prime}a}(\bm{q},\bm{p})\chi^{j}_{nm,aa^{\prime}}(\bm{p},\bm{q})V_{ij}(\bm{p}-\bm{q}+\bm{K}_{a}-\bm{K}_{a^{\prime}})n_{ma^{\prime}\sigma}(\bm{q})n_{na\sigma}(\bm{p}).$
(S7)

Put this together with the intravalley contribution Eq.(11), we can calculate
the total $E_{\mathrm{ex}}/S$ and then $\Delta E/S$. As shown in Fig. S2(d),
$\Delta E/S$ of the FM state is lower than that of the NM state, manifesting
that the FM state is more favorable and stable at equilibrium than the NM
state. Notice that the energy differences between the two states are tiny,
which is due to that the intervalley scattering contribution to the exchange
energy $E^{\mathrm{inter}}_{\mathrm{ex}}$ is two or three orders of magnitudes
smaller than $E^{\mathrm{intra}}_{\mathrm{ex}}$ contributed by intravalley
scattering. The results given in Fig. S2(d) are slightly deviated from that
shown in Fig. S2(b) without $E^{\mathrm{inter}}_{\mathrm{ex}}$ contribution,
making the differences between the FM and NM states insignificant.

##  Appendix C Derivation of the spin supercurrent

We give here a detailed derivation of the spin supercurrent in the spin SC
heterojunction shown in Fig. 5, following the general non-equilibrium Green’s
function method. The equilibrium current flowing from lead $\beta$ to the
central region is given by

$I_{\beta}(t)=-e\langle\dot{\hat{N}}_{\beta}\rangle=\frac{ie}{\hbar}\langle[\hat{N}_{\beta},\bar{H}](t)\rangle,$
(S8)

where
$\hat{N}_{\beta}=\sum_{\bm{k},\eta}\alpha^{\dagger}_{\beta\bm{k}+\eta\uparrow}\alpha_{\beta\bm{k}+\eta\uparrow}+\alpha^{\dagger}_{\beta\bm{k}-\eta\downarrow}\alpha_{\beta\bm{k}-\eta\downarrow}$
is particle number operator of the lead, and $\bar{H}=H-\mu\hat{N}$ is the
thermodynamic potential of the junction with $\hat{N}$ the total particle
number operator. The commutator in the above equation is calculated by
$[\hat{N}_{\beta},\bar{H}]=\sum_{\bm{k},\eta}t_{\beta}\alpha^{\dagger}_{\beta\bm{k}+\eta\uparrow}d_{\uparrow}+t_{\beta}\alpha^{\dagger}_{\beta\bm{k}-\eta\downarrow}d_{\downarrow}-t^{\ast}_{\beta}d^{\dagger}_{\uparrow}\alpha_{\beta\bm{k}+\eta\uparrow}-t^{\ast}_{\beta}d^{\dagger}_{\downarrow}\alpha_{\beta\bm{k}-\eta\downarrow}$.
Then the current $I_{\beta}$ can be recast to

$\begin{split}I_{\beta}(t)=&\frac{ie}{\hbar}\sum_{\bm{k},\eta}t_{\beta}\langle\alpha^{\dagger}_{\beta\bm{k}+\eta\uparrow}(t)d_{\uparrow}(t)\rangle-t^{\ast}_{\beta}\langle
d^{\dagger}_{\uparrow}(t)\alpha_{\beta\bm{k}+\eta\uparrow}(t)\rangle+t_{\beta}\langle\alpha^{\dagger}_{\beta\bm{k}-\eta\downarrow}(t)d_{\downarrow}(t)\rangle-t^{\ast}_{\beta}\langle
d^{\dagger}_{\downarrow}(t)\alpha_{\beta\bm{k}-\eta\downarrow}(t)\rangle\\\
=&\frac{2e}{\hbar}\mathrm{Re}\sum_{\bm{k},\eta}t_{\beta}i\langle\alpha^{\dagger}_{\beta\bm{k}+\eta\uparrow}(t)d_{\uparrow}(t)\rangle+t_{\beta}i\langle\alpha^{\dagger}_{\beta\bm{k}-\eta\downarrow}(t)d_{\downarrow}(t)\rangle\\\
=&\frac{2e}{\hbar}\mathrm{Re}\sum_{\bm{k},\eta}t_{\beta}G^{<}_{d\beta\bm{k}\eta,11}(t,t)+t_{\beta}G^{<}_{d\beta\bm{k}\eta,22}(t,t).\end{split}$
(S9)

We can extract the spin-resolved currents from above equation, which read

$\begin{split}I_{\beta\uparrow}(t)=&\mathrm{Re}\frac{2et_{\beta}}{\hbar}\sum_{\bm{k},\eta}G^{<}_{d\beta\bm{k}\eta,11}(t,t);\\\
I_{\beta\downarrow}(t)=&\mathrm{Re}\frac{2et_{\beta}}{\hbar}\sum_{\bm{k},\eta}G^{<}_{d\beta\bm{k}\eta,22}(t,t),\end{split}$
(S10)

where the lesser Green’s functions are given as the diagonal elements of the
Green’s function matrix defined in the Nambu representation

$G^{<}_{d\beta\bm{k}\eta}(t,t^{\prime})=i\left(\begin{array}[]{cc}\langle\alpha^{\dagger}_{\beta\bm{k}+\eta\uparrow}(t^{\prime})d_{\uparrow}(t)\rangle&\langle\alpha^{\dagger}_{\beta\bm{k}-\eta\downarrow}(t^{\prime})d_{\uparrow}(t)\rangle\\\
\langle\alpha^{\dagger}_{\beta\bm{k}+\eta\uparrow}(t^{\prime})d_{\downarrow}(t)\rangle&\langle\alpha^{\dagger}_{\beta\bm{k}-\eta\downarrow}(t^{\prime})d_{\downarrow}(t)\rangle\\\
\end{array}\right).$ (S11)

Furthermore, we also define the lesser and advanced Green’s functions in the
Nambu basis for the lead $\beta$ without coupling to the quantum dot by

$\begin{split}g^{<}_{\beta\bm{k}\eta}(t,t^{\prime})=&i\left(\begin{array}[]{cc}\langle\alpha^{\dagger}_{\beta\bm{k}+\eta\uparrow}(t^{\prime})\alpha_{\beta\bm{k}+\eta\uparrow}(t)\rangle_{0}&\langle\alpha^{\dagger}_{\beta\bm{k}-\eta\downarrow}(t^{\prime})\alpha_{\beta\bm{k}+\eta\uparrow}(t)\rangle_{0}\\\
\langle\alpha^{\dagger}_{\beta\bm{k}+\eta\uparrow}(t^{\prime})\alpha_{\beta\bm{k}-\eta\downarrow}(t)\rangle_{0}&\langle\alpha^{\dagger}_{\beta\bm{k}-\eta\downarrow}(t^{\prime})\alpha_{\beta\bm{k}-\eta\downarrow}(t)\rangle_{0}\\\
\end{array}\right);\\\
g^{a}_{\beta\bm{k}\eta}(t,t^{\prime})=&i\theta(t^{\prime}-t)\left(\begin{array}[]{cc}\langle\\{\alpha_{\beta\bm{k}+\eta\uparrow}(t),\alpha^{\dagger}_{\beta\bm{k}+\eta\uparrow}(t^{\prime})\\}\rangle_{0}&\langle\\{\alpha_{\beta\bm{k}+\eta\uparrow}(t),\alpha^{\dagger}_{\beta\bm{k}-\eta\downarrow}(t^{\prime})\\}\rangle_{0}\\\
\langle\\{\alpha_{\beta\bm{k}-\eta\downarrow}(t),\alpha^{\dagger}_{\beta\bm{k}+\eta\uparrow}(t^{\prime})\\}\rangle_{0}&\langle\\{\alpha_{\beta\bm{k}-\eta\downarrow}(t),\alpha^{\dagger}_{\beta\bm{k}-\eta\downarrow}(t^{\prime})\\}\rangle_{0}\\\
\end{array}\right),\end{split}$ (S12)

and the lesser and retarded/advanced Green’s functions for the quantum dot
with coupling to the leads by

$\begin{split}G^{<}_{dd}(t,&t^{\prime})=i\left(\begin{array}[]{cc}\langle
d^{\dagger}_{\uparrow}(t^{\prime})d_{\uparrow}(t)\rangle&\langle
d^{\dagger}_{\downarrow}(t^{\prime})d_{\uparrow}(t)\rangle\\\ \langle
d^{\dagger}_{\uparrow}(t^{\prime})d_{\downarrow}(t)\rangle&\langle
d^{\dagger}_{\downarrow}(t^{\prime})d_{\downarrow}(t)\rangle\\\
\end{array}\right);\\\ G^{r/a}_{dd}(t,&t^{\prime})=\mp i\theta(\pm t\mp
t^{\prime})\left(\begin{array}[]{cc}\langle\\{d_{\uparrow}(t),d^{\dagger}_{\uparrow}(t^{\prime})\\}\rangle&\langle\\{d_{\uparrow}(t),d^{\dagger}_{\downarrow}(t^{\prime})\\}\rangle\\\
\langle\\{d_{\downarrow}(t),d^{\dagger}_{\uparrow}(t^{\prime})\\}\rangle&\langle\\{d_{\downarrow}(t),d^{\dagger}_{\downarrow}(t^{\prime})\\}\rangle\\\
\end{array}\right).\end{split}$ (S13)

We next solve the diagonal lesser Green’s functions in Eq. (S10) from the
Dyson’s equation:

$G^{<}=g^{<}+G^{r}\Sigma^{r}g^{<}+G^{r}\Sigma^{<}g^{a}+G^{<}\Sigma^{a}g^{a}.$
(S14)

Its non-diagonal matrix can be written as

$G^{<}_{d\beta\bm{k}\eta}(t,t^{\prime})=\frac{1}{\hbar}\int
dt_{1}~{}G^{r}_{dd}(t,t_{1})\Sigma^{r}_{d\beta\bm{k}\eta}g^{<}_{\beta\bm{k}\eta}(t_{1},t^{\prime})+G^{<}_{dd}(t,t_{1})\Sigma^{a}_{d\beta\bm{k}\eta}g^{a}_{\beta\bm{k}\eta}(t_{1},t^{\prime}),$
(S15)

where the retarded, advanced and lesser self-energies are given by

$\Sigma^{r,a}_{d\beta\bm{k}\eta}=\left(\begin{array}[]{cc}t^{\ast}_{\beta}&0\\\
0&t^{\ast}_{\beta}\\\
\end{array}\right),~{}~{}~{}\Sigma^{<}_{d\beta\bm{k}\eta}=\left(\begin{array}[]{cc}0&0\\\
0&0\\\ \end{array}\right).$ (S16)

Therefore, we can obtain the diagonal lesser Green’s functions from Eq. (S15)
as

$\begin{split}G^{<}_{d\beta\bm{k}\eta,11}(t,t^{\prime})=\frac{1}{\hbar}\int
dt_{1}&~{}G^{r}_{dd,11}(t,t_{1})\Sigma^{r}_{d\beta\bm{k}\eta,11}g^{<}_{\beta\bm{k}\eta,11}(t_{1},t^{\prime})+G^{r}_{dd,12}(t,t_{1})\Sigma^{r}_{d\beta\bm{k}\eta,22}g^{<}_{\beta\bm{k}\eta,21}(t_{1},t^{\prime})+\\\
&G^{<}_{dd,11}(t,t_{1})\Sigma^{a}_{d\beta\bm{k}\eta,11}g^{a}_{\beta\bm{k}\eta,11}(t_{1},t^{\prime})+G^{<}_{dd,12}(t,t_{1})\Sigma^{a}_{d\beta\bm{k}\eta,22}g^{a}_{\beta\bm{k}\eta,21}(t_{1},t^{\prime});\\\
G^{<}_{d\beta\bm{k}\eta,22}(t,t^{\prime})=\frac{1}{\hbar}\int
dt_{1}&~{}G^{r}_{dd,22}(t,t_{1})\Sigma^{r}_{d\beta\bm{k}\eta,22}g^{<}_{\beta\bm{k}\eta,22}(t_{1},t^{\prime})+G^{r}_{dd,21}(t,t_{1})\Sigma^{r}_{d\beta\bm{k}\eta,11}g^{<}_{\beta\bm{k}\eta,12}(t_{1},t^{\prime})+\\\
&G^{<}_{dd,22}(t,t_{1})\Sigma^{a}_{d\beta\bm{k}\eta,22}g^{a}_{\beta\bm{k}\eta,22}(t_{1},t^{\prime})+G^{<}_{dd,21}(t,t_{1})\Sigma^{a}_{d\beta\bm{k}\eta,11}g^{a}_{\beta\bm{k}\eta,12}(t_{1},t^{\prime}).\end{split}$
(S17)

Put above into Eq. (S10), we get the spin-resolved currents

$\begin{split}I_{\beta\uparrow}(t)=\frac{2e}{\hbar^{2}}\int
dt_{1}~{}\mathrm{Re}\sum_{\bm{k},\eta}|t_{\beta}|^{2}[&G^{r}_{dd,11}(t,t_{1})g^{<}_{\beta\bm{k}\eta,11}(t_{1},t)+G^{<}_{dd,11}(t,t_{1})g^{a}_{\beta\bm{k}\eta,11}(t_{1},t)+\\\
&G^{r}_{dd,12}(t,t_{1})g^{<}_{\beta\bm{k}\eta,21}(t_{1},t)+G^{<}_{dd,12}(t,t_{1})g^{a}_{\beta\bm{k}\eta,21}(t_{1},t)];\\\
I_{\beta\downarrow}(t)=\frac{2e}{\hbar^{2}}\int
dt_{1}~{}\mathrm{Re}\sum_{\bm{k},\eta}|t_{\beta}|^{2}[&G^{r}_{dd,22}(t,t_{1})g^{<}_{\beta\bm{k}\eta,22}(t_{1},t)+G^{<}_{dd,22}(t,t_{1})g^{a}_{\beta\bm{k}\eta,22}(t_{1},t)+\\\
&G^{r}_{dd,21}(t,t_{1})g^{<}_{\beta\bm{k}\eta,12}(t_{1},t)+G^{<}_{dd,21}(t,t_{1})g^{a}_{\beta\bm{k}\eta,12}(t_{1},t)]\end{split}$
(S18)

For the advanced and lesser Green’s functions of the spin SC lead without
coupling to the quantum dot, they can be solved under the wide-band
approximation, which are given by

$\begin{split}&\sum_{\bm{k}}g^{a}_{\beta\bm{k}\eta}(t,t^{\prime})\approx
i\theta(t^{\prime}-t)\rho_{N}(0)\int
d\varepsilon~{}\frac{e^{-i\varepsilon(t-t^{\prime})/\hbar}}{\sqrt{\varepsilon^{2}-\Delta^{2}}}\left(\begin{array}[]{cc}|\varepsilon|&\Delta
e^{i\phi_{\beta}}\frac{|\varepsilon|}{\varepsilon}\\\ \Delta
e^{-i\phi_{\beta}}\frac{|\varepsilon|}{\varepsilon}&|\varepsilon|\\\
\end{array}\right);\\\
&\sum_{\bm{k}}g^{<}_{\beta\bm{k}\eta}(t,t^{\prime})\approx i\rho_{N}(0)\int
d\varepsilon~{}f(\varepsilon)\frac{|\varepsilon|}{\sqrt{\varepsilon^{2}-\Delta^{2}}}e^{-i\varepsilon(t-t^{\prime})/\hbar}\left(\begin{array}[]{cc}1&\Delta
e^{i\phi_{\beta}}\frac{1}{\varepsilon}\\\ \Delta
e^{-i\phi_{\beta}}\frac{1}{\varepsilon}&1\\\ \end{array}\right),\end{split}$
(S19)

where $\rho^{N}(0)$ is the density of states of a single lead in the norm
state ($\Delta=0$) at the Fermi energy $E_{f}=0$. Insert this equation into
Eq. (S18), take the Fourier transformation of the Green’s functions from time
domain to frequency domain and use the relation
$G^{<}_{dd}(\varepsilon)=f(\varepsilon)[G^{a}_{dd}(\varepsilon)-G^{r}_{dd}(\varepsilon)]$,
we obtain

$\begin{split}I_{\beta\uparrow}(t)\approx&-\frac{4e}{\hbar^{2}}\mathrm{Im}\int\frac{d\varepsilon}{2\pi}~{}\Gamma\frac{|\varepsilon|}{2\sqrt{\varepsilon^{2}-\Delta^{2}}}f(\varepsilon)\left\\{[G^{a}_{dd,11}(\varepsilon)+G^{r}_{dd,11}(\varepsilon)]+\frac{\Delta
e^{-i\phi_{\beta}}}{\varepsilon}[G^{a}_{dd,12}(\varepsilon)+G^{r}_{dd,12}(\varepsilon)]\right\\};\\\
I_{\beta\downarrow}(t)\approx&-\frac{4e}{\hbar^{2}}\mathrm{Im}\int\frac{d\varepsilon}{2\pi}~{}\Gamma\frac{|\varepsilon|}{2\sqrt{\varepsilon^{2}-\Delta^{2}}}f(\varepsilon)\left\\{[G^{a}_{dd,22}(\varepsilon)+G^{r}_{dd,22}(\varepsilon)]+\frac{\Delta
e^{+i\phi_{\beta}}}{\varepsilon}[G^{a}_{dd,21}(\varepsilon)+G^{r}_{dd,21}(\varepsilon)]\right\\},\end{split}$
(S20)

where $\Gamma=2\pi|t_{\beta}|^{2}\rho^{N}(0)$ is the linewidth function of the
lead. When $\Delta=0$, the first bracket in the above equations represents the
current tunneling from normal metal to norm metal at equilibrium, which should
be zero as indeed. For $\Delta\neq 0$, this term contributes small and can be
dropped. Therefore, the spin-resolved currents reduce to

$\begin{split}&I_{\beta\uparrow}(t)\approx-\frac{4e}{\hbar^{2}}\mathrm{Im}\int\frac{d\varepsilon}{2\pi}~{}\Gamma^{\beta}\frac{\Delta
e^{-i\phi_{\beta}}}{2\sqrt{\varepsilon^{2}-\Delta^{2}}}\frac{|\varepsilon|}{\varepsilon}f(\varepsilon)[G^{a}_{dd,12}(\varepsilon)+G^{r}_{dd,12}(\varepsilon)];\\\
&I_{\beta\downarrow}(t)\approx-\frac{4e}{\hbar^{2}}\mathrm{Im}\int\frac{d\varepsilon}{2\pi}~{}\Gamma^{\beta}\frac{\Delta
e^{+i\phi_{\beta}}}{2\sqrt{\varepsilon^{2}-\Delta^{2}}}\frac{|\varepsilon|}{\varepsilon}f(\varepsilon)[G^{a}_{dd,21}(\varepsilon)+G^{r}_{dd,21}(\varepsilon)],\end{split}$
(S21)

To evaluate the retarded Green’s function $G^{r}_{dd}(\varepsilon)$ of the
quantum dot, let us first calculate the Green’s function
$g^{r}_{dd}(\varepsilon)$ of the dot without coupling to the leads. It can be
given by

$g^{r}_{dd}(\varepsilon)=\left[(\varepsilon+i0^{\dagger})I-\left(\begin{array}[]{cc}\varepsilon_{0\uparrow}&0\\\
0&\varepsilon_{0\downarrow}\\\
\end{array}\right)\right]^{-1}=\left(\begin{array}[]{cc}\frac{1}{\varepsilon-\varepsilon_{0\uparrow}+i0^{\dagger}}&0\\\
0&\frac{1}{\varepsilon-\varepsilon_{0\downarrow}+i0^{\dagger}}\\\
\end{array}\right).$ (S22)

And the retarded self-energies of the quantum dot contributed by the leads are

$\Sigma^{r}(\varepsilon)=\sum_{\beta}\Sigma^{r}_{\beta}(\varepsilon)=\sum_{\beta}|t_{\beta}|^{2}\sum_{\bm{k},a}g^{r}_{\beta\bm{k}a\beta\bm{k}a}(\varepsilon)=-i\sum_{\beta}\frac{\Gamma}{\sqrt{\varepsilon^{2}-\Delta^{2}}}\left(\begin{array}[]{cc}|\varepsilon|&\Delta
e^{i\phi_{\beta}}\frac{|\varepsilon|}{\varepsilon}\\\ \Delta
e^{-i\phi_{\beta}}\frac{|\varepsilon|}{\varepsilon}&|\varepsilon|\\\
\end{array}\right).$ (S23)

Then, the Green’s function of the quantum dot with coupling to the leads can
be given by

$G^{r}_{dd}(\varepsilon)=\left[(\varepsilon+i0^{\dagger})I-\left(\begin{array}[]{cc}\varepsilon_{0\uparrow}&0\\\
0&\varepsilon_{0\downarrow}\\\
\end{array}\right)-\Sigma^{r}(\varepsilon)\right]^{-1}=\frac{1}{B(\varepsilon)}\left(\begin{array}[]{cc}g^{r,-1}_{00,22}-\Sigma^{r}_{22}&\Sigma^{r}_{12}\\\
\Sigma^{r}_{21}&g^{r,-1}_{00,11}-\Sigma^{r}_{11}\\\ \end{array}\right),$ (S24)

where
$B(\varepsilon)=(g^{r,-1}_{00,11}-\Sigma^{r}_{11})(g^{r,-1}_{00,22}-\Sigma^{r}_{22})-\Sigma^{r}_{12}\Sigma^{r}_{21}$.
Bring $G^{r}_{dd}(\varepsilon)$ into Eq. (S21) and use
$G^{a}_{dd}(\varepsilon)=[G^{r}_{dd}(\varepsilon)]^{\dagger}$, we arrive at
the final form of the equilibrium spin-resolved currents, flowing in lead
$\beta=L$

$\begin{split}&I_{L\uparrow}\approx\frac{4e}{\hbar}\int\frac{d\varepsilon}{2\pi}~{}\frac{\Gamma^{2}\Delta^{2}}{\varepsilon^{2}-\Delta^{2}}f(\varepsilon)\mathrm{Im}\frac{1}{B^{\ast}(\varepsilon)}\mathrm{sin}(\phi_{R}-\phi_{L});\\\
&I_{L\downarrow}\approx\frac{4e}{\hbar}\int\frac{d\varepsilon}{2\pi}~{}\frac{\Gamma^{2}\Delta^{2}}{\varepsilon^{2}-\Delta^{2}}f(\varepsilon)\mathrm{Im}\frac{1}{B^{\ast}(\varepsilon)}\mathrm{sin}(\phi_{L}-\phi_{R}),\end{split}$
(S25)

which is just Eq. (32) presented in the main text.

## References

  * van Delft and Kes [2010] D. van Delft and P. Kes, Phys. Today 63, 38 (2010).
  * Bardeen  _et al._ [1957] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. Rev. 108, 1175 (1957).
  * Josephson [1962] B. Josephson, Phys. Lett. 1, 251 (1962).
  * Tsuei and Kirtley [2000] C. C. Tsuei and J. R. Kirtley, Rev. Mod. Phys. 72, 969 (2000).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Sato and Ando [2017] M. Sato and Y. Ando, Rep. Prog. Phys. 80, 076501 (2017).
  * Zhou  _et al._ [2021] X. Zhou, W.-S. Lee, M. Imada, N. Trivedi, P. Phillips, H.-Y. Kee, P. Törmä, and M. Eremets, Nat. Rev. Phys. 3, 462 (2021).
  * Denholme [2023] S. Denholme, Nat. Rev. Phys. 5, 557 (2023).
  * Hahn  _et al._ [2019] S. Hahn, K. Kim, K. Kim, X. Hu, T. Painter, I. Dixon, S. Kim, K. R. Bhattarai, S. Noguchi, J. Jaroszynski, and D. C. Larbalestier, Nature 570, 496 (2019).
  * MacManus-Driscoll and Wimbush [2021] J. L. MacManus-Driscoll and S. C. Wimbush, Nat. Rev. Mater. 6, 587 (2021).
  * Eschrig [2011] M. Eschrig, Phys. Today 64, 43 (2011).
  * Eschrig [2015] M. Eschrig, Rep. Prog. Phys. 78, 104501 (2015).
  * Linder and Robinson [2015] J. Linder and J. W. A. Robinson, Nat. Phys. 11, 307 (2015).
  * Fert [2008] A. Fert, Rev. Mod. Phys. 80, 1517 (2008).
  * Lin  _et al._ [2019] X. Lin, W. Yang, K. L. Wang, and W. Zhao, Nat. Electron. 2, 274 (2019).
  * Buzdin [2005] A. I. Buzdin, Rev. Mod. Phys. 77, 935 (2005).
  * Bergeret  _et al._ [2005] F. S. Bergeret, A. F. Volkov, and K. B. Efetov, Rev. Mod. Phys. 77, 1321 (2005).
  * Aoki  _et al._ [2001] D. Aoki, A. Huxley, E. Ressouche, D. Braithwaite, J. Flouquet, J.-P. Brison, E. Lhotel, and C. Paulsen, Nature 413, 613 (2001).
  * Huy  _et al._ [2007] N. T. Huy, A. Gasparini, D. E. de Nijs, Y. Huang, J. C. P. Klaasse, T. Gortenmulder, A. de Visser, A. Hamann, T. Görlach, and H. v. Löhneysen, Phys. Rev. Lett. 99, 067006 (2007).
  * Ran  _et al._ [2019] S. Ran, C. Eckberg, Q.-P. Ding, Y. Furukawa, T. Metz, S. R. Saha, I.-L. Liu, M. Zic, H. Kim, J. Paglione, and N. P. Butch, Science 365, 684 (2019).
  * Mackenzie and Maeno [2003] A. P. Mackenzie and Y. Maeno, Rev. Mod. Phys. 75, 657 (2003).
  * Kallin and Berlinsky [2009] C. Kallin and A. J. Berlinsky, J. Phys.: Condens. Mat. 21, 164210 (2009).
  * Yang  _et al._ [2021] J. Yang, J. Luo, C. Yi, Y. Shi, Y. Zhou, and G. qing Zheng, Sci. Advan. 7, eabl4432 (2021).
  * Sun  _et al._ [2011] Q.-f. Sun, Z.-t. Jiang, Y. Yu, and X. C. Xie, Phys. Rev. B 84, 214501 (2011).
  * Liu  _et al._ [2012] H. Liu, H. Jiang, X. C. Xie, and Q.-f. Sun, Phys. Rev. B 86, 085441 (2012).
  * Sun and Xie [2013] Q.-f. Sun and X. C. Xie, Phys. Rev. B 87, 245427 (2013).
  * Jiang  _et al._ [2020] Z. Jiang, W. Lou, Y. Liu, Y. Li, H. Song, K. Chang, W. Duan, and S. Zhang, Phys. Rev. Lett. 124, 166401 (2020).
  * Wang  _et al._ [2013] Z.-b. Wang, Q.-f. Sun, and X. Xie, Eur. Phys. J. B 86, 496 (2013).
  * Bao  _et al._ [2013] Z.-q. Bao, X. Xie, and Q.-f. Sun, Nat. Commun. 4, 2951 (2013).
  * Sethi  _et al._ [2021] G. Sethi, Y. Zhou, L. Zhu, L. Yang, and F. Liu, Phys. Rev. Lett. 126, 196403 (2021).
  * Sethi  _et al._ [2023] G. Sethi, M. Cuma, and F. Liu, Phys. Rev. Lett. 130, 186401 (2023).
  * Yuan  _et al._ [2018] W. Yuan, Q. Zhu, T. Su, Y. Yao, W. Xing, Y. Chen, Y. Ma, X. Lin, J. Shi, R. Shindou, X. C. Xie, and W. Han, Sci. Adv. 4, eaat1098 (2018).
  * Afonso and Kuneš [2017] J. F. Afonso and J. Kuneš, Phys. Rev. B 95, 115131 (2017).
  * Ikeda  _et al._ [2023] A. Ikeda, Y. H. Matsuda, K. Sato, Y. Ishii, H. Sawabe, D. Nakamura, S. Takeyama, and J. Nasu, Nat. Commun. 14, 1744 (2023).
  * Varsano  _et al._ [2020] D. Varsano, M. Palummo, E. Molinari, and M. Rontani, Nat. Nanotechnol. 15, 367 (2020).
  * Yang  _et al._ [2024] H. Yang, J. Zeng, Y. Shao, Y. Xu, X. Dai, and X.-Z. Li, Phys. Rev. B 109, 075167 (2024).
  * Wu  _et al._ [2022] Y. Wu, H. Jiang, H. Chen, H. Liu, J. Liu, and X. C. Xie, Phys. Rev. Lett. 128, 106804 (2022).
  * Wu  _et al._ [2023] Y. Wu, J. Liu, and X. Xie, Sci. China Phys., Mech. & Astr. 66, 267004 (2023).
  * Deveaud  _et al._ [1991] B. Deveaud, F. Clérot, N. Roy, K. Satzke, B. Sermage, and D. S. Katzer, Phys. Rev. Lett. 67, 2355 (1991).
  * Palummo  _et al._ [2015] M. Palummo, M. Bernardi, and J. C. Grossman, Nano Lett. 15, 2794 (2015).
  * Haugen  _et al._ [2008] H. Haugen, D. Huertas-Hernando, and A. Brataas, Phys. Rev. B 77, 115406 (2008).
  * Kerelsky  _et al._ [2021] A. Kerelsky, C. Rubio-Verdú, L. Xian, D. M. Kennes, D. Halbertal, N. Finney, L. Song, S. Turkel, L. Wang, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, D. N. Basov, A. Rubio, and A. N. Pasupathy, Proc. Natl. Acad. Sci. 118, e2017366118 (2021).
  * Liu  _et al._ [2024] K. Liu, J. Zheng, Y. Sha, B. Lyu, F. Li, Y. Park, Y. Ren, K. Watanabe, T. Taniguchi, J. Jia, W. Luo, Z. Shi, J. Jung, and G. Chen, Nat. Nanotechnol. 19, 188 (2024).
  * Sha  _et al._ [2024] Y. Sha, J. Zheng, K. Liu, H. Du, K. Watanabe, T. Taniguchi, J. Jia, Z. Shi, R. Zhong, and G. Chen, Science 384, 414 (2024).
  * Myhro  _et al._ [2018] K. Myhro, S. Che, Y. Shi, Y. Lee, K. Thilahar, K. Bleich, D. Smirnov, and C. N. Lau, 2D Mater. 5, 045013 (2018).
  * Nilsson  _et al._ [2006] J. Nilsson, A. H. Castro Neto, N. M. R. Peres, and F. Guinea, Phys. Rev. B 73, 214418 (2006).
  * Castro Neto  _et al._ [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  * Dutreix and Katsnelson [2016] C. Dutreix and M. I. Katsnelson, Phys. Rev. B 93, 035413 (2016).
  * Das Sarma  _et al._ [2011] S. Das Sarma, S. Adam, E. H. Hwang, and E. Rossi, Rev. Mod. Phys. 83, 407 (2011).
  * van Gelderen  _et al._ [2011] R. van Gelderen, L.-K. Lim, and C. M. Smith, Phys. Rev. B 84, 155446 (2011).
  * Kotov  _et al._ [2012] V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea, and A. H. Castro Neto, Rev. Mod. Phys. 84, 1067 (2012).
  * Min  _et al._ [2012] H. Min, E. H. Hwang, and S. Das Sarma, Phys. Rev. B 86, 081402 (2012).
  * Peres  _et al._ [2005] N. M. R. Peres, F. Guinea, and A. H. Castro Neto, Phys. Rev. B 72, 174406 (2005).
  * Stoner [1947] E. C. Stoner, Rep. Prog. Phys. 11, 43 (1947).
  * Stepanov  _et al._ [2020] P. Stepanov, I. Das, X. Lu, A. Fahimniya, K. Watanabe, T. Taniguchi, F. H. L. Koppens, J. Lischner, L. Levitov, and D. K. Efetov, Nature 583, 375 (2020).
  * Liu  _et al._ [2021] X. Liu, Z. Wang, K. Watanabe, T. Taniguchi, O. Vafek, and J. I. A. Li, Science 371, 1261 (2021).
  * Fulde and Ferrell [1964] P. Fulde and R. A. Ferrell, Phys. Rev. 135, A550 (1964).
  * Kumagai  _et al._ [2006] K. Kumagai, M. Saitoh, T. Oyaizu, Y. Furukawa, S. Takashima, M. Nohara, H. Takagi, and Y. Matsuda, Phys. Rev. Lett. 97, 227002 (2006).
  * Matsuda and Shimahara [2007] Y. Matsuda and H. Shimahara, J. Phys. Soc. Jpn 76, 051005 (2007).
  * Kinnunen  _et al._ [2018] J. J. Kinnunen, J. E. Baarsma, J.-P. Martikainen, and P. Törmä, Rep. Prog. Phys. 81, 046401 (2018).
  * Kasahara  _et al._ [2020] S. Kasahara, Y. Sato, S. Licciardello, M. Čulo, S. Arsenijević, T. Ottenbros, T. Tominaga, J. Böker, I. Eremin, T. Shibauchi, J. Wosnitza, N. E. Hussey, and Y. Matsuda, Phys. Rev. Lett. 124, 107001 (2020).
  * Wehling  _et al._ [2011] T. O. Wehling, E. Şaşıoğlu, C. Friedrich, A. I. Lichtenstein, M. I. Katsnelson, and S. Blügel, Phys. Rev. Lett. 106, 236805 (2011).