# Stabilizing fluctuating spin-triplet superconductivity in graphene via
induced spin-orbit coupling

Jonathan B. Curtis<EMAIL_ADDRESS>College of Letters and Science,
University of California, Los Angeles, CA 90095, USA  Department of Physics,
Harvard University, Cambridge, MA 02138, USA  Nicholas R. Poniatowski
Yonglong Xie  Amir Yacoby  Department of Physics, Harvard University,
Cambridge, MA 02138, USA  Eugene Demler  Institute for Theoretical Physics,
ETH Zürich, 8093 Zürich, Switzerland  Prineha Narang<EMAIL_ADDRESS>College
of Letters and Science, University of California, Los Angeles, CA 90095, USA

###### Abstract

A recent experiment showed that proximity induced Ising spin-orbit coupling
enhances the spin-triplet superconductivity in Bernal bilayer graphene. Here,
we show that, due to the nearly perfect spin rotation symmetry of graphene,
the fluctuations of the spin orientation of the triplet order parameter
suppress the superconducting transition to nearly zero temperature. Our
analysis shows that both Ising spin-orbit coupling and in-plane magnetic field
can eliminate these low-lying fluctuations and can greatly enhance the
transition temperature, consistent with the recent experiment. Our model also
suggests the possible existence of a phase at small anisotropy and magnetic
field which exhibits quasi-long-range ordered spin-singlet charge 4e
superconductivity, even while the triplet 2e superconducting order only
exhibits short-ranged correlations. Finally, we discuss relevant experimental
signatures.

Introduction —Graphene based two-dimensional materials have emerged as a
highly tunable platform for studying spin-triplet superconductivity. Notable
examples include rhombohedral trilayer graphene (RTG) [1], moiré graphene [2,
3], and Bernal bilayer graphene (BBG) [4, 5]. In the case of BBG,
superconductivity was initially found to emerge exclusively in the presence of
an in-plane magnetic field that exceeds the paramagnetic limit, consistent
with a spin triplet order parameter [4]. A very recent experiment showed that
introducing a tungsten diselenide (WSe2) layer —a semiconducting material with
strong spin-orbit coupling —on top of BBG, as schematically depicted in Fig.
1(a), can stabilize superconductivity at strictly zero magnetic field and can
enhance the critical temperature by an order of magnitude [5].

An important task prompted by this recent experiment is to identify the role
of the proximity-induced spin-orbit coupling in promoting spin-triplet
superconductivity in BBG. One possibility is that spin-orbit coupling selects
a particular symmetry-broken parent state from which the superconductivity
emerges, and previous work [6, 7, 8] shows that the intervalley coherent
ground state is conducive to spin triplet pairing. Alternatively, the presence
of WSe2 can introduce virtual tunneling process that enhances pair binding
energies and hence boosts the critical temperature [9].

In this work, we present a novel mechanism by which the proximity-induced
spin-orbit coupling promotes the spin-triplet superconductivity observed in
graphene. Specifically, we show that the presence of Ising spin-orbit coupling
breaks the nearly perfect spin-rotation symmetry in graphene and thus
suppresses the order parameter fluctuations caused by the otherwise gapless
Goldstone modes corresponding to fluctuations of the orientation (in spin
space) of the triplet order parameter (for a precise definition of Ising-type
spin-orbit coupling see the discussion below or Supplemental Material). While
triplet superconductivity is formally impossible in intrinsic graphene on
account of the Mermin-Wagner theorem, it is allowed via a BKT-type transition
once spin-rotation symmetry is broken by the WSe2. A similar proposal was
recently put forth in Ref. [10], motivated by the potential observations of
triplet pairing in graphene.

To some extent, this phenomenology is relevant to all reports of spin-triplet
superconductivity in graphene and may offer a degree of universality once one
takes superconducting fluctuations in to account. In particular, we expect
that all perturbations which break spin-rotation symmetry, such as in-plane
magnetic field or induced Ising spin-orbit coupling, can help stabilize
fluctuating spin-triplet superconductivity; a magnetic field stabilizes a
spin-polarized condensate [10], whereas Ising spin-orbit coupling stabilizes a
unitary (spin nematic) triplet state. By suppressing order-parameter
fluctuations, these perturbations help restore the mean-field pairing
temperature. Our theory is based entirely on general symmetry arguments and is
therefore broadly applicable to many graphene-based triplet superconductors
including BBG, RTG, and moiré systems and highlights the important role of
collective modes in designing spin-triplet superconductors with higher
critical temperatures.

Figure 1:  Schematic overview of system considered. (a) Illustration of Bernal
bilayer graphene (BBG) on top of WSe2. The WSe2 proximity induces an Ising
spin-orbit coupling in the graphene by virtual tunneling processes. (b)
Pairing in the graphene Brillouin zone, for the case of $f$-wave triplet
pairing. Pairing is between $K$ and $K^{\prime}$ so that electrons are in
Kramers pairs with total momentum zero while Ising spin-orbit coupling acts
like a Zeeman field for spins pointing up at $K$ and down at $K^{\prime}$. (c)
Illustration of how spin-orbit coupling governs the scale of the d-vector
pinning $\Gamma_{\perp}$. For small $\Gamma_{\perp}$ the fluctuations of the
orientation of $\mathbf{d}$ are soft, whereas for large $\Gamma_{\perp}$ the
$\mathbf{d}$-vector is strongly pinned along the $z$-axis. (d) BCS calculation
of dimensionless pinning energy
$\tilde{\Gamma}_{\perp}=\Gamma_{\perp}/N(E_{F})$ as a function of the induced
Ising spin-orbit coupling $\lambda$, relative to the BCS pairing temperature.
One can see that $\tilde{\Gamma}_{\perp}$ increases by over two orders of
magnitude.

Microscopic Model—We begin by considering a minimal microscopic model for BBG,
consisting of two spin-degenerate parabolic bands at valleys $\mathbf{K}$ and
$\mathbf{K}^{\prime}$ in the presence of an induced Ising spin orbit coupling
of strength $\lambda$. This is depicted in Fig. 1(b) which shows the
intervalley pairing and induced spin-orbit coupling in the Brillouin zone.

Introducing the four-component spinor $\psi_{\mathbf{k}}$ and Pauli matrices
$\sigma_{i}$ ($\rho_{i}$) which act on the spin (valley) degrees of freedom,
we can write the normal state Hamiltonian as
$H=\sum_{\mathbf{k}}\psi_{\mathbf{k}}^{\dagger}h_{\mathbf{k}}\psi_{\mathbf{k}}$
with
$h_{\mathbf{k}}=\xi_{\mathbf{k}}\rho_{0}\sigma_{0}+\frac{1}{2}\lambda\rho_{3}\sigma_{3}$
where $\xi_{\mathbf{k}}$ is the quadratic dispersion of each band. To include
the effects of triplet pairing, we express the Hamiltonian in terms of the
eight-component Nambu spinor
$\Psi_{\mathbf{k}}=\left(\psi_{\mathbf{k}}\;\bar{\psi}_{-\mathbf{k}}i\sigma_{2}\right)^{T}$
as
$H=\frac{1}{2}\sum_{\mathbf{k}}\Psi^{\dagger}_{\mathbf{k}}\hat{\mathcal{H}}_{\mathbf{k}}\Psi_{\mathbf{k}}$
with

$\hat{\mathcal{H}}_{\mathbf{k}}=\begin{pmatrix}\xi_{\mathbf{k}}+\frac{1}{2}\lambda\rho_{3}\sigma_{3}&\mathbf{d}\cdot\bm{\sigma}i\rho_{2}\\\
-\mathbf{d}^{\star}\cdot\bm{\sigma}i\rho_{2}&-\xi_{\mathbf{k}}+\frac{1}{2}\lambda\rho_{3}\sigma_{3}\end{pmatrix}\,.$
(1)

This model features vector order parameter
$\mathbf{d}=g_{f}\sum_{\mathbf{k}}\langle\psi_{\mathbf{k}}^{T}\left(\rho_{2}\sigma_{2}\bm{\sigma}\right)\psi_{-\mathbf{k}}\rangle$
which describes an isotropically gapped $f$-wave, valley-singlet, spin-triplet
state such that the gap takes the opposite sign in the $\mathbf{K}$ and
$\mathbf{K}^{\prime}$ valleys 111The pairing is assumed isotropic in
$\mathbf{k}$ with respect to the $\mathbf{K,K^{\prime}}$ points, but is valley
singlet meaning that under $C_{3}$ rotations the valleys exchange and the
wavefunction is overall odd, making it $f$-wave..

This spin-singlet, valley-triplet order parameter has recently been argued to
generally correspond to the most energetically favorable state in graphene-
based systems [12, 13, 8], but the analysis to follow should generalize to
other triplet order parameters. Our crucial assumption is that in the absence
of spin-orbit coupling the system has perfect $SU(2)$ symmetry, and that all
subdominant pairing interactions are negligible such that, in particular,
there are no valley collective modes (which unlike spin are not guaranteed to
exist by symmetry). We leave the case of specific broken symmetry normal
states to future work. This assumption allows us to project the spin-orbit
coupling onto the triplet manifold while neglecting the mixing between singlet
and triplet subspaces.

Using this Hamiltonian, we can carry out the standard procedure for deriving
the Ginzburg-Landau free-energy for the triplet order parameter $\mathbf{d}$
near the critical point (see Supplemental Material for details). We work
within the weak-coupling BCS limit 222In fact, it is likely that BBG is not in
the weak coupling regime and this may lead to additional interesting effects,
though the argument we advance is based on symmetry and is still a relevant
consideration in such case., and find the free-energy density

$f=K|\nabla\mathbf{d}|^{2}+r|\mathbf{d}|^{2}+\Gamma_{\perp}|\mathbf{d}_{\perp}|^{2}+\frac{1}{2}u|\mathbf{d}|^{4}+\frac{1}{2}v(i\mathbf{d}\times\mathbf{d}^{\star})^{2}.$
(2)

The first term with $K$ gives the superconducting phase stiffness, the second
term is the bare pairing susceptibility, with $r\sim\log(T/T_{\rm BCS})$
changing sign at the mean-field BCS temperature, and $\Gamma_{\perp}$ is the
spin-orbit induced pinning energy for the out-of-plane orientations, with
$|\mathbf{d}_{\perp}|^{2}=|d_{x}|^{2}+|d_{y}|^{2}$ (see discussion below). The
nonlinearity $u$ respects an enlarged $SU(3)$ symmetry while the nonlinearity
$v=u$ further penalizes the non-unitary equal-spin pairing states, and breaks
the symmetry down to $SO(3)$.

Crucially, for nonzero symmetry-breaking spin-orbit coupling $\lambda$ (see
Hamiltonian in Eq. 1 for definition) we find a finite $\mathbf{d}$-vector
pinning energy $\Gamma_{\perp}$ which selects the $z$-axis as the favored
orientation of the $\mathbf{d}$-vector, depicted in Fig. 1(c). This pinning
energy, $\Gamma_{\perp}$ is plotted as a function of $\lambda$ in Fig. 1(d).
For small $\lambda/T_{\rm BCS}$ we find $\Gamma_{\perp}\sim\lambda^{2}/T_{\rm
BCS}^{2}$, while for larger values the growth slows to a logarithmic scaling
$\sim\log(\lambda)$. We note that the value of the induced spin-orbit coupling
$\lambda\sim 0.7$ meV inferred from experiments [5] is likely to be one of the
larger energy scales in the problem and may possibly rival the Fermi energy of
the small pockets, pointing to the need for a more detailed treatment in
future work. For small $\Gamma_{\perp}$, the actual transition temperature
$T_{c}$ will be suppressed due to strong fluctuations, as we now show.

Goldstone Modes—We now explore the phase diagram of the Ginzburg-Landau
functional in Eq. (2) as a function of temperature $T$ and the pinning energy
$\Gamma_{\perp}$. In fact, similar models have already been studied
extensively in the context of BKT transitions in ultracold spinor condensates
[15, 16, 17], where Monte Carlo studies have been conducted for a range of
possible anisotropy and interaction parameters [17]. In this work, we will
focus on the qualitative effect that Gaussian fluctuations have on the phase
diagram, rather than quantitatively determine the phase diagram since there is
anyways still great debate about the model and parameters, which can also vary
from system to system.

We treat this system using a large-$N$ expansion in the number of spin
components, as detailed in the Supplemental Material. Since this is a Gaussian
approximation, it does not capture topological defects such as vortices and is
expected to underestimate the role of fluctuations. This approach does not
allow us to distinguish between long-range ordered phases and quasi-long-range
ordered (QLRO) phases. However, it is known that in two-dimensions at finite
temperature fluctuations of the superconducting phase do not permit true long-
range order, and thus when we discuss the symmetry breaking phases with
$\langle d_{z}\rangle\neq 0$ it is understood that these are in reality QLRO.
Nevertheless, this approach allows us to study the impact of spin-orbit
coupling on fluctuations in a controlled way by finding the renormalized
critical temperature $T_{c}$ where the $\langle d_{z}\rangle$ QLRO condensate
develops. We also find a small region of $4e$ QLRO superconductivity at small
anisotropy which will be discussed below.

The key quantity of interest is the condensate $\psi\sim\langle d_{z}\rangle$.
Within mean-field theory, this is simply the solution to the saddle-point
equation and is given by $|\psi|=\sqrt{-r/u}$ when $T<T_{\rm BCS}$ and zero
above. Once the order-parameter fluctuations are taken in to account [18, 19],
we find the renormalized critical temperature $T_{c}$ given by

$\log\left(\frac{T_{c}}{T_{\rm BCS}}\right)=-{\rm
Gi}_{(2)}\log\left(\frac{C}{\tilde{\Gamma}_{\perp}}\right).$ (3)

Here $\tilde{\Gamma}_{\perp}=\Gamma_{\perp}/N(E_{F})$ is the unitless pinning
energy relative to the density-of-states at the Fermi level, $C$ is a unitless
ultraviolet cutoff on fluctuations given by $C=Kq_{\rm max}^{2}/N(E_{F})\sim
E_{F}^{2}/T_{\rm BCS}^{2}\sim 1/\textrm{Gi}_{(2)}^{2}$. The important
parameter here is the Ginzburg-Levanyuk [18, 19] number

$\textrm{Gi}_{(2)}=T_{\rm BCS}\frac{u}{2\pi KN(E_{F})}=\frac{T_{\rm BCS}}{\pi
v_{F}^{2}N(E_{F})}=\frac{T_{\rm BCS}}{E_{F}},$ (4)

which is a measure of the strength of superconducting fluctuations.

Solving this for $T_{c}$ we find a power-law behavior in $\Gamma_{\perp}$
scaling as $T_{c}\sim\Gamma_{\perp}^{\textrm{Gi}_{(2)}}$. The Ginzburg-
Levanyuk parameter is not known well experimentally in this system since it
requires a knowledge of the “pseudogap” temperature characterized by a
spectral gap appearing at $T_{\rm BCS}$, and spectroscopic experiments have
not yet been performed on this material. In other graphene-based
superconductors however, it has been argued that this may be quite large with
values of $T_{\rm BCS}/T_{F}\sim 0.1$ [20]. Tunneling spectroscopy
measurements have also found evidence for a substantial quasiparticle
pseudogap [21, 22]. Furthermore, superconductivity in BBG occurs proximal to
different Fermi surface topologies which may feature small Fermi surfaces with
small Fermi energies, which may only be as small as $0.6$ meV (for details see
Supplemental Material). We will roughly estimate a pairing temperature of
$T_{\rm BCS}\sim 0.5$ K 333While this may be somewhat optimistic, closer
examination of the data in Ref. [5], and in particular the inset of Fig. 2(d),
reveals a possible Azlamazov-Larkin downturn in $R_{xx}$ around $T\sim 500$
mK, indicating this may be a serious possibility., and a Fermi temperature of
$T_{F}\sim 0.6\;\textrm{meV}\sim 7$ K, in line with strong-coupling, with
$\textrm{Gi}_{(2)}\sim 0.07$. Intrinsic spin-orbit coupling in graphene is
expected to be no more than 40 $\mu$eV or so [24], and once placed on top of
WSe2 is estimated to be of order $0.7$ meV, corresponding to at least a
20-fold increase in $\lambda$, meaning an increase in ${\Gamma}_{\perp}$ of
order $60$-fold is feasible 444This is a lower bound on the increase. In
particular, it is thought that the intrinsic spin-orbit coupling in graphene
takes the form of $H_{\rm soc}=\lambda\sigma_{3}\rho_{3}s_{3}$ where $s_{a}$
are the sublattice Pauli matrices. Thus, in AB stacked graphene the average
spin-orbit coupling still vanishes and the intrinsic value is likely much less
than 40 $\mathrm{\SIUnitSymbolMicro eV}$ to begin with.. This could increase
$T_{c}$ by a factor of 33%, however since we have only gone to leading order
in $1/N$ (recall $N=3$), our approximations likely underestimate the effect of
fluctuations and give a conservative estimate.

This behavior is illustrated schematically in Fig. 2. We briefly mention that
this is the phase diagram before including the topological vortex defects,
which are responsible for the BKT transition known to occur in all two-
dimensional superfluids. We qualitatively expect that this will lead to a
further renormalization of the boundary of the superfluid phase [17] with the
relevant BKT temperature based on the topological classification of the
defects. Our estimates therefore indicate that at least some part of the
observed spin-orbit induced enhancement of $T_{c}$ may be due to suppression
of superconducting fluctuations.

Figure 2:  Phase diagram of triplet superconductivity in graphene in the
$T-\Gamma_{\perp}$ plane. For temperatures below the mean-field BCS transition
$T_{\rm BCS}$ (dashed line) there are strong superconducting fluctuations
which are especially pronounced in the isotropic triplet system. Ising spin-
orbit coupling $\Gamma_{\perp}$ pins the $\mathbf{d}$ vector and induces a gap
for the unitary Goldstone modes, suppressing fluctuations and enhancing the
renormalized transition into a unitary triplet superconducting state with
$\langle d_{z}\rangle$ pairing. This leads to a quasi-long-range ordered
unitary triplet superconducting phase (blue region) below the renormalized
critical temperature (dashed line C-D) via BKT transition. At small spin-orbit
coupling the Goldstone modes strongly suppress the triplet transition, but
allow for a quasi-long-range ordered singlet $4e$ superconductor (dashed line
A-C), also via a BKT transition. In principle, at lower temperatures in the
$4e$ phase there may be a second Ising transition which spontaneously breaks
the residual $\mathbb{Z}_{2}$ spin-symmetry of $d_{z}\to-d_{z}$ (solid line
B-C), terminating in a tricritical point. In the presence of a substrate with
large spin-orbit coupling, such as the WSe2, we reside at the star shown,
where the transition is essentially restored to $T_{\rm BCS}$. The size of the
QLRO $4e$ phase is over-exaggerated for clarity in this diagram.

Magnetic Field—In addition to proximity inducing spin-orbit coupling, we
briefly consider breaking the spin-rotation symmetry by applying an in-plane
magnetic field, as was done in the Santa Barbara experiments which first saw
putative triplet pairing in BBG [4]. The relevant term in the Ginzburg-Landau
expansion is not the $\Gamma|\mathbf{d}_{\perp}|^{2}$ term but rather the term

$f_{\rm Zeeman}=-\mu_{\rm
eff}\mathbf{H}_{\parallel}\cdot(i\mathbf{d}\times\mathbf{d}^{\star}),$ (5)

where $\mathbf{H}_{\parallel}=H_{\rm ext}\mathbf{e}_{x}$ is the applied in-
plane magnetic field that couples to the non-unitary condensate moment, which
has an effective magnetic moment of $\mu_{\rm eff}$.

While we find that within the BCS model and quasiclassical approximation this
coupling vanishes, we expect that in a strongly-coupled superconductor this
coupling will in general appear with a coefficient proportional to the degree
of particle-hole symmetry breaking; in the Supplemental Material we find
$\mu_{\rm eff}\sim 2\mu_{B}N^{\prime}(E_{F})/\tilde{g}_{f}$ where
$N^{\prime}(E_{F})$ is a measure of the particle-hole asymmetry and
$\tilde{g}_{f}=N(E_{F})g_{f}$ is the dimensionless BCS coupling constant in
the relevant channel. We estimate that for a $100$ mT field, $\mu_{\rm
eff}H_{\rm ext}\sim\textrm{Gi}_{(2)}\times\frac{2\mu_{B}H_{\rm ext}}{T_{\rm
BCS}}\times\frac{N(E_{F})}{\tilde{g}_{f}}\sim 0.2\,N(E_{F})$, assuming a BCS
coupling constant of $\tilde{g}_{f}\sim 0.1$ 555This is to be distinguished
from the ratio of $T_{\rm BCS}/E_{F}$, though in general it is a conceptual
question how to characterize the strength of the pairing interaction in the
regime of large Ginzburg number..

At the level of mean-field theory, this term has two effects. The first is
that it actually slightly raises $T_{c}$, and the second is that it causes the
transition to now enter into a fully spin-polarized condensate, analogous to
the A1 phase of 3He [27]. To see this, we can change basis to the Zeeman
sublevels $m=1,0,-1$ of the pair spin along the direction of the in-plane
field, with condensates $d_{m}$. The free energy up to quadratic order is then
characterized by the potential

$V(d_{+1})=(r+H_{\parallel})|d_{-1}|^{2}+r|d_{0}|^{2}+(r-H_{\parallel})|d_{+1}|^{2}.$
(6)

Evidently, the system first orders into the $m=1$ condensate once
$r-H_{\parallel}$ becomes zero. This gain in energy is linear in the
condensate density, whereas the energy penalty for non-unitary pairing due to
the $v$ coefficient is only quadratic and thus near $T_{c}$, when the
condensate density is small, it will always enter into a fully-polarized
condensate.

In the polarized phase the magnetic field is expected to also suppress
fluctuations since it breaks spin-rotation symmetry. One can check that there
is only one remaining Goldstone mode, which is the overall $U(1)$ gauge
symmetry; the remaining modes are pinned by the field with a gap of order
$|\mu_{\rm eff}H_{\rm ext}|$. While a quantitative treatment of this effect is
left to future studies, one would then expect a similar behavior as to the
case of the Ising spin-orbit coupling, with the replacement of
$\Gamma_{\perp}\to\mu_{\rm eff}H_{\rm ext}$. Comparing the pinning energy of
the lowest-lying Goldstone mode between the case of magnetic field and spin-
orbit coupling, we see that the spin-orbit interaction, with
$\Gamma_{\perp}\sim 30\,N(E_{F})$, is roughly $150$ times stronger than the
pinning energy of the 100 mT field, which gives $gH_{\rm ext}\sim
0.2\,N(E_{F})$. It is therefore plausible that there will still be a
significant decrease in the fluctuations when passing from a 100 mT magnetic
field to the 0.7 meV spin-orbit interaction.

4e-Condensate—We now return to the existence of $4e$-pairing. In the ansatz
used to solve the large-$N$ theory, this corresponds to a complex order
parameter $\chi$, which couples to $\mathbf{d}^{2}$. This parameter allows for
a breaking of the $U(1)$ gauge symmetry while still leaving the $SO(3)$ spin-
rotation symmetry unbroken, with a condensate of
$\langle\mathbf{d}^{2}\rangle\neq 0$. Therefore, this $4e$ condensate is more
robust against fluctuations and can in principle appear at small spin-orbit
coupling $\Gamma_{\perp}$.

This is seen in the phase-diagram (greatly exaggerated relative to the mean-
field condensate) in Fig. 2. The boundary of this phase is also determined in
the Supplemental Material by minimizing the variational ansatz. We find the
$4e$ condensate appears at a temperature $T^{\star}$ found by solving

$\log\left(\frac{T^{\star}}{T_{\rm
BCS}}\right)=\frac{1}{2}\textrm{Gi}_{(2)}-\tilde{\Gamma}_{\perp}-\textrm{Gi}_{(2)}\log\left(\frac{2C}{\textrm{Gi}_{(2)}}\right).$
(7)

Thus, the phase boundary decreases exponentially with increasing
$\Gamma_{\perp}$ before eventually colliding with the phase boundary of $2e$
condensate at a multi-critical point, with critical spin-orbit coupling
$\Gamma_{\perp c}=N(E_{F})\textrm{Gi}_{(2)}$.

This phase can be thought of as a fluctuating superfluid of triplet pairs
which exhibit correlations in the phase of their fluctuations while remaining
randomly oriented in their spin orientation. $4e$ superconductivity has been
proposed recentlyin graphene based systems [28, 29, 30] as well as quantum
gases [15, 17], and would be accompanied by a characteristic flux quantum of
$h/4e$ —half the usual value.

Further lowering the temperature, we expect that the $4e$ phase will exhibit a
second phase transition into the unitary (spin nematic) phase characterized by
a condensate in $\langle d_{z}\rangle$, shown in Fig. 2 as the blue region at
lower temperatures below the red region. Because the $4e$ condensate has half
the standard flux quantum, it only fixes the phase of the condensate modulo
$\pi$ —alternatively, this is understood as a Josephson coupling between the
$4e$ condensate $\chi$ and the $2e$ condensate of the form $f_{\rm
int}\sim\overline{\chi}\psi^{2}+\textrm{c.c.}$ and thus the phase of $\psi$
still has a remnant $\mathbb{Z}_{2}$ symmetry which is broken at the second
transition. It is worth mentioning that this phase is in reality quite subtle
[17] and may be overwhelmed by competing orders such as the intervalley
coherent order [31]

Experimental Signatures It is likely possible to tune the strength of the
induced Ising spin-orbit coupling, and hence the pinning energy
$\Gamma_{\perp}$, by adding spacer layers of an inert material like hexagonal
boron-nitride or even an oxide barrier in order to determine the functional
dependence of $T_{c}$ on $\Gamma_{\perp}$. This may also be possible to vary
by tuning the twist-angle between the BBG and WSe2 [9]. Our predictions also
should carry over to other semiconductors besides WSe2 so long as they feature
large Ising spin-orbit coupling and interface well with graphene.

As to signatures of the fluctuations, a good starting point may be
measurements of the Nernst-effect, which is sensitive to mobile vortices
forming above $T_{c}$ [32, 33, 34, 35]. It may also be possible to image the
spin-fluctuations directly using scanning SQUID [36] or NV magnetometry [37,
38], though this is complicated somewhat by the unitary nature of the
fluctuations, which have no net magnetic dipole moment.

In general, we expect that the presence of strong triplet fluctuations will
lead to a strongly non-BCS character of the spin-susceptibility in the normal
state just above $T_{c}$. This is explored in more detail in the Supplemental
Material, but neglecting the $4e$ pairing, we expect that the zero-field in-
plane magnetic susceptibility (which probes only the spin susceptbility) has a
characteristic behavior of $\chi_{\parallel}\sim\textrm{Gi}_{(2)}/(r_{\rm
eff}+\Gamma_{\perp})$ and for small spin-orbit coupling this has an apparent
divergence as the critical point $r_{\rm eff}=0$ is approached, though this
may be complicated by the emergence of $4e$-singlet or other competing orders
at small spin-orbit coupling. We note this is in contrast to the standard case
of a BCS triplet superconductor, which has no anomalous behavior at the
critical point in the transverse spin-susceptibility.

Finally, spectroscopic measurements would be useful in establishing the
existence of a pseudogap phase above the actual transition, which in turn
would indicate the presence of strong-fluctuations. This has already been
accomplished in the case of twisted bilayer and trilayer graphene, which have
indicated at least strong-coupling physics [20, 21]. Potentially, there may
also be spectroscopic signatures of the collective modes themselves which
would be interesting to probe [39, 40, 41, 42]. In particular, signatures of
dipolar condensate waves [39], or degenerate Bardasis-Schrieffer modes [40]
should point towards triplet pairing in this system.

Conclusion—In conclusion, we have presented evidence which supports the claim
that graphene multilayers, and in particular BBG, may support triplet
superconductivity. In its intrinsic form, graphene has almost perfect spin-
rotation symmetry and this, combined with the two-dimensional nature of
graphene, leads to remarkably strong fluctuation effects which profoundly
affect the phase diagram. Applying either Ising spin-orbit coupling or an in-
plane Zeeman field would suppress these fluctuations and lead to an
enhancement of the superconducting transition temperature, as have been
experimentally observed.

While our work is to some degree universally applicable to triplet pairing in
graphene, there remain a number of qualitative complications due to the
particular broken symmetry states which are proximal to the triplet pairing.
It therefore remains important to understand what role the fermiology plays
and how this couples to the fluctuation physics we have outlined here. In
particular, it is interesting to consider the interplay between the
fluctuating triplet order and possible charge-density wave orders, such as the
intervalley coherent order [31, 4, 5]. Treating the joint fluctuations of
these two order parameters is a very interesting problem which may help
explain certain aspects of the system at low magnetic field. One can imagine
that by expanding the symmetry group from $SO(3)\times U(1)$ to include other
symmetries, one can obtain a unified [43, 44, 45] model of both the isospin
ordered and triplet paired states, like in Refs. [46, 47]. Similarly, recent
proposals for a spontaneous Ising spin-orbit coupling [48] may also naturally
fall out of a higher-symmetry parent phase which spontaneously breaks $SO(3)$
before $U(1)$ (and thus should also have soft Goldstone modes associated to
the precession of the spontaneous Ising axis).

In the future, more direct evidence of triplet pairing, and of fluctuating
superconductivity with a pseudogap, is still needed. Theoretically, it will be
important in the future to account for the role of broken symmetry normal
states and their potentially competing order parameters in our model.
Similarly, more careful consideration of the role of vortices and BKT physics
is also needed, especially to make contact with experiments. We believe that
generally, the study of fluctuation effects in low-dimensional, clean, highly-
symmetric graphene systems provides a unique opportunity to study
superconductivity beyond mean-field, and in the presence of strong
correlations.

###### Acknowledgements.

The authors would like to acknowledge invaluable discussions with Étienne
Lantagne-Hurtubise, Andrea Young, Ludwig Holleis, Pavel Dolgirev, Ilya
Esterlis, and Andrew Zimmerman, and Olivia Liebman. J.B.C., N.R.P., Y.X.,
A.Y., and P.N. are supported by the Quantum Science Center (QSC), a National
Quantum Information Science Research Center of the U.S. Department of Energy
(DOE). J.B.C. and Y.X. are HQI Prize Postdoctoral Fellows and gratefully
acknowledge support from the Harvard Quantum Initiative. N.R.P. is supported
by the Army Research Office through an NDSEG fellowship. E.D. acknowledges
support from the ETH Zürich Institute for Theoretical Physics. P.N. is a Moore
Inventor Fellow and gratefully acknowledges support through Grant GBMF8048
from the Gordon and Betty Moore Foundation.

## References

  * Zhou  _et al._ [2021] H. Zhou, T. Xie, T. Taniguchi, K. Watanabe, and A. F. Young, Superconductivity in rhombohedral trilayer graphene, Nature 598, 434 (2021).
  * Cao  _et al._ [2021] Y. Cao, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Pauli-limit violation and re-entrant superconductivity in moiré graphene, Nature 595, 526 (2021).
  * Park  _et al._ [2022] J. M. Park, Y. Cao, L.-Q. Xia, S. Sun, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Robust superconductivity in magic-angle multilayer graphene family, Nature Materials 10.1038/s41563-022-01287-1 (2022).
  * Zhou  _et al._ [2022] H. Zhou, L. Holleis, Y. Saito, L. Cohen, W. Huynh, C. L. Patterson, F. Yang, T. Taniguchi, K. Watanabe, and A. F. Young, Isospin magnetism and spin-polarized superconductivity in Bernal bilayer graphene, Science 375, 774 (2022), 2110.11317 .
  * Zhang  _et al._ [2022] Y. Zhang, R. Polski, A. Thomson, É. Lantagne-Hurtubise, C. Lewandowski, H. Zhou, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, Spin-Orbit Enhanced Superconductivity in Bernal Bilayer Graphene, arXiv (2022), 2205.05087 .
  * Qin  _et al._ [2022] W. Qin, C. Huang, T. Wolf, N. Wei, I. Blinov, and A. H. MacDonald, Functional Renormalization Group Study of Superconductivity in Rhombohedral Trilayer Graphene (2022), arXiv:2203.09083 .
  * You and Vishwanath [2022] Y. Z. You and A. Vishwanath, Kohn-Luttinger superconductivity and intervalley coherence in rhombohedral trilayer graphene, Physical Review B 105, 1 (2022), arXiv:2109.04669 .
  * Lu  _et al._ [2022] D.-C. Lu, T. Wang, S. Chatterjee, and Y.-Z. You, Correlated metals and unconventional superconductivity in rhombohedral trilayer graphene: a renormalization group analysis (2022), arXiv:2206.01213 .
  * Chou  _et al._ [2022a] Y.-Z. Chou, F. Wu, and S. D. Sarma, Enhanced superconductivity through virtual tunneling in Bernal bilayer graphene coupled to WSe$_2$, arXiv (2022a), 2206.09922 .
  * Cornfeld  _et al._ [2021] E. Cornfeld, M. S. Rudner, and E. Berg, Spin-polarized superconductivity: Order parameter topology, current dissipation, and multiple-period Josephson effect, Physical Review Research 3, 013051 (2021), 2006.10073 .
  * Note [1] The pairing is assumed isotropic in $\mathbf{k}$ with respect to the $\mathbf{K,K^{\prime}}$ points, but is valley singlet meaning that under $C_{3}$ rotations the valleys exchange and the wavefunction is overall odd, making it $f$-wave.
  * Chou  _et al._ [2021] Y.-Z. Chou, F. Wu, J. D. Sau, and S. Das Sarma, Correlation-induced triplet pairing superconductivity in graphene-based moiré systems, Phys. Rev. Lett. 127, 217001 (2021).
  * Chou  _et al._ [2022b] Y.-Z. Chou, F. Wu, J. D. Sau, and S. Das Sarma, Acoustic-phonon-mediated superconductivity in moiréless graphene multilayers, Phys. Rev. B 106, 024507 (2022b).
  * Note [2] In fact, it is likely that BBG is not in the weak coupling regime and this may lead to additional interesting effects, though the argument we advance is based on symmetry and is still a relevant consideration in such case.
  * Mukerjee  _et al._ [2006] S. Mukerjee, C. Xu, and J. E. Moore, Topological Defects and the Superfluid Transition of the s=1 Spinor Condensate in Two Dimensions, Physical Review Letters 97, 120406 (2006), cond-mat/0605102 .
  * Podolsky  _et al._ [2009] D. Podolsky, S. Chandrasekharan, and A. Vishwanath, Phase transitions of S=1 spinor condensates in an optical lattice, Physical Review B 80, 214513 (2009), 0707.0695 .
  * Kobayashi [2019] M. Kobayashi, Berezinskii–Kosterlitz–Thouless Transition of Spin-1 Spinor Bose Gases in the Presence of the Quadratic Zeeman Effect, Journal of the Physical Society of Japan 88, 094001 (2019), 1902.00210 .
  * Larkin and Varlamov [2005] A. Larkin and A. Varlamov,  _Theory of Fluctuations in Superconductors_ (Oxford University Press, 2005).
  * Varlamov  _et al._ [1999] A. A. Varlamov, G. Balestrino, E. Milani, and D. V. Livanov, The role of density of states fluctuations in the normal state properties of high T c superconductors, Advances in Physics 48, 655 (1999).
  * Park  _et al._ [2021] J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Tunable strongly coupled superconductivity in magic-angle twisted trilayer graphene, Nature 590, 249 (2021), 2012.01434 .
  * Oh  _et al._ [2021] M. Oh, K. P. Nuckolls, D. Wong, R. L. Lee, X. Liu, K. Watanabe, T. Taniguchi, and A. Yazdani, Evidence for unconventional superconductivity in twisted bilayer graphene, Nature 600, 240 (2021), 2109.13944 .
  * Kim  _et al._ [2022] H. Kim, Y. Choi, C. Lewandowski, A. Thomson, Y. Zhang, R. Polski, K. Watanabe, T. Tanigushi, J. Alicea, and S. Nadj-Perge, Evidence for unconventional super- conductivity in twisted trilayer graphene, Nature 606, 494 (2022).
  * Note [3] While this may be somewhat optimistic, closer examination of the data in Ref. [5], and in particular the inset of Fig. 2(d), reveals a possible Azlamazov-Larkin downturn in $R_{xx}$ around $T\sim 500$ mK, indicating this may be a serious possibility.
  * Sichau  _et al._ [2019] J. Sichau, M. Prada, T. Anlauf, T. J. Lyon, B. Bosnjak, L. Tiemann, and R. H. Blick, Resonance Microwave Measurements of an Intrinsic Spin-Orbit Coupling Gap in Graphene: A Possible Indication of a Topological State, Physical Review Letters 122, 046403 (2019), 1709.05705 .
  * Note [4] This is a lower bound on the increase. In particular, it is thought that the intrinsic spin-orbit coupling in graphene takes the form of $H_{\rm soc}=\lambda\sigma_{3}\rho_{3}s_{3}$ where $s_{a}$ are the sublattice Pauli matrices. Thus, in AB stacked graphene the average spin-orbit coupling still vanishes and the intrinsic value is likely much less than 40 $\mathrm{\SIUnitSymbolMicro eV}$ to begin with.
  * Note [5] This is to be distinguished from the ratio of $T_{\rm BCS}/E_{F}$, though in general it is a conceptual question how to characterize the strength of the pairing interaction in the regime of large Ginzburg number.
  * Kojima and Ishimoto [2008] H. Kojima and H. Ishimoto, Spin Polarized Superfluid 3He A1, Journal of the Physical Society of Japan 77, 111001 (2008).
  * Maccari  _et al._ [2022] I. Maccari, J. Carlström, and E. Babaev, Possible time-reversal-symmetry-breaking fermionic quadrupling condensate in twisted bilayer graphene, arXiv (2022), 2206.02698 .
  * Khalaf  _et al._ [2022] E. Khalaf, P. Ledwith, and A. Vishwanath, Symmetry constraints on superconductivity in twisted bilayer graphene: Fractional vortices, 4e condensates, or nonunitary pairing, Physical Review B 105, 224508 (2022).
  * Fernandes and Fu [2021] R. M. Fernandes and L. Fu, Charge-4e Superconductivity from Multicomponent Nematic Pairing: Application to Twisted Bilayer Graphene, Physical Review Letters 127, 047001 (2021), 2101.07943 .
  * Chatterjee  _et al._ [2021] S. Chatterjee, T. Wang, E. Berg, and M. P. Zaletel, Inter-valley coherent order and isospin fluctuation mediated superconductivity in rhombohedral trilayer graphene, arXiv (2021), 2109.00002 .
  * Ussishkin  _et al._ [2002] I. Ussishkin, S. L. Sondhi, and D. A. Huse, Gaussian Superconducting Fluctuations, Thermal Transport, and the Nernst Effect, Physical Review Letters 89, 287001 (2002), cond-mat/0204484 .
  * Cyr-Choinière  _et al._ [2018] O. Cyr-Choinière, R. Daou, F. Laliberté, C. Collignon, S. Badoux, D. LeBoeuf, J. Chang, B. J. Ramshaw, D. A. Bonn, W. N. Hardy, R. Liang, J.-Q. Yan, J.-G. Cheng, J.-S. Zhou, J. B. Goodenough, S. Pyon, T. Takayama, H. Takagi, N. Doiron-Leyraud, and L. Taillefer, Pseudogap temperature T* of cuprate superconductors from the Nernst effect, Physical Review B 97, 064502 (2018), 1703.06927 .
  * Levchenko  _et al._ [2011] A. Levchenko, M. R. Norman, and A. A. Varlamov, Nernst effect from fluctuating pairs in the pseudogap phase of the cuprates, Physical Review B 83, 020506 (2011), 1009.2213 .
  * Jotzu  _et al._ [2021] G. Jotzu, G. Meier, A. Cantaluppi, A. Cavalleri, D. Pontiroli, M. Riccò, A. Ardavan, and M.-S. Nam, Superconducting fluctuations observed far above Tc in the isotropic superconductor K3C60, arXiv (2021), 2109.08679 .
  * Vasyukov  _et al._ [2013] D. Vasyukov, Y. Anahory, L. Embon, D. Halbertal, J. Cuppens, L. Neeman, A. Finkler, Y. Segev, Y. Myasoedov, M. L. Rappaport, M. E. Huber, and E. Zeldov, A scanning superconducting quantum interference device with single electron spin sensitivity, Nature Nanotechnology 8, 639 (2013).
  * Dolgirev  _et al._ [2022] P. E. Dolgirev, S. Chatterjee, I. Esterlis, A. A. Zibrov, M. D. Lukin, N. Y. Yao, and E. Demler, Characterizing two-dimensional superconductivity via nanoscale noise magnetometry with single-spin qubits, Physical Review B 105, 024507 (2022), 2106.05283 .
  * Chatterjee  _et al._ [2022] S. Chatterjee, P. E. Dolgirev, I. Esterlis, A. A. Zibrov, M. D. Lukin, N. Y. Yao, and E. Demler, Single-spin qubit magnetic spectroscopy of two-dimensional superconductivity, Physical Review Research 4, L012001 (2022), 2106.03859 .
  * Poniatowski  _et al._ [2021] N. R. Poniatowski, J. B. Curtis, C. G. L. Bøttcher, V. M. Galitski, A. Yacoby, P. Narang, and E. Demler, Surface Cooper pair spin waves in triplet superconductors, arXiv (2021), 2112.12146 .
  * Curtis  _et al._ [2022] J. B. Curtis, N. R. Poniatowski, A. Yacoby, and P. Narang, Proximity-induced collective modes in an unconventional superconductor heterostructure, arXiv (2022), 2201.04635 .
  * Poniatowski  _et al._ [2022] N. R. Poniatowski, J. B. Curtis, A. Yacoby, and P. Narang, Spectroscopic signatures of time-reversal symmetry breaking superconductivity, Communications Physics 5, 44 (2022), 2103.05641 .
  * Dolgirev  _et al._ [2021] P. E. Dolgirev, A. Zong, M. H. Michael, J. B. Curtis, D. Podolsky, A. Cavalleri, and E. Demler, Periodic dynamics in superconductors induced by an impulsive optical quench, arXiv (2021), 2104.07181 .
  * Zhang [1997] S.-C. Zhang, A Unified Theory Based on $SO(5)$ Symmetry of Superconductivity and Antiferromagnetism, Science 275, 1079 (1997).
  * Podolsky  _et al._ [2004] D. Podolsky, E. Altman, T. Rostunov, and E. Demler, SO(4) Theory of Antiferromagnetism and Superconductivity in Bechgaard Salts, Phys. Rev. Lett. 93, 246402 (2004).
  * Murakami  _et al._ [1999] S. Murakami, N. Nagaosa, and M. Sigrist, SO(5) Model of $p$-Wave Superconductivity and Ferromagnetism, Phys. Rev. Lett. 82, 2939 (1999).
  * Khalaf  _et al._ [2021] E. Khalaf, S. Chatterjee, N. Bultinck, M. P. Zaletel, and A. Vishwanath, Charged skyrmions and topological origin of superconductivity in magic-angle graphene, Science Advances 7, eabf5299 (2021), 2004.00638 .
  * Christos  _et al._ [2020] M. Christos, S. Sachdev, and M. S. Scheurer, Superconductivity, correlated insulators, and Wess–Zumino–Witten terms in twisted bilayer graphene, Proceedings of the National Academy of Sciences 117, 29543 (2020), 2007.00007 .
  * Lake  _et al._ [2022] E. Lake, A. S. Patri, and T. Senthil, Pairing symmetry of twisted bilayer graphene: a phenomenological synthesis, arXiv (2022), 2204.12579 .
  * de la Barrera  _et al._ [2022] S. C. de la Barrera, S. Aronson, Z. Zheng, K. Watanabe, T. Taniguchi, Q. Ma, P. Jarillo-Herrero, and R. Ashoori, Cascade of isospin phase transitions in Bernal-stacked bilayer graphene at zero magnetic field, Nature Physics 18, 771 (2022).

*

##  Appendix A Microscopic BCS Model

In this section, we derive the Ginzburg-Landau free energy functional from the
microscopic weak coupling theory. As introduced in the main text, we consider
a two-dimensional system subject to an Ising spin-orbit coupling of strength
$\lambda$, and an attractive interaction in the spin-triplet valley-singlet
channel, with the coupling constant $g_{f}$. We treat the system within the
functional integral formalism in the the Matsubara representation, and
introduce the four-momenta $q=(i\omega_{m},\mathbf{q})$ and
$k=(i\epsilon_{n},\mathbf{k})$ where $i\omega_{m}$ ($i\epsilon_{n}$) is a
bosonic (fermionic) Matsubara frequency. The action governing the system is
then written as

$S=\sum_{k}\psi_{k}^{\dagger}\left(-i\epsilon_{n}+\xi_{\mathbf{k}}+\frac{\lambda}{2}\,\sigma_{3}\rho_{3}\right)\psi_{k}-\frac{g_{f}}{2}\sum_{q}\mathbf{P}_{q}^{\dagger}\mathbf{P}_{q}$
(8)

where $\psi_{k}$ is a four-component spinor in spin $\otimes$ valley space and
$\sigma_{i}$ ($\rho_{i}$) are the Pauli matrices acting on the spin (valley)
degrees of freedom. We take the single-particle dispersion $\xi_{\mathbf{k}}$
to be parabolic in each of the valleys, and have defined the pair bilinear

$\mathbf{P}_{q}^{\dagger}=\sum_{k}\bar{\psi}^{T}_{k+q/2}\left(i\bm{\sigma}\sigma_{2}\otimes
i\rho_{2}\right)\bar{\psi}_{-k+q/2}\,.$ (9)

We introduce the vector superconducting order parameter $\mathbf{d}_{q}$ via a
Hubbard-Stratonovich transformation. We express the action by doubling our
basis and defining the eight-component Nambu spinor in spin $\otimes$ valley
$\otimes$ particle-hole space,
$\Psi_{k}=\left(\psi_{k}\,,\bar{\psi}_{-k}i\sigma_{2}\ \right)^{T}$, so that
the action takes the form

$S=\frac{1}{2g_{f}}\sum_{q}\bar{\mathbf{d}}_{q}\cdot\mathbf{d}_{q}-\frac{1}{2}\sum_{kq}\Psi_{k+q/2}^{\dagger}\left[\check{\mathbb{G}}^{-1}_{k}-\check{\Delta}_{q}\right]\Psi_{k-q/2}$
(10)

where the normal-state Green function and pairing vertex are

$\displaystyle\check{\mathbb{G}}^{-1}_{k}$
$\displaystyle=i\epsilon_{n}-\xi_{\mathbf{k}}\tau_{3}-\frac{\lambda}{2}\sigma_{3}\rho_{3}\,,$
(11) $\displaystyle\check{\Delta}_{q}$
$\displaystyle=\mathbf{d}_{q}\cdot\bm{\sigma}(i\rho_{2})\,\tau^{+}+\bar{\mathbf{d}}_{-q}\cdot\bm{\sigma}(-i\rho_{2})\,\tau^{-}\,.$
(12)

In the above, we have introduced the Pauli matrices $\tau_{i}$ in Nambu space,
as well as the combinations $\tau^{\pm}=\frac{1}{2}\left(\tau_{1}\pm
i\tau_{2}\right)$. We may now perform the Gaussian integral over the fermion
fields to furnish an effective action for the order parameter field,

$S=\frac{1}{2g_{f}}\sum_{q}\bar{\mathbf{d}}_{q}\cdot\mathbf{d}_{q}-\frac{1}{2}\mathbf{tr}\log\left[\check{\mathbb{G}}^{-1}_{k}-\check{\Delta}_{q}\right]\,.$
(13)

In the vicinity of the superconducting transition where the order parameter is
small, we can expand the action in powers of $\mathbf{d}$ which, in the static
$i\omega_{n}=0$ limit, amounts to a derivation of the Ginzburg-Landau
functional studied in the main text.

To begin, we expand the functional logarithm to quadratic order, yielding

$S^{(2)}=\frac{1}{2g_{f}}\sum_{q}\bar{\mathbf{d}}_{q}\cdot\mathbf{d}_{q}+\frac{T}{4}\sum_{kq}\mathbf{tr}\left[\check{\mathbb{G}}_{k+q/2}\check{\Delta}_{q}\check{\mathbb{G}}_{k-q/2}\check{\Delta}_{-q}\right]\,.$
(14)

The Ising spin-orbit coupling breaks spin rotation symmetry and leads to
different propagators for the $d_{z}$ and $d_{x},d_{y}$ fluctuations. The
quadratic action can then be expressed as

$S^{(2)}=\frac{1}{2}\sum_{q}\left[\mathcal{L}^{-1}_{zz}(0,\mathbf{q})\bar{d}^{z}_{q}d^{z}_{q}+\mathcal{L}^{-1}_{\perp}(0,\mathbf{q})\left(\bar{d}^{x}_{q}d_{q}^{x}+\bar{d}^{y}_{q}d_{q}^{y}\right)\right]\,,$
(15)

The fluctuation propogators can be further expanded for small $\mathbf{q}$ as

$\displaystyle\mathcal{L}_{zz}^{-1}(0,\mathbf{q})$
$\displaystyle=r+K_{zz}\mathbf{q}^{2}+\dots$ (16)
$\displaystyle\mathcal{L}_{\perp}^{-1}(0,\mathbf{q})$
$\displaystyle=r+\Gamma+K_{\perp}\mathbf{q}^{2}+\dots\,.$ (17)

which allows us to compute the quadratic coefficients appearing in the
Ginzburg-Landau functional studied in the main text. In particular, we find
that the $r$ coefficient reflects the singular piece of the fluctuation
propogator that determines the mean field transition temperature
$T_{BCS}=(4\mathrm{e}^{\gamma}/\pi)\,\mathrm{e}^{-1/N(E_{F})g_{f}}$,

$r=\mathcal{L}^{-1}_{zz}(0,0)=\frac{1}{g_{f}}-2T\sum_{k}\left[\frac{1}{\epsilon_{n}^{2}+(\xi_{\mathbf{k}}+\lambda/2)^{2}}+\frac{1}{\epsilon_{n}^{2}+(\xi_{\mathbf{k}}-\lambda/2)^{2}}\right]=\nu\log\left(\frac{T}{T_{BCS}}\right)\,.$
(18)

Meanwhile, the $\mathbf{d}$-vector pinning is

$\Gamma=\mathcal{L}^{-1}_{\perp}(0,0)-\mathcal{L}_{zz}^{-1}(0,0)=2T\sum_{k}\frac{\lambda^{2}}{[\epsilon_{n}^{2}+(\xi_{\mathbf{k}}+\lambda/2)^{2}][\epsilon_{n}^{2}+(\xi_{\mathbf{k}}-\lambda/2)^{2}]}\,.$
(19)

This quantity is plotted as function of the spin orbit strength $\lambda$ in
Fig. 1 of the main text, evaluated in the approximation that the density of
states is constant over the energy range relevant for superconductivity, i.e.
$N(\xi_{\mathbf{k}})\approx N(E_{F})$.

For simplicity, we evaluate the coefficients of the gradient terms
$K_{zz}=K_{\perp}\equiv K$ to lowest order in $\lambda$. We also expand the
action to fourth order for $\mathbf{q}=0$, and again evaluate the resulting
expansion coefficients to leading order in $\lambda$. The results for these
coefficients are presented in Table 1.

GL Parameter | Weak Coupling Value
---|---
$r$ | $N(E_{F})\log(T/T_{\rm BCS})$
$u$ | $7\zeta(3)N(E_{F})/(8\pi^{2}T_{\rm BCS}^{2})$
$v$ | $7\zeta(3)N(E_{F})/(8\pi^{2}T_{\rm BCS}^{2})$
$\Gamma_{\perp}$ | $7\zeta(3)N(E_{F})\lambda^{2}/(4\pi^{2}T_{\rm BCS}^{2})$
$K$ | $7\zeta(3)N(E_{F})v_{F}^{2}/(16\pi^{2}T_{\rm BCS}^{2})$
Table 1: Coefficients appearing in Ginzburg-Landau functional and their weak
coupling values, evaluated to leading order in the induced spin orbit
coupling, $\lambda$.

##  Appendix B Zeeman field

In this section, we consider the Ginzburg-Landau functional for a system
subject to an in-plane magnetic field, $H_{\parallel}$, in the absence of
spin-orbit coupling. The technical details of the calculation are identical to
the previous section, with the replacement of the normal-state Green function
by

$\check{\mathbb{G}}^{-1}_{k}=i\epsilon_{n}-\xi_{\mathbf{k}}\tau_{3}-H_{\parallel}\sigma_{3}\tau_{3}\,.$
(20)

By computing the quadratic action, a new term appears which corresponds to a
linear coupling at $\mathbf{q}=0$ between the magnetic field and the
condensate polarization,

$S_{\text{lin}}=-\mu_{\text{eff}}\sum_{q}H_{\parallel}\left(i\mathbf{d}_{q}\times\bar{\mathbf{d}}_{q}\right)_{z}\,,$
(21)

where the coupling constant is

$\mu_{\text{eff}}=8T\sum_{k}\frac{\xi_{\bf
k}}{[(i\epsilon_{n}-\xi_{\mathbf{k}})^{2}-H_{\parallel}^{2}][(i\epsilon_{n}+\xi_{\mathbf{k}})^{2}-H_{\parallel}^{2}]}\,.$
(22)

To capture the linear response to weak fields, we can evaluate this coupling
constant for $H_{\parallel}=0$, and express the sum over momenta as an
integral over energy weighted by the density of states $N(\xi)$, so that

$\mu_{\text{eff}}=8T\sum_{\epsilon_{n}}\int\mathrm{d}\xi\;N(\xi)\,\frac{\xi}{(\epsilon_{n}^{2}+\xi^{2})^{2}}\,.$
(23)

This expression is odd under $\xi\to-\xi$, and thus vanishes in the usual
quasiclassical approximation $N(\xi)\approx N(E_{F})$ which treats the
approximate particle-hole symmetry of the weak coupling theory as exact. Thus,
we see that the coupling between the condensate polarization and field is
proportional to the degree of particle-hole asymmetry. Taking this to be
small, we may then expand the density of states around its value at the Fermi
level and write $N(\xi)\approx N(E_{F})+N^{\prime}(E_{F})\,\xi$, which
generates the leading order contribution to $\mu_{\text{eff}}$,

$\mu_{\text{eff}}=4\pi
TN^{\prime}(E_{F})\sum_{\epsilon_{n}}\frac{1}{\epsilon_{n}}\approx
4N^{\prime}(E_{F})\left[\frac{1}{N(E_{F})g_{f}}+\log\left(\frac{T_{BCS}}{T}\right)\right]\,.$
(24)

##  Appendix C Large-$N$ Expansion

Here we approach the problem of order-parameter fluctuations via a $1/N$
expansion in the number of spin components, which in reality is $N=3$. But as
we all know, $3=\infty$. Of course, in reality this means that the effect of
fluctuations is likely underestimated in this approach, and furthermore only
Gaussian fluctuations are accounted for so that topological textures like
vortices and half-vortices are also unaccounted for.

This approach will only be applied to the case of the Ising spin-orbit
coupling, since in this case it is clear how to generalize to the
$N$-component spinor, but in the case of magnetic field it is less clear. The
appropriate generalization to the $N$-component unitary triplet pairing is
found by writing
$(i\mathbf{d}^{*}\times\mathbf{d})^{2}=|\mathbf{d}|^{4}-|\mathbf{d}\cdot\mathbf{d}|^{2}$,
which can be generalized beyond the vector representation.

We then have

$f=K|\nabla\mathbf{d}|^{2}+r|\mathbf{d}|^{2}+\Gamma_{\perp}|\mathbf{d}_{\perp}|^{2}+\frac{(u+v)}{2N}|\mathbf{d}|^{4}-\frac{v}{2N}|\mathbf{d}\cdot\mathbf{d}|^{2},$
(25)

where the separation into $\perp$ and $\parallel$ spin components is selected
by the Ising spin-orbit $\Gamma_{\perp}$.

We now solve this through a Gaussian variational ansatz for the free energy,
evaluated in the large $N$ limit. The ansatz we use is based on the Feynman-
Bogoliubov-Gibbs inequality such that

$\mathcal{F}\leq\mathcal{F}_{\rm
eff}=\mathcal{F}_{0}+\langle\Delta\mathcal{F}\rangle_{\mathcal{F}_{0}}$ (26)

where $\mathcal{F}_{0}$ is the free energy of the ansatz and
$\Delta\mathcal{F}=\mathcal{F}-\mathcal{F}_{0}$. If we choose
$\mathcal{F}_{0}$ to be Gaussian then this can be evaluated easily. By
minimizing over our ansatz parameters we can best approximate the true free
energy $\mathcal{F}$.

To this end, we make the choice $\mathcal{F}_{0}=\int d^{2}rf_{0}$ with

$f_{0}=K|\nabla\mathbf{d}|^{2}+(r_{\rm
eff}+\Gamma_{\perp})|\mathbf{d}_{\perp}|^{2}+r_{\rm
eff}|d_{z}-\sqrt{N}\psi|^{2}.$ (27)

The parameters here are $r_{\rm eff}$ and $\psi$ which characterize the
renormalized critical temperature and condensate respectively.

We first evaluate the free energy of this ansatz. To leading order in $N$ this
is simply

$\mathcal{F}_{0}=NT\int_{\bf q}\log\left(K\mathbf{q}^{2}+r_{\rm
eff}+\Gamma_{\perp}\right).$ (28)

This is independent of the condensate by construction.

We also need to calculate the various expectation values using this ansatz. We
find the relevant correlation functions are

$\langle d_{z}\rangle=\sqrt{N}\psi$ (29)

and

$\langle|\mathbf{d}_{\perp}|^{2}\rangle=NT\int_{\bf
q}\frac{1}{K\mathbf{q}^{2}+r_{\rm eff}+\Gamma_{\perp}}\equiv NI$ (30)

where $I$ is the integral, which gives the normal fluid component.

We must also evaluate $\langle\Delta\mathcal{F}\rangle$ in this ansatz, which
gives to leading order in $N$

$\langle\Delta\mathcal{F}\rangle=(r-r_{\rm
eff})\langle|\mathbf{d}_{\perp}|^{2}\rangle+r\langle|d_{z}|^{2}\rangle+\frac{u+v}{2N}\langle|\mathbf{d}|^{4}\rangle-\frac{v}{2N}\langle|\mathbf{d}^{2}|^{2}\rangle.$
(31)

The first two terms give $N(r-r_{\rm eff})I+Nr|\psi|^{2}$. To leading order in
$N$ the interactions give

$\frac{1}{N}\langle|\mathbf{d}|^{4}\rangle=N\left[|\psi|^{4}+2|\psi|^{2}I+I^{2}\right]$
(32)

and

$\frac{1}{N}\langle|\mathbf{d}^{2}|^{2}\rangle=N|\psi|^{4}.$ (33)

We then obtain the effective free energy functional of

$\mathcal{F}=N\left[T\int_{\bf q}\log\left(K\mathbf{q}^{2}+r_{\rm
eff}+\Gamma_{\perp}\right)+(r-r_{\rm
eff})I+(r+(u+v)I+\frac{u}{2}|\psi|^{2})|\psi|^{2}+\frac{u+v}{2}I^{2}\right].$
(34)

We now minimize over $\psi$ and $r_{\rm eff}$ to obtain the phase diagram.
This yields two equations

$\displaystyle\left[r+(u+v)I+u|\psi|^{2}\right]\psi=0$ (35a) $\displaystyle
I+(r-r_{\rm eff})\frac{\partial I}{\partial r_{\rm eff}}-I+(u+v)\frac{\partial
I}{\partial r_{\rm eff}}|\psi|^{2}+(u+v)I\frac{\partial I}{\partial r_{\rm
eff}}=0.$ (35b)

Since $I$ is in general dependent on $r_{\rm eff}$ this last equation can be
simplified to read

$r+(u+v)(|\psi|^{2}+I)=r_{\rm eff}.$ (36)

We see the condensate density follows $(r_{\rm eff}-v|\psi|^{2})\psi=0$ and
thus the transition to non-zero $|\psi|$ occurs when $r_{\rm eff}=0$. This
gives the phase boundary of

$r+(u+v)I=0\Rightarrow r+\frac{u+v}{4\pi K}T\log(\frac{C}{\Gamma_{\perp}})=0.$
(37)

For simplicity, we neglect the temperature dependence of the fluctuations and
set $T=T_{\rm BCS}$, which simplifies the system and is formally justified
provided the renormalization is not to too low temperature. Defining the
Ginzburg-Levanyuk number as $u/(4\pi KN(E_{F}))T_{\rm
BCS}=\frac{1}{2}\textrm{Gi}_{(2)}$ and taking $C/N(E_{F})\sim(E_{F}/T_{\rm
BCS})^{2}$ we find the renormalized critical temperature of

$T_{c}=T_{\rm BCS}\left(\frac{\tilde{\Gamma}_{\perp}T_{\rm
BCS}}{E_{F}}\right)^{\textrm{Gi}_{(2)}}.$ (38)

We have defined untiless pinning energy
$\tilde{\Gamma}_{\perp}=\Gamma_{\perp}/N(E_{F})$.

This goes to zero as $\Gamma_{\perp}^{\textrm{Gi}_{(2)}}$ and for large
$\Gamma_{\perp}$ will ultimately be expected to saturate to $T_{\rm BCS}$.
Above the renormalized $T_{c}$ we have appreciable pair fluctuations but no
condensate. We now include the possibility for $4e$ pairing in order to
complete the picture of the phase diagram.

###  C.1 4$e$ Pairing

In order to accommodate $4e$ pairing, which we expect to reside at small
$\Gamma_{\perp}$, we must modify the ansatz we use to evaluate the correlation
functions. In particular, we include an off-diagonal symmetric coupling via

$f_{0}=K|\nabla\mathbf{d}|^{2}+(r_{\rm
eff}+\Gamma_{\perp})|\mathbf{d}_{\perp}|^{2}+r_{\rm
eff}|d_{z}-\sqrt{N}\psi|^{2}-\frac{1}{2}\left[\overline{\chi}(\mathbf{d}_{\perp}^{2}+(d_{z}-\sqrt{N}\psi)^{2})+\chi(\overline{\mathbf{d}}_{\perp}^{2}+(\overline{d}_{z}-\sqrt{N}\overline{\psi})^{2})\right].$
(39)

The parameters here are $r_{\rm eff}$ and $\psi$ which characterize the
renormalized critical temperature and condensate respectively. The free energy
of this ansatz is

$\mathcal{F}_{0}=N\frac{T}{2}\int_{\bf
q}\left[\log\left(K\mathbf{q}^{2}+r_{\rm
eff}+\Gamma_{\perp}+|\chi|\right)+\log\left(K\mathbf{q}^{2}+r_{\rm
eff}+\Gamma_{\perp}-|\chi|\right)\right].$ (40)

We have condensate expectation value

$\displaystyle\langle d_{z}\rangle=\sqrt{N}\psi$ (41a)
$\displaystyle\langle|\mathbf{d}_{\perp}|^{2}\rangle=NI$ (41b)
$\displaystyle\langle\mathbf{d}_{\perp}^{2}\rangle=NF$ (41c)
$\displaystyle\langle\overline{\mathbf{d}}_{\perp}^{2}\rangle=NF^{*}.$ (41d)

The integrals are in turn given by

$\displaystyle I=\frac{T}{2}\int_{\bf q}\frac{1}{K\mathbf{q}^{2}+r_{\rm
eff}+\Gamma_{\perp}+|\chi|}+\frac{1}{K\mathbf{q}^{2}+r_{\rm
eff}+\Gamma_{\perp}-|\chi|}=\frac{1}{4\pi
K}\frac{T}{2}\left[\log\left(\frac{C}{r_{\rm
eff}+\Gamma_{\perp}+|\chi|}\right)+\log\left(\frac{C}{r_{\rm
eff}+\Gamma_{\perp}-|\chi|}\right)\right]$ (42a) $\displaystyle
F=-\frac{T}{2}\frac{\chi}{|\chi|}\int_{\bf q}\frac{1}{K\mathbf{q}^{2}+r_{\rm
eff}+\Gamma_{\perp}+|\chi|}-\frac{1}{K\mathbf{q}^{2}+r_{\rm
eff}+\Gamma_{\perp}-|\chi|}=\frac{T}{2}\frac{\chi}{4\pi
K|\chi|}\log\left(\frac{r_{\rm eff}+\Gamma_{\perp}+|\chi|}{r_{\rm
eff}+\Gamma_{\perp}-|\chi|}\right).$ (42b)

The evaluation of $\langle\Delta\mathcal{F}\rangle$ gives

$\langle\Delta\mathcal{F}\rangle=(r-r_{\rm
eff})\langle|\mathbf{d}_{\perp}|^{2}\rangle+r\langle|d_{z}|^{2}\rangle+\frac{u+v}{2N}\langle|\mathbf{d}|^{4}\rangle-\frac{v}{2N}\langle|\mathbf{d}^{2}|^{2}\rangle+\frac{1}{2}\left[\overline{\chi}\langle(\mathbf{d}_{\perp}^{2}+(d_{z}-\sqrt{N}\psi)^{2})\rangle+\chi\langle(\overline{\mathbf{d}}_{\perp}^{2}+(\overline{d}_{z}-\sqrt{N}\overline{\psi})^{2})\rangle\right].$
(43)

The nonlinear terms are now evaluated as

$\frac{1}{N}\langle|\mathbf{d}|^{4}\rangle=N\left[|\psi|^{4}+2|\psi|^{2}I+I^{2}\right]$
(44)

as before, and

$\frac{1}{N}\langle|\mathbf{d}^{2}|^{2}\rangle=N\left[|\psi|^{4}+\psi^{2}\overline{F}+F\overline{\psi}^{2}+|F|^{2}\right].$
(45)

We then get the free energy functional of

$\frac{1}{N}\mathcal{F}=\frac{1}{2}T\int_{\bf
q}\log\left(K\mathbf{q}^{2}+r_{\rm
eff}+\Gamma_{\perp}+|\chi|\right)+\log\left(K\mathbf{q}^{2}+r_{\rm
eff}+\Gamma_{\perp}-|\chi|\right)\\\ +(r-r_{\rm
eff})I+(r+(u+v)I+\frac{u}{2}|\psi|^{2})|\psi|^{2}+\frac{u+v}{2}I^{2}-\frac{v}{2}\left(\psi^{2}\overline{F}+F\overline{\psi}^{2}+|F|^{2}\right)+\frac{1}{2}(\overline{\chi}F+\overline{F}\chi).$
(46)

The variational procedure gives a system of equations

$\displaystyle\left(r+(u+v)I+u|\psi|^{2}\right)\psi-vF\overline{\psi}=0$ (47a)
$\displaystyle r_{\rm eff}=r+(u+v)(I+|\psi|^{2})$ (47b)
$\displaystyle\chi=v(\psi^{2}+F).$ (47c)

With the integrals in turns being functions of $\chi,\overline{\chi},r_{\rm
eff}$.

We have four cases; normal state, $2e$ condensate, $4e$ condensate, and
coexistence. Starting from the normal state, where $\psi=\chi=0$ we imagine
lowering the temperature and find which phase we enter into first. For small
$\Gamma$ we expect to first enter into the $4e$ condensate, which has a phase
boundary determined by

$\displaystyle r_{\rm eff}-r=(u+v)I=(u+v)\frac{T}{4\pi
K}\frac{1}{2}\left[\log\left(\frac{C}{r_{\rm
eff}+\Gamma_{\perp}+|\chi|}\right)+\log\left(\frac{C}{r_{\rm
eff}+\Gamma_{\perp}-|\chi|}\right)\right]$ (48a)
$\displaystyle\chi=vF=v\frac{T}{2}\frac{\chi}{4\pi
K|\chi|}\log\left(\frac{r_{\rm eff}+\Gamma_{\perp}+|\chi|}{r_{\rm
eff}+\Gamma_{\perp}-|\chi|}\right).$ (48b)

We again approximate these by replacing $T=T_{\rm BCS}$ and
$\tilde{\chi}=|\chi|/N(E_{F})$ to get

$2\tilde{\chi}=\textrm{Gi}_{(2)}\textrm{artanh}\left(\frac{\tilde{\chi}}{\tilde{r}_{\rm
eff}+\tilde{\Gamma}_{\perp}}\right).$ (49)

The first equation gives

$\tilde{r}_{\rm
eff}=\tilde{r}+\textrm{Gi}_{(2)}\log\left(\frac{C}{\sqrt{(\tilde{r}_{\rm
eff}+\tilde{\Gamma}_{\perp})^{2}-\tilde{\chi}^{2}}}\right).$ (50)

We find the critical temperature of the transition by noting that a nontrivial
solution for $\tilde{\chi}$ only appears once

$\tilde{r}_{\rm
eff}^{\star}=\frac{1}{2}\textrm{Gi}_{(2)}-\tilde{\Gamma}_{\perp}.$ (51)

At this point $\tilde{\chi}=0$ and we can solve for the physical temperature
$T^{\star}=T_{\rm BCS}e^{\tilde{r}^{\star}}$ to get

$\tilde{r}^{\star}=\frac{1}{2}\textrm{Gi}_{(2)}-\tilde{\Gamma}_{\perp}-\textrm{Gi}_{(2)}\log\left(2C/\textrm{Gi}_{(2)}\right).$
(52)

We see the critical temperature decays exponentially with the induced spin-
orbit coupling.

We see that the two phases boundaries coincide when

$-\textrm{Gi}_{(2)}\log(C/\tilde{\Gamma}_{\perp})=\frac{1}{2}\textrm{Gi}_{(2)}-\tilde{\Gamma}_{\perp}-\textrm{Gi}_{(2)}\log\left(2C/\textrm{Gi}_{(2)}\right).$
(53)

We find solution for the critical $\tilde{\Gamma}_{c}$ of

$\tilde{\Gamma}_{c}+\textrm{Gi}_{(2)}\log(2\tilde{\Gamma}_{c}/\textrm{Gi}_{(2)})=\frac{1}{2}\textrm{Gi}_{(2)}.$
(54)

In fact, this equation can be solved easily to give
$\tilde{\Gamma}_{c}=\frac{1}{2}\textrm{Gi}_{(2)}$. We therefore will have an
intersection of the two different condensate phases and we must consider their
coexistence.

At $\Gamma_{\perp}=0$ and estimateing $1/\sqrt{C}=\textrm{Gi}_{(2)}\sim.1$ we
find $T^{*}$ is of order $.5T_{\rm BCS}$ at zero spin-orbit coupling, which is
remarkably large, atleast before the BKT physics is accounted for.

We start by considering the appearance of the $d_{z}$ condensate out of the
$4e$ condensate, which proceeds by a $\mathbb{Z}_{2}$ spin-nematic transition
associated to the phase of the $d_{z}$ condensate with respect to the phase of
the quasi-long range ordered $4e$ condensate. We find the boundary of this
transition by finding when the finite condensate appears. The relevant
equation is

$\left(r+(u+v)I+u|\psi|^{2}\right)\psi-vF\overline{\psi}=0.$ (55)

We can take the $4e$ order parameter $F$ to be real, in which case we find the
condensate locks into phase with $4e$ condensate, so as to minimize the
Josephson coupling due to the $v$ nonlinearity. Alternatively, we can recall
that this nonlinearity essentially penalized the nonunitary fluctuations and
therefore the lower energy state comes from the condensate being in phase so
that the pairing remains unitary. Thus, the transition is spin nematic but
preserves time-reversal symmetry, rather than the other case of a spin
polarized state emerging which happens when the Josephson coupling between the
$4e$ and $d_{z}$ condensates is negative.

Armed with this, we then proceed to linearize this equation around small
$\psi$ real. We get $\left(r+(u+v)I\right)\psi-vF\psi=0$. This gives a
critical point of

$r_{\rm eff}-vF=0.$ (56)

Here we have noted that $r_{\rm eff}=r+(u+v)I$. The integrals in turn require
solving

$\chi=r_{\rm eff}.$ (57)

Together these give equation

$\tilde{r}_{\rm eff}=\frac{1}{4}\textrm{Gi}_{(2)}\log(\frac{2\tilde{r}_{\rm
eff}+\tilde{\Gamma}_{\perp}}{\tilde{\Gamma}_{\perp}}).$ (58)

This always has a solution for positive $r_{\rm eff}$ provided
$\textrm{Gi}_{(2)}/2\tilde{\Gamma}_{\perp}>1$. At the point where
$\tilde{\Gamma}_{\perp}=\frac{1}{2}\textrm{Gi}_{(2)}$ we see that
$\tilde{r}_{\rm eff}$ goes to zero and thus $\tilde{\chi}$ goes to zero,
indicating indeed this curve terminates at $\tilde{\Gamma}_{c}$ where the $2e$
condensate takes over from the $4e$ condensate. We therefore see that there is
a second transition which occurs within the $4e$ condensed phase characterized
by the spin nematic transition, and this is a continuous second-order
transition persisting up to the point where the $4e$ condensate disappears in
a lambda point.

##  Appendix D Spin-Susceptibility

We now briefly consider the susceptibility of the condensate to an in-plane
magnetic field in the pseudogap region, which may serve as an experimental
signature to help confirm this model. The spin-susceptibility is evaluated as

$\chi_{\parallel}=-\frac{\delta^{2}\mathcal{F}}{\delta\mathbf{H_{\parallel}}^{2}}$
(59)

where the coupling to magnetic field is governed by the coupling
$\mathbf{H}_{\parallel}\cdot(i\overline{\mathbf{d}}\times\mathbf{d})$, and
evaluated at zero external field. We thus get

$\chi_{\parallel}=\frac{\mu_{\rm
eff}^{2}}{T}\langle(i\overline{\mathbf{d}}\times\mathbf{d})(x)(i\overline{\mathbf{d}}\times\mathbf{d})(x^{\prime})\rangle.$
(60)

The correlation functions are evaluated in the disordered phase and to leading
order in $N$ are given by (ignoring the $4e$ condensate)

$\langle{d}_{\alpha}(x)\overline{{d}}_{\beta}(x^{\prime})\rangle=-\int_{\bf
q}e^{i\mathbf{q}\cdot(\mathbf{x-x^{\prime}})}\frac{\mu_{\rm
eff}^{2}T}{K\mathbf{q}^{2}+r_{\rm eff}+\Gamma_{\perp}}\delta_{\alpha\beta}.$
(61)

We then get by Wick’s theorem

$\chi_{\parallel}^{ml}=-\frac{\mu_{\rm
eff}^{2}}{T}\epsilon^{abm}\epsilon^{jkl}\langle\overline{d}_{a}(x)d_{k}(x^{\prime})\rangle\langle
d_{b}(x)\overline{d}_{j}(x^{\prime})\rangle.$ (62)

Contracting indices and picking out the dominant term in $1/N$ we find for the
response to homogeneous field

$\chi_{\parallel}^{lm}(\mathbf{Q}\to 0)=\delta_{lm}N\mu_{\rm
eff}^{2}T\int_{\bf q}\left(\frac{1}{K\mathbf{q}^{2}+r_{\rm
eff}+\Gamma_{\perp}}\right)^{2}.$ (63)

This integral gives the fluctuation contribution of

$\chi_{\parallel}^{lm}(\mathbf{Q}\to 0)=\delta_{lm}\frac{N\mu_{\rm
eff}^{2}T}{4\pi K}\frac{1}{r_{\rm eff}+\Gamma_{\perp}}.$ (64)

We find that due to the superconducting spin fluctuations above $T_{c}$ we
find a divergent paramagnetic response, similar to a Curie-Weiss behavior.
This contribution will be proportional to $\textrm{Gi}_{(2)}$ and is cutoff
for finite spin-orbit coupling at the critical point as

$\chi_{\parallel}\sim\frac{\textrm{Gi}_{(2)}}{\tilde{\Gamma}_{\perp}},$ (65)

while the divergence at zero spin-orbit coupling has form

$\chi_{\parallel}\sim\frac{\textrm{Gi}_{(2)}}{r_{\rm eff}}.$ (66)

This will have the same characteristic critical exponent as the one associated
with $r_{\rm eff}$, which also governs the order parameter susceptibility. We
note that this change in the spin susceptibility is purely due to the soft
triplet order parameter fluctuations, as evidenced by it being proportional to
$\textrm{Gi}_{(2)}$, and in a mean-field transition into a triplet phase it is
expected to vanish. For simplicity, we have neglected the role of the $4e$
pairing although it would be interesting to consider them in the future since
it may have a competing effect since it is essentially spin-singlet. This may
then manifest as an apparent divergence of the spin susceptibility which
abruptly turns over into a singlet-like response at $T^{*}$. Depending on the
strength of fluctuations, this apparent initial upturn may still be visible
and would still be anomalous in comparison to a BCS-like singlet
superconductor.

##  Appendix E Estimates of density of states and the Fermi energy

In this section, we discuss our procedure for estimating the density of states
at the Fermi level and the Fermi energy that are required for calculating the
Ginzburg-Landau parameters.

In Ref[4], inverse compressibility measurements were used to identify the
broken isospin phases in BBG at higher temperatures before the spin-triplet
superconductivity emerges. Inverting the measured inverse compressibility data
directly yields the many-body density of states in the presence of strong
correlation. The blue trace in Fig. 3 shows the density of state at
$D=1.02$V/nm obtained from Ref.[4]. The grey shaded region marks the density
range within which superconductivity is observed. We therefore estimate the
density of state to be 7
$\textrm{meV}^{-1}$\mathrm{\SIUnitSymbolAngstrom}$^{-2}$.

Integrating the same inverse compressibility data allows us to extract the
electronic chemical potential as shown in Fig. 3 (red trace). The black dotted
line denotes the valence band edge and thus the origin of the Fermi energy.
The chemical potential reaches approximately 23 meV at the optimal doping. To
interpret the extracted chemical potential, we note that the system undergoes
a cascade of isospin phase transitions that lead to one large and one small
Fermi surfaces at the optimal doping [4, 49, 5]. These two Fermi surfaces are
each two-fold degenerate and characterized by two different Fermi energies.
Although it is challenging to precisely disentangle their contributions to the
change in chemical potential, one rough estimate is to assume that the large
Fermi surface dominates the total change in chemical potential counting from
the charge neutrality, whereas the small Fermi surface is mainly responsible
for the changes counting from the phase transition point (see green dotted
line in the inset of Fig. 3). Under these assumptions, we find that the Fermi
energy associated with the small Fermi surface is roughly 0.6 meV. Since the
system has gone through already one phase transition that resets the chemical
potential, we estimate the Fermi energy associated with the large Fermi
surface to be roughly half of the total change in chemical potential, that is
approximately 10 meV. While we do not know which Fermi surface is more
responsible for the spin-triplet pairing, we note that Pauli-limit violation
is also observed in twisted graphene multilayers[2, 3] and the Fermi energies
in these systems ($\sim$ 1 meV) are very similar to that of the small Fermi
surface in BBG. Since our theory is motivated solely based on general symmetry
arguments and is broadly applicable to all graphene-based spin-triplet
superconductors, we use 0.6 meV for the Fermi energy in our calculations.

Figure 3:  Density of states and chemical potential for Bernal bilayer
graphene Extracted from Ref. [4].. The black dotted line marks the edge of the
valence band of BBG in the presence of electric field. The grey shaded region
denotes the density range in which spin triplet superconductivity is observed
under small in-plane magnetic field. The inset shows a zoom in near the grey
shaded region. The green dotted line marks the isospin phase transition near
the superconducting dome edge.