Multiple intertwined pairing states and temperature-sensitive gap anisotropy for superconductivity at a nematic quantum-critical point

The proximity of many strongly correlated superconductors to density-wave or nematic order has led to an extensive search for fingerprints of pairing mediated by dynamical quantum-critical (QC) fluctuations of the corresponding order parameter. Here we study anisotropic s-wave superconductivity induced by anisotropic QC dynamical nematic fluctuations. We solve the non-linear gap equation for the pairing gap Δ(θ,ωm)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta (\theta ,{\omega }_{m})$$\end{document} and show that its angular dependence strongly varies below Tc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{\rm{c}}}$$\end{document}. We show that this variation is a signature of QC pairing and comes about because there are multiple s-wave pairing instabilities with closely spaced transition temperatures Tc,n\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{\rm{c}},n}$$\end{document}. Taken alone, each instability would produce a gap Δ(θ,ωm)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta (\theta ,{\omega }_{m})$$\end{document} that changes sign 8n\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$8n$$\end{document} times along the Fermi surface. We show that the equilibrium gap Δ(θ,ωm)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta (\theta ,{\omega }_{m})$$\end{document} is a superposition of multiple components that are nonlinearly induced below the actual Tc=Tc,0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{{\rm{c}}}={T}_{{\rm{c}},0}$$\end{document}, and get resonantly enhanced at T=Tc,n<Tc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T={T}_{{\rm{c}},n}\ <\ {T}_{{\rm{c}}}$$\end{document}. This gives rise to strong temperature variation of the angular dependence of Δ(θ,ωm)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta (\theta ,{\omega }_{m})$$\end{document}. This variation progressively disappears away from a QC point.


INTRODUCTION
Recent observations of superconductivity (SC) near a nematic quantum-critical point (QCP) in Fe-based superconductors, such as FeSe, [1][2][3][4][5][6][7] renewed interest in studies of SC in proximity to a densitywave or a nematic order. 8,9 Several researchers hinted 2,10-14 that the SC dome observed in several Fe-SCs may be the consequence of the pairing mediated by nematic fluctuations. However, a SC dome can appear for other reasons as well, e.g., owing to the "fight" for the Fermi surface (FS) between SC and density-wave orders, even when each is described within BCS/mean-field theory. [15][16][17][18][19][20] The question we address here is whether there are unique features of SC near a nematic QCP, encoded in the SC gap structure on the FS, Δðθ; ω m Þ (θ is the angle along the FS and ω m is Matsubara frequency). Previous studies have focused on signatures of QC pairing in the frequency dependence of the gap function. [21][22][23] Such emphasis stems from the understanding that the pairing kernel in the QC regime is a singular function of frequency. In contrast, the angular dependence of the gap along the FS was assumed to be set either by the non-s-wave pairing symmetry (e.g., d-wave in the cuprates 24,25 ), or, for s-wave, by some material specific non-singular angular dependencies of interactions and band structures, as in the Fe-based superconductors. [26][27][28][29][30] In either case the angular variation of the gap Δðθ; ω m Þ was expected to be set at T c and not vary strongly in the SC state.
In this paper, we demonstrate that QC pairing can give rise to a strong temperature evolution of Δðθ; ω m Þ below T c at any given ω m . Specifically, we argue that this is the case for s-wave SC near a nematic transition in 2D, e.g., a transition into a state with spontaneously broken symmetry between d xz and d yz orbitals in FeSe. 3,4,7,11,26,27,29,[31][32][33][34][35] Several previous works, including by some of us, 36,37 compared the values of T c for s-wave and d-wave pairings owing to anisotropic dynamical nematic fluctuations, which give rise to an attraction in both channels, 36 and analyzed s-wave and d-wave gap structures right at T c by solving the linearized gap equations for Δðθ; ω m Þ. These studies found that T c for s-wave is larger. The s-wave gap is anisotropic, with four maxima along the FS. This gap anisotropy reflects the anisotropy of the pairing interaction and by itself is not a signature of quantum criticality. Here, we argue that the gap evolution below T c is the signature of quantum criticality. We show that the mechanism, driving the evolution of Δðθ; ω m Þ below T c , is the existence of multiple s-wave pairing states with closely spaced transition temperatures T c;n . These solutions yield Δ n ðθ; ω m Þ of the same symmetry, but with a different number (¼ 8n) of sign changes of the gap along the FS. The closeness of T c;n stems from the long range nature of the interaction near a QCP, which reduces the energy cost of gap oscillations. Taken separately, the state with the largest condensation energy is a conventional swave state with sign-preserving gap Δ 0 ðθ; ω m Þ, which emerges below the highest T c;0 ¼ T c . However, all other Δ n are induced below T c , and each gets resonantly enhanced below T c;n . As a result, the actual Δðθ; ω m Þ coincides with Δ 0 ðθ; ω m Þ only near T c , while at smaller T the gap structure at any frequency ω m is a mixture of different Δ n ðθ; ω m Þ.
We emphasize that all Δ n have the same symmetry (s-wave). In this respect, the gap evolution below T c , which we consider here, is different from the one induced by a change of a gap symmetry from, e.g., s-wave to d-wave, or owing to the emergence of a mixed s-d state. We argue that multiple closely spaced solutions for T c;n exist only in proximity to the nematic QCP. Away from the QCP, other solutions shift to smaller T and progressively disappear as the nematic correlation length gets smaller. An experimental observation of a strong variation of the shape of Δðθ; ω m Þ below T c would then be conclusive evidence for QC pairing.

RESULTS
We begin with a qualitative explanation for the existence of multiple pairing states and the temperature evolution of the gap. We consider s-wave SC of 2D fermions, minimally coupled to order parameter fluctuations near a nematic QCP, at which the system spontaneously breaks C 4 lattice symmetry. [38][39][40][41][42][43][44][45][46][47][48][49] The strength of fermion-boson interaction is described by a dimensionless coupling λ (defined below), which we treat as a small parameter to keep calculations under control. The nematic order develops at q ¼ 0, so the pairing interaction is peaked at zero momentum transfer and involves fermions at any angle θ on the FS. Still, the pairing interaction does depend on θ via the square of the nematic form-factor f 2 ðθÞ ¼ cos 2 2θ. As a result, the pairing interaction is larger in 'hot' regions near θ ¼ π 2 m, and is smaller in 'lukewarm' regions near θ ¼ π 4 þ m π 2 (see Fig. 1). 36,37 This separation is not sharp as the width of both lukewarm and hot regions is of order one.
At the QCP, the pairing boson is massless, and its dynamics cannot be neglected. The scale of bosonic dynamics (the Landau damping) is set by the same interaction, which gives rise to the pairing. As a result, the strongest pairing occurs for fermions within a small angular separation of order δθ $ λ, which is small compared with the width of a hot region. To leading order in λ, the gap equation then becomes local and allows a continuous set of solutions ΔðθÞ / δðθ À νÞ with T c ðνÞ / f 4 ðνÞ, where ν is a continuous parameter. The physical T c is the highest temperature in the set, and it corresponds to ν ¼ πm=2, where the form-factor f ðνÞ ¼ 1. The actual gap structure is determined at the next approximation, when one properly accounts for weaker interactions at angle transfers above δθ. The result is that in each octant, e.g., at 0 < jθj < π=4, the actual gap magnitude is of order Δðθ ¼ 0Þ for jθj < θ h $ λ 1=3 (δθ ( θ h ( 1) and rapidly drops as ðθ h =θÞ 4 at larger θ. This can be interpreted as if a Cooper pair is 'trapped' in a potential well within the range θ h . Figure 2 depicts the interplay of the different scales in the problem. The correlation function between fermions in a pair (the gap function) can vary inside the trap at the cost of a small kinetic energy $ λ=λ 1=3 ¼ λ 2=3 . In this situation, the continuous set of T 0 c s transforms into a discrete set T c;n , each corresponds to the solution of a Schroedinger-like equation for the gap function in a box with the width 2θ h . The discrete T c;n differ from T c;0 by multiples of the kinetic energy cost The solutions exist up to n max $ 1=λ ) 1. The corresponding gap functions Δ n ðθÞ are all fourfold periodic, i.e., s-wave, but have n nodes in the first octant (8n total along the FS). We depict several such oscillating solutions in Fig. 4. The gap oscillations can be interpreted semiclassically as back-and-forth motion of a localized Cooper pair wave-packet within the gap's "trapping" potential. In our approximate analytical treatment they appear as the solutions of the Airy equation, which naturally emerges when we linearize the 'trapping' potential (see Methods section). At T < T c , the physical gap is a superposition of Δ 0 and Δ n>0 , which are all induced by Δ 0 , because all Δ n have the same symmetry. This superposition causes destructive interference in the trap region, where oscillations occur, and constructive interference outside it. As a result, the gap width increases strongly with decreasing temperature, as more oscillating Δ n>0 are superimposed on the non-oscillating Δ 0 . Figure 3 shows the numerical solution of the non-linear gap equation. A strong increase of the gap width with decreasing temperature is clearly visible.
Multiple pairing states of s-wave symmetry have been found previously by Yang and Sondhi (YS) 50 in their analysis of the pairing owing to a near-local static interaction. In their case, Δ n is angle-independent, but oscillates n times as a function of the distance to the FS, k À k F . YS argued that additional solutions affect superconducting stiffness and reduce T c to a smaller value, proportional to the size of their short-range potential. Our dynamical theory is different in two aspects. First, the gap oscillates on the FS, i.e., at k ¼ k F . Second, the potential range is self-consistently set by the dynamics, and in this self-consistent framework the reduction of superfluid stiffness (and, hence, of T c ) is at most Oð1Þ.

Model and gap equation
We consider 2D fermions with Fermi energy E F , minimally coupled to a nematic order parameter field Δ nem ðqÞ by where f ðkÞ is a form-factor, which has d-wave symmetry with respect to C 4 lattice rotations. We assume that the static susceptibility χðqÞ of the nematic field is peaked at q ¼ 0: where ξ 0 is the bare correlation length. We define the dimensionless coupling to be λ ¼ g 2 χ 0 =4E F . We assume for simplicity a circular FS, but our results are readily generalized to other C 4 -symmetric FSs. Because relevant fermions are near the FS, we can approximate f ðkÞ by f ðθÞ ¼ cos 2θ (see Fig. 1). Below we focus on the octant 0 < θ < π=4.  A bosonic excitation with momentum q connects two fermions on the FS with angles θ; θ þ ϕ such that q ¼ 2k F sin jϕ=2j % k F jϕj. The effective pairing interaction is then the function of both ϕ (via χðqÞ) and θ (via the form-factor f ðθÞ). This interaction modifies both bosonic and fermionic variables. Fermionic self-energy scales as Σðω m Þ $ ðλ 2 E F Þ 1=3 jω m j 2=3 sgnðω m Þ, modulo logarithmic corrections from higher-order planar diagrams, 43,46,51,52 which we neglect here, and singular contributions from thermal fluctuations at T > 0. The thermal contribution cancels out between the selfenergy and the pairing vertex 53-55 and we eliminate it in Σðω m Þ and in the gap Eq. (4) below. The bosonic self-energy renormalizes the bare correlation length ξ 0 into the true ξ, which vanishes at the QCP, and also generates a dynamical Landau damping term. The dressed bosonic susceptibility at a QCP is To obtain T c and the gap function we need to analyze the equation for the pairing vertex Φðθ; ω m Þ. The pairing gap Δ is related to Φ in a usual way, as Φðθ; The discrete set of solutions of the linearized gap equation We first analyze the gap equation for infinitesimally small Δðθ; ω m Þ. We explicitly integrate out fermionic momenta transverse to the FS in the gap equation and obtain the equation for the gap function Δðθ; ω m Þ on the FS: Equation (4) allows solutions with different gap symmetry. Earlier works found 36,37,56 that the s-wave solution has the largest T c , so we focus on s-wave gap. Because χðθ; ϕ; ω m 0 À ω m Þ is strongly peaked at ϕ ¼ 0, to first approximation we may set This yields a local equation for Δðθ; ω m Þ ¼ Δðθ; mÞ, with θ acting as a parameter. 38,47,55,57,58 At the QCP we have Equation (5) has a continuous set of solutions Δðθ; mÞ / δðθ À νÞ with arbitrary ν from the interval 0 < ν < π=4, and T c ðνÞ $ λ 2 E F f 4 ðνÞ. The maximum T c ð0Þ % 0:022λ 2 E F f 4 ð0Þ corresponds to ν ¼ 0. To determine the actual structure of Δðθ; mÞ and the correct number of solutions, we need to go beyond the leading approximation and keep the dependence on ϕ in the numerator in (4). The problem is analytically tractable if we use the fact that typical ω m ; ω m 0 $ T, i.e., typical m; m 0 ¼ Oð1Þ and typical Z ¼ Oð1Þ, and simplify the gap equation by neglecting the frequency dependence of Δðθ; mÞ and setting Z ¼ 1. In this approximation, Eq. (4) becomes an effective 1D integral equation over the angle. We expand in small angles near θ ¼ 0 and obtain Δðθ þ ϕÞ À ΔðθÞ (See the Methods section for the detailed derivation.) Here, ηðTÞ $ ððT c ð0Þ=TÞ 1=3 À 1Þ=θ 2 h and ðθ; ϕÞ ¼ θ À1 h ðθ; ϕÞ, where θ h $ λ 1=3 sets the width of the gap function ΔðθÞ. Transforming to a Fourier representation ΔðθÞ ¼ 2π ð Þ À1 R dx expðixθÞΔðxÞ, we find from (6) ηðTÞΔðxÞ ¼ À∂ 2 x ΔðxÞ þ jxjΔðxÞ: For even ΔðxÞ, which we consider, this reduces to the Airy equation. It has orthogonal solutions η n ¼ x n ; Δ n ðxÞ / Aiðjxj À x n Þ; where Àx n is the nth zero of the derivative of the Airy function Ai 0 ðxÞ. For later convenience, we choose Δ n ðθÞ to be orthonormal: R dθΔ n ðθÞΔ m ðθÞ ¼ δ mn . As we described earlier in this Section, the appearance of the Airy equation is attributed to the fact that the Cooper pair can oscillate in the 'trapping potential' set by θ h .
The smallest eigenvalue η n is for the even solution with n ¼ 0. For ΔðθÞ in the original coordinates, it yields a non-oscillating gap Δ 0 ðθÞ, peaked at θ ¼ 0, with the width Oðθ h Þ. The corresponding T c;0 is the actual T c for the pairing at a nematic QCP. It differs from T c ð0Þ from Eq. (5) by a numerical factor of order λ 2=3 . For other solutions, Δ n ðθÞ changes sign n times in the first octant. Extracting T c;n from Eq. (8), we obtain T c;n ¼ T c;0 1 À OðnλÞ 2=3 , up to n max $ 1=λ ) 1. Figure 4 depicts the first few even solutions. Although the detailed form of these solutions depends on λ and f ðθÞ, their appearance is the universal feature of QC pairing by nematic fluctuations.
Non-linear gap equation and evolution of θ h At T only slightly below T c;0 ¼ T c , the angular dependence of the pairing gap coincides with Δ 0 ðθÞ, just the overall gap magnitude increases with decreasing T. To obtain the form of ΔðθÞ at lower T, we solve the full non-linear gap equation (see Methods) in the temperature range down to T ( T c . Figure 3 depicts the evolution of the amplitude and the width of the gap θ h ðtÞ as a function of the reduced temperature t ¼ 1 À T=T c . We see from Fig. 3 that θ h ðtÞ exhibits strong temperature evolution, i.e., ΔðθÞ substantially broadens below T c . We are not aware of other models that show such behavior of the gap function in the SC state (see Discussion section). We now argue that this drastic evolution of the shape of ΔðθÞ comes about because the width of ΔðθÞ is strong influenced by orthogonal components Δ n ðθÞ. These components have the same symmetry as Δ 0 ðθÞ and are all induced by Δ 0 ðθÞ immediately below T c . They are initially small, but the Δ 0 ðθÞ induces the orthogonal gap components with n < 0, and the nth component Δ n ðθÞ gets resonantly enhanced at t $ ðnλÞ 2=3 , and its magnitude becomes comparable to Δ 0 ðθÞ. These additional gap components interfere destructively in the region jθj < θ h , where they oscillate, but constructively for θ > θ h .
To demonstrate this, we express the gap function as ΔðθÞ ¼ P n!0 Δ n Δ n ðθÞ and derive the Landau free energy expansion for the coefficients Δ n . We obtain where t n ¼ 1 À T c;n =T c (t À t n ¼ ðT c;n À TÞ=T c ), a nk / R dθΔ 4Àk n ðθÞΔ k 0 ðθÞ, and stand for less relevant terms in F, e.g., the coupling between two states with n > 0.
We stress that the solution for ΔðθÞ, obtained by minimizing Eq. (9), is equivalent to the solution of the non-linear gap equation. At small θ, Δ n ðθÞ / ðÀ1Þ n , because the sign at the origin depends on the number of oscillations in each octant. At large θ ) θ h all Δ n ðθÞ decay as 1=θ 4 (see Methods section). Solving the saddle-point equations, we find that immediately below T c , when 0 < t ( 1, Δ 2 0 % t=a 0 and jΔ n j % jðΔ 0 Þ 3 a n3 =t n j ( Δ 0 . In the opposite case, when t is finite and t À t n % t, we have, using the numerical results for a n and a ni (i ¼ 1 À 3) (see Methods section), Δ n $ ðÀ1Þ nþ1 Δ 0 , i.e., jΔ n j is of the same order as Δ 0 . Figure 5 illustrates the evolution of Δ n with temperature. Because at small θ, Δ n ðθÞ / ðÀ1Þ n Δ 0 and Δ n / ðÀ1Þ nþ1 Δ 0 , all n > 0 contributions to ΔðθÞ (the terms Δ n Δ n ðθÞ) are of opposite sign compared with Δ 0 Δ 0 ðθÞ. The original and the induced gap components then interfere destructively, and the total Δðθ % 0Þ decreases with decreasing T. On the other hand, at θ ) θ h , Δ n ðθÞ / Δ 0 , hence Δ n Δ n ðθÞ / ðÀ1Þ nþ1 oscillates in sign between even and odd n. Because Δ n Δ n ðθÞ decrease in magnitude with increasing n, and Δ 1 Δ 1 ðθÞ has the same sign as Δ 0 Δ 0 ðθÞ, the original and the induced gap components then interfere mostly constructively. As the consequence, the effective 'width' of the gap, θ h ðtÞ / R Δðt; θÞdθ, increases with decreasing T. To leading order at t ( 1, θ h obeys We emphasize that (i) the temperature variation is strong, particularly at small λ, and (ii) it is a non-linear effect owing to the existence of multiple sign-changing Δ n ðθÞ, all having the same s-wave symmetry. The inset of Fig. 3 shows that θ h ðtÞ indeed increases as t α with α % 1, in agreement with Eq. (10). We also verified numerically (see Fig. 7 in the Methods section) that the scaling form of Eq. (10) holds even when λ is of order one and θ h is comparable to the width of the hot region.
Away from the QCP Our results are readily generalized (see Methods section) to the case when the correlation length ξ for nematic fluctuations is large but finite, i.e., the system is at a finite distance from the QCP. For ðk F ξÞ À1 < λ, the effect of ξ is a uniform reduction of all T c;n . However, for ðk F ξÞ À1 > λ, we find that t n ¼ 1 À T c;n =T c $ ðnλÞ 2=3 ððk F ξÞ À1 =λÞ 5=3 gets larger. As a result, fewer gap components get resonantly enhanced, and the temperature variation of the shape of ΔðθÞ weakens. It becomes undetectable at ðk F ξÞ À1 λ 3=5 . This clearly points out that the variation of the gap width θ h with t is a fingerprint of QC pairing.

DISCUSSION
We analyzed s-wave SC, induced by QC nematic fluctuations, near a nematic QCP in 2D. The corresponding gap function ΔðθÞ is peaked at θ ¼ πm=2, where the pairing interaction is the strongest, and has the width θ h $ λ 1=3 , where λ is a dimensionless coupling. Away from a QCP, the structure of ΔðθÞ is set at T c and varies only weakly below T c . We showed that near a QCP the situation is different, and ΔðθÞ has a distinctive temperature evolution within the superconducting state. Namely, the width of the peak strongly grows with decreasing temperature. The source of this evolution is the existence of multiple solutions for the pairing gap Δ n ðθÞ of the same s-wave symmetry, with closely spaced T c;n . These solutions can be thought of as oscillating excited states in an effective trapping potential within the range θ h . The actual T c coincides with T c;0 , and below this temperature the non-oscillating solution Δ 0 ðθÞ develops. However, other Δ n ðθÞ Different colored dots correspond to T c;n , and dashed lines with the same color are the corresponding Δ n . The contribution of Δ n Δ n ðθÞ to ΔðθÞ is resonantly enhanced for T < T c;n . Observe that the sign of Δ n oscillates between even and odd n: Δ n / ðÀ1Þ nþ1 Δ 0 . Because at small θ, Δ n ðθÞ / ðÀ1Þ n Δ 0 ðθÞ, the original (n ¼ 0) and the induced (n > 0) gap components interfere destructively are induced by Δ 0 ðθÞ and each gets resonantly enhanced below T c;n ¼ T c ð1 À OððnλÞ 2=3 ÞÞ. Interference effects from these resonantly induced components modify θ h and make it temperaturedependent. The effect is strong near the nematic QCP and rapidly disappears away from it. We obtained this behavior in a controllable analytical calculation at small λ, but verified numerically that it holds at λ ¼ Oð1Þ. This allows us to estimate temperature ranges where our theory applies. For E F $ 20 À 30 meV, as in QC FeSe 1Àx S x with x $ 0:2 59 we obtain T c $ 5 À 8K for λ ¼ 1, consistent with the experimental T c $ 8K (we caution that this is only an order-of magnitude estimate for T c because of the complex multi-band electronic structure of FeSe 1Àx S x ). We now briefly discuss other mechanisms of anisotropic SC and argue that they do not lead to strong temperature evolution of the gap width. First, an anisotropic gap ΔðθÞ naturally emerges if the pairing is driven by order parameter fluctuations at a finite q, such as spin-or charge-d-wave fluctuations. However, in this case, the gap width θ h is comparable to the width of the hot region already for small coupling λ, and the requirement for closely spaced T c;n does not hold. Second, lattice effects do give rise to some variation of θ h ðtÞ, but this variation is (i) small in T c =ω D and also small in the electron-phonon coupling λ ph , and (ii) does not critically depends on the closeness to a QCP. It was argued 60 that the coupling to the lattice may prevent an electronic system from fully reaching a nematic QCP because the phonon that softens at a nematic QCP only does so in the cold regions. This will not affect our results if λ 2 ) λ ph ω D =E F .
We therefore believe that the observation of the temperature evolution of the gap width will be a "smoking gun" proof of QC SC in a proximity to a nematic QCP.

METHODS
In this section, we give a detailed description of our analytic calculations, both for the linearized and the non-linear gap equation. In the final section, we briefly describe our numerical procedure. In what follows, we first consider the pairing at a QCP, and then away from it. For the latter case, we introduce a finite correlation length ξ ) k À1 F .

The linearized gap equation at a QCP
We consider a generic d-wave form-factor f ðθÞ and focus on the region Àπ=2 < θ < π=2. We use as an input previous works, 36,37 in which the Eliashberg gap equation has been derived. It has the form (Eq. (4) of the main text) Zðθ; ω n ÞΔðθ; ω n Þ ¼ T where the bosonic propagator is As we explained in the main text, the key effect of fermionic Zðθ; ω n Þ is to cancel the term with ω m ¼ ω n in the frequency sum in the r.h.s. of (11). We eliminate this term and thereafter set Zðθ; ω n Þ ¼ 1. This simplifies the analytical consideration. We keep the full Zðθ; ω n Þ contribution in numerical calculations. In order to make manifest the different roles played by the strong local fluctuations at the smallest angular transfers ϕ and the weak tail of the interaction at larger ϕ, it is convenient to split Eq. (11) into two parts. To do so, we subtract Δðθ; ω m Þ from Δðθ þ ϕ; ω m Þ in the numerator of the RHS of Eq. (11) and add it as a separate term. We then obtain Δðθ þ ϕ; ω m Þ À Δðθ; ω m Þ ð Þ f 2 θ þ ϕ 2 À Á jω m j Dðθ; ϕ; ω m À ω n Þ: HereΛ 0 takes care of strong local fluctuations, Dðθ; ϕ; ω m À ω n Þ: (14) Because the integration over ϕ is confined to small jϕj ( 1, we can approximate f ðθ þ ϕ=2Þ by f ðθÞ, and extend the limits of the ϕ integration to ± 1. Integrating over ϕ we then obtain which is Eq. (5) of the main text. The local gap equation Δðθ; ω n Þ ¼Λ 0 ðθ; ω n ÞΔðθ; ω n Þ has a continuous set of solutions Δðθ; ω n Þ ¼ Δ ν ðω n Þδðθ À νÞ (16) which progressively emerge at T c ðνÞ ¼ T c ð0Þð1 À a ν ν 2 þ Þ, where a ν is a constant of order one. The largest T c ð0Þ is obtained by solving (15) with f ðθÞ ¼ f ð0Þ. By order of magnitude, at a QCP, T c ðνÞ $ λ 2 E F f 4 ðνÞ. To get the exact prefactor, we note the gap equation, Eq. (15), with fermionic Z-factor re-introduced in the l.h.s., is equivalent to that for the QC γ model with γ ¼ 1=3 58 . Using the results for the γ model, we find the exact expression The existence of a continuous set of solutions is a consequence of neglecting the second, non-local term in the r.h.s. of Eq. (13). To obtain the actual gap function we need to account for this non-local term. We proceed with an analytic treatment by making two approximations. First, we use the fact that the frequency sum converges, i.e., typical ω m ¼ OðTÞ.
We then neglect the frequency dependence in the gap equation by taking ω n and ω m to be the smallest non-equal Matsubara frequencies, i.e., set ω n ¼ πT; ω m ¼ ÀπT. Second, we expand f ðθÞ and f ðθ þ ϕ=2Þ in small angles near θ ¼ 0; ϕ ¼ 0, i.e., around the center of the hot region. We assume and then verify that the expansion is justified everywhere in the hot region, since we will see that the width of Δðθ; nÞ, viewed as a function of θ, is small, θ h $ λ 1=3 , as long as λ ( 1. Using these approximations, we re-write Eq. (13) as Δðθ þ ϕÞ À ΔðθÞ where we used f ðθÞ ¼ f ð0Þð1 À a θ θ 2 Þ and a λ is a constant of order one. A quick study of Eq. (18) shows that there are two relevant scales for θ: a smaller (lower) scale θ l ¼ OðλÞ, below which ΔðθÞ % Δðθ ¼ 0Þ and a larger (higher) scale θ h ¼ Oðλ 1=3 Þ, which is set by balancing θ 2 ΔðθÞ in the l.h.s. of Eq. (18) and λ R dϕ π ΔðθþϕÞÀΔðθÞ ϕ 2 $ λΔðθÞ=jθj in the r.h.s. We will be interested in the gap function at θ $ θ h ) λ. In this region, one can drop the λ 3 =jϕj term in the r.h.s. of (18) and re-write this equation as Δðθ þ ϕÞ À ΔðθÞ ϕ 2 Rescaling θ by θ h and choosing the prefactor in θ h / λ 1=3 to eliminate the numerical factor between θ 2 and integral terms in the r.h.s. of Eq. (19), we re-write Eq. (19) as Δðθ þ ϕÞ À ΔðθÞ where In original variables, Eq. (20) shows that ΔðθÞ rapidly drops once θ exceeds θ h . Inside the gapped region, at θ < 1, we can transform to Fourier components ΔðθÞ ¼ ð2πÞ À1 R dxe ixθ ΔðxÞ and reduce Eq. (6) to the Airy equation ηΔðxÞ % À∂ 2 x ΔðxÞ þ jxjΔðxÞ This is Eq. (7) from the main text. The boundary condition for the even gap function is Δ 0 ð0Þ ¼ 0. Another boundary condition is Δðx ) 1Þ ! 0. Using asymptotic expressions for the Airy functions we then obtain a discrete set of solutions, specified by integer numbers. The solutions are Δ n ðxÞ ¼ d n Aiðjxj À x n Þ; η n ¼ x n ; Equation (8) of the main text, where Àx n is the position of the nth zero of the derivative of the Airy function Ai 0 ðxÞ, and d n is a constant. This implies that there exists a discrete set of solutions for the gap Δ n ðxÞ with T c;n set by ηðTÞ ¼ η n . The corresponding Δ n ðθÞ changes sign n times in the first octant (8n times over the whole FS). For n ) 1, such that The highest T c ¼ T c;0 is for the solution with n ¼ 0. The corresponding Δ 0 ðθÞ does not change sign.

Non-linear gap equation
In this section, we solve the non-linear gap equation. We extend Eq. (13) by changing the frequency term in the denominator of the RHS to For Δ=T c ( 1, it is enough expand to 3rd order in Δ and keep only the contribution to this order from the operator Δ 0 of Eq. (13). We then repeat the steps to give us an effective equation for angles θ; ϕ like Eq. (6), but with additional terms that we will show below.
Equivalently, we can write a Ginzburg Landau free energy and compute its saddle-point equations. We write the gap function as an expansion in the states Δ n ðθÞ that are solutions of the linearized equation, ΔðθÞ ¼ P n Δ n Δ n ðθÞ. For convenience, we normalize them to be orthonormal and such that Δ n ðθ ¼ 0Þ / ðÀ1Þ n , i.e., non-oscillating Δ 0 is positive for Δðθ ¼ 0Þ, then Δ 1 ðθ ¼ 0Þ is negative, etc. The free energy has the form, where t n ¼ ðT c À T c;n Þ=T c ; t ¼ ðT c À TÞ=T c , and A is a constant of order one. We simplify the analysis of Eq. (24) by keeping only those cross-terms (terms with at least two different Δ n ) that have a power of Δ 0 . The justification for this is that the numerical solution of the gap equation shows no oscillations, implying Δ 0 ) Δ n>0 . In addition, the numerical solution is real, indicating that the coefficients Δ n are real. The result appears in Eq. (9) of the main text, which we reproduce here, Qualitatively, we expect that ja n3 j; ja n1 j ( ja n2 jtja n j. This is because a 2n is an integral over a positive-definite quantity, whereas the integrals that determine a n3 ; a n1 oscillate. In addition, we expect that a n3 ; a n1 should decay rapidly with increasing n because of the oscillations. Finally, we expect that a n3 / ðÀ1Þ n . This is because Δ 3 0 ðθÞ is peaked in the region near θ ¼ 0, and so the overall sign should go as the sign of Δ n ðθ ¼ 0Þ, as long as the integral in Eq. (26) is determined by the peak region θ < 1. Numerically, we find that the overlap integral in Eq. (26) These numbers are consistent with our qualitative arguments. As a n1 $ a n3 and jΔ n j ( jΔ 0 j, we can neglect the a n1 terms in F 4 . We neglect the contribution of the induced states Δ n > 0 on the non-oscillating Δ 0 , in which case the saddle-point equation for Δ 0 immediately yields Δ 0 % ffiffiffiffiffiffiffiffiffi t=a 0 p . The saddle-point equation of Eq. (24) for n > 0 is The results quoted in the main text (after Eq. (9)) correspond to the limits a n3 Δ 3 0 ) a n Δ 3 n and a n3 Δ 3 0 ( a n Δ 3 n . Note, that in the latter limit there are two possible solutions to the equation, with opposite signs. The correct sign is determined by the behavior for small Δ n , i.e., by the sign of Àa n3 . Figure 6 depicts the numerical solution of Eq. (28). Using the solution of Eq. (28) we can compute the t dependence of the width of the hot region, We can simplify Eq. (29) using the following analytic properties of Airy functions: AiðÀx n Þ / ðÀ1Þ n =n 1=4 ; R dxAi 2 ðjxj À x n Þ / ffiffiffiffi ffi x n p $ n 1=3 ; R dxAi ðjxj À x n Þ % const. These in turn imply Δðθ ¼ 0Þ / ðÀ1Þ n =n 1=6 ; Δðx ¼ 0Þ / 1=n 1=6þ1=4¼5=12 . Then we have Δn Δ0 where a; b are constants. The second sum dominates over the first sum (which oscillates rapidly), so the overall trend of θ h ðtÞ is positive. For t ( 1, as long as is not too small, jΔ n =Δ 0 j % a n3 =a 0 ðt=t n Þ $ t=ðnλÞ 2=3 . Then the sum converges rapidly since a n3 decay rapidly. For the smallest λ ¼ 0:03 that we studied numerically in this work, t 1 $ 0:3 and our assumptions are justified. This yields, where c is a constant, in agreement with Eq. (10) of the main text. In Fig. 7, we depict the evolution of θ h ðtÞ À θ h ð0Þ for λ ¼ 0:03; 0:15; 0:3 extracted from the numerical solution of the full Eliashberg equations. The data show good agreement with the predicted scaling. In obtaining Eq. (31) we assumed that the evolution of θ h is a consequence only of the existence of multiple solutions, and neglected the temperature evolution of the characteristic width θ 0 of each independent solution by itself. To further verify that the strong evolution of θ h is a result of multiple solutions, we investigate the evolution of θ 0 ðtÞ for t > 0. The width θ 0 can be extracted numerically from the width of Δ 0 , the non- Fig. 6 An illustration of the numerical solution of Eq. (28). We chose the parameters t n ¼ 0:25; a 0 ¼ 4; a n ¼ 3:14; a 2n ¼ 0:59; a 3n ¼ À0:42, corresponding to n ¼ 1 for λ ¼ 0:03. We chose A ¼ 1 for simplicity oscillating solution of the linearized gap equation, for T < T c . It is straightforward to show, from Eq. (18), that the expected behavior is i.e., θ 0 ðtÞ=θ 0 ð0Þ depends only weakly on temperature, does not depend on λ, and actually goes down, not up. Figure 8 depicts the behavior of θ 0 ðtÞ for several different λ, showing excellent agreement with Eq. (32). This implies that the strong growth of θ h ðtÞ cannot be attributed to just the non-oscillating Δ 0 but comes from multiple solutions.
Away from a critical point The analysis of the discrete set of solutions of the gap equation can be readily extended to a finite nematic correlation length ξ. The bosonic propagator at a finite ξ is Dðθ; ϕ; ΩÞ ¼ λ The computational steps are the same as before. The local gap equation is given by (15) and still has an infinite number of solutions. To study the effect of a finite correlation length, it is useful to consider first the perturbative regime ðk F ξÞ À1 ( λ and then go to the opposite regime ðk F ξÞ À1 ) λ. The operatorΛ 0 ðθ; ω n Þ is given by (14). For ðk F ξÞ À1 ( λ we havê where a ξ ¼ Oð1Þ. In the non-local part of the equation, the finite ξ has no role since the behavior at small ϕ is smooth. Consequently, the gap equation with the non-local term has the same form as Eq. (6), but with a new scaling for η, The new onset temperatures are, T c;n % T c;0 1 À a ξ ðkFξÞ À2 λ 2 1 þ bðλnÞ 2=3 where k F $ k F and λ $ λ. We see that, at ðk F ξÞ À1 ( λ, the effect of the finite correlation length is a uniform reduction of all onset temperatures, but multiple solutions survive and the width of the hot region θ h does not depend on ξ. In the opposite limit ðk F ξÞ À1 ) λ, Λ 0 has a BCS form, Equation (37) shows that the frequency sum is no longer limited to a few first Matsubara frequencies, and instead gives a logarithm and determines an onset temperature T c ð0Þ $ ðE F ðk F ξÞ À3 =λÞ exp Àðk F ξÞ À1 =λ À Á : This logarithm then appears in the non-local part of the gap equation as well, and affects the width of the hotspot θ h . The gap equation with the non-local term becomes ΔðθÞ T c ð0Þ À T T c ð0Þ λ 2ðk F ξÞ À1 / Àa θ θ 2 ΔðθÞ þ ðk F ξÞ À1 Z dϕ Δðθ þ ϕÞ À ΔðθÞ ϕ 2 ; Equation (39) can be recast into the same form as the dimensionless Eq. (6), but now and Equation (41) implies that in order for T c;n to be close in temperature the system must obey ðk F ξÞ À5=3 λ ( 1: Equation (42) demonstrates that the strong evolution of θ h with temperature is a signature of QC pairing. When ðk F ξÞ À1 $ λ 3=5 , the solutions with n > 0 do not develop at finite T=T c , and the SC is of a conventional nature, although it remains strongly anisotropic till k F ξ $ 1.

Details of the numerical solutions
We solved both the linearized and the non-linear gap equation numerically. The details of the numerical solution of the linearized equation have appeared previously in ref. 37 . Figure 4 in this work was obtained for a value of λ ¼ 0:025, with a numerical mesh of 512 points in the range Àπ=2 < θ < π=2, and 48 Matsubara frequencies (half negative and half positive). For Fig. 8, the number of Matsubara freuencies was 24.
Regarding the non-linear gap equation, all the figures and numerical values that appear in this paper are for values of λ ¼ 0:03; 0:15; 0:3. They were obtained with a mesh of 1000 points in the range Àπ < θ < π and 101 Matsubara frequencies, namely we take ω n ¼ ð2n þ 1ÞπT with n ranging from −50 to 50. The numerical solution was obtained by iteration. Both linear and non-linear equations were solved in MATLAB (various versions). To produce Fig. 6 in the Methods section as well as Fig. 5 in the main text, we employed Mathematica 11 and its implementation of the Airy function, with λ ¼ 0:03. We created Fig. 5 using the same parameters as in Fig. 6, except for T c;n , which we estimated from solutions of the Airy equation.