Vestigial singlet pairing in a fluctuating magnetic triplet superconductor and its implications for graphene superlattices

Stacking and twisting graphene layers allows to create and control a two-dimensional electron liquid with strong correlations. Experiments indicate that these systems exhibit strong tendencies towards both magnetism and triplet superconductivity. Motivated by this phenomenology, we study a 2D model of fluctuating triplet pairing and spin magnetism. Individually, their respective order parameters, d and N, cannot order at finite temperature. Nonetheless, the model exhibits a variety of vestigial phases, including charge-4e superconductivity and broken time-reversal symmetry. Our main focus is on a phase characterized by finite d ⋅ N, which has the same symmetries as the BCS state, a Meissner effect, and metastable supercurrents, yet rather different spectral properties: most notably, the suppression of the electronic density of states at the Fermi level can resemble that of either a fully gapped or nodal superconductor, depending on parameters. This provides a possible explanation for recent tunneling experiments in the superconducting phase of graphene moiré systems.

Strongly correlated systems often exhibit complex phase diagrams with multiple phases, characterized by long-range or quasi-long-range order (QLRO) of different order parameters.Aside from phase competition as a possible origin, a rich set of phases might also be understood as different manifestations of an underlying primary order-a concept often referred to as "intertwined orders" [1].For instance, thermal or quantum fluctuations can disorder a primary order parameter, while higher-order composite order parameters can still survive.An example of such a "vestigial phase" [2,3], is the charge-4e superconducting state that emerges when a charge-2e pair density wave order parameter, ∆ Q , itself vanishes, yet ∆ Q ∆ −Q does not [4]; this and other forms of charge-4e superconductivity have attracted a lot of attention [5][6][7][8][9][10][11][12][13][14][15][16][17][18], in particular, as a result of recent experiments [19,20].
Another exciting recent development is the emergence of twisted graphene moiré superlattices as versatile playgrounds for strongly correlated physics [21,22].These systems display a variety of different phases such as nematic [23][24][25] and density-wave order [26][27][28], different forms of magnetism [29][30][31][32][33], and, possibly unconventional [34,35], superconductivity [36]; magnetism and superconductivity appear in the same density range [34,35,[37][38][39][40][41] and recent experiments [33,42] demonstrate that they can coexist microscopically.Motivated by these observations, we here study the case of two primary order parameters: a fully gapped spin-triplet superconductor (d) and, in line with the conclusions of [41,43], magnetic order (N ) with antiparallel spins in the two valleys.At finite temperature, T > 0, it must hold d = N = 0, in two-dimensions (2D).However, there are several different vestigial phases, see Fig. 1(a), characterized by the composite order parameters φ dd = d • d, φ dN = d•N , and φ ddN = i(d † ×d)•N .These include not only a charge-4e superconductor [44,45], see Fig. 1(b), but also a charge-2e state, which has the same symme-tries as and is, hence, adiabatically connected to the BCS state.However, it should primarily be thought as a condensate of three electrons and a hole, see Fig. 1(c), or, more formally, QLRO of φ dN .We develop a theory for this state and study its spectral properties at finite T , which are rather different from those of the BCS state.Depending on T and φ dN , we obtain a low-energy suppression of the density of states (DOS) similar to a fully gaped or nodal state.This could provide an alternative explanation [43,46,47] to the tunneling data of [34,35], which does not require any momentum dependence in the superconducting order parameter.Matsubara frequencies and 2D momentum, they couple as Note that N couples anti-ferromagnetically in the two valleys; while the ferromagnetic case-with τ 0 instead of τ z in the first term of S c above-can be studied similarly, we focus on antiferromagnetism not only for concreteness here but also because recent microwave experiments [41] and a systematic analysis [43] of multiple other experiments on graphene moiré systems favor this scenario.
The bare dynamics of d and N is governed by We take the susceptibilitites to be χ µ (q) = χ 0 µ /(r µ +Ω 2 n + v 2 µ q 2 ), µ = N, d, where q = (iΩ n , q) and Ω n are bosonic Matsubara frequencies.The nature of the phase realized in the system depends crucially on the interactions between the bosonic fields.Up to quartic order, the local terms allowed by the symmetries listed in Table I can be written as Finally, the bare electronic action is given by where we already used that the band structures in the two valleys are related by time-reversal Θ.
Mean-field and possible phases.-Toprobe the possible phases, we start with a mean-field analysis with respect to d and N .Absorbing the impact of the coupling to the electrons [48] by a redefinition of the parameters of V , we obtain the four distinct zero-temperature phases labeled (A), (B 1,2 ), and (C) in Fig. 1(a), where we assumed that both d and N are non-zero and homogeneous.Using ê1,2,3 ∈ R 3 to denote orthogonal unit vectors, we have N = N 0 ê1 and d = d 0 e iα ê2 in phase (A), which breaks SO(3) completely, while Θ is preserved (in any gaugeinvariant observable); as for any phase with N = 0, C 2z is broken.In phase (B 1 ), N and d are aligned; we, thus, obtain a residual spin-rotation symmetry SO(2) along that direction and Θ is preserved too.Beyond a critical value of b 2 , an additional component with relative phase π/2 emerges in d, defining phase (B 2 ) where N = N 0 ê1 and d = d 0 e iα (ê 1 + iηê 2 ), with 0 < η < 1; this is a distinct phase as η = 0 breaks both the residual SO(2) spin symmetry and Θ.Finally, phase (C) is characterized by N = N 0 ê1 and d = d 0 e iα (ê 2 + iê 3 ).Consequently, Θ is also broken but a residual SO(2) spin-symmetry remains.
Importantly, d , N = 0 is only possible and, thus, our discussion of symmetries is only valid for T = 0 in TABLE I: Relevant symmetries g and their action on the field operators.Here Rϕ is the orthogonal matrix obeying e −iϕ•s se iϕ•s = R(ϕ)s.All symmetries are linear except for Θ which is anti-linear.
To analyze the resulting vestigial phases at finite T , where SO(3) spin-rotation symmetry is preserved and d = N = 0, it is convenient to define the following composite order parameters • N , with symmetry properties listed in Table I.Crucially, all of them transform trivially under SO(3) spin-rotations and, hence, can exhibit long-range (in case of the last one) or QLRO (in case of the former two) at finite T .We indicate this in Fig. 1(a) for the different phases.This immediately tells us that, in spite of d = 0, phase (A) transitions for finite T into state where φ dd has QLRO and, thus, constitutes a charge-4e superconductor (as φ dN = 0), which does not break C 2z or Θ (as φ ddN = 0); intuitively, one can think of this state as a condensate of four electrons forming a spinsinglet out of two triplets, see Fig. 1(b).At finite T , (B 1 ) and (B 2 ) will both preserve all normal-state symmetries and become the same phase, which we denote by (B) in the following.It is characterized by QLRO not only in φ dd but also in φ dN ; as the latter has charge 2e, it is a charge-2e superconductor and adiabatically connected to the conventional BCS state.Nonetheless, in our current description, this state should rather be thought of as the condensation of three electrons and a hole, see Fig. 1(c), consisting of a pair of electrons in a triplet state forming a singlet with a spin-1 particle-hole excitation.In fact, we will see below that it exhibits spectral properties rather different from those of the BCS state at finite T .Finally, while phase (C) does not exhibit any vestigial pairing at T > 0, it will have long-range order in φ ddN and, as such, continues to break both C 2z and Θ.
Theory for phase (B).-As c 1 < 0 is found when the coefficients in V are computed by integrating out electrons [48], we next focus on phase (B).To obtain an efficient description of this phase that properly captures the preserved SO(3) symmetry at finite temperature, we first decouple the four terms in V using four Hubbard-Stratonovich fields, ψ d for d † d, ψ N for N 2 , φ d for d • d, and φ dN for d • N .We treat them on the saddle-point level, which becomes exact in the limit where the number of components of d and N is taken to be infinitely large [49].The saddle point values of ψ d and ψ N will in general be non-zero, which we absorb into a redefinition of r d,N .Then, the effective action for phase (B) becomes S eff = S χ + S e + S c + S φ where While generically both saddle point values φ 0 dN and φ 0 dd are expected to be non-zero simultaneously in phase (B), we take φ 0 dd → 0 and φ 0 dN ≡ φ 0 = 0 for the following explicit calculations.Setting φ 0 dd = 0 does not change any symmetries of the phase, allows for a more compact discussion of the results, and can formally be seen as the large b 2 limit of the theory where φ 0 dd is suppressed [cf.Fig. 1(a)].More generally than the derivation of S eff via Hubbard-Stratonovich transformations, it can also be thought of as the simplest field theory capturing the key aspects of phase (B) in Fig. 1(a) at finite T .Electronic self energy.-Tocompute the spectral properties of the electrons within this model, we employ a large-N technique similar to [50,51]: we add extra indices to the electrons and bosons, c k,τ,s → c k,τ,s,a , d ab → d ab and similarly for N , where a, b = 1, 2, . . ., N , which are contracted in all terms of S eff so as to ensure O(N ) symmetry.In the limit N → ∞, the electronic self-energy Σ is given by the "rainbow diagrams" [50,51] shown in Fig. 2(a).In our case, however, Σ involves both normal and anomalous contributions as a result of the anomalous bosonic term ∝ φ 0 in Eq. (1).To make this more explicit, we integrate out the bosons, yielding the effec-tive fermionic interactions S int = S 1 + S 2 with where k+q sis y τ y c † −k .The two terms in S 1 describe spin and superconducting triplet fluctuations, respectively; their associated self-energy contributions are normal in the sense that U (1) symmetry is preserved, with leading terms represented by the first two diagrams Σ 1,2 in Fig. 2(b).Conversely, S 2 breaks U (1) symmetry, when φ 0 attains a mean-field value, and results in an anomalous contribution to the self-energy, with leading term given by the last diagram Σ 3 in Fig. 2(b).
To represent the diagrams algebraically, we shift to the Bogoliubov-de Gennes basis (c q,+ , is y c † −q,− ) T , with Pauli matrices γ i acting on this space.In this basis, the free Green's function is G 0 (iω, ) = iω − γ z .Up to first order in λ 2 , the spin-spin self energy term can be written as After performing a gauge transformation to make φ 0 real, the anomalous term from the spin-triplet interaction is given by For concreteness and since spin fluctuations are believed to occur already at higher energies than superconducting fluctuations in graphene moiré systems [37,38], we focus on r d > r N ; we will use Density of states.-Figure2(c) shows the effect of the normal contributions of the self energy Σ 1,2 on the DOS of a 2D parabolic band.The effect of Σ 1 is to push the peak of the free spectral function at energy away from ω = 0.This results in the opening of a gap (which can be soft depending on the parameter regime), very similar to fluctuating anti-ferromagnetism discussed in the cuprates [52][53][54].Σ 2 on the other hand has the opposite effect, where it pushes states towards ω = 0.This is because Σ 1 and Σ 2 have the exact same functional form with one key difference: k+q of Σ 1 is replaced by − k+q in Σ 2 .The effect of the total normal self energy Σ 1 + Σ 2 is to enhance the DOS in the vicinity of the Fermi level, see Fig. 2(c).The anomalous contribution Σ 3 does not interfere in these effects since it occurs in the γ y channel.The role of Σ 1 + Σ 2 can, thus, be intuitively thought of as providing a renormalized DOS in the normal state on top of which the anomalous Σ 3 opens up FIG.3: Spectral weight as a function of ω with (blue) and without (purple) Σ3 (a) close to k = 0 and (b) including a larger energy range; in both cases, we focus on the q = 0 contribution (see text).(c) The effect of all three self energy contributions Σ1 + Σ2 + Σ3 (including the momentum integration) on the DOS.For small φ0, there is suppression of the DOS at ω = 0 which resembles the V-shaped DOS of a nodal state.For large φ0, the gap resembles a hard BCS gap.
a gap.We have checked [48] by numerically solving the self-consistency equation for the self-energy [Fig.2(a)] in the limit (of large v µ ) where only the q = 0 term of the momentum sum contributes that higher-order corrections do not change our results qualitatively for small φ 0 .For instance, Fig. 2(d) shows the numerical solution for the Green's function G = iω − (iω)γ z + ∆(iω)γ y in Matsubara space upon including the effect of S 2 ; the difference to the first-order result is small.To gain intuition for the impact of Σ 3 on the DOS, we first focus again on the q = 0 term of the momentum sum in Eq. (3).In this limit, one can easily see [48] that Σ 3 vanishes linearly in k for small energies.Since Σ 3 is in the γ y channel, the effect of any non-zero value is to generically open a gap.As a result of the linear behavior, the states exactly at zero energy are unaffected, but slightly away from it, the states get pushed away to higher energy; this is clearly visible in Fig. 3(a).In contrast, for large energies, Σ 3 is readily seen to tend to zero.The spectral function, thus, remains asymptotically unaffected, as can be seen in Fig. 3(b).Taken together, we expect the DOS to be reduced (but not fully suppressed for small φ 0 ) in an energy range around the Fermi level, exhibiting an enhancement with respect to its normal-state value at intermediate energies, and then approaching the normal-state limit at larger energies.
To demonstrate this explicitly beyond the simple q = 0 limit, we approximate k+q k +v F q +q 2 /(2m), where q is the component of q along k, and numerically evaluate the momentum integrals to find the total self energy for concreteness, Fig. 3(c) shows the resulting DOS.As expected, we see that there is a suppression of the DOS.However, for small values of φ 0 , the resulting DOS has a V-shaped behavior, which is typically only seen in nodal states (with either nodal lines or points).Recall that the superconducting phase in our model is symmetryequivalent to a conventional BCS state and that the triplet superconductor that arises at T = 0 in phase (B) will be fully gapped.For larger φ 0 , the gap at ω = 0 increases, and resembles a hard BCS gap.The suppression of the DOS ρ F at ω = 0 can be estimated analytically at finite temperature by again taking the limit (of large v µ ) where the integration over q can be replaced by an evaluation at q = 0; we find Note that φ 2 0 is bounded above by r d r N , at which point the bosonic fields would condense and continuous symmetries would be broken, which cannot happen at finite T .As φ 0 increases, α increases the suppression of the DOS, and near the instability point of φ 2 0 = r d r N , there are no states near the Fermi surface.
To complement this analysis, we have also studied the Hamiltonian associated with setting q = 0 in Eq. (2b) within self-consistent Hartree-Fock, only allowing for spin-rotation invariant operators to condense [48].For small α, one also finds only a partial suppression of the low-energy spectral weight, akin to Eq. ( 4); including higher-order corrections leads to a hard gap for α ≥ 1.
Electromagnetic response.-Wewill finally demonstrate that the superconducting phase characterized by φ 0 = 0 has the same electromagnetic phenomenology as BCS superconductors, despite the unusual electronic spectral properties.
First focusing on the electrons, we show that c to leading (first) order in φ 0 ; as non-zero Ψ F to linear order in φ 0 implies that it cannot vanish identically for generic φ 0 , this is sufficient to show the presence of ODLRO.We find the "macroscopic wave function" to be a singlet, Ψ F (x) = is y ψ F (x), as expected since spin-rotation symmetry is preserved at finite T , with ψ F (x) shown in Fig. 4(a).Alternatively, one can demonstrate ODLRO to arbitrary order in φ 0 , by focusing on the bosons: to zeroth order in λ, we find 4(b) along with an analytic asymptotic form for large x; in [48], we show that this leads to the same constraints as the conventional form of bosonic ODLRO [55,56].Finally, the connection to the textbook theory of superconductivity can be made more explicit by deriving the analogue of the time-dependent Ginzburg-Landay theory: we reinstate fluctuations via φ 0 → φ(x, τ ) and integrate out all other degrees of freedom yielding to leading order in φ and gauge-covariant derivatives (D τ , D) µ = ∂ µ − i2eA µ .For demonstration purposes, we evaluated the coefficients in S GL to leading (zeroth) order in S c and find ρ, v φ > 0 and r φ < 0 for low T [see Fig. 4(c,d)]; the state with QLRO in φ 0 thus corresponds, as usual, to the Higgs phase, with Meissner effect and massive Higgs mode, but without Goldstone modes.
Conclusion.-We have studied the finite-T vestigial phases, see Fig. 1(a), associated with two primary order parameters, d and N , describing a fully gapped triplet superconductor and spin magnetism, respectively.A crucial result is the DOS of phase (B 1,2 ) in Fig. 3(c): varying φ 0 changes the low-energy DOS from partial suppression, akin to that of a nodal superconducting state, to a hard gap.As φ 0 is expected to change with electron filling, this could explain the tunneling data in [34,35].We finally point out that the suppression of N would immediately also suppress φ 0 in our model and could, therefore, explain why superconductivity is connected to the reset behavior in trilayer graphene [34,35,39,40].
can change the sign of b 2 to positive values [45].

Appendix B: Evaluation of the self-energies at leading order
In this section, we show the evaluation of the self energies up to first order in perturbation theory.We first evaluate the anomalous part of the self energy, Σ 3 in Fig. 2(b), which is contributed by the anomalous term of the action given by In the following, we work in the c q,+ is y c † −q,− T Bogoliubov-de Gennes basis, with the Pauli matrices γ i acting on it.The free Green's function then reads as Choosing φ 0 to be real, we have where , and g µ = r µ + v 2 µ q 2 .Thus, The Matsubara sum can be evaluated using where n f /B ( ) = 1 e β ±1 .Thus we get, where we performed a partial fraction decomposition of The normal part of the self energy, Σ 1,2 in Fig. 2(b), is contributed by the following term of the action Defining γ ± = 1 2 (γ x ± iγ y ) , the corresponding contribution to the self energy is given by As a result, if we consider the self energies as function of iω and k+q , we find that G 0 (iω, − k+q ).This allows us to argue the effect of Σ 2 pushing high energy states towards the vicinity of ω = 0, while Σ 1 pushes states away from ω = 0.
To perform the Matsubara sums, we define with K(iω, , E) as defined in (B8).In terms of these functions, the self energy is given by We can expand the total self energy Σ = Σ 1 + Σ 2 + Σ 3 in terms of Pauli matrices in Nambu space, where Appendix C: Suppression of DOS at ω = 0 In this section, we derive a compact approximate analytical expression for the suppression of the density of states (DOS) as a result of the anomalous term Σ 3 = Σ γy γ y .To this end, we focus on the limit of large bosonic velocities v µ in χ µ and replace the q integral in Eq. (B20) with the value of the integrand at q = 0, Note that we would first need to re-parametrize the integral in terms of q = q √ r N /v N and then set q = 0.This approximation would then be valid in the large v d /v N limit with this re-scaling.We then Taylor expand f (z, ) with respect to , ω, at a non-zero finite T (satisfying T 4 r d r N − φ 2 0 ).In this limit, we find the self energy to be This expression is in agreement with the result in the main text [Fig.3(a)] which shows that as → 0, the contribution of Σ y vanishes.With such a self-energy, the spectral function is given by A simple way to look at this, is that the band structure is simply renormalized as k → √ 1 + α 2 k .This reduces the effective band mass, and thus the DOS is suppressed by a factor of √ 1 + α 2 , as stated in the main text.
-  -  6: The first order solution to ε(iω) (red) and the self consistent solution (green) after including the effects of Σ1 (left column) and Σ2 (right column) separately.Same parameters as in Fig. 5.
not get renormalized since Σ 3 acquires a γ z term only if the Green's function has a γ y term.Such a γ y term does not exist in the normal state about which we perform perturbation theory.As φ 0 increases, we find that the self consistent solution is lower in magnitude that the first order solution.
In Fig. 6, we show the corrections in ε(iω n ) after including the effects of Σ 1 (left column) and Σ 2 (right column).As expected and argued in the main text, we find that Σ 1 and Σ 2 have qualitatively the opposite effects on the renormalization of ε(iω n ).In both the cases, we find that the magnitude of the self consistent solution is higher than the perturbative corrections.However, since the fermionic Matsubara frequencies do not contain 0, we cannot directly say what this implies for the solution on the real axis.The magnitude of φ 0 has little effect on the solution since the effect of spin and triplet fluctuations are controlled by g N and g d , respectively, which we keep constant.
In Fig. 7, we plot the corrections in ε(iω n ) and ∆(iω n ) after including the effects of all the self energies Σ = Σ 1 + Σ 2 + Σ 3 .We find that the inclusion Σ 1 and Σ 2 together reduces the difference between the self consistent and perturbative solution (refer to the plot near φ 0 ∼ 0).As we increase φ 0 , the difference between the self consistent and perturbative solution increases due to the effect of Σ 3 which is controlled by φ 0 .The first order solution to ε(iω), ∆(iω) (red) and the self consistent solution (green) after including the effects of all the terms of the self energy Σ1 + Σ2 + Σ3.Same parameters as in Fig. 5.
Taken together, we see that the inclusion of second-and higher-order diagrams that contribute in the large-N limit defined in the main text yields qualitatively similar behavior on the imaginary axis compared to the first-order diagrams.We therefore expect that the qualitative picture that S 1 renormalizes the DOS close to the Fermi level on top of which S 2 reduces the low-energy spectral weight still applies.Since the impact of S 2 is controlled by small φ 0 and good quantitative agreement is found for φ 0 up to 0.6r N , we expect that Fig. 3(c) would look similar when higher-order corrections were included.
Since the superconducting pairing takes place only between electrons between opposite valleys, we will have only τ 2 = −τ 1 giving non-zero correlators.Without loss of generality we chose τ 1 = +, τ 2 = −.Up to first order in φ 0 , we have where ... is the average with respect to the interacting and ... 0 with respect to the non-interacting ground state.We define to be the Green's function in the fermionic basis (assuming k = −k ).Equation (F9) can then be evaluated as, We continue by calculating the Matsubara sum over iΩ n and over iω n [see Eq. (F8)], For simplicity, we here focus on the limit where the remaining sum over q in Eq. (F9) is determined by its q = 0 component.With E ± ≡ E ± (q = 0) and v N , r N = 1, we have we can then finally write In the second line, we assumed k = 2 (k 2 − k 2 F )/2m.Using this expression, we calculate the spatial profile of the fermionic ODLRO wavefunction numerically for various values of F ≡ k F in Fig. 4(a).Unlike the case of the bosonic ODLRO (which was exponentially decaying), the fermionic ODLRO has an oscillating component superimposed on an exponentially decaying envelope.

Appendix G: Ginzburg-Landau theory
We here calculate the Landau-Ginzburg theory for the bosonic superfluid condensate parameter to leading (zeroth) order in the fermion-boson coupling λ.To tis end, we assume that φ 0 is now spatially and temporally varying.This results in non-zero Fourier modes φ q for q, iΩ = 0.
In momentum space, the bosonic action is generalized according to So after integrating out d and N , the effective action for φ reads as where To derive the Ginzburg-Landau theory for φ, we expand T r ln G −1 upto second order in Γ, which is equivalent to second order in φ.This gives us Because of the diagonal structure of G 0 , and the off diagonal structure of Γ, the linear term T rG 0 Γ is 0. The quadratic term becomes q ,q T rG 0 (q + q)Γ(q + q, q )G 0 (q )Γ(q , q + q) = q ,q T r 0 χ N (q + q)φ q χ d (q + q)φ * −q 0 0 χ N (q )φ −q χ d (q )φ * q 0 (G9) = q ,q χ N (q + q)χ d (q )φ q φ * q + χ N (q )χ d (q + q)φ −q φ * −q (G10) = q ,q (χ N (q + q)χ d (q ) + χ N (q )χ d (q − q)) φ q φ * q (G11) = q ,q (χ N (q + q)χ d (q ) + χ N (q + q)χ d (q )) φ q φ * q (G12) = 2 q ,q χ N (q + q)χ d (q )φ q φ * q (G13) We need to evaluate By expanding the above expression up to second order in iΩ, q, we find the effective action for the φ field to be where the coefficients are given by with g µ = r µ + v 2 µ q 2 .We numerically calculate the quantities r φ , ρ, v 2 and plot it in Fig. 4(c,d) of the main text.
while the free Hamiltonian is given by We consider only the effective Hamiltonian at time t = 0, which is why there are no Matsuabra indices.We perform a Hartree-Fock decomposition of V , which gives us where , choosing a gauge with real φ 0 ; further take φ 0 to be positive such that c > 0. Note that this correlator is related to the Green's function G by C k = T iωn G(k).Note that all the Hartree terms vanish since we do not allow for spontaneous breaking of spin-rotation invariance (recall we study finite T in 2D).The effective 2−particle Hamiltonian is given by where ˜ k , ∆k are the self consistent band structure and gap.Making connection with the diagrammatic self consistency relationship to be discussed below, we can foresee that the resulting self consistent equation we get will be the same as (H18) but with ˜ , ∆ replaced with the corresponding iω n averaged value, and the whole equation itself will be iω n averaged.
The correlators in terms of ˜ , ∆ are given by where E k = ˜ 2 k + ∆2 k > 0. Thus, using (H7), the self consistency equations become and first assume β k < 1, which always holds as long as T > c/4.The self consistency equations can then be rearranged as Using this, we find Note, however, that β k also depends on E k and, thus, this relation should be thought of as a self consistency equation, to be solved for β k or E k .
Equations (H12) allow to derive asymptotic relations.In the limit k → 0, we then have E k → 0 and β k → c 4T , ensuring the self-consistent solutions are well controlled in the k → 0 regime that we are interested in.Near k = 0 and for large T c (β k 1), the renormalized spectrum is given by E k = √ 2 .Figure 8 illustrates the behavior obtained by numerical solution of the self-consistency equations.

Zero energy-momentum transfer
In this section, we consider the limit where the bosonic fields N , d do not transfer any momentum or Matsubara frequency in the interaction (q = 0 in S c ).Additionally, we consider only the effect of S 2 on the self energy to study the effect of the anomalous contribution.In this limit, we would like to analyze the self consistent solution of the Green's function up to all orders in λ within the large-N theory of the main text.The ansatz of the full Green's function is given by G −1 = iω n − ˜ k γ z − ∆k γ y , since Σ 3 renormalizes only the anomalous term ∆k and the spectrum ˜ k .We have Thus the self-consistent analogue of Σ 3 in Eq. ( 3) becomes (where we have replaced the integration over q by the q = 0 value of the integrand, and φ0 = φ 0 λ 2 r N /v 2 N ) From the self-energy equation we get where c = 6 φ0 r N r d −φ 2 0 .Right at the Fermi surface, k = 0, the self consistency equations reduce to There are two possible solutions to Eqs. (H19) and (H20).The first is ˜ k = ∆k = 0; this is exactly what we find within perturbation theory.For a solution with ˜ k = 0 to exist, it must hold (assuming ˜ k , ∆k ∈ R as expected in the gauge that we use) Thus, a non-zero solution only exists if T < c/π 2 ∼ c/9.As compared to Hartree-Fock, the critical temperature for a non-perturbative solution is lower.

ΦFIG. 1 :
FIG. 1: (a) Mean-field phase diagram for r d = rN , b3 = b1, c2 = 0, where we indicate the symmetries at T = 0 (blue), those of the resulting vestigial phases at T > 0 (red), and which composite order parameters are finite.Solid (dashed) orange lines are phase transitions at T = 0 and T > 0 (become a crossover at T > 0).(b,c) illustrate the finite-T pairing in phases (A) and (B) schematically.

and set χ 0 µ = 1
by rescaling of the fields.

FIG. 5 :
FIG.5:The first order solution to ε(iω), ∆(iω) (red) and the self consistent solution (green) for the self energy in Matsubara space.Note the offset by 0.1 in the y axis in the left column.We chose k = 0.1, r d = 9, rN = 1, T = 1 β = 0.2, λ = 1 and measured all energies in units of √ rN .
FIG. 7:The first order solution to ε(iω), ∆(iω) (red) and the self consistent solution (green) after including the effects of all the terms of the self energy Σ1 + Σ2 + Σ3.Same parameters as in Fig.5.

, 6 FIG. 8 :4T 1 − 1 12 E 2 kT 2 .
FIG.8:The self consistent solution for ˜ k , ∆k and E k as a function of k for various temperatures.At T /c = 1/4, the self-consistent solutions become non-analytic having an infinite slope at k = 0, and a gap opens up as the temperature decreases.There is a discontinuity in ˜ k , ∆k at k = 0, where the gap value has different signs for k → 0 − , 0 + .