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 protected] Jonathan Ruhman Department of Physics, Bar Ilan University, Ramat Gan 5290002, Israel [email protected] Avraham Klein Physics Department, Ariel University, Ariel 40700, Israel

(May 1, 2024)

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].

Refer to caption — 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.

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}$ oder $\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.

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 Sr₂RuO₄ 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-TaS₂, 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 UTe₂, 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).

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.

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$ .

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$ oder $\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).

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.

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.