# Singlet, triplet, and mixed all-to-all pairing states emerging from
incoherent fermions

Jagannath Sutradhar  Physics Department, Ariel University, Ariel 40700, Israel
Department of Physics, Bar Ilan University, Ramat Gan 5290002, Israel
<EMAIL_ADDRESS>Jonathan Ruhman  Department of Physics, Bar Ilan University,
Ramat Gan 5290002, Israel<EMAIL_ADDRESS>Avraham Klein  Physics
Department, Ariel University, Ariel 40700, Israel

###### Abstract

The electron-electron and electron-phonon coupling in complex materials can be
more complicated than simple density-density interactions, involving
intertwined dynamics of spin, charge, and spatial symmetries. This motivates
studying universal models with complex interactions, and studying whether in
this case BCS-type singlet pairing is still the “natural” fate of the system.
To this end, we construct a Yukawa-SYK model with nonlocal couplings in both
spin and charge channels. Furthermore, we provide for time-reversal-symmetry
breaking dynamics by averaging over the Gaussian Unitary ensemble rather than
the Orthogonal ensemble. We find that the ground state of the system can be an
orbitally nonlocal superconducting state arising from incoherent fermions with
no BCS-like analog. The superconductivity has an equal tendency to triplet and
singlet pairing states separated by a non-Fermi liquid phase. We further study
the fate of the system within the superconducting phase and find that the
expected ground state, away from the critical point, is a mixed
singlet/triplet state. Finally, we find that while at $T_{c}$ the triplet and
singlet transitions are dual to one another, below $T_{c}$ the duality is
broken, with the triplet state more susceptible to orbital fluctuations just
by virtue of its symmetry. Our results indicate that such fluctuation-induced
mixed states may be an inherent feature of strongly correlated materials.

††preprint: APS/123-QED

Introduction– One of the best-known avenues to analyze superconductivity
beyond the classic Bardeen-Cooper-Schrieffer (BCS) paradigm of weakly coupled
Fermi liquids, is studying toy models of fermions coupled to a soft boson [1,
2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. The boson encodes both the
dynamics and the symmetry properties of the interaction. Of particular focus
has been the fate of such models when approaching a quantum critical point
(QCP), such that the diverging correlation length creates both long-ranged and
strong pairing interactions, and drives the fermions incoherent. An
overwhelming majority of studies consider _spin-singlet_ pairing [11, 12, 13,
14, 15]. This is probably because most known superconductors are singlets, and
the most common pairing mediators with the simplest symmetry properties (e.g.
phonons) tend to prefer singlet pairing [16].

Figure 1: (a) Schematic of the model Hamiltonian [Eq. (1)]: $N$ spinful
Fermionic orbitals (circles) coupled to $M$ bosons (wiggly line) via random
coupling $g_{ij,k}$, scattering electrons between orbitals $i$ and $j$. The
scattering occurs in both density ($\sigma^{0}$) and spin ($\sigma^{z}$)
channels with variances $g_{0}^{2}$ and $g_{z}^{2}$. (b) Earlier studies
averaged $g_{ij,k}$ over the GOE, leading to intra-orbital singlet pairing.
(c) In our study, $g_{ij,k}$ is drawn from the GUE, resulting in an “all-to-
all” inter-orbital pairing between all possible up-down spin pairs. (d) The
phase diagram showing dual triplet and singlet SC phases, terminating at QCPs
and separated by an NFL phase. At $T<T_{c}$, both singlet and triplet phases
can induce a coexisting singlet/triplet state, with the SC ground state
transiting to a mixed pairing state. The dots represent a numerical solution
of the gap equation and the lines depict a fit to our analytical results, see
the text for further details.

In truth, though, sufficiently complex materials should host interactions with
intertwined attraction and repulsion, involving both charge and spin, as well
as fluctuations that break lattice, time-reversal, and inversion symmetries.
Such interactions have been considered for a host of candidate materials [17,
18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]. In addition, a variety of
superconducting systems appear to evince phenomena that are neither clearly
singlet-like nor triplet-like, raising the question of how likely mixed
singlet/triplet states are in these systems [19, 31, 32, 33, 34]. This
motivates finding a universal model of strongly correlated pairing that is
“agnostic” as regards spin singlet and triplet pairing, as well as various
spacetime and lattice symmetries. With such a model one can study the
properties of strongly correlated superconductors, investigate whether spin
singlet and triplet pairing have equal footing or not, and look for exotic
pairing states that might be missing in “garden-variety” models.

In this Letter, we develop and study such a model. Our starting point is a
Yukawa-Sachdev-Ye-Kitaev (Y-SYK) model [30, 35, 13], where $N$ fermions are
randomly coupled to $M$ bosons in 0+1D, as shown in Fig. 1(a). This model is
analytically solvable in the large-$N$ large-$M$ limit and exhibits both non-
Fermi liquid (NFL) behavior and high-temperature superconductivity. We
consider both spin and charge coupling, with independent coupling constants,
$g_{z}$ and $g_{0}$, respectively. Moreover, we allow for time-reversal-
symmetry (TRS) breaking fluctuations by drawing the coupling matrix from the
Gaussian Unitary Ensemble (GUE), rather than the more usual Orthogonal
Ensemble (GOE).

We find that upon lowering the temperature the bosons become critical and the
NFL becomes unstable to an _inter-orbital_ pairing state $\Phi\sim\langle
c_{i}c_{j}\rangle,i\neq j$ [Fig. 1(c)], where $i,j$ denote orbitals and we
suppressed spin indices. Such pairing is qualitatively distinct from the
“standard” BCS type $\Phi\sim\langle c_{i}c_{i}\rangle$ that arises in the GOE
averaging [Fig. 1(b)]. Furthermore, this nonlocal state is tuned between
triplet and singlet pairing depending on the relative coupling strengths as
shown in the phase diagram Fig. 1(d). The two states are separated by an NFL
phase. We also study the pairing state below $T_{c}$. We find that the low-$T$
pairing state is generically a nonlocal triplet-singlet mixture. Our results
show that the previously deemed NFL phase can become unstable to pairing once
the non-local state is taken into account. This includes models where the
random couplings are a sum of two random numbers drawn from the GOE and GUE
[13].

Thus, systems with both charge and spin fluctuations have a rich landscape of
superconducting phases that go beyond the usual paradigms of unconventional
superconductivity and can involve not just broken symmetries such as time
reversal, but also spatially nonlocal correlations even at the mean-field
level. In these systems, an introduction of nonlinear contributions plays a
crucial role in shaping the phase diagram, even in the absence of static
symmetry breaking. We describe our model and results in detail below.

Model– We consider $M$ bosonic fields, $\phi_{k}$, randomly coupled to $N$
fermions $c_{i\alpha}$ as schematically depicted in Fig. 1(a) and expressed as

$\displaystyle H$ $\displaystyle=H_{z}+H_{0}\
+\sum_{k}(\pi_{k}^{2}+m^{2}_{0}\phi_{k}^{2})\,,$ (1) $\displaystyle H_{a}$
$\displaystyle=1/\sqrt{MN}\sum_{ijk\alpha\beta}g^{a}_{ij,k}\phi_{k}c_{i\alpha}^{\dagger}\sigma^{a}_{\alpha\beta}c_{j\beta}\quad(a=z,0),$

where $\pi_{k}$ are the momenta fields conjugate to $\phi_{k}$, and
$\alpha,\beta$ are spin indices. The random couplings $g^{a}_{ij,k}$ are
completely uncorrelated between the spin and charge components, denoted by
$a=z,0$ respectively. The interaction form is chosen to allow for universal
spin and charge fluctuations, while avoiding stability issues that arise from
SU(2) spin fluctuations [36, 37]. The couplings are drawn from the GUE, such
that the average of the real and imaginary parts:
$\overline{\mathrm{Re}(g^{a}_{ij,k})\mathrm{Re}(g^{a}_{i^{\prime}j^{\prime},k^{\prime}}})=g_{a}^{2}\delta_{kk^{\prime}}(\delta_{ii^{\prime}}\delta_{jj^{\prime}}+\delta_{ij^{\prime}}\delta_{ji^{\prime}})$
and
$\overline{\mathrm{Im}(g^{a}_{ijk})\mathrm{Im}(g^{a}_{i^{\prime}j^{\prime}k^{\prime}})}=g_{a}^{2}\delta_{kk^{\prime}}(\delta_{ii^{\prime}}\delta_{jj^{\prime}}-\delta_{ij^{\prime}}\delta_{ji^{\prime}})$,
respectively and $g_{a}$ is the coupling strength. Note that the complex
$g^{a}_{ij,k}$ manifestly breaks TRS, which is restored upon disorder-
averaging. Previous studies considered a similar model with TRS by drawing the
coupling from the GOE. In that case, there is a large-$N$ instability towards
an intra-orbital singlet pairing phase characterized by the non-vanishing
expectation value $\sum_{i}\langle
c_{i\uparrow}^{\dagger}c_{i\downarrow}^{\dagger}\rangle\neq 0$. In contrast,
the GUE-averaged action avoids the intra-orbital pairing and opens a path to
more exotic pairing states. The distinct pairing forms are illustrated in
Figs. 1(b) and (c).

The triplet normal state and pairing instability– Before jumping into the full
Hamiltonian in Eq. (1), we start by focusing on the spin fluctuation $H_{z}$,
which governs the triplet pairing and set $g_{0}=0$. We then show that the
complimentary singlet pairing scenario where $g_{z}=0,g_{0}\neq 0$ can be
obtained directly from the triplet case. Finally we study the interplay
between both spin and charge couplings.

As usual, we disorder-average the partition function associated with the
Hamiltonian in Eq. (1) with $H_{0}=0$ using the replica trick and consider
only the replica-diagonal term in the action [30]. The form of the interacting
part of the action is

$\mathcal{S}_{I}\sim-g_{z}^{2}\sum_{k,i,..}\int
d\tau\phi_{k}c_{i\alpha}^{\dagger}\sigma^{z}_{\alpha\beta}c_{j\beta}\int
d\tau^{\prime}\phi_{k}c_{j\alpha^{\prime}}^{\dagger}\sigma^{z}_{\alpha^{\prime}\beta^{\prime}}c_{i\beta^{\prime}}\,,$
(2)

(see the supplemental materials [SM] for details), and the corresponding
Feynman diagram is shown in Fig. 2(a). Superficially, Eq. (2) appears similar
to “standard” spin-fermion interactions whose QC dynamics have been studied
extensively [11, 12, 13, 14, 15, 38, 39]. However, the interaction in $H_{z}$
of Eq. (1) involves both spin-flips _and_ nonlocality, very different from the
usual onsite magnetic interactions, giving rise to some nontrivial physics. In
what follows we study the fate of the fermions in the presence of this
interaction.

Figure 2:  Feynman diagrams: (a) Schematic of the interacting four-fermionic
term mediated by the bosonic propagator $D$ (wavy line) after disorder average
represented by the dashed line. (b) The linearized gap equation. Note that the
exchange in $ij$ indices in $\hat{\Phi}$ on the right side makes one of the
triplet channels attractive.

Connecting different legs in Fig. 2(a) we obtain the self-energies in the
large $N,M$ limit. The Green’s function is diagonal both in orbital and spin
basis [SM] and it is defined as
$G(\omega_{n})=1/[\mathrm{i}\omega_{n}+\mathrm{i}\Sigma(\omega_{n})]$. At
finite $T$, the electronic self-energy is given by

$\Sigma(\omega_{n})=\mathrm{i}g^{2}_{z}T\sum_{m}D(\omega_{n}-\omega_{m})\dfrac{1}{2}\operatorname{Tr}\big{[}\sigma^{z}G(\omega_{m})\sigma^{z}\big{]}\,,$
(3)

where $\omega_{n}=(2n+1)\pi T$, $n\in\mathbb{Z}$ are fermionic Matsubara
frequencies and Tr represents trace over the spin indices.
$D(\Omega_{n})=[\Omega_{n}^{2}+m_{0}^{2}+\Pi(\Omega_{n})]^{-1}$ is the bosonic
propagator. The bosonic self-energy is given by

$\Pi(\Omega_{n})=2g^{2}_{z}T\dfrac{N}{M}\sum_{m}\operatorname{Tr}\big{[}\sigma^{z}G(\omega_{m}+\Omega_{n})\sigma^{z}G(\omega_{m})\big{]}\,,$
(4)

where $\Omega_{n}=2n\pi T$ are the bosonic Matsubara frequencies. Similar to
what occurs in previously studied models [35, 13], Eqs. (3) and (4) give rise
to NFL behavior, with
$\Sigma(\omega)\sim\mathrm{sgn}(\omega)|\omega|^{(1-\eta)/2}$ and
$\Pi(\Omega)-\Pi(0)\sim g^{2}|\Omega|^{\eta}$ at low $T,\omega$, and self-
tuned criticality that renormalizes the boson mass to exactly zero. The
exponent $\eta$ varies from $0$ to $1$ as a function of $N/M$ [SM].

The form of the interaction naturally leads us to consider an inter-orbital
pairing state
$\Phi_{ij}^{\alpha\beta}=\braket{c_{i\alpha}^{\dagger}c_{j\beta}^{\dagger}}$,
which includes all possible orbital pairs. A somewhat similar type of pairing
was considered in Ref. [35]. The linearized equation for the inter-orbital
order parameter, which corresponds to Fig. 2(b), is given by

$\displaystyle\hat{\Phi}^{\alpha\beta}_{ij}(\omega_{n})=\dfrac{g^{2}_{z}}{N}T\sum_{m,\alpha^{\prime},\beta^{\prime}}D(\omega_{n}-\omega_{m})\sigma^{z}_{\alpha\alpha^{\prime}}\qquad\qquad$
$\displaystyle
G(\omega_{m})\hat{\Phi}^{\alpha^{\prime}\beta^{\prime}}_{ji}(\omega_{m})G(-\omega_{m})\sigma^{z}_{\beta\beta^{\prime}}\,.$
(5)

Note the orbital exchange $i\leftrightarrow j$ between the left and right-hand
side of the equation, which is a consequence of the nonlocality discussed
above. The antisymmetry of the pairing function is enforced by the orbital
exchange, i.e.,
$\hat{\Phi}_{ij}^{\alpha\beta}=-\hat{\Phi}_{ji}^{\alpha\beta}$. Consequently,
Eq. (5) indicates attraction only when
$\hat{\Phi}^{\alpha\beta}_{ij}=\hat{\Phi}_{ij}\sigma^{x}_{\alpha\beta}$, which
signifies a triplet pairing state. Note that this is the _only_ attractive
channel, with the other two spin channels, as well as the singlet channel,
being repulsive. This is in stark contrast to what occurs in quantum critical
models with a local spin-fermion interaction, where each magnetic channel is
attractive in two spin channels and repulsive in the third [36]. To simplify
Eq. (5), we neglect the orbital phase fluctuation and consider
$\hat{\Phi}_{ij}=\Phi_{t}A_{ij}$, where $\Phi_{t}$ is the amplitude of the
pairing function dependent only on frequency, and $A_{ij}$ is a skew-symmetric
matrix, such that all elements $A_{i>j}=1,A_{i<j}=-1,A_{i=j}=0$. Choosing a
different antisymmetric matrix $A_{ij}$ yields the same results [SM]. We will
discuss the origin and consequence of this degeneracy later on.

At $T=0$ Eq. (5) has a SC instability in the regime $\sqrt{2M}>N$ and in the
large $M$ and $N$ limit [35, 13]. In contrast to the BCS theory, for
$\sqrt{2M}<N$ even strong attraction does not lead to SC. At finite $T$, the
bosons develop a mass correction, and the pairing equation resembles that in
Ref. [13]. The critical temperature ($T_{c}$) for the phase transition from
the NFL to SC follows $T_{c}\sim g^{2}_{z}$ [SM].

The singlet channel and mixed case– We next consider the pure singlet case,
$g_{z}=0,\ g_{0}\neq 0$, the gap equation exhibits attraction in the singlet
channel, characterized by the pairing function
$\hat{\Phi}_{ij}^{\alpha\beta}=\mathrm{i}\Phi_{s}S_{ij}\sigma^{y}_{\alpha\beta}$,
where $S$ is a symmetric traceless matrix. As before, we neglect the
fluctuations and set all the off-diagonal elements of the matrix to be $1$,
$S_{ij}=1-\delta_{ij}$ (for further discussion see SM). The normal state and
SC gap equations remain the same as in the triplet case. Consequently, the
phase diagram resembles that in the earlier case, replacing the triplet SC
with a singlet SC. Thus, above the critical $\sqrt{M}/N$ ratio, the ground
state for our model with either of the couplings is always SC arising from
within the incoherent NFL state. This is again in contrast with the case of
GOE averaging, where the BCS-type local pairing sets in at the same energy
scale as NFL physics [30, 13]. For the singlet case, a BCS-type local pairing
$\hat{\Phi}^{\alpha\beta}_{ij}\propto\delta_{ij}i\sigma^{y}_{\alpha\beta}$ is
degenerate with the nonlocal pairing. For simplicity we consider henceforth
only the nonlocal state, since it keeps the physics qualitatively the same
[SM].

We are now ready to consider the full Hamiltonian of Eq. (1). As
$g^{0}_{ij,k}$ and $g^{z}_{ij,k}$ are uncorrelated, the interacting part of
the disorder averaged action is a sum of two terms similar to Fig. 2(a): one
for $\sigma^{z}$ as earlier and another for $\sigma^{0}$, with their
respective strengths $g_{z}$ and $g_{0}$. The normal state solution remains
NFL. Within the linearized gap equation, there is a critical strength
$\lambda_{c}$ of $\lambda=|g_{z}^{2}-g_{0}^{2}|$ beyond which the SC phase
emerges. The state is either singlet (for $g_{0}>g_{z}$) or triplet (for
$g_{z}>g_{0}$), with no mixing [see Fig. 1(d)]. The critical temperature
follows $T_{c}\sim\mathrm{exp}(-\pi/\sqrt{\lambda-\lambda_{c}})$ [13][SM].
Between the SC phases, the system remains an NFL down to $T=0$.

Solution of the gap equation below $T_{c}$– To explore the phase diagram at
$T<T_{c}$ we must turn to the nonlinear gap equation. To begin the analysis,
let us consider the case with only one finite coupling, $\sigma^{z}$ or
$\sigma^{0}$. Obtaining an analytic solution to the full non-linear equations
is challenging due to the form of the matrices $A$ and $S$, which couple all
orbitals. However, near $T_{c}$, we can study the first non-linear correction
to the linearized gap equation by dressing one of the Green’s functions in the
diagram Fig. 2(b) with a correction of order
$\sim|\Phi_{t(s)}(\omega_{m})|^{2}\Phi_{t(s)}(\omega_{m})$, as depicted in
Fig. 3(a) [SM]. The nonlinear term has an opposite sign compared to the linear
term [Eq. (6)] and gives rise to the usual $\Phi_{a}\sim\sqrt{T_{c}-T}$
behavior.

We now turn to the case involving both $\sigma^{z}$ and $\sigma^{0}$
couplings. We introduce a _mixed_ pairing function
$\hat{\Phi}=\Phi_{t}A\sigma^{x}+\mathrm{i}\Phi_{s}S\sigma^{y}$. As mentioned
earlier, in the linear approximation, $\Phi_{t}$ and $\Phi_{s}$ do not
interact. However, higher-order terms introduce an interaction between
$\Phi_{t}$ and $\Phi_{s}$. The equation for $\Phi_{t(s)}$ can be written in
simplified form as,

$\Phi_{t(s)}(\omega_{n})\approx\big{[}g_{z(0)}^{2}-g_{0(z)}^{2}\big{]}\dfrac{T}{N}\sum_{m}\dfrac{D(\omega_{n}-\omega_{m})}{|\omega_{m}+\Sigma(\omega_{m})|^{2}}\Bigg{[}\Phi_{t(s)}(\omega_{m})-\dfrac{\Phi_{t(s)}(\omega_{m})}{|\omega_{m}+\Sigma(\omega_{m})|^{2}}\Big{\\{}\zeta_{t(s)}\Phi_{t(s)}^{2}(\omega_{m})+\dfrac{2}{3}\Phi^{2}_{s(t)}(\omega_{m})\Big{\\}}\Bigg{]}\,,$
(6)

where we rescaled $\Phi_{t(s)}\rightarrow\Phi_{t(s)}/N$ and neglected
$1/N^{2}$ terms. The factor $\zeta_{t}=1/3$ and $\zeta_{s}=1$. We also took
both fields to be real, since their phase turns out to be locked by the
nonlinear terms (see SM for the full equation and analysis).

The non-linear term implies the emergence of a mixed pairing state at lower
temperatures. For example, let us consider $g_{z}^{2}>g_{0}^{2}$ and $T$ is
just below $T_{c}$, which results in $\Phi_{t}\neq 0$ and $\Phi_{s}=0$, i.e.,
the SC phase consists purely of triplet pairing. Now, in the equation for
$\Phi_{s}$ the nonlinear contribution $\propto\Phi_{t}^{2}\Phi_{s}$ is
attractive although the linear term is repulsive. Consequently, as we lower
the temperature and $|\Phi_{t}|$ increases, the attraction at low frequencies
dominates over the linear repulsive term, triggering a second instability
below some critical $T$. Thus, the SC phase becomes a superposition of singlet
and triplet pairing. This mechanism is reminiscent of the competition between
repulsive and attractive interactions in the Anderson-Morel model [40]. The
arguments for $g_{0}^{2}>g_{z}^{2}$ are exactly analogous, and this gives rise
to the phase diagram plotted in Fig. 1(d) for $\eta=0.68$ and $N=1$. We
confirm our analysis by numerically determining the critical temperature
$T_{t(s)\rightarrow m}$ [SM]. An analysis of the nonlinear gap equation within
the SC phase reveals that the mixed state is stable all the way to the
dominant order’s critical point at $\lambda=\lambda_{c}$, and our numerics
confirm this analysis. This contradicts the decoupling of the two states that
occurs at the level of the linearized gap equation. Thus, the QCP appears to
be a singular point with different behavior depending on how it is approached
in the phase diagram. This is probably because the Ginzburg-Landau expansion
has coefficients that diverge at low temperatures. To conclusively determine
whether the QCPs to the single state and mixed state converge or not requires
an analysis of the full gap equation, which is beyond the scope of this work.

Figure 3: Nonlinearity in the gap equation: (a) The Feynman diagram represents
the first nonlinear term in the expansion of the gap equation. (b) Variation
of the critical temperature $T_{t\rightarrow m}$ to the mixed state as a
function of $g_{z}^{2}-g_{0}^{2}$. This figure zooms at the segment of Fig.
1(d) depicting the transition from the triplet to the mixed state. The squares
denote the numerical data points and the dashed line represents a fitting to
the curve $T_{t\rightarrow
m}\sim\mathrm{exp}\big{(}-\pi(4/(3+\eta))/\sqrt{\lambda-\lambda_{c}}\big{)}$
[SM] for $g_{0}^{2}=1$.

The phase diagram Fig. 1(d) implies a duality between the singlet and triplet
states. However, an interesting distinction exists between them. To understand
it, we recall that the $T_{c}$ obtained in the triplet state does not depend
on the choice of the antisymmetric matrix $A$. The existence of a degeneracy
is a manifestation of an $S_{N}$ permutation symmetry of the disorder-averaged
action Eq. (2). On the contrary, no such degeneracy exists in the singlet
case, such that $S$ is unique. Indeed, the matrix $S$ forms a _one-
dimensional_ irreducible representation (irrep) of $S_{N}$, while $A$ belongs
to a multidimensional irrep, with $N(N-1)/2$ different elements. Nonetheless,
it should be noted that in any case the gap between the different irreps
vanishes in the $N\to\infty$ limit.

Discussion– We have studied the Yukawa-SYK model for a random coupling between
$M$ bosons and $N$ fermions drawn from the GUE, which corresponds to the limit
where time-reversal symmetry is absent. We couple the bosons to both spin and
charge. The intra-orbital pairing term vanishes upon averaging the random
coupling (diagonal replica symmetry breaking). In the absence of this term,
the leading pairing instability is a highly non-local inter-orbital pairing
state with an equal superposition on all orbitals. By tuning the relative
strength between the charge and spin couplings the ground state is tuned
through a sequence of QCPs, separating between singlet, NFL and triplet
phases. Going deeper into each one of these phases, there is an additional
transition to a singlet-triplet mixed state, as shown in the phase diagram
Fig. 1(d).

It is interesting to compare the inter- and intra-orbital pairing order
parameters with the standard BCS pairing wave function in a finite-dimensional
system with translational invariance (i.e. with a Fermi surface). In the BCS
wave function, the pair correlations are between $k$ and $-k$ and are
therefore equivalent to intra-orbital pairing. On the other hand, the inter-
orbital order parameter we presented here would imply pairing correlations
between all momentum states. Thus, the wavefunction corresponding to the
inter-orbital pairing state is inconsistent with the picture of momentum-space
quasiparticles. In that sense, one may think of a GOE-averaged Y-SYK system as
the quasiparticle weak coupling limit, and the GUE averaged system as a strong
coupling limit. Moreover, the existence of an instability in the entire phase
diagram implies that the NFL phase driven by TRS breaking is generically
unstable to pairing in the Y-SYK model.

It is also interesting that the triplet order parameter belongs to a _multi-
dimensional_ irreducible representation of the permutation group, while the
singlet belongs to a one-dimensional irrep. Consequently, we expect the
triplet state to be more susceptible to quantum fluctuations, and thus the
phase diagram in Fig. 1(d) will not retain the symmetry between the
$g_{z}>g_{0}$ and $g_{z}<g_{0}$ limits where the transition into the triplet
state will likely be suppressed.

Finally, the Y-SYK action has a 1+1 AdS dual [41, 42, 43, 44], where $T_{c}$
is dual to an event horizon (indicating a black hole formation). Looking
forward it is interesting to understand the gravity dual of the inter-orbital
state found here. In particular, it may give access to black holes with finite
angular momentum.

###### Acknowledgements.

We thank J. Schmalian, A.V. Chubukov, S.K. Saha, and H. Yerzhakov for
interesting discussions. A.K. and J.R. acknowledge support by the Israel
Science Foundation (ISF), and the Israeli Directorate for Defense Research and
Development (DDR&D) under grant No. 3467/21.

## References

  * Millis [1993] A. J. Millis, Effect of a nonzero temperature on quantum critical points in itinerant fermion systems, Phys. Rev. B 48, 7183 (1993).
  * Abanov and Chubukov [2000] A. Abanov and A. V. Chubukov, Spin-fermion model near the quantum critical point: One-loop renormalization group results, Physical Review Letters 84, 5608–5611 (2000).
  * Wang  _et al._ [2001] Z. Wang, W. Mao, and K. Bedell, Superconductivity near itinerant ferromagnetic quantum criticality, Phys. Rev. Lett. 87, 257001 (2001).
  * Abanov  _et al._ [2001] A. Abanov, A. V. Chubukov, and J. Schmalian, Quantum-critical superconductivity in underdoped cuprates, Europhysics Letters (EPL) 55, 369–375 (2001).
  * Abanov  _et al._ [2003] A. Abanov, A. V. Chubukov, and J. Schmalian, Quantum-critical theory of the spin-fermion model and its application to cuprates: Normal state analysis, Advances in Physics 52, 119 (2003).
  * Chubukov and Schmalian [2005] A. V. Chubukov and J. Schmalian, Superconductivity due to massless boson exchange in the strong-coupling limit, Phys. Rev. B 72, 174520 (2005).
  * Wang  _et al._ [2016] Y. Wang, A. Abanov, B. L. Altshuler, E. A. Yuzbashyan, and A. V. Chubukov, Superconductivity near a quantum-critical point: The special role of the first Matsubara frequency, Phys. Rev. Lett. 117, 157001 (2016).
  * Metlitski  _et al._ [2015] M. A. Metlitski, D. F. Mross, S. Sachdev, and T. Senthil, Cooper pairing in non-fermi liquids, Phys. Rev. B 91, 115111 (2015).
  * Raghu  _et al._ [2015] S. Raghu, G. Torroba, and H. Wang, Metallic quantum critical points with finite bcs couplings, Phys. Rev. B 92, 205104 (2015).
  * Lederer  _et al._ [2015] S. Lederer, Y. Schattner, E. Berg, and S. A. Kivelson, Superconductivity and non-fermi liquid behavior near a nematic quantum critical point, Phys. Rev. B 92, 205104 (2015).
  * Pan  _et al._ [2021] G. Pan, W. Wang, A. Davis, Y. Wang, and Z. Y. Meng, Yukawa-SYK model and self-tuned quantum criticality, Phys. Rev. Res. 3, 013250 (2021).
  * Choi  _et al._ [2022] W. Choi, O. Tavakol, and Y. B. Kim, Pairing instabilities of the Yukawa-SYK models with controlled fermion incoherence, SciPost Phys. 12, 151 (2022).
  * Classen and Chubukov [2021] L. Classen and A. Chubukov, Superconductivity of incoherent electrons in the Yukawa Sachdev-Ye-Kitaev model, Phys. Rev. B 104, 125120 (2021).
  * Wang  _et al._ [2021] W. Wang, A. Davis, G. Pan, Y. Wang, and Z. Y. Meng, Phase diagram of the spin-$\frac{1}{2}$ Yukawa–Sachdev-Ye-Kitaev model: Non-fermi liquid, insulator, and superconductor, Phys. Rev. B 103, 195108 (2021).
  * Valentinis  _et al._ [2023] D. Valentinis, G. A. Inkof, and J. Schmalian, BCS to incoherent superconductivity crossover in the Yukawa-Sachdev-Ye-Kitaev model on a lattice, Phys. Rev. B 108, L140501 (2023).
  * Brydon  _et al._ [2014] P. M. R. Brydon, S. Das Sarma, H.-Y. Hui, and J. D. Sau, Odd-parity superconductivity from phonon-mediated pairing: Application to ${\mathrm{cu}}_{x}{\mathrm{bi}}_{2}{\mathrm{se}}_{3}$, Phys. Rev. B 90, 184512 (2014).
  * Gor’kov and Rashba [2001] L. P. Gor’kov and E. I. Rashba, Superconducting 2d system with lifted spin degeneracy: Mixed singlet-triplet state, Phys. Rev. Lett. 87, 037004 (2001).
  * Bauer  _et al._ [2004] E. Bauer, G. Hilscher, H. Michor, C. Paul, E. W. Scheidt, A. Gribanov, Y. Seropegin, H. Noël, M. Sigrist, and P. Rogl, Heavy fermion superconductivity and magnetic order in noncentrosymmetric ${\mathrm{c}\mathrm{e}\mathrm{p}\mathrm{t}}_{3}\mathrm{S}\mathrm{i}$, Phys. Rev. Lett. 92, 027003 (2004).
  * Frigeri  _et al._ [2004] P. A. Frigeri, D. F. Agterberg, A. Koga, and M. Sigrist, Superconductivity without inversion symmetry: Mnsi versus ${\mathrm{c}\mathrm{e}\mathrm{p}\mathrm{t}}_{3}\mathrm{S}\mathrm{i}$, Phys. Rev. Lett. 92, 097001 (2004).
  * Vorontsov  _et al._ [2009] A. B. Vorontsov, M. G. Vavilov, and A. V. Chubukov, Interplay between magnetism and superconductivity in the iron pnictides, Phys. Rev. B 79, 060508 (2009).
  * Chubukov  _et al._ [2003a] A. V. Chubukov, A. M. Finkel’stein, R. Haslinger, and D. K. Morr, First-order superconducting transition near a ferromagnetic quantum critical point, Phys. Rev. Lett. 90, 077002 (2003a).
  * Rech  _et al._ [2006] J. Rech, C. Pépin, and A. V. Chubukov, Quantum critical behavior in itinerant electron systems: Eliashberg theory and instability of a ferromagnetic quantum critical point, Phys. Rev. B 74, 195126 (2006).
  * Bert  _et al._ [2011] J. A. Bert, B. Kalisky, C. Bell, M. Kim, Y. Hikita, H. Y. Hwang, and K. A. Moler, Direct imaging of the coexistence of ferromagnetism and superconductivity at the LaAlO3/SrTiO3 interface, Nature Physics 7, 767 (2011).
  * Chubukov [2012] A. Chubukov, Pairing mechanism in Fe-based superconductors, Annual Review of Condensed Matter Physics 3, 57 (2012).
  * Scalapino [2012] D. J. Scalapino, A common thread: The pairing interaction for unconventional superconductors, Rev. Mod. Phys. 84, 1383 (2012).
  * Hinojosa  _et al._ [2014] A. Hinojosa, R. M. Fernandes, and A. V. Chubukov, Time-reversal symmetry breaking superconductivity in the coexistence phase with magnetism in Fe Pnictides, Phys. Rev. Lett. 113, 167001 (2014).
  * Klein and Chubukov [2018] A. Klein and A. Chubukov, Superconductivity near a nematic quantum critical point: Interplay between hot and lukewarm regions, Phys. Rev. B 98, 220501 (2018).
  * Pimenov and Chubukov [2022] D. Pimenov and A. Chubukov, Quantum phase transition in a clean superconductor with repulsive dynamical interaction, npj Quantum Mater 7, 45 (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).
  * Esterlis and Schmalian [2019] I. Esterlis and J. Schmalian, Cooper pairing of incoherent electrons: An electron-phonon version of the Sachdev-Ye-Kitaev model, Physical Review B 100, 10.1103/physrevb.100.115132 (2019).
  * Pustogow  _et al._ [2019] A. Pustogow, Y. Luo, A. Chronister, Y.-S. Su, D. Sokolov, F. Jerzembeck, A. P. Mackenzie, C. W. Hicks, N. Kikugawa, S. Raghu, _et al._ , Constraints on the superconducting order parameter in Sr2RuO4 from oxygen-17 nuclear magnetic resonance, Nature 574, 72 (2019).
  * Ribak  _et al._ [2020] A. Ribak, R. M. Skiff, M. Mograbi, P. Rout, M. Fischer, J. Ruhman, K. Chashka, Y. Dagan, and A. Kanigel, Chiral superconductivity in the alternate stacking compound 4Hb-TaS2, Science advances 6, eaax9480 (2020).
  * Salmani-Rezaie  _et al._ [2020] S. Salmani-Rezaie, K. Ahadi, and S. Stemmer, Polar nanodomains in a ferroelectric superconductor, Nano Letters 20, 6542 (2020).
  * Hayes  _et al._ [2021] I. Hayes, T. Metz, S. Saha, J. Collini, N. Butch, D. Agterberg, A. Kapitulnik, and J. Paglione, Multicomponent superconducting order parameter in UTe2, Science 373, 797 (2021).
  * Wang [2020] Y. Wang, Solvable Strong-Coupling Quantum-Dot Model with a Non-Fermi-Liquid Pairing Transition, Phys. Rev. Lett. 124, 017002 (2020).
  * Chubukov  _et al._ [2003b] A. V. Chubukov, A. M. Finkel’stein, R. Haslinger, and D. K. Morr, First-order superconducting transition near a ferromagnetic quantum critical point, Phys. Rev. Lett. 90, 077002 (2003b).
  * Chubukov and Maslov [2009] A. V. Chubukov and D. L. Maslov, Spin conservation and fermi liquid near a ferromagnetic quantum critical point, Phys. Rev. Lett. 103, 216401 (2009).
  * Wu  _et al._ [2020] Y.-M. Wu, A. Abanov, Y. Wang, and A. V. Chubukov, Interplay between superconductivity and non-fermi liquid at a quantum critical point in a metal. II. the $\gamma$ model at a finite $t$ for $0<\gamma<1$, Phys. Rev. B 102, 024525 (2020).
  * Wu  _et al._ [2021] Y.-M. Wu, S.-S. Zhang, A. Abanov, and A. V. Chubukov, Interplay between superconductivity and non-fermi liquid behavior at a quantum-critical point in a metal. V. The $\gamma$ model and its phase diagram: The case $\gamma=2$, Phys. Rev. B 103, 184508 (2021).
  * Morel and Anderson [1962] P. Morel and P. W. Anderson, Calculation of the superconducting state parameters with retarded electron-phonon interaction, Phys. Rev. 125, 1263 (1962).
  * Hartnoll  _et al._ [2008] S. A. Hartnoll, C. P. Herzog, and G. T. Horowitz, Building a holographic superconductor, Phys. Rev. Lett. 101, 031601 (2008).
  * Sachdev [2015] S. Sachdev, Bekenstein-Hawking entropy and strange metals, Phys. Rev. X 5, 041025 (2015).
  * Schmalian [2022] J. Schmalian, Holographic superconductivity of a critical fermi surface, arXiv 2209.00474 (2022).
  * Inkof  _et al._ [2022] G.-A. Inkof, K. Schalm, and J. Schmalian, Quantum critical Eliashberg theory, the Sachdev-Ye-Kitaev superconductor and their holographic duals, npj Quantum Materials 7, 56 (2022).

Supplementary Material

In these supplementary material notes, we provide details on the calculations
we referred to in the main text. We start with technical details regarding the
action and the evaluation of the normal state properties. Most of these
details can be found in previous works, but we include them for completeness.
Then we give details on our solutions of both linear and nonlinear gap
equations at finite $T$, and with both $g_{z}$ and $g_{0}$.

##  S1 Evaluation of the action

The partition function is defined as $Z=e^{-H/T}$ and the Hamiltonian $H$ is
given in Eq. (1) of the main text. We begin with deriving the action only for
the spin channel, i.e. for $g_{0}=0$. The disorder-averaged action is obtained
using the replica trick. To do so, we first calculate
$\braket{Z^{R}}_{\text{GUE}}$, where $R$ is an integer and
$\braket{..}_{\text{GUE}}$ denotes average over the GUE disorder
$g^{z}_{ij,k}$. The interacting part of the action is expressed as follows:

$\mathcal{S}_{I}=-g_{z}^{2}\int d\tau
d\tau^{\prime}\sum_{a,b=1}^{R}\sum_{k}\phi_{k,a}(\tau)\phi_{k,b}(\tau^{\prime})\sum_{ij\alpha\beta\alpha^{\prime}\beta^{\prime}}c_{i\alpha,a}^{\dagger}(\tau)\sigma^{z}_{\alpha\beta}c_{j\beta,a}(\tau)\
c^{\dagger}_{j\alpha^{\prime},b}(\tau^{\prime})\sigma^{z}_{\alpha^{\prime}\beta^{\prime}}c_{i\beta^{\prime},b}(\tau^{\prime})\,.$
(S1)

Here $a,b=1,...,R$ signifies the replica indices. For simplification and
further analysis, we assume the contribution in $\mathcal{S}_{I}$ arising only
from the diagonal elements in the replica basis, allowing us to eliminate the
replica structure. Consequently, the interacting part is given by

$\mathcal{S}_{I}=-g_{z}^{2}\int d\tau
d\tau^{\prime}\sum_{k}\phi_{k}(\tau)\phi_{k}(\tau^{\prime})\sum_{ij\alpha\beta\alpha^{\prime}\beta^{\prime}}c_{i\alpha}^{\dagger}(\tau)\sigma^{z}_{\alpha\beta}c_{j\beta}(\tau)\
c^{\dagger}_{j\alpha^{\prime}}(\tau^{\prime})\sigma^{z}_{\alpha^{\prime}\beta^{\prime}}c_{i\beta^{\prime}}(\tau^{\prime})\,,$
(S2)

which is represented using Feynman diagram in Fig. 2(a). In the presence of
both spin and charge couplings, the interaction is expressed as

$\mathcal{S}_{I}=-\int d\tau
d\tau^{\prime}\sum_{k}\phi_{k}(\tau)\phi_{k}(\tau^{\prime})\sum_{ij\alpha\beta\alpha^{\prime}\beta^{\prime}}\sum_{a=0,z}g_{a}^{2}\
c_{i\alpha}^{\dagger}(\tau)\sigma^{a}_{\alpha\beta}c_{j\beta}(\tau)\
c^{\dagger}_{j\alpha^{\prime}}(\tau^{\prime})\sigma^{a}_{\alpha^{\prime}\beta^{\prime}}c_{i\beta^{\prime}}(\tau^{\prime})\,.$
(S3)

##  S2 Normal state solution

The calculations shown in this section are only in the presence of the spin
channel. From Fig. 2(a), it is evident that the self-energy is diagonal in the
orbital basis. To obtain its spin structure let us consider a general self-
energy $\mathbf{\Sigma}_{f}=\Sigma\sigma^{0}+\vec{\Lambda}\cdot\vec{\sigma}$
incorporating all the spin components. The self-consistent self-energy
equation is then given by

$\begin{split}\mathbf{\Sigma}_{f}(\omega_{n})&=\mathrm{i}g^{2}_{z}T\sum_{m}D(\omega_{n}-\omega_{m})\sigma^{z}[\omega_{m}\sigma^{0}+\mathbf{\Sigma}_{f}(\omega_{m})]^{-1}\sigma^{z}\,,\\\
&=\mathrm{i}g^{2}_{z}T\sum_{m}D(\omega_{n}-\omega_{m})\sigma^{z}\dfrac{[\omega_{m}+\Sigma(\omega_{m})]\sigma^{0}-\vec{\Lambda}(\omega_{m})\cdot\vec{\sigma}}{[\omega_{m}+\Sigma(\omega_{m})]^{2}-|\vec{\Lambda}(\omega_{m})|^{2}}\sigma^{z}\,.\end{split}$

Due to the presence of the bare term $\omega_{m}$ in the $\sigma^{0}$
component, $\Sigma$ must be nonzero. This implies

$(\omega_{m}+\Sigma(\omega_{n}))^{2}>|\vec{\Lambda}|^{2},$ (S4)

and consequently $\Lambda_{z}=0$, as there is a relative sign between the left
and right side of the equation for the $z$-component. In the limit
$\omega\rightarrow 0$, there are two possibilities for $\Lambda_{x}$ and
$\Lambda_{y}$. One possibility is that they decay with with a power-law form
that is subleading compared to $\Sigma$, which we can neglect. Alternatively
we consider the following ansatz, $(\Lambda_{x},\Lambda_{y})=\vec{Q}\Sigma$,
where $|\vec{Q}|<1$ is a constant due to the constraint of Eq. (S4). In this
case,

$\displaystyle\vec{Q}\Sigma(\omega_{n})$
$\displaystyle=-\mathrm{i}g^{2}_{z}T\sum_{m}D(\omega_{n}-\omega_{m})\dfrac{\vec{Q}\Sigma(\omega_{m})}{[\omega_{m}+\Sigma(\omega_{m})]^{2}-|\vec{\Lambda}(\omega_{m})|^{2}}\,,$
(S5) $\displaystyle\Rightarrow\vec{Q}\Sigma(\omega_{n})$
$\displaystyle=\vec{Q}\Sigma(\omega_{n})+\mathrm{i}g^{2}_{z}T\sum_{m}D(\omega_{n}-\omega_{m})\dfrac{\vec{Q}\omega_{m}}{[\omega_{m}+\Sigma(\omega_{m})]^{2}-|\vec{\Lambda}(\omega_{m})|^{2}}\,.$
(S6)

Therefore, $\vec{Q}=0$, leaving the self-energy contribution only from
$\Sigma\sigma^{0}$.

In the normal state $\Phi_{t}=0$. Here, we calculate the $T=0$ case. The self-
consistent Eqs. (3) and (4) are obtained using the Feynman diagram illustrated
in Fig. S1(a).

Figure S1:  Linear gap equation: (a) Feynman diagrams for evaluating the
electronic self-energy (top) and bosonic self-energy (bottom). The bold lines
signify the electronic Green’s functions. (b) Phase diagram at $T=0$. For
large $M$ and $N$ when $\sqrt{2M}/N>1$, the system transits from NFL to the
triplet SC phase. (c) Phase diagram at finite temperature $T$. The transition
temperature to the SC phase $T_{c}\sim g_{z}^{2}$.

Following the approach in Ref. [13], we consider a power-law form for the
electronic self-energy, namely

$\Sigma(\omega)=\mathrm{sgn}(\omega)a_{\eta}^{(1+\eta)/2}|\omega|^{(1-\eta)/2}\,,$
(S7)

i.e., the propagator is dominated by the self-energy at low energy,
$\omega\rightarrow 0$. Consequently, the bosonic self-energy at low energy is
given by

$\begin{split}\delta\Pi(\omega)&=\Pi(\omega)-\Pi(0)\,,\\\
&=-2g_{z}^{2}\dfrac{N}{M}\int\dfrac{d\omega^{\prime}}{2\pi}\Bigg{[}\dfrac{1}{\Sigma(\omega^{\prime}+\omega/2)\Sigma(\omega^{\prime}-\omega/2)}-\dfrac{1}{\Sigma(\omega^{\prime})\Sigma(\omega^{\prime})}\Bigg{]}\,,\\\
&=-2g_{z}^{2}\dfrac{N}{M}\dfrac{\Gamma^{2}(-x)}{2\Gamma(-2x)}\dfrac{1+\sec(\pi
x)}{1/x-2}\dfrac{|\omega|^{\eta}}{2\pi a_{\eta}^{1+\eta}}\,.\end{split}$ (S8)

Here $x=(1-\eta)/2$. The condition required for critical boson,
$\Pi(0)=-m_{0}^{2}$, yields

$a_{\eta}=\dfrac{2}{\pi}\Big{(}\dfrac{g_{z}}{\omega_{0}}\Big{)}^{2}\dfrac{N}{M}\Gamma\Bigg{(}\dfrac{2\eta}{1+\eta}\Bigg{)}\Gamma\Bigg{(}\dfrac{3+\eta}{1+\eta}\Bigg{)}\,.$

The evaluation of $\Pi(0)$ requires the bare terms in $G(\omega)$ to avoid
divergence in the integration [13]. We use the above expressions to calculate
the electronic self-energy and verify the self-consistency:

$\begin{split}&\Sigma(\omega)=g_{z}^{2}\int\dfrac{d\omega^{\prime}}{2\pi}\dfrac{D(\omega-\omega^{\prime})}{\omega+\Sigma(\omega^{\prime})}\,,\\\
\Rightarrow\Sigma(\omega)&=g_{z}^{2}\int\dfrac{d\omega^{\prime}}{2\pi}\dfrac{1}{\delta\Pi(\omega-\omega^{\prime})\Sigma(\omega^{\prime})}\,,\
\ \text{for }\Sigma(\omega)\ \text{and}\ \delta\Pi(\omega)\gg\omega,\\\
&=\dfrac{M}{2N}\dfrac{1/x-2}{1+\sec(\pi
x)}a_{\eta}^{(1+\eta)/2}\mathrm{sgn}(\omega)|\omega|^{(1-\eta)/2}\,.\\\
\end{split}$

Self consistency of $\Sigma(\omega)$ relates $\eta$ with $N/M$ as

$\dfrac{2N}{M}=\dfrac{1/x-2}{1+\sec(\pi
x)}=\dfrac{2\eta}{1-\eta}\dfrac{\tan(\pi\eta/2)}{\tan[\pi(1+\eta)/4]}\,,$ (S9)

and plotted in Fig. S2. Using Eqs. (S8) and (S9), we write
$\delta\Pi(\omega)=g_{z}^{2}b_{\eta}|\omega|^{\eta}/a_{\eta}^{1+\eta}$, where

$b_{\eta}=-\dfrac{1}{4\pi}\dfrac{\Gamma^{2}[(\eta-1)/2]}{\Gamma(\eta-1)}.$

The expression of $b_{\eta}$ given above differs from that in Ref. [13]
because of a typo present in their expression.

Figure S2: $\eta$, which determines the power law of the self-energies, is
shown as a function of $N/M$ using Eq. (S9).

The solution at $T>0$ is a more complicated version of the above, and we omit
it for brevity as the details can be found in Ref. [13]. We note that the
normal state of our model has other regimes of behaviour, namely, a regime of
free fermions and an impurity regime, [13] apart from the SYK regime
considered here, depending on temperature, coupling strength, and the ratio of
$N/M$. However, we focus only on the SYK regime at a very low-temperature
limit, where the self-energies exhibit the power-law forms.

##  S3 Linearized gap solution for $T>0$

The linearized gap equation appears in the main text, Eq. (5). Rewriting this
equation we have,

$\hat{\Phi}^{\alpha\beta}_{ij}(\omega_{n})=\dfrac{g^{2}_{z}}{N}T\sum_{m,\alpha^{\prime},\beta^{\prime}}D(\omega_{n}-\omega_{m})\sigma^{z}_{\alpha\alpha^{\prime}}G(\omega_{m})\hat{\Phi}^{\alpha^{\prime}\beta^{\prime}}_{ji}(\omega_{m})G(-\omega_{m})\sigma^{z}_{\beta\beta^{\prime}}\,.$
(S10)

The only attractive channel is $\sigma_{x}=-i\sigma_{y}\sigma_{z}$, as can be
verified by substituting various Pauli matrices into $\hat{\Phi}$ and summing
over spin indices. Using the relationships
$\sigma^{z}\sigma^{x}\sigma^{z}=-\sigma^{x}$,
$\hat{\Phi}^{\alpha\beta}_{ij}=-\hat{\Phi}^{\alpha\beta}_{ji}$, and
$\hat{\Phi}=\sigma^{x}\Phi_{t}\hat{A}$ for the triplet pairing, where
$\Phi_{t}$ and $\hat{A}$ are defined in the main text, we obtain

$\Phi_{t}(\omega_{n})=\dfrac{1}{N}g_{z}^{2}T\sum_{m}\dfrac{\Phi_{t}(\omega_{m})}{(\omega_{m}+\Sigma(\omega_{m}))^{2}}D(\omega_{n}-\omega_{m})\,.$

We see that the choice of $\hat{A}$ is inconsequential in determining the
above equation from Eq. (S10). This leads us to the conclusion that $T_{c}$
for the triplet SC is the same irrespective of $\hat{A}$. For convenience, we
define the SC gap function as
$\Delta_{t}(\omega_{n})=\omega_{n}\Phi(\omega_{n})/[\omega_{n}+\Sigma(\omega_{n})]$,
and using the above equation we derive:

$\Delta_{t}(\omega_{n})=g_{z}^{2}T\sum_{m}\dfrac{D(\omega_{n}-\omega_{m})}{\omega_{n}+\Sigma(\omega_{n})}\Bigg{[}\dfrac{1}{N}\dfrac{\Delta_{t}(\omega_{m})}{\omega_{m}}-\dfrac{\Delta_{t}(\omega_{n})}{\omega_{n}}\Bigg{]}\,.$
(S11)

Comparing this equation with Eq. (10) in Ref. [13], we see an equivalence
between the parameter $1/N$ and $(1-\alpha)$, where $\alpha$ in Ref. [13]
determines the strength of time-reversal symmetry breaking disorder. For
$\alpha=0$, indicating time-reversal symmetry, random couplings are drawn from
the GOE. Conversely, a nonzero $\alpha$ breaks time-reversal symmetry,
manifesting as an imaginary part proportional to $\alpha$ in the coupling.
Specifically, the random coupling matrix reflects the GUE for $\alpha=1$. We
stress however, that while the gap _equation_ is equivalent for our model and
that of Ref. [13], the gap _structures_ are completely different.

A solution of the gap equation has been shown to exist only for $\alpha$ below
a critical value $\alpha_{c}$ dependent on $\eta$. This suggests a critical
value $N=N_{c}$ in our case as well, below which the SC phase emerges. $N=1$,
indicating $\alpha=0$, corresponds to the largest critical temperature.
Moreover, the rigorous validation of our theory requires accessibility of a
large $N$ limit, hence a large $N_{c}$ for the SC phase. As
$N_{c}\rightarrow\infty$, corresponding to $\alpha\rightarrow 1$, it
associates with $\eta\rightarrow 0$, which is achievable if $M\gg N$,
indicating the presence of a large number of bosons for each fermion. This is
also consistent with the constraint $\sqrt{2M}>N$ for the SC phase to occur.

Eq. (S11) can be formulated as an eigenvalue problem having multiple
eigensolutions. The eigenvectors are denoted as $\Delta_{t}^{(i)}$, where $i$
indicates the eigenvalue index arranged in descending order.
$\Delta_{t}^{(i)}$ associated with the eigenvalue denoted in Eq. (S11) is the
gap solution and these solutions correspond to distinct critical temperature
$T_{c}^{(i)}$. In the main text, the highest critical temperature
$T_{c}^{(0)}$ is referred to as $T_{c}$, which follows $T_{c}\sim g_{z}^{2}$.
Several solutions are shown in Fig. S3. Values of $T_{c}^{(i)}$ are given in
the unit of $m_{0}^{2}/g_{z}^{2}$. In the context of singlet
superconductivity, the significance of solutions with $i>0$ has been discussed
previously in Refs. [13, 39] and references within. They are not relevant for
this work and we do not consider them further here.

All numerical calculations, both here and below, have been performed using the
parameters $N=1$ and $\eta=0.68$. Note that, $\eta\approx 0.68$ corresponds to
$\alpha_{c}\approx 0.63$ as cited in Ref. [13], resulting in $N_{c}=\lfloor
2.7\rfloor=2$.

Figure S3:  The solutions $\Delta_{t}^{(i)}$ as a function of Matsubara
frequency $\omega_{n}$ are determined using Eq. S11 for
$T_{c}^{(0)}=5.372\times 10^{-2}$, $T_{c}^{(1)}=0.891\times 10^{-3}$,
$T_{c}^{(2)}=1.59\times 10^{-5}$ in (a), (b) and (c), respectively. The
remaining parameters, namely $\eta=0.68$, $N=1$, are held constant. For the
computation of $\Delta_{t}^{(i)}$ across a wide range of $\omega_{n}$, the
hybrid frequency technique detailed in Ref. [38] is employed.

##  S4 Nonlinear correction in the gap equation for $g_{0}=0$

We derive the nonlinear equation in the presence of a single coupling term.
First, let us consider the triplet case, where $g_{z}\neq 0$ and $g_{0}=0$. To
visualize the structure of the self-energy and the pairing equation, we
express the entire Green’s function $\hat{\mathcal{G}}$ in the Nambu space. In
this space, the basis is represented by the following set of operators:

$[c_{1\uparrow}(\omega),c_{1\downarrow}(\omega),..,c_{N\uparrow}(\omega),c_{N\downarrow}(\omega);c_{1\uparrow}^{\dagger}(-\omega),c_{1\downarrow}^{\dagger}(-\omega),..,c_{N\uparrow}^{\dagger}(-\omega),c_{N\downarrow}^{\dagger}(-\omega)].$

The matrix $\hat{\mathcal{G}}$ is then expressed as

$\hat{\mathcal{G}}^{-1}=\begin{pmatrix}\hat{G}^{-1}(\omega)&\hat{\Phi}\\\
\hat{\Phi}^{\dagger}&-\hat{G}^{-1}(-\omega)\\\ \end{pmatrix},$ (S12)

where $\hat{G}^{-1}(\omega)=G^{-1}(\omega)\mathds{1}_{2N\times 2N}$,
$G^{-1}(\omega)=i(\omega+\Sigma(\omega))$, and $\hat{\Phi}$ represents the
matrix for all-to-all pairing. $\hat{\mathcal{G}}$ can be further expressed in
terms of the particle-hole ($G_{ph}$) and particle-particle ($G_{pp}$)
components:

$\hat{\mathcal{G}}=\begin{pmatrix}\hat{G}_{ph}&\hat{G}_{pp}\\\
\hat{G}_{hh}&\hat{G}_{hp}\\\ \end{pmatrix}\,.$

First, we evaluate the equations for triplet pairing and the derivation for
singlet pairing readily follows. The self-energy is given by

$\hat{\Sigma}(\omega_{n})=\dfrac{g_{z}^{2}}{N}T\sum_{m}D(\omega_{n}-\omega_{m})\hat{G}_{ph}(\omega_{m})\,.$

As shown above, the normal state Green’s function obeys
$\hat{\Sigma}=\mathrm{i}\Sigma\mathds{1}_{2N\times 2N}$ and
$\hat{G}_{ph}=G\mathds{1}_{2N\times 2N}$. The spin indices are not shown
explicitly in the following discussion.

The pairing equation is given by

$\hat{\Phi}_{ij}(\omega_{n})=\dfrac{g_{z}^{2}}{N}T\sum_{m}D(\omega_{n}-\omega_{m})\sigma^{z}[\hat{G}_{pp}(\omega_{m})]_{ji}\sigma^{z}\,,$
(S13)

where $\hat{G}_{pp}$ is evaluated using blockwise inversion in Eq. (S12):

$\hat{G}_{pp}=\hat{G}(\omega)\hat{\Phi}\hat{G}_{f}(-\omega),$

where we define

$\hat{G}_{f}(-\omega)=[\hat{G}^{-1}(-\omega)+\hat{\Phi}^{\dagger}\hat{G}(\omega)\hat{\Phi}]^{-1}.$

Expanding $\hat{G}_{f}$ for small $\hat{\Phi}$, we get

$\hat{G}_{f}(-\omega)=\hat{G}(-\omega)-\hat{G}(-\omega)\hat{\Phi}^{\dagger}\hat{G}(\omega)\hat{\Phi}\hat{G}(-\omega)+[\hat{G}(-\omega)\hat{\Phi}^{\dagger}\hat{G}(\omega)\hat{\Phi}]^{2}\hat{G}(-\omega)-...$

The approximation $\hat{G}_{f}\approx\hat{G}(-\omega)$ considered in Eq. (S13)
reproduces the linear gap equation. The gap equation, including the higher
order term in $\hat{G}_{f}$, is given as

$\hat{\Phi}_{ij}(\omega_{n})=\dfrac{g_{z}^{2}}{N}T\sum_{m}D(\omega_{n}-\omega_{m})\sigma^{z}\Big{[}G(\omega_{m})\hat{\Phi}G(-\omega_{m})-G(\omega_{m})\hat{\Phi}G(-\omega_{m})\hat{\Phi}^{\dagger}G(\omega_{m})\hat{\Phi}G(-\omega_{m})+...\Big{]}_{ji}\sigma^{z}\,.$
(S14)

Using $\hat{\Phi}=\sigma^{x}\Phi_{t}\hat{A}$ for the triplet pairing, where
$\Phi$ and $\hat{A}$ are defined in the main text, the above equation yields

$\begin{split}\Phi_{t}\hat{A}&=-\dfrac{g_{z}^{2}}{N}T\sum_{m}D(\omega_{n}-\omega_{m})\Big{[}|G(\omega_{m})|^{2}\Phi_{t}\hat{A}-|G(\omega_{m})|^{4}|\Phi_{t}|^{2}\Phi_{t}\hat{A}\hat{A}^{T}\hat{A}+...\Big{]}^{\mathcal{T}_{o}}\,,\\\
&=\dfrac{g_{z}^{2}}{N}T\sum_{m}D(\omega_{n}-\omega_{m})\Big{[}|G(\omega_{m})|^{2}\Phi_{t}\hat{A}-|G(\omega_{m})|^{4}|\Phi_{t}|^{2}\Phi_{t}\hat{A}\hat{A}^{T}\hat{A}+...\Big{]}\,.\\\
\end{split}$

$\mathcal{T}_{o}$ represents the transpose operator in the orbital basis
($ij\rightarrow ji$) and does not affect the spin, hence
$A^{\mathcal{T}_{o}}=-A$. It is noteworthy that for the singlet case, $S$ does
not flip sign under $\mathcal{T}_{o}$. Consequently, the gap shows attraction
in the $\mathrm{i}\sigma^{y}$ channel for the $\sigma^{o}$ coupling in the
Yukawa term. Multiplying both sides of the above equation with $\hat{A}^{T}$
and taking trace, we obtain

$\begin{split}\Phi_{t}&\approx\dfrac{g_{z}^{2}}{N}T\sum_{m}D(\omega_{n}-\omega_{m})\Big{[}|G(\omega_{m})|^{2}\Phi_{t}-|G(\omega_{m})|^{4}|\Phi_{t}|^{2}\Phi_{t}\underbrace{\dfrac{\operatorname{Tr}(\hat{A}^{T}\hat{A}\hat{A}^{T}\hat{A})}{\operatorname{Tr}(\hat{A}^{T}\hat{A})}}_{\approx
N^{2}/3}\Big{]}\,,\\\ \Rightarrow\Phi_{t}(\omega_{n})&\approx
g_{z}^{2}T\sum_{\omega_{m}}D(\omega_{n}-\omega_{m})\Big{[}\dfrac{1}{N}\dfrac{1}{|\omega_{m}+\Sigma(\omega_{m})|^{2}}-\dfrac{N}{3}\dfrac{|\Phi_{t}(\omega_{m})|^{2}}{|\omega_{m}+\Sigma(\omega_{m})|^{4}}\Big{]}\Phi_{t}(\omega_{m})\,.\end{split}$

For the linear term, $\operatorname{Tr}(\hat{A}^{T}\hat{A})$ cancels on both
sides of the equation. However,
$\operatorname{Tr}(\hat{A}^{T}\hat{A}\hat{A}^{T}\hat{A})/\operatorname{Tr}(\hat{A}^{T}\hat{A})$
depends on $A$, and the consequences of this dependence here and in the mixed
state are not studied in this work. For the results shown here, we consider
$A$ as given in the main text.

$\Delta_{t}$ from the above equation follows,

$\Delta_{t}(\omega_{n})=g_{z}^{2}T\sum_{m}\dfrac{D(\omega_{n}-\omega_{m})}{\omega_{m}+\Sigma(\omega_{m})}\Bigg{[}\dfrac{1}{N}\dfrac{\Delta_{t}(\omega_{m})}{\omega_{m}}-\dfrac{N}{3}\dfrac{\Delta^{3}_{t}(\omega_{m})}{\omega^{3}_{m}}-\dfrac{\Delta_{t}(\omega_{n})}{\omega_{n}}\Bigg{]}\,.$
(S15)

A similar equation can be obtained for the singlet gap function, where
$g_{z}=0$ and $g_{0}\neq 0$, and the equation is given by

$\Delta_{s}(\omega_{n})=g_{0}^{2}T\sum_{m}\dfrac{D(\omega_{n}-\omega_{m})}{\omega_{m}+\Sigma(\omega_{m})}\Bigg{[}\dfrac{1}{N}\dfrac{\Delta_{s}(\omega_{m})}{\omega_{m}}-N\dfrac{\Delta^{3}_{s}(\omega_{m})}{\omega^{3}_{m}}-\dfrac{\Delta_{s}(\omega_{n})}{\omega_{n}}\Bigg{]}\,.$
(S16)

The triplet and singlet gap equations differ by a factor of $1/3$ in the
second term on the right-hand side, arising from
$\operatorname{Tr}(\hat{A}^{T}\hat{A}\hat{A}^{T}\hat{A})\approx\operatorname{Tr}(\hat{S}^{4})/3\sim
N^{4}$ for large $N$. Consequently, as the temperature decreases below
$T_{c}$, the triplet gap function grows faster than the singlet one. Solutions
of the two equations are related as

$\Delta_{s}(\omega_{n})=\dfrac{1}{\sqrt{3}}\Delta_{t}(\omega_{n})\,,$ (S17)

with the corresponding interchange between $g_{0}^{2}$ and $g_{z}^{2}$.

##  S5 The gap equation in the presence of both $\sigma^{z}$ and $\sigma^{0}$

Here, we derive the gap equation in the presence of both coupling terms.
Following the derivation in Sec. S4, the gap equation is expressed as

$\hat{\Phi}=\sum_{m}\dfrac{T}{N}D(\omega_{n}-\omega_{m})|G(\omega_{m})|^{-2}\sum_{a=z,0}g_{\lambda}^{2}\sigma^{a}\Big{[}\hat{\Phi}-|G(\omega_{m})|^{-2}\hat{\Phi}\hat{\Phi}^{\dagger}\hat{\Phi}\Big{]}^{T_{o}}\sigma^{a}\,,$
(S18)

where
$\hat{\Phi}=\Phi_{t}\hat{A}\sigma^{x}+\Phi_{s}\hat{S}\mathrm{i}\sigma^{y}$, as
described in the main text. We define
$\tilde{\Sigma}(\omega_{m})=\omega_{m}+\Sigma(\omega_{m})$. To obtain the
equation for the triplet pairing function, we take trace after multiplying
both sides of the above equation by $\hat{A}^{T}$. This yields

$\Phi_{t}(\omega_{n})=(g_{z}^{2}-g_{0}^{2})T\sum_{m}D(\omega_{n}-\omega_{m})\Bigg{[}\dfrac{1}{N}\dfrac{\Phi_{t}}{|\tilde{\Sigma}(\omega_{m})|^{2}}-\dfrac{1}{|\tilde{\Sigma}(\omega_{m})|^{4}}\Big{(}\dfrac{N}{3}|\Phi_{t}|^{2}\Phi_{t}+\dfrac{2N}{3}|\Phi_{s}|^{2}\Phi_{t}+\dfrac{2}{3}(\Phi_{s})^{2}\Phi^{*}_{t}\Big{)}\Bigg{]}\,.$
(S19)

Similarly, for the singlet gap, we multiply both sides of Eq. (S18) by
$\hat{S}^{T}$ and take trace, resulting in

$\Phi_{s}(\omega_{n})=(g_{0}^{2}-g_{z}^{2})T\sum_{m}D(\omega_{n}-\omega_{m})\Bigg{[}\dfrac{1}{N}\dfrac{\Phi_{s}}{|\tilde{\Sigma}(\omega_{m})|^{2}}-\dfrac{1}{|\tilde{\Sigma}(\omega_{m})|^{4}}\Big{(}N|\Phi_{s}|^{2}\Phi_{s}+\dfrac{2N}{3}|\Phi_{t}|^{2}\Phi_{s}+\dfrac{2}{3}(\Phi_{t})^{2}\Phi^{*}_{s}\Big{)}\Bigg{]}\,.$
(S20)

In the derivation, we use the relations:
$\operatorname{Tr}[\hat{A}^{T}\hat{A}]\sim N^{2}$,
$\operatorname{Tr}[\hat{A}^{T}\hat{A}\hat{A}^{T}\hat{A}]\approx N^{4}/3$,
$\operatorname{Tr}[\hat{A}^{T}\hat{A}\hat{S}\hat{S}]\approx N^{4}/3$,
$\operatorname{Tr}[\hat{A}^{T}\hat{S}\hat{A}^{T}\hat{S}]\approx 2N^{3}/3$,
$\operatorname{Tr}[\hat{S}^{T}\hat{S}]\approx N^{2}$,
$\operatorname{Tr}[\hat{S}\hat{A}]=0$, and
$\operatorname{Tr}[\hat{S}^{4}]\approx N^{4}$ for the large-$N$ approximation.

It is evident that the singlet and triplet gap equations are equivalent,
except for the prefactor, where the difference in disorder strengths exhibits
opposite signs in the two cases, and the factor of $1/3$ present in the first
nonlinear term only in the triplet equation. Upon rescaling
$\Phi_{s(t)}\rightarrow\Phi_{s(t)}/N$ in Eqs. (S19) and (S20), we recover Eq.
(6) in the main text.

The last terms in Eqs. (S19) and (S20) serve to fix the relative phase between
$\Phi_{s}$ and $\Phi_{t}$ when both pairings are present in the SC state. The
phase difference between $\Phi_{t}$ and $\Phi_{s}$ is either $0$ or $\pi$. To
see this, consider the system in a singlet pairing state, i.e.
$g_{0}^{2}>g_{z}^{2}$, with temperature just below to $T_{c}$. Here,
$\Phi_{s}$ can be assumed to be real. Upon lowering the temperature, a nonzero
triplet pairing $\Phi_{t}$ can be induced. The last term in Eq. (S19) acts as
an attractive potential in the triplet channel if $(\Phi_{s})^{2}\Phi^{*}_{t}$
remains positive. Since $\Phi_{s}$ is real, this implies $\Phi_{t}$ is real as
well, fixing the relative phase up to $\pi$. Note that the induced phase
$\Phi_{t}$ also has an out-of-phase feedback effect on $\Phi_{s}$, but this is
a higher-order effect that we neglect here. In principle, deep in the mixed
phase, the feedback and induced terms may generate a nontrivial phase relation
between $\Phi_{s}$ and $\Phi_{t}$.

###  S5.1 Linearized equation and evaluation of $T_{c}$

While at the linear approximation, the two gap equations are completely
decoupled, the results in the presence of either of the couplings [Sec. S3]
cannot be extrapolated here just by replacing, e.g., $g_{z}^{2}$ with
$g_{z}^{2}-g_{0}^{2}$ for the triplet gap function. To see this, we consider
Eq. (S19) which after linearization yields

$\Phi_{t}(\omega_{n})=(g_{z}^{2}-g_{0}^{2})\dfrac{1}{N}T\sum_{m}D(\omega_{n}-\omega_{m})\dfrac{\Phi_{t}}{|\tilde{\Sigma}(\omega_{m})|^{2}}\,.$
(S21)

The effective coupling in the gap equation is $\lambda=g_{z}^{2}-g_{0}^{2}$,
whereas that for the self-energies is
$\lambda_{{}_{(+)}}=g_{z}^{2}+g_{0}^{2}$. Consequently, the corresponding
equation for $\Delta_{t}$ is expressed as

$\Delta_{t}(\omega_{n})=\lambda_{{}_{(+)}}T\sum_{m}\dfrac{D(\omega_{n}-\omega_{m})}{\omega_{m}+\Sigma(\omega_{m})}\Bigg{[}\dfrac{\lambda}{\lambda_{{}_{(+)}}}\dfrac{1}{N}\dfrac{\Delta_{t}(\omega_{m})}{\omega_{m}}-\dfrac{\Delta_{t}(\omega_{n})}{\omega_{n}}\Bigg{]}\,.$
(S22)

We define

$A_{nm}=\dfrac{\lambda_{{}_{(+)}}}{\pi}\dfrac{D(\omega_{n}-\omega_{m})}{\tilde{\Sigma}(\omega_{m})}\,,$

rescale $T\rightarrow Tm_{0}^{2}/\lambda_{{}_{(+)}}$ [13], and rewrite Eq.
(S22) as

$\Rightarrow.\sum_{m}\dfrac{1}{(1+Q_{n})}\dfrac{A_{nm}}{2m+1}\Delta_{m}=N\dfrac{\lambda_{{}_{(+)}}}{\lambda}\Delta_{n}\,,$
(S23)

where, $Q_{n}=\sum_{m}A_{nm}/(2n+1)$. We vary $T$ to satisfy the above
equation for different values of $0<(\lambda/\lambda_{{}_{(+)}})<1$ and the
temperature corresponds to $T_{c}$, while $N$ is held at a constant value of
1. Without loss of generality, we consider $\lambda_{{}_{(+)}}=1$ and
$0<\lambda<1$. As discussed previously, a comparison with Ref. [13] ensures a
critical value $\lambda_{c}$ for a given $N$ and $\eta$, above which the SC
phase emerges.

Far away from $\lambda_{c}$ the difference between $\lambda\sim g_{z}^{2}$ and
$\lambda_{{}_{(+)}}\sim\lambda$ can be neglected, and $T_{c}\sim\lambda\sim
g_{z}^{2}$. Near the transition, the critical behavior is
$T_{c}\sim\mathrm{exp}\big{[}-\pi/\sqrt{\lambda-\lambda_{c}}\big{]}$ [13]. For
the parameter $N=1$ and $\eta=0.68$, we find $\lambda_{c}\approx 0.37$.

Note that at the linear level, the triplet and singlet gap equations are dual
to each other. Therefore, a similar dependency of $T_{c}$ for the singlet
phase holds true as a function of $|\lambda|=(g_{0}^{2}-g_{z}^{2})$ as well.

###  S5.2 Gap function below $T_{c}$

We aim to obtain $\Delta_{t}$ for temperatures lesser than but close to
$T_{c}$ and $g_{z}^{2}>g_{0}^{2}$. In this region, we assume the singlet gap
function $\Phi_{s}=\Delta_{s}=0$ and verify this assumption later. By using
Eq. (S19), we derive the nonlinear equation for $\Delta_{t}$ as

$\tilde{\Delta}_{t}(\omega_{n})=\dfrac{1}{(1+Q_{n})}\dfrac{\lambda}{\lambda_{{}_{(+)}}}\sum_{m}\Big{[}\dfrac{1}{N}\dfrac{A_{nm}}{2m+1}\tilde{\Delta}_{t}(\omega_{m})-\dfrac{N}{3\pi^{2}}\sum_{m}\dfrac{A_{nm}}{(2m+1)^{3}}\tilde{\Delta}_{t}^{3}(\omega_{m})\Big{]}\,,$
(S24)

where $\tilde{\Delta}_{t}=\Delta_{t}/T$ and we set $\lambda_{{}_{(+)}}=1$ as
before. To solve Eq. (S24), we use the linear solution
$\Delta_{t}^{(0)}(\omega_{n})$ as the initial input and iterate the equation
until convergence is achieved. Several gap solutions for different $T$ and
$\lambda$ are shown in Fig. S4. The solution of $\Delta_{s}$ for
$g^{2}_{0}>g^{2}_{z}$ can be inferred using Eq. (S17).

Figure S4: The solutions of $\Delta_{t}(\omega_{n})$ obtained from Eq. (S24)
for several temperatures $T<T_{c}$ and coupling $\lambda$.  Figure S5: (a) and
(b): Various solutions corresponding to different temperatures shown in Figs.
S4(b) and (c) are collapsed onto their respective curves corresponding to
$T=0.1T_{c}$ using $f(T)$ as a fitting parameter. (c) The empty markers show
$f(T)$ obtained from the collapse. The dashed lines represent the fit of
$f(T)\simeq p_{1}T/T_{c}+p_{2}$ for low temperatures with $p_{1}$ and $p_{2}$
given in Tab. 1.

Our Ginzburg-Landau expression is formally valid only for $|T_{c}-T|\ll T$.
However, in order to obtain the qualitative form of the phase diagram, we
solve it for low $T$ as well. To derive an expression for $\Delta_{t}$ at low
$T$, the finite temperature solutions, such as those represented in Fig. S4
are scaled as $\bar{\Delta}_{t}(T)=\Delta_{t}/\sqrt{T_{c}-T}$ and then
collapsed onto the lowest temperature curve obtained numerically, denoted as
$\bar{\Delta}_{ref}=\bar{\Delta}_{t}(T=0.1Tc)$ for each value of $\lambda$.
The expression used is $\bar{\Delta}_{t}(\omega_{n})\sim
f(T,\lambda)\bar{\Delta}_{ref}$, where $f(T,\lambda)$ is a fitting parameter.
The fittings are depicted in Figs. S5(a) and (b) for $\lambda=0.80$ and
$0.70$, respectively. Figure S5(c) illustrates the functional form of
$f(T,\lambda)$, which agrees well with the following analytic form:

$f(T,\lambda)\simeq p_{1}(\lambda)\dfrac{T}{T_{c}}+p_{2}(\lambda)\,,$ (S25)

for low temperatures with $p_{1}$ and $p_{2}$ as the fitting parameters.
Values of $p_{1}$ and $p_{2}$ are given in Tab. 1. An examination of Eq. (S25)
as well as Figs. S4 and S5 reveals that the gap _decreases_ at low
temperatures at the lowest $\omega_{n}$. This is apparently an artifact of the
singular form of the expansion at low frequencies, see e.g. Eq. (S15). Hence,
the quantitative form of the phase diagram should not be taken too seriously.

$\lambda$ |  $p_{1}$ |  $p_{2}$
---|---|---
0.55 |  1.24 |  0.91
0.60 |  1.03 |  0.91
0.65 |  1.15 |  0.89
0.70 |  1.30 |  0.87
0.75 |  1.46 |  0.85
0.80 |  1.64 |  0.83
0.85 |  1.84 |  0.80
0.90 |  2.08 |  0.77
0.95 |  2.34 |  0.74
1.0 |  2.64 |  0.71
Table 1: Fitting parameters for $f(T,\lambda)$ expressed in Eq. (S25) at
several $\lambda$.

##  S6 Calculation of the critical temperature for the mixed state

Let us consider the case mentioned in the main text, i.e.,
$g_{0}^{2}<g_{z}^{2}$ and $T$ is just below $T_{c}$, which results in
$\Phi_{t}\neq 0$ and $\Phi_{s}=0$. The SC phase consists purely of triplet
pairing. Now consider lowering the temperature further. We anticipate that
$\Phi_{s}$ will be small near the onset of the mixed state and neglect the
$\sim\Phi_{s}^{3}$ term in Eq. (S20), which leads to a linear equation for
$\Phi_{s}$:

$\Phi_{s}(\omega_{n})=-\lambda\dfrac{T}{N}\sum_{m}D(\omega_{n}-\omega_{m})\Bigg{[}\dfrac{\Phi_{s}}{|\tilde{\Sigma}(\omega_{m})|^{2}}-\dfrac{1}{|\tilde{\Sigma}(\omega_{m})|^{4}}\dfrac{2}{3}|\Phi_{t}|^{2}\Phi_{s}\Bigg{]}\,.$
(S26)

Using the definition of $\Delta_{s}$, we obtain

$\displaystyle\tilde{\Sigma}(\omega_{n})\Delta_{s}(\omega_{n})=\omega_{n}\Phi_{s}(\omega_{n})\,,$
$\displaystyle\Rightarrow\Delta_{s}(\omega_{n})=\lambda_{{}_{(+)}}T\sum_{m}\dfrac{D(\omega_{n}-\omega_{m})}{\tilde{\Sigma}(\omega_{m})}\Bigg{[}\dfrac{\lambda}{\lambda_{{}_{(+)}}}\dfrac{1}{N}\Bigg{(}\dfrac{2}{3}\dfrac{\Delta_{t}^{2}(\omega_{m})}{\omega_{m}^{2}}-1\Bigg{)}\dfrac{\Delta_{s}(\omega_{m})}{\omega_{m}}-\dfrac{\Delta_{s}(\omega_{n})}{\omega_{n}}\Bigg{]}\,.$
(S27)

This can be written as

$\dfrac{\lambda}{\lambda_{{}_{(+)}}}\dfrac{1}{N}\sum_{m}\dfrac{A_{nm}}{(1+Q_{n})(2m+1)}\Bigg{(}\dfrac{2}{3}\dfrac{\Delta_{t}^{2}(\omega_{m})}{\omega_{m}^{2}}-1\Bigg{)}\Delta_{s}(\omega_{m})=\Delta_{s}(\omega_{n})\,,$
$\Rightarrow\dfrac{\lambda}{\lambda_{{}_{(+)}}}\dfrac{1}{N}\sum_{m}\dfrac{A_{nm}}{(1+Q_{n})(2m+1)}\Bigg{(}\dfrac{2}{3}\dfrac{\tilde{\Delta}_{t}^{2}(\omega_{m})}{(2m+1)^{2}\pi^{2}}-1\Bigg{)}\tilde{\Delta}_{s}(\omega_{m})=\tilde{\Delta}_{s}(\omega_{n})\,.$
(S28)

Using the finite temperature form of $\Delta_{t}$ discussed in the previous
section, we solve Eq. (S28) as an eigenvalue problem. The critical temperature
$T_{t\rightarrow m}$ for the mixed state is obtained as a function of
$\lambda$, as shown in Fig. S6(b). We can estimate the critical temperature
using Eq. (S30). At low energy $\Delta_{t}\sim T^{(1+\eta)/2}\Phi_{t}$ and
assuming the kernel to be zero at the critical temperature, we obtain

$T_{t\rightarrow
m}\sim\mathrm{exp}\Big{[}-\dfrac{\pi}{\sqrt{\lambda-\lambda_{m_{1}}}}\dfrac{4}{3+\eta}\Big{]}\,,$
(S29)

where $\Delta_{t}(T=0)\sim T_{c}$. The numerical data points fit well with the
above form as shown in Fig. S6(b) for $\lambda_{m_{1}}\approx\lambda_{c}$.

Similarly, if we begin from a pure singlet state, i.e., $g_{0}^{2}>g_{z}^{2}$,
and lower the temperature, due to the slower growth of the singlet gap
function [Eq. (S17)], the critical temperature for the transition to the mixed
state $T_{s\rightarrow m}$ will be smaller than $T_{t\rightarrow m}$. This is
verified numerically and shown in Fig. S6(a), where a similar functional form
as that for $T_{t\rightarrow m}$ is used to fit $T_{s\rightarrow m}$. The
critical value of $g_{0}^{2}-g_{z}^{2}$ is approximately
$\lambda_{m_{2}}\approx-\lambda_{m_{1}}$ for the mixed state to appear [Fig.
S6(a)]. Our numerical analysis and the value of
$\lambda_{m_{1}}\approx-\lambda_{m_{2}}\approx\lambda_{c}$ suggest that the
ground state of the superconductor always consists of mixed pairings. We now
discuss this issue and show that the reality is a bit more complex.

Figure S6: Numerically evaluated critical temperatures $T_{s\rightarrow m}$
and $T_{t\rightarrow m}$, denoted by squares, are shown in (a) and (b),
respectively. We see that $T_{s\rightarrow m}<T_{t\rightarrow m}$. The solid
lines in (a) and (b) illustrate fittings to a curve of the form shown in Eq.
(S29).

Consider the behavior of the kernel in Eq. (S26) deep in the SC phase and away
from $\lambda_{c}$. We replace the sum with integration with the UV cutoff
$\omega_{0}$, where $\omega_{0}\sim\Sigma(\omega_{0})$. At $\omega_{n}=\pi T$,
this gives

$\Phi_{s}(T)\sim
A\lambda\int_{T}^{\omega_{0}}d\omega\Bigg{[}\dfrac{1}{B\omega^{(3-\eta)/2}}|\Phi_{t}|^{2}-\dfrac{1}{\omega}\Bigg{]}\Phi_{s}\,,$
(S30)

where all $\eta$-dependent constants are incorporated into $A$ and $B$.
Assuming $|\Phi_{t}|^{2}$ to be a constant, the kernel simplifies to

$\sim\dfrac{2|\Phi_{t}|^{2}}{(1-\eta)BT^{(1-\eta)/2}}-\ln\dfrac{\omega_{0}}{T},$

which diverges as $T\rightarrow 0$. The implication is that for low enough
temperatures, the nonlinear term always overtakes the linear term. As a
result, in the SC phase far away from the QCP, at sufficiently low
temperatures, the pure pairing state always transits to a mixed state
(although formally the Ginzburg-Landau treatment also breaks down at the same
scale). However, near $\lambda_{c}$, the divergence at $T\rightarrow 0$
renders our entire Ginzburg-Landau approach invalid, as can be seen by e.g.
the fact that $\Phi_{t}$ vanishes at low frequencies in the $T\to 0$ limit,
see Fig. S5. Along the $T=0$ axis, since the pure singlet/triplet pairing
phase undergoes a second-order phase transition at $\lambda_{c}$, the gap
function, being infinitesimally small close to $\lambda_{c}$, cannot dominate
over the repulsive linear term for the induced phase. Consequently, the
induced gap function remains repulsive and becomes attractive only when at
higher coupling strength the nonlinear term dominates, thereby suggesting the
existence of the second QCP. To pinpoint $\lambda$ corresponding to this QCP,
the full nonlinear equation needs to be solved, which is beyond the scope of
the current endeavor.

###  S6.1 Onsite pairing

Up to this point, we have not considered intra-orbital pairing in our
discussion. To incorporate it, we extend the pairing function as follows:

$\hat{\Phi}=\Phi_{t}\hat{A}\sigma^{x}+[\Phi_{nl}\hat{S}+\Phi_{l}\hat{\mathds{I}}]\mathrm{i}\sigma^{y}\,,$
(S31)

where diag$(S)=\\{0\\}$ as before and $\hat{\mathds{I}}$ represents the
identity matrix. To distinguish between the local and nonlocal singlet
pairing, we introduce $\Phi_{l}$ for the local pairing and replace the
notation $\Phi_{s}$ with $\Phi_{nl}$ for the nonlocal pairing.

Following the previously mentioned method, we evaluate the equations for the
gap functions as:

$\displaystyle\Phi_{t}(\omega_{n})=(g_{z}^{2}-g_{0}^{2})\dfrac{T}{N}\sum_{m}\dfrac{D(\omega_{n}-\omega_{m})}{|\tilde{\Sigma}(\omega_{m})|^{2}}\Bigg{[}\Phi_{t}-\dfrac{\Phi_{t}}{|\tilde{\Sigma}(\omega_{m})|^{2}}\Big{\\{}\dfrac{1}{3}|\Phi_{t}|^{2}+\dfrac{2}{3}|\Phi_{nl}|^{2}+2|\Phi_{l}|^{2}+\dfrac{2}{3}(\Phi_{nl}^{*}\Phi_{l}+\Phi_{nl}\Phi_{l}^{*})\Big{\\}}$
$\displaystyle-\dfrac{\Phi_{t}^{*}}{|\tilde{\Sigma}(\omega_{m})|^{2}}\Big{\\{}\dfrac{2}{3N}(\Phi_{nl})^{2}-(\Phi_{l})^{2}-\dfrac{2}{3}\Phi_{nl}\Phi_{l}\Big{\\}}\Big{)}\Bigg{]}\,$
(S32)
$\displaystyle\Phi_{nl}(\omega_{n})=(g_{0}^{2}-g_{z}^{2})\dfrac{T}{N}\sum_{m}\dfrac{D(\omega_{n}-\omega_{m})}{|\tilde{\Sigma}(\omega_{m})|^{2}}\Bigg{[}\Phi_{nl}-\dfrac{\Phi_{nl}}{|\tilde{\Sigma}(\omega_{m})|^{2}}\Big{\\{}|\Phi_{nl}|^{2}+\dfrac{2}{3}|\Phi_{t}|^{2}+2|\Phi_{l}|^{2}+\Phi_{nl}^{*}\Phi_{l}+\Phi_{nl}\Phi_{l}^{*}\Big{\\}}$
$\displaystyle-\dfrac{\Phi_{nl}^{*}}{|\tilde{\Sigma}(\omega_{m})|^{2}}\Big{\\{}\dfrac{2}{3N}(\Phi_{t})^{2}+(\Phi_{l})^{2}+\Phi_{nl}\Phi_{l}\Big{\\}}\Bigg{]}\,$
(S33)

and

$\displaystyle\Phi_{l}(\omega_{n})=(g_{0}^{2}-g_{z}^{2})\dfrac{T}{N}\sum_{m}\dfrac{D(\omega_{n}-\omega_{m})}{|\tilde{\Sigma}(\omega_{m})|^{2}}\Bigg{[}\Phi_{l}-\dfrac{\Phi_{l}}{|\tilde{\Sigma}(\omega_{m})|^{2}}\Big{\\{}|\Phi_{l}|^{2}+\dfrac{2}{N}|\Phi_{t}|^{2}+\dfrac{2}{N}|\Phi_{nl}|^{2}\Big{\\}}$
$\displaystyle-\dfrac{\Phi_{l}^{*}}{|\tilde{\Sigma}(\omega_{m})|^{2}}\Big{\\{}-\dfrac{1}{N}(\Phi_{t})^{2}+\dfrac{1}{N}(\Phi_{nl})^{2}\Big{\\}}\Bigg{]}.$
(S34)

In these equations, we rescale the triplet and nonlocal singlet pairing
functions: $\Phi_{t}\rightarrow\Phi_{t}/N$, $\Phi_{nl}\rightarrow\Phi_{nl}/N$.
For complex order parameters, the equations suggest complicated relative
phases among them, which are left for future exploration.