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 density-wave or a nematic order.8,9 Several researchers hinted2,10,11,12,13,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, \(\Delta (\theta ,{\omega }_{m})\) (\(\theta\) is the angle along the FS and \({\omega }_{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 cuprates24,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 \(\Delta (\theta ,{\omega }_{m})\) was expected to be set at \({T}_{{\rm{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 \(\Delta (\theta ,{\omega }_{m})\) below \({T}_{{\rm{c}}}\) at any given \({\omega }_{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}_{{\rm{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}_{{\rm{c}}}\) by solving the linearized gap equations for \(\Delta (\theta ,{\omega }_{m})\). These studies found that \({T}_{{\rm{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}_{{\rm{c}}}\) is the signature of quantum criticality. We show that the mechanism, driving the evolution of \(\Delta (\theta ,{\omega }_{m})\) below \({T}_{{\rm{c}}}\), is the existence of multiple s-wave pairing states with closely spaced transition temperatures \({T}_{{\rm{c}},n}\). These solutions yield \({\Delta }_{n}(\theta ,{\omega }_{m})\) of the same symmetry, but with a different number (\(=8n\)) of sign changes of the gap along the FS. The closeness of \({T}_{{\rm{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 s-wave state with sign-preserving gap \({\Delta }_{0}(\theta ,{\omega }_{m})\), which emerges below the highest \({T}_{{\rm{c}},0}={T}_{{\rm{c}}}\). However, all other \({\Delta }_{n}\) are induced below \({T}_{{\rm{c}}}\), and each gets resonantly enhanced below \({T}_{{\rm{c}},n}\). As a result, the actual \(\Delta (\theta ,{\omega }_{m})\) coincides with \({\Delta }_{0}(\theta ,{\omega }_{m})\) only near \({T}_{{\rm{c}}}\), while at smaller \(T\) the gap structure at any frequency \({\omega }_{m}\) is a mixture of different \({\Delta }_{n}(\theta ,{\omega }_{m})\).

We emphasize that all \({\Delta }_{n}\) have the same symmetry (s-wave). In this respect, the gap evolution below \({T}_{{\rm{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}_{{\rm{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 \(\Delta (\theta ,{\omega }_{m})\) below \({T}_{{\rm{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 \(\lambda\) (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 \(\theta\) on the FS. Still, the pairing interaction does depend on \(\theta\) via the square of the nematic form-factor \({f}^{2}(\theta )={\cos }^{2}2\theta\). As a result, the pairing interaction is larger in ‘hot’ regions near \(\theta =\frac{\pi }{2}m\), and is smaller in ‘lukewarm’ regions near \(\theta =\frac{\pi }{4}+m\frac{\pi }{2}\) (see Fig. 1).36,37 This separation is not sharp as the width of both lukewarm and hot regions is of order one.

Fig. 1
figure 1

The s-wave gap function around the Fermi surface for a superconducting state near a nematic QCP. The gap is sharply peaked in four hot regions surrounding the points \(\theta =\pi m/2\), where the interaction mediated by QC nematic fluctuations is the strongest, and is strongly suppressed around \(\theta =\pi /4+\pi m/2\), where the interaction is the weakest

Full size image

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 \(\delta \theta \sim \lambda\), which is small compared with the width of a hot region. To leading order in \(\lambda\), the gap equation then becomes local and allows a continuous set of solutions \(\Delta (\theta )\propto \delta (\theta -\nu )\) with \({T}_{{\rm{c}}}(\nu )\propto {f}^{4}(\nu )\), where \(\nu\) is a continuous parameter. The physical \({T}_{{\rm{c}}}\) is the highest temperature in the set, and it corresponds to \(\nu =\pi m/2\), where the form-factor \(f(\nu )=1\). The actual gap structure is determined at the next approximation, when one properly accounts for weaker interactions at angle transfers above \(\delta \theta\). The result is that in each octant, e.g., at \(0\ <\ | \theta | \ <\ \pi /4\), the actual gap magnitude is of order \(\Delta (\theta =0)\) for \(| \theta | \ <\ {\theta }_{{\rm{h}}} \sim {\lambda }^{1/3}\) (\(\delta \theta \ll {\theta }_{{\rm{h}}}\ll 1\)) and rapidly drops as \({({\theta }_{{\rm{h}}}/\theta )}^{4}\) at larger \(\theta\). This can be interpreted as if a Cooper pair is ‘trapped’ in a potential well within the range \({\theta }_{{\rm{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 \(\sim \lambda /{\lambda }^{1/3}={\lambda }^{2/3}\). In this situation, the continuous set of \(T_{c}^{\prime} s\) transforms into a discrete set \({T}_{{\rm{c}},n}\), each corresponds to the solution of a Schroedinger-like equation for the gap function in a box with the width \(2{\theta }_{{\rm{h}}}\). The discrete \({T}_{{\rm{c}},n}\) differ from \({T}_{{\rm{c}},0}\) by multiples of the kinetic energy cost

$${T}_{{\rm{c}},n}={T}_{{\rm{c}}}\left(1-O{(n\lambda )}^{2/3}\right),{T}_{{\rm{c}}}={T}_{{\rm{c}},0}$$
(1)

The solutions exist up to \({n}_{\rm{max}} \sim 1/\lambda \gg 1\). The corresponding gap functions \({\Delta }_{n}(\theta )\) 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).

Fig. 2
figure 2

The typical scales in the formation of the SC gap on the FS. In proximity to the QCP, pairing fluctuations (represented by the sharp gray peak) are almost local on the FS, limited to a small angle \(\lambda\), on order of the effective dimensionless coupling. The anisotropy changes on a much larger scale, of order one. The gap forms a region with an intermediate width \(\lambda \ll {\theta }_{{\rm{h}}} \sim {\lambda }^{1/3}\ll 1\). Fermions on the FS bind into Cooper pairs that are localized on the scale \(\lambda\), but can oscillate in the gapped region at the cost of a small kinetic energy. This results in a series of distinct pairing states with closely spaced instability temperatures

Full size image

At \(T\ <\ {T}_{{\rm{c}}}\), the physical gap is a superposition of \({\Delta }_{0}\) and \(\Delta _{n> 0}\), which are all induced by \({\Delta }_{0}\), because all \({\Delta }_{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 \(\Delta _{n> 0}\) are superimposed on the non-oscillating \({\Delta }_{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.

Fig. 3
figure 3

Evolution of the gap anisotropy with temperature. a The 3D plot depicts the pairing gap \(\Delta (t,\theta )\) at \({\omega }_{m}=\pi T\), as a function of the angle along the FS, counted from \(\theta =0\), and reduced temperature \(t=({T}_{{\rm{c}}}-T)/{T}_{{\rm{c}}}\). The gap function has been obtained by numerically solving the full non-linear Eliashberg gap equation at a nematic QCP with coupling \(\lambda =0.03\). \(\Delta\) has been normalized to its maximum value at \(T=0\). b The temperature-dependent width of the gap \({\theta }_{{\rm{h}}}(t)\propto \int \Delta (t,\theta )d\theta\). The gray dashed line is a linear fit, in agreement with Eq. (10)

Full size image

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, \({\Delta }_{n}\) is angle-independent, but oscillates n times as a function of the distance to the FS, \(k-{k}_{{\rm{F}}}\). YS argued that additional solutions affect superconducting stiffness and reduce \({T}_{{\rm{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}_{{\rm{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}_{{\rm{c}}}\)) is at most \(O(1)\).

Model and gap equation

We consider 2D fermions with Fermi energy \({E}_{{\rm{F}}}\), minimally coupled to a nematic order parameter field \({\Delta }^{{\rm{nem}}}(q)\) by

$${H}_{I}=g\sum _{{\rm{k}},{\rm{q}},\sigma }{\Delta }^{{\rm{nem}}}({\rm{q}})f({\rm{k}}){\psi }_{\sigma }^{\dagger }\left({\rm{k}}+\frac{{\rm{q}}}{2}\right){\psi }_{\sigma }\left({\rm{k}}-\frac{{\rm{q}}}{2}\right),$$
(2)

where \(f({\rm{k}})\) is a form-factor, which has d-wave symmetry with respect to \({C}_{4}\) lattice rotations. We assume that the static susceptibility \(\chi (q)\) of the nematic field is peaked at \(q=0\): \(\chi (q)={\chi }_{0}/({\xi }_{0}^{-2}+{q}^{2})\), where \({\xi }_{0}\) is the bare correlation length. We define the dimensionless coupling to be \(\lambda ={g}^{2}{\chi }_{0}/4{E}_{{\rm{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({\rm{k}})\) by \(f(\theta )=\cos 2\theta\) (see Fig. 1). Below we focus on the octant \(0\ <\ \theta \ <\ \pi /4\).

A bosonic excitation with momentum \({\it{q}}\) connects two fermions on the FS with angles \(\theta ,\theta +\phi\) such that \(q=2{k}_{{\rm{F}}}\sin | \phi /2| \approx {k}_{{\rm{F}}}| \phi |\). The effective pairing interaction is then the function of both \(\phi\) (via \(\chi (q)\)) and \(\theta\) (via the form-factor \(f(\theta )\)). This interaction modifies both bosonic and fermionic variables. Fermionic self-energy scales as \(\Sigma ({\omega }_{m}) \sim {({\lambda }^{2}{E}_{{\rm{F}}})}^{1/3}| {\omega }_{m}{| }^{2/3}{\rm{sgn}}({\omega }_{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 self-energy and the pairing vertex53,54,55 and we eliminate it in \(\Sigma ({\omega }_{m})\) and in the gap Eq. (4) below. The bosonic self-energy renormalizes the bare correlation length \({\xi }_{0}\) into the true \(\xi\), which vanishes at the QCP, and also generates a dynamical Landau damping term. The dressed bosonic susceptibility at a QCP is

$$\chi {(\theta ,\phi ,{\Omega }_{m})}^{-1}=\frac{{k}_{F}^{2}}{{\chi }_{0}}\left(| \phi {| }^{2}+\frac{\lambda }{\pi {E}_{{\rm{F}}}}{f}^{2}\left(\theta +\frac{\phi }{2}\right)\left|\frac{{\Omega }_{m}}{\phi }\right|\right).$$
(3)

To obtain \({T}_{{\rm{c}}}\) and the gap function we need to analyze the equation for the pairing vertex \(\Phi (\theta ,{\omega }_{m})\). The pairing gap \(\Delta\) is related to \(\Phi\) in a usual way, as \(\Phi (\theta ,{\omega }_{m})=Z(\theta ,{\omega }_{m})\Delta (\theta ,{\omega }_{m})\), where \(Z(\theta ,{\omega }_{m})=1+\Sigma (\theta /{\omega }_{m})/{\omega }_{m}\) is is the inverse quasiparticle residue.

The discrete set of solutions of the linearized gap equation

We first analyze the gap equation for infinitesimally small \(\Delta (\theta ,{\omega }_{m})\). We explicitly integrate out fermionic momenta transverse to the FS in the gap equation and obtain the equation for the gap function \(\Delta (\theta ,{\omega }_{m})\) on the FS:

$$\begin{array}{rcl}Z(\theta ,{\omega }_{m})\Delta (\theta ,{\omega }_{m})=\frac{\lambda {k}_{F}^{2}}{{\chi }_{0}}T\sum _{m^{\prime} \ne m}\int \frac{d\phi }{2\pi }\frac{\Delta (\theta +\phi ,{\omega }_{m^{\prime} })}{| {\omega }_{m^{\prime} }| }{S}_{m-m^{\prime} }(\theta ,\phi )\\ {S}_{m-m^{\prime} }(\theta ,\phi )={f}^{2}\left(\theta +\frac{\phi }{2}\right)\chi (\theta ,\phi ,{\omega }_{m^{\prime} }-{\omega }_{m}).\end{array}$$
(4)

Equation (4) allows solutions with different gap symmetry. Earlier works found36,37,56 that the s-wave solution has the largest \({T}_{{\rm{c}}}\), so we focus on s-wave gap. Because \(\chi (\theta ,\phi ,{\omega }_{m^{\prime} }-{\omega }_{m})\) is strongly peaked at \(\phi =0\), to first approximation we may set \(\phi =0\) in \(\Delta (\theta +\phi ,{\omega }_{m}){f}^{2}\left(\theta +\frac{\phi }{2}\right)\). This yields a local equation for \(\Delta (\theta ,{\omega }_{m})=\Delta (\theta ,m)\), with \(\theta\) acting as a parameter.38,47,55,57,58 At the QCP we have

$$\Delta (\theta ,m)\approx \frac{{({\lambda }^{2}{E}_{{\rm{F}}}{f}^{4}(\theta ))}^{1/3}}{{(2{T}_{{\rm{c}}})}^{1/3}{3}^{3/2}\pi Z(\theta ,m)}\sum _{m^{\prime} }\frac{\Delta (\theta ,m^{\prime} )}{| m^{\prime} +1/2| }\frac{1}{| m^{\prime} -m{| }^{1/3}}.$$
(5)

Equation (5) has a continuous set of solutions \(\Delta (\theta ,m)\propto \delta (\theta -\nu )\) with arbitrary \(\nu\) from the interval \(0\ <\ \nu \ <\ \pi /4\), and \({T}_{{\rm{c}}}(\nu ) \sim {\lambda }^{2}{E}_{{\rm{F}}}{f}^{4}(\nu )\). The maximum \({T}_{{\rm{c}}}(0)\approx 0.022{\lambda }^{2}{E}_{{\rm{F}}}{f}^{4}(0)\) corresponds to \(\nu =0\). To determine the actual structure of \(\Delta (\theta ,m)\) and the correct number of solutions, we need to go beyond the leading approximation and keep the dependence on \(\phi\) in the numerator in (4). The problem is analytically tractable if we use the fact that typical \({\omega }_{m},{\omega}_{m^{\prime}} \sim T\), i.e., typical \(m,m^{\prime} =O(1)\) and typical \(Z=O(1)\), and simplify the gap equation by neglecting the frequency dependence of \(\Delta (\theta ,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 \(\theta =0\) and obtain

$$\eta (T)\Delta (\bar{\theta })\approx {\bar{\theta }}^{2}\Delta (\bar{\theta })-{\int_{\hskip-2pt-\infty }^{\infty }}\frac{d\bar{\phi }}{\pi }\frac{\Delta (\bar{\theta }+\bar{\phi })-\Delta (\bar{\theta })}{{\bar{\phi }}^{2}}.$$
(6)

(See the Methods section for the detailed derivation.) Here, \(\eta (T) \sim ({({T}_{{\rm{c}}}(0)/T)}^{1/3}-1)/{\theta }_{{\rm{h}}}^{2}\) and \((\bar{\theta },\bar{\phi })={\theta }_{{\rm{h}}}^{-1}(\theta ,\phi )\), where \({\theta }_{{\rm{h}}} \sim {\lambda }^{1/3}\) sets the width of the gap function \(\Delta (\theta )\). Transforming to a Fourier representation \(\Delta (\bar{\theta })={\left(2\pi \right)}^{-1}\int dx\exp ({\rm{i}}x\bar{\theta })\Delta (x)\), we find from (6)

$$\eta (T)\Delta (x)=-{\partial }_{x}^{2}\Delta (x)+| x| \Delta (x).$$
(7)

For even \(\Delta (x)\), which we consider, this reduces to the Airy equation. It has orthogonal solutions

$${\eta }_{n}={x}_{n},\,{\Delta }_{n}(x)\propto {\rm{Ai}}(| x| -{x}_{n}),$$
(8)

where \(-{x}_{n}\) is the nth zero of the derivative of the Airy function \({\rm{Ai}}^{\prime} (x)\). For later convenience, we choose \({\Delta }_{n}(\theta )\) to be orthonormal: \(\int d\theta {\Delta }_{n}(\theta ){\Delta }_{m}(\theta )={\delta }_{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 \({\theta }_{{\rm{h}}}\).

The smallest eigenvalue \({\eta }_{n}\) is for the even solution with \(n=0\). For \(\Delta (\theta )\) in the original coordinates, it yields a non-oscillating gap \({\Delta }_{0}(\theta )\), peaked at \(\theta =0\), with the width \(O({\theta }_{{\rm{h}}})\). The corresponding \({T}_{{\rm{c}},0}\) is the actual \({T}_{{\rm{c}}}\) for the pairing at a nematic QCP. It differs from \({T}_{{\rm{c}}}(0)\) from Eq. (5) by a numerical factor of order \({\lambda }^{2/3}\). For other solutions, \({\Delta }_{n}(\theta )\) changes sign \(n\) times in the first octant. Extracting \({T}_{{\rm{c}},n}\) from Eq. (8), we obtain \({T}_{{\rm{c}},n}={T}_{{\rm{c}},0}\left(1-O{(n\lambda )}^{2/3}\right)\), up to \({n}_{\rm{max}} \sim 1/\lambda \gg 1\). Figure 4 depicts the first few even solutions. Although the detailed form of these solutions depends on \(\lambda\) and \(f(\theta )\), their appearance is the universal feature of QC pairing by nematic fluctuations.

Fig. 4
figure 4

Orthogonal solutions of the linearized gap equation. We used \(\lambda =0.025\) to produce this figure

Full size image

Non-linear gap equation and evolution of \({\theta }_{{\rm{h}}}\)

At \(T\) only slightly below \({T}_{{\rm{c}},0}={T}_{{\rm{c}}}\), the angular dependence of the pairing gap coincides with \({\Delta }_{0}(\theta )\), just the overall gap magnitude increases with decreasing \(T\). To obtain the form of \(\Delta (\theta )\) at lower \(T\), we solve the full non-linear gap equation (see Methods) in the temperature range down to \(T\ll {T}_{{\rm{c}}}\). Figure 3 depicts the evolution of the amplitude and the width of the gap \({\theta }_{{\rm{h}}}(t)\) as a function of the reduced temperature \(t=1-T/{T}_{{\rm{c}}}\). We see from Fig. 3 that \({\theta }_{{\rm{h}}}(t)\) exhibits strong temperature evolution, i.e., \(\Delta (\theta )\) substantially broadens below \({T}_{{\rm{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 \(\Delta (\theta )\) comes about because the width of \(\Delta (\theta )\) is strong influenced by orthogonal components \({\Delta }_{n}(\theta )\). These components have the same symmetry as \({\Delta }_{0}(\theta )\) and are all induced by \({\Delta }_{0}(\theta )\) immediately below \({T}_{{\rm{c}}}\). They are initially small, but the \({\Delta }_{0}(\theta )\) induces the orthogonal gap components with n < 0, and the nth component \({\Delta }_{n}(\theta )\) gets resonantly enhanced at \(t \sim {(n\lambda )}^{2/3}\), and its magnitude becomes comparable to \({\Delta }_{0}(\theta )\). These additional gap components interfere destructively in the region \(| \theta | \ <\ {\theta }_{{\rm{h}}}\), where they oscillate, but constructively for \(\theta\, {>} \,{\theta }_{{\rm{h}}}\). To demonstrate this, we express the gap function as \(\Delta (\theta )={\sum }_{n\ge 0}{\bar{\Delta }}_{n}{\Delta }_{n}(\theta )\) and derive the Landau free energy expansion for the coefficients \({\bar{\Delta }}_{n}\). We obtain

$$\begin{array}{ll}F&=-\frac{1}{2}\sum _{n}(t-{t}_{n}){\bar{\Delta }}_{n}^{2}+\sum _{n}\frac{{a}_{n}}{4}{\bar{\Delta }}_{n}^{4}\\ & + \sum _{n>0}\left[{a}_{n1}{\bar{\Delta }}_{n}^{3}{\bar{\Delta }}_{0}+{a}_{n2}{\bar{\Delta }}_{0}^{2}{\bar{\Delta }}_{n}^{2}+{a}_{n3}{\bar{\Delta }}_{0}^{3}{\bar{\Delta }}_{n}\right]+\ldots \end{array}$$
(9)

where \({t}_{n}=1-{T}_{{\rm{c}},n}/{T}_{{\rm{c}}}\) (\(t-{t}_{n}=({T}_{{\rm{c}},n}-T)/{T}_{{\rm{c}}}\)), \({a}_{nk}\propto \int d\theta {\Delta }_{n}^{4-k}(\theta ){\Delta }_{0}^{k}(\theta )\), and \(\ldots\) stand for less relevant terms in \(F\), e.g., the coupling between two states with \(n\, {>} \,0\).

We stress that the solution for \(\Delta (\theta )\), obtained by minimizing Eq. (9), is equivalent to the solution of the non-linear gap equation. At small \(\theta\), \({\Delta }_{n}(\theta )\propto {(-1)}^{n}\), because the sign at the origin depends on the number of oscillations in each octant. At large \(\theta \gg {\theta }_{{\rm{h}}}\) all \({\Delta }_{n}(\theta )\) decay as \(1/{\theta }^{4}\) (see Methods section). Solving the saddle-point equations, we find that immediately below \({T}_{{\rm{c}}}\), when \(0\ <\ t\ll 1\), \({\bar{\Delta }}_{0}^{2}\approx t/{a}_{0}\) and \(| {\bar{\Delta }}_{n}| \approx | {({\bar{\Delta }}_{0})}^{3}{a}_{n3}/{t}_{n}| \ll {\bar{\Delta }}_{0}\). In the opposite case, when \(t\) is finite and \(t-{t}_{n}\approx t\), we have, using the numerical results for \({a}_{n}\) and \({a}_{ni}\) (\(i=1-3\)) (see Methods section), \({\bar{\Delta }}_{n} \sim {(-1)}^{n+1}{\bar{\Delta }}_{0}\), i.e., \(| {\bar{\Delta }}_{n}|\) is of the same order as \({\bar{\Delta }}_{0}\). Figure 5 illustrates the evolution of \({\bar{\Delta }}_{n}\) with temperature. Because at small \(\theta\), \({\Delta }_{n}(\theta )\propto {(-1)}^{n}{\Delta }_{0}\) and \({\bar{\Delta }}_{n}\propto {(-1)}^{n+1}{\bar{\Delta }}_{0}\), all \(n\, {>} \,0\) contributions to \(\Delta (\theta )\) (the terms \({\bar{\Delta }}_{n}{\Delta }_{n}(\theta )\)) are of opposite sign compared with \({\bar{\Delta }}_{0}{\Delta }_{0}(\theta )\). The original and the induced gap components then interfere destructively, and the total \(\Delta (\theta \approx 0)\) decreases with decreasing \(T\). On the other hand, at \(\theta \gg {\theta }_{{\rm{h}}}\), \({\Delta }_{n}(\theta )\propto {\Delta }_{0}\), hence \({\bar{\Delta }}_{n}{\Delta }_{n}(\theta )\propto {(-1)}^{n+1}\) oscillates in sign between even and odd \(n\). Because \({\bar{\Delta }}_{n}{\Delta }_{n}(\theta )\) decrease in magnitude with increasing \(n\), and \({\bar{\Delta }}_{1}{\Delta }_{1}(\theta )\) has the same sign as \({\bar{\Delta }}_{0}{\Delta }_{0}(\theta )\), the original and the induced gap components then interfere mostly constructively. As the consequence, the effective ‘width’ of the gap, \({\theta }_{{\rm{h}}}(t)\propto \int \Delta (t,\theta )d\theta\), increases with decreasing \(T\). To leading order at \(t\ll 1\), \({\theta }_{{\rm{h}}}\) obeys

$$\frac{{\theta }_{{\rm{h}}}(t)-{\theta }_{{\rm{h}}}(0)}{{\theta }_{{\rm{h}}}(0)}\propto \frac{t}{{\lambda }^{2/3}}$$
(10)

We emphasize that (i) the temperature variation is strong, particularly at small \(\lambda\), and (ii) it is a non-linear effect owing to the existence of multiple sign-changing \({\Delta }_{n}(\theta )\), all having the same s-wave symmetry. The inset of Fig. 3 shows that \({\theta }_{{\rm{h}}}(t)\) indeed increases as \({t}^{\alpha }\) with \(\alpha \approx 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 \(\lambda\) is of order one and \({\theta }_{{\rm{h}}}\) is comparable to the width of the hot region.

Fig. 5
figure 5

Temperature evolution of the induced pairing states. The figure depicts an illustration of the temperature evolution of the induced \({\bar{\Delta }}_{n}\) of \(\Delta (\theta )={\bar{\Delta }}_{0}{\Delta }_{0}(\theta )+{\bar{\Delta }}_{1}{\Delta }_{1}(\theta )+\cdots ,\) where \({\Delta_n}(\theta )\) are orthonormal eigenfunctions of the linearized gap equation. Different colored dots correspond to \({T}_{{\text{c}},n}\), and dashed lines with the same color are the corresponding \({\bar{\Delta }}_{n}\). The contribution of \({\bar{\Delta }}_{n}{\Delta }_{n}(\theta )\) to \(\Delta (\theta )\) is resonantly enhanced for \(T\, {<}\, {T}_{{\text{c}},n}\). Observe that the sign of \({\bar{\Delta }}_{n}\) oscillates between even and odd \(n\): \({\bar{\Delta }}_{n}\propto {(-1)}^{n+1}{\bar{\Delta }}_{0}\). Because at small \(\theta\), \({\Delta }_{n}(\theta )\propto {(-1)}^{n}{\Delta }_{0}(\theta )\), the original (\(n=0\)) and the induced (\(n\, {>} \,0\)) gap components interfere destructively

Full size image

Away from the QCP

Our results are readily generalized (see Methods section) to the case when the correlation length \(\xi\) for nematic fluctuations is large but finite, i.e., the system is at a finite distance from the QCP. For \({({k}_{{\rm{F}}}\xi )}^{-1}\,<\,\lambda\), the effect of \(\xi\) is a uniform reduction of all \({T}_{{\rm{c}},n}\). However, for \({({k}_{{\rm{F}}}\xi )}^{-1}\,>\,\lambda\), we find that \({t}_{n}=1-{T}_{{\rm{c}},n}/{T}_{{\rm{c}}} \sim {(n\lambda )}^{2/3}{({({k}_{{\rm{F}}}\xi )}^{-1}/\lambda )}^{5/3}\) gets larger. As a result, fewer gap components get resonantly enhanced, and the temperature variation of the shape of \(\Delta (\theta )\) weakens. It becomes undetectable at \({({k}_{{\rm{F}}}\xi )}^{-1}\le {\lambda }^{3/5}\). This clearly points out that the variation of the gap width \({\theta }_{{\rm{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 \(\Delta (\theta )\) is peaked at \(\theta =\pi m/2\), where the pairing interaction is the strongest, and has the width \({\theta }_{{\rm{h}}} \sim {\lambda }^{1/3}\), where \(\lambda\) is a dimensionless coupling. Away from a QCP, the structure of \(\Delta (\theta )\) is set at \({T}_{{\rm{c}}}\) and varies only weakly below \({T}_{{\rm{c}}}\). We showed that near a QCP the situation is different, and \(\Delta (\theta )\) 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 \({\Delta }_{n}(\theta )\) of the same s-wave symmetry, with closely spaced \({T}_{{\rm{c}},n}\). These solutions can be thought of as oscillating excited states in an effective trapping potential within the range \({\theta }_{{\rm{h}}}\). The actual \({T}_{{\rm{c}}}\) coincides with \({T}_{{\rm{c}},0}\), and below this temperature the non-oscillating solution \({\Delta }_{0}(\theta )\) develops. However, other \({\Delta }_{n}(\theta )\) are induced by \({\Delta }_{0}(\theta )\) and each gets resonantly enhanced below \({T}_{{\rm{c}},n}={T}_{{\rm{c}}}(1-O({(n\lambda )}^{2/3}))\). Interference effects from these resonantly induced components modify \({\theta }_{{\rm{h}}}\) and make it temperature-dependent. 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 \(\lambda\), but verified numerically that it holds at \(\lambda =O(1)\). This allows us to estimate temperature ranges where our theory applies. For \({E}_{{\rm{F}}} \sim 20-30\,{\rm{meV}}\), as in QC FeSe\({}_{1-x}\)S\({}_{x}\) with \(x \sim 0.2\)59 we obtain \({T}_{{\rm{c}}} \sim 5-8K\) for \(\lambda =1\), consistent with the experimental \({T}_{{\rm{c}}} \sim 8K\) (we caution that this is only an order-of magnitude estimate for \({T}_{{\rm{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 \(\Delta (\theta )\) 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 \({\theta }_{{\rm{h}}}\) is comparable to the width of the hot region already for small coupling \(\lambda\), and the requirement for closely spaced \({T}_{{\rm{c}},n}\) does not hold. Second, lattice effects do give rise to some variation of \({\theta }_{{\rm{h}}}(t)\), but this variation is (i) small in \({T}_{{\rm{c}}}/{\omega }_{{\rm{D}}}\) and also small in the electron-phonon coupling \({\lambda }_{{\rm{ph}}}\), and (ii) does not critically depends on the closeness to a QCP. It was argued60 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 \({\lambda }^{2}\gg {\lambda }_{{\rm{ph}}}{\omega }_{{\rm{D}}}/{E}_{{\rm{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 \(\xi \gg {k}_{{\rm{F}}}^{-1}\).

The linearized gap equation at a QCP

We consider a generic d-wave form-factor \(f(\theta )\) and focus on the region \(-\pi /2\ <\ \theta \ <\ \pi /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(\theta ,{\omega }_{n})\Delta (\theta ,{\omega }_{n})=T\sum _{{\omega }_{m}}\int \frac{d\phi }{2\pi }\frac{\Delta (\theta +\phi ,{\omega }_{m}){f}^{2}\left(\theta +\frac{\phi }{2}\right)}{| {\omega }_{m}| }D(\theta ,\phi ,{\omega }_{m}-{\omega }_{n}),$$
(11)

where the bosonic propagator is

$$D(\theta ,\phi ,\Omega )=\frac{\lambda }{{\phi }^{2}+{f}^{2}\left(\theta +\frac{\phi }{2}\right)\frac{\lambda | \Omega | }{{E}_{{\rm{F}}}\pi | \phi | }}.$$
(12)

As we explained in the main text, the key effect of fermionic \(Z(\theta ,{\omega }_{n})\) is to cancel the term with \({\omega }_{m}={\omega }_{n}\) in the frequency sum in the r.h.s. of (11). We eliminate this term and thereafter set \(Z(\theta ,{\omega }_{n})=1\). This simplifies the analytical consideration. We keep the full \(Z(\theta ,{\omega }_{n})\) contribution in numerical calculations.

In order to make manifest the different roles played by the strong local fluctuations at the smallest angular transfers \(\phi\) and the weak tail of the interaction at larger \(\phi\), it is convenient to split Eq. (11) into two parts. To do so, we subtract \(\Delta (\theta ,{\omega }_{m})\) from \(\Delta (\theta +\phi ,{\omega }_{m})\) in the numerator of the RHS of Eq. (11) and add it as a separate term. We then obtain

$$\Delta (\theta ,{\omega }_{n})=\,{\hat{\Lambda }}_{0}\Delta (\theta ,{\omega }_{n})+T\sum _{{\omega }_{m}\ne {\omega }_{n}}\int \frac{d\phi }{2\pi }\frac{\left(\Delta (\theta +\phi ,{\omega }_{m})-\Delta (\theta ,{\omega }_{m})\right){f}^{2}\left(\theta +\frac{\phi }{2}\right)}{| {\omega }_{m}| }D(\theta ,\phi ,{\omega }_{m}-{\omega }_{n}).$$
(13)

Here \({\hat{\Lambda }}_{0}\) takes care of strong local fluctuations,

$${\hat{\Lambda }}_{0}\Delta (\theta ,{\omega }_{n})=T\sum _{{\omega }_{m}\ne {\omega }_{n}}\frac{\Delta (\theta ,{\omega }_{m})}{| {\omega }_{m}| }\times \int \frac{d\phi }{2\pi }{f}^{2}\left(\theta +\frac{\phi }{2}\right)D(\theta ,\phi ,{\omega }_{m}-{\omega }_{n}).$$
(14)

Because the integration over \(\phi\) is confined to small \(| \phi | \ll 1\), we can approximate \(f(\theta +\phi /2)\) by \(f(\theta )\), and extend the limits of the \(\phi\) integration to \(\pm \infty\). Integrating over \(\phi\) we then obtain

$${\hat{\Lambda }}_{0}\Delta (\theta ,{\omega }_{n})\approx \frac{2}{{3}^{3/2}}T\sum _{m\ne n}\frac{\Delta (\theta ,{\omega }_{m})}{| {\omega }_{m}| }\frac{1}{{\left(\frac{| {\omega }_{m}-{\omega }_{n}| }{\pi {\lambda }^{2}{E}_{{\rm{F}}}{f}^{4}(\theta )}\right)}^{1/3}},$$

which is Eq. (5) of the main text.

The local gap equation

$$\Delta (\theta ,{\omega }_{n})={\hat{\Lambda }}_{0}(\theta ,{\omega }_{n})\Delta (\theta ,{\omega }_{n})$$
(15)

has a continuous set of solutions

$$\Delta (\theta ,{\omega }_{n})={\Delta }_{\nu }({\omega }_{n})\delta (\theta -\nu )$$
(16)

which progressively emerge at \({T}_{{\rm{c}}}(\nu )={T}_{{\rm{c}}}(0)(1-{a}_{\nu }{\nu }^{2}+\ldots )\), where \({a}_{\nu }\) is a constant of order one. The largest \({T}_{{\rm{c}}}(0)\) is obtained by solving (15) with \(f(\theta )=f(0)\). By order of magnitude, at a QCP, \({T}_{{\rm{c}}}(\nu ) \sim {\lambda }^{2}{E}_{{\rm{F}}}{f}^{4}(\nu )\). 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 \(\gamma\) model with \(\gamma =1/3\)58. Using the results for the \(\gamma\) model, we find the exact expression

$${T}_{{\rm{c}}}(0)=0.022{\lambda }^{2}{E}_{{\rm{F}}}{[f(0)]}^{4}$$
(17)

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 \({\omega }_{m}=O(T)\). We then neglect the frequency dependence in the gap equation by taking \({\omega }_{n}\) and \({\omega }_{m}\) to be the smallest non-equal Matsubara frequencies, i.e., set \({\omega }_{n}=\pi T,{\omega }_{m}=-\pi T\). Second, we expand \(f(\theta )\) and \(f(\theta +\phi /2)\) in small angles near \(\theta =0,\phi =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 \(\Delta (\theta ,n)\), viewed as a function of \(\theta\), is small, \({\theta }_{{\rm{h}}} \sim {\lambda }^{1/3}\), as long as \(\lambda \ll 1\).

Using these approximations, we re-write Eq. (13) as

$$\Delta (\theta )\left[1-{\left(1-{a}_{\theta }{\theta }^{2}\right)}^{4/3}{\left(\frac{{T}_{{\rm{c}}}(0)}{T}\right)}^{1/3}\right]=\frac{\lambda {f}^{2}(0)}{\pi }\int \frac{d\phi }{2\pi }\frac{\Delta (\theta +\phi )-\Delta (\theta )}{{\phi }^{2}+{a}_{\lambda }\frac{{\lambda }^{3}}{| \phi | }}$$
(18)

where we used \(f(\theta )=f(0)(1-{a}_{\theta }{\theta }^{2})\) and \({a}_{\lambda }\) is a constant of order one.

A quick study of Eq. (18) shows that there are two relevant scales for \(\theta\): a smaller (lower) scale \({\theta }_{l}=O(\lambda )\), below which \(\Delta (\theta )\approx \Delta (\theta =0)\) and a larger (higher) scale \({\theta }_{{\rm{h}}}=O({\lambda }^{1/3})\), which is set by balancing \({\theta }^{2}\Delta (\theta )\) in the l.h.s. of Eq. (18) and \(\lambda \int \frac{d\phi }{\pi }\frac{\Delta (\theta +\phi )-\Delta (\theta )}{{\phi }^{2}} \sim \lambda \Delta (\theta )/| \theta |\) in the r.h.s. We will be interested in the gap function at \(\theta \sim {\theta }_{{\rm{h}}}\gg \lambda\). In this region, one can drop the \({\lambda }^{3}/| \phi |\) term in the r.h.s. of (18) and re-write this equation as

$$\Delta (\theta )\left[1-{\left(\frac{{T}_{{\rm{c}}}(0)}{T}\right)}^{1/3}\right]=-\frac{4}{3}{a}_{\theta }{\theta }^{2}\Delta (\theta )+\frac{\lambda {f}^{2}(0)}{\pi }\int \frac{d\phi }{2\pi }\frac{\Delta (\theta +\phi )-\Delta (\theta )}{{\phi }^{2}}$$
(19)

Rescaling \(\theta\) by \({\theta }_{{\rm{h}}}\) and choosing the prefactor in \({\theta }_{{\rm{h}}}\propto {\lambda }^{1/3}\) to eliminate the numerical factor between \({\theta }^{2}\) and integral terms in the r.h.s. of Eq. (19), we re-write Eq. (19) as

$$\eta (T)\Delta (\bar{\theta })={\bar{\theta }}^{2}\Delta (\bar{\theta })-{\int }_{-\infty }^{\infty }\frac{d\bar{\phi }}{\pi }\frac{\Delta (\bar{\theta }+\bar{\phi })-\Delta (\bar{\theta })}{{\bar{\phi }}^{2}}$$

where \(\bar{\theta },\bar{\phi }\) are rescaled angles and we again extended the integration to \(\pm \infty\). This is Eq. (6) of the main text. The width of the \(\Delta\) in Eq. (6) (defined as \(| \theta | \le {\theta }_{{\rm{h}}}\)) is of order one in the rescaled units. The parameter \(\eta \sim ({({T}_{{\rm{c}}}(0)/T)}^{1/3}-1)/{\lambda }^{2/3}\). For \(| T-{T}_{{\rm{c}}}(0)| \ll {T}_{{\rm{c}}}(0)\), \(\eta \sim ({T}_{{\rm{c}}}(0)-T)/({T}_{{\rm{c}}}(0){\lambda }^{2/3})\).

To solve Eq. (6) we treat separately the behavior in and out of the hot region, i.e., at \(\bar{\theta }\ll 1\) and \(\bar{\theta }\gg 1\). For \(\bar{\theta }\gg 1\), we may neglect the \(\Delta (\bar{\theta })\) term in the integrand. The largest contribution to the integral comes from the region where the gap is large, i.e., from \(\bar{\phi } \sim -\bar{\theta }\), with the width of order \(\bar{\phi }=O(1)\), so to leading order we have

$$\Delta (\bar{\theta }\gg 1)\approx \frac{\bar{\Delta }}{{\bar{\theta }}^{4}},$$
(20)

where

$$\bar{\Delta }=\int \frac{dx}{\pi }\Delta (x)$$
(21)

In original variables, Eq. (20) shows that \(\Delta (\theta )\) rapidly drops once \(\theta\) exceeds \({\theta }_{{\rm{h}}}\). Inside the gapped region, at \(\bar{\theta }\ <\ 1\), we can transform to Fourier components \(\Delta (\bar{\theta })={(2\pi )}^{-1}\int dx{{\rm{e}}}^{{\rm{i}}x\bar{\theta }}\Delta (x)\) and reduce Eq. (6) to the Airy equation

$$\eta \Delta (x)\approx -{\partial }_{x}^{2}\Delta (x)+| x| \Delta (x)$$

This is Eq. (7) from the main text. The boundary condition for the even gap function is \(\Delta ^{\prime} (0)=0\). Another boundary condition is \(\Delta (x\gg 1)\to 0\). Using asymptotic expressions for the Airy functions we then obtain a discrete set of solutions, specified by integer numbers. The solutions are

$${\Delta }_{n}(x)={d}_{n}{\rm{Ai}}(| x| -{x}_{n}),{\eta }_{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 \({\rm{Ai}}^{\prime} (x)\), and \({d}_{n}\) is a constant. This implies that there exists a discrete set of solutions for the gap \({\Delta }_{n}(x)\) with \({T}_{{\rm{c}},n}\) set by \(\eta (T)={\eta }_{n}\). The corresponding \({\Delta }_{n}(\bar{\theta })\) changes sign \(n\) times in the first octant (\(8n\) times over the whole FS). For \(n\gg 1\),

$$\frac{2}{3}{\eta }_{n}^{3/2}\approx \frac{\pi }{4}+n\pi ,$$
(22)

such that

$${T}_{{\rm{c}},n}={T}_{{\rm{c}}}(0)\left(1-O{(\lambda n)}^{2/3}\right)$$
(23)

The highest \({T}_{{\rm{c}}}={T}_{{\rm{c}},0}\) is for the solution with \(n=0\). The corresponding \({\Delta }_{0}(\bar{\theta })\) 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 \(| {\omega }_{m}| \to \sqrt{{\omega }_{m}^{2}+{\Delta }^{2}(\theta +\phi ,{\omega }_{m})}\). For \(\Delta /{T}_{{\rm{c}}}\ll 1\), it is enough expand to 3rd order in \(\Delta\) and keep only the contribution to this order from the operator \({\hat{\Delta }}_{0}\) of Eq. (13). We then repeat the steps to give us an effective equation for angles \(\bar{\theta },\bar{\phi }\) 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 \({\Delta }_{n}(\bar{\theta })\) that are solutions of the linearized equation, \(\Delta (\bar{\theta })={\sum }_{n}{\bar{\Delta }}_{n}{\Delta }_{n}(\bar{\theta })\). For convenience, we normalize them to be orthonormal and such that \({\Delta }_{n}(\bar{\theta }=0)\propto {(-1)}^{n}\), i.e., non-oscillating \({\Delta }_{0}\) is positive for \(\Delta (\bar{\theta }=0)\), then \({\Delta }_{1}(\bar{\theta }=0)\) is negative, etc. The free energy has the form,

$$F\propto \int -\frac{1}{2}\sum _{n}(t-{t}_{n}){\bar{\Delta }}_{n}^{2}+{F}_{4}+\cdots$$
(24)

where \({t}_{n}=({T}_{{\rm{c}}}-{T}_{{\rm{c}},n})/{T}_{{\rm{c}}},t=({T}_{{\rm{c}}}-T)/{T}_{{\rm{c}}}\),

$${F}_{4}=A\sum _{ijkl}{\bar{\Delta }}_{i}{\bar{\Delta }}_{j}{\bar{\Delta }}_{k}{\bar{\Delta }}_{l}\int d\bar{\theta }{\Delta }_{i}(\bar{\theta }){\Delta }_{j}(\bar{\theta }){\Delta }_{k}(\bar{\theta }){\Delta }_{l}(\bar{\theta }),$$
(25)

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 \({\bar{\Delta }}_{n}\)) that have a power of \({\bar{\Delta }}_{0}\). The justification for this is that the numerical solution of the gap equation shows no oscillations, implying \({\bar{\Delta }}_{0}\gg {\bar{\Delta }}_{n{>}0}\). In addition, the numerical solution is real, indicating that the coefficients \({\bar{\Delta }}_{n}\) are real. The result appears in Eq. (9) of the main text, which we reproduce here,

$${F}_{4}\approx \sum _{n}\frac{{a}_{n}}{4}{\bar{\Delta }}_{n}^{4}+\sum _{n {>} 0}\left[{a}_{n1}{\bar{\Delta }}_{n}^{3}{\bar{\Delta }}_{0}+{a}_{n2}{\bar{\Delta }}_{0}^{2}{\bar{\Delta }}_{n}^{2}+{a}_{n3}{\bar{\Delta }}_{0}^{3}{\bar{\Delta }}_{n}\right].$$
(26)

Qualitatively, we expect that \(| {a}_{n3}| ,| {a}_{n1}| \ll | {a}_{n2}| \lesssim | {a}_{n}|\). 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}\propto {(-1)}^{n}\). This is because \({\Delta }_{0}^{3}(\bar{\theta })\) is peaked in the region near \(\bar{\theta }=0\), and so the overall sign should go as the sign of \({\Delta }_{n}(\bar{\theta }=0)\), as long as the integral in Eq. (26) is determined by the peak region \(\bar{\theta }\ <\ 1\). Numerically, we find that the overlap integral in Eq. (26) yields (up to the constant \(A\)), for \(\lambda =0.03\),

$$\begin{array}{rcl}{a}_{n}&=&4,3.14,2.77,2.56,\ldots \qquad (n\;\ge\;0)\\ {a}_{n1}&=&-0.1,-0.04,-0.02,\ldots \qquad (n\,>\,0)\\ {a}_{n2}&=&0.59,0.49,0.43,\ldots \qquad (n\,>\,0)\\ {a}_{n3}&=&-0.42,0.06,-0.03,\ldots \qquad (n\,>\,0)\end{array}$$
(27)

These numbers are consistent with our qualitative arguments. As \({a}_{n1} \sim {a}_{n3}\) and \(| {\bar{\Delta }}_{n}| \ll | {\bar{\Delta }}_{0}|\), we can neglect the \({a}_{n1}\) terms in \({F}_{4}\). We neglect the contribution of the induced states \({\Delta }_n\,{>}\,0\) on the non-oscillating \({\Delta }_{0}\), in which case the saddle-point equation for \({\bar{\Delta }}_{0}\) immediately yields \({\bar{\Delta }}_{0}\approx \sqrt{t/{a}_{0}}\). The saddle-point equation of Eq. (24) for \(n\, {>} \,0\) is

$${\bar{\Delta }}_{n}\approx -\frac{{a}_{n3}{\bar{\Delta }}_{0}^{3}+{a}_{n}{\bar{\Delta }}_{n}^{3}}{{t}_{n}+2{a}_{n2}{\bar{\Delta }}_{0}^{2}-t}.$$
(28)

The results quoted in the main text (after Eq. (9)) correspond to the limits \({a}_{n3}{\bar{\Delta }}_{0}^{3}\gg {a}_{n}{\bar{\Delta }}_{n}^{3}\) and \({a}_{n3}{\bar{\Delta }}_{0}^{3}\ll {a}_{n}{\bar{\Delta }}_{n}^{3}\). 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 \({\bar{\Delta }}_{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,

$$\begin{array}{rcl}{\theta }_{{\rm{h}}}(t)\propto \int d\bar{\phi }\frac{\Delta (\bar{\phi })}{\Delta (\bar{\theta }=0)}\approx \frac{{\sum }_{n}{\bar{\Delta }}_{n}\int d\bar{\phi }{\Delta }_{n}(\bar{\phi })}{{\sum }_{n}{\bar{\Delta }}_{n}{\Delta }_{n}(\bar{\theta }=0)}\\ \approx \frac{{\sum }_{n}{\bar{\Delta }}_{n}{\Delta }_{n}(x=0)}{{\sum }_{n}{\bar{\Delta }}_{n}{\Delta }_{n}(\bar{\theta }=0)}\end{array}$$
(29)

We can simplify Eq. (29) using the following analytic properties of Airy functions: \({\rm{Ai}}(-{x}_{n})\propto {(-1)}^{n}/{n}^{1/4},\int dx{{\rm{Ai}}}^{2}(| x| -{x}_{n})\propto \sqrt{{x}_{n}} \sim {n}^{1/3},\int dx{\rm{Ai}}(| x| -{x}_{n})\approx {\rm{const}}\). These in turn imply \(\Delta (\bar{\theta }=0)\propto {(-1)}^{n}/{n}^{1/6},\Delta (x=0)\propto 1/{n}^{1/6+1/4=5/12}\). Then we have

$${\theta }_{{\rm{h}}}(t)\propto \frac{1-a{\sum }_{n>0}\left|\frac{{\bar{\Delta }}_{n}}{{\bar{\Delta }}_{0}}\right|\frac{{(-1)}^{n+1}}{{n}^{5/12}}}{1-b{\sum }_{n>0}\left|\frac{{\bar{\Delta }}_{n}}{{\bar{\Delta }}_{0}}\right|\frac{1}{{n}^{1/6}}}\approx 1-a\sum _{n>0}\left|\frac{{\bar{\Delta }}_{n}}{{\bar{\Delta }}_{0}}\right|\frac{{(-1)}^{n+1}}{{n}^{5/12}}+b\sum _{n>0}\left|\frac{{\bar{\Delta }}_{n}}{{\bar{\Delta }}_{0}}\right|\frac{1}{{n}^{1/6}}$$
(30)

where \(a,b\) are constants. The second sum dominates over the first sum (which oscillates rapidly), so the overall trend of \({\theta }_{{\rm{h}}}(t)\) is positive. For \(t\ll 1\), as long as \({t}_{1} \sim ({\eta }_{1}-{\eta }_{0}){\lambda }^{2/3}\) is not too small, \(| {\bar{\Delta }}_{n}/{\bar{\Delta }}_{0}| \approx {a}_{n3}/{a}_{0}(t/{t}_{n}) \sim t/{(n\lambda )}^{2/3}\). Then the sum converges rapidly since \({a}_{n3}\) decay rapidly. For the smallest \(\lambda =0.03\) that we studied numerically in this work, \({t}_{1} \sim 0.3\) and our assumptions are justified. This yields,

$${\theta }_{{\rm{h}}}(t)\propto 1+\frac{c}{{\lambda }^{2/3}}t,$$
(31)

where \(c\) is a constant, in agreement with Eq. (10) of the main text. In Fig. 7, we depict the evolution of \({\theta }_{{\rm{h}}}(t)-{\theta }_{{\rm{h}}}(0)\) for \(\lambda =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.

Fig. 6
figure 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 \(\lambda =0.03\). We chose \(A=1\) for simplicity

Full size image
Fig. 7
figure 7

Evolution and scaling of \({\theta }_{{\rm{h}}}(t)\) for different \(\lambda\)

Full size image

In obtaining Eq. (31) we assumed that the evolution of \({\theta }_{{\rm{h}}}\) is a consequence only of the existence of multiple solutions, and neglected the temperature evolution of the characteristic width \({\theta }_{0}\) of each independent solution by itself. To further verify that the strong evolution of \({\theta }_{{\rm{h}}}\) is a result of multiple solutions, we investigate the evolution of \({\theta }_{0}(t)\) for t > 0. The width \({\theta }_{0}\) can be extracted numerically from the width of \({\Delta }_{0}\), the non-oscillating solution of the linearized gap equation, for \(T\ <\ {T}_{{\rm{c}}}\). It is straightforward to show, from Eq. (18), that the expected behavior is

$${\theta }_{0}(T)/{\theta }_{0}(0)={(T/{T}_{{\rm{c}}})}^{1/9},$$
(32)

i.e., \({\theta }_{0}(t)/{\theta }_{0}(0)\) depends only weakly on temperature, does not depend on \(\lambda\), and actually goes down, not up. Figure 8 depicts the behavior of \({\theta }_{0}(t)\) for several different \(\lambda\), showing excellent agreement with Eq. (32). This implies that the strong growth of \({\theta }_{{\rm{h}}}(t)\) cannot be attributed to just the non-oscillating \({\Delta }_{0}\) but comes from multiple solutions.

Fig. 8
figure 8

Evolution of the width of the non-oscillating solution of the linearized gap equation with temperature. The dashed line is the prediction of Eq. (32)

Full size image

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 \(\xi\). The bosonic propagator at a finite \(\xi\) is

$$D(\theta ,\phi ,\Omega )=\frac{\lambda }{{\phi }^{2}+{({k}_{{\rm{F}}}\xi )}^{-2}+{f}^{2}\left(\theta +\frac{\phi }{2}\right)\frac{\lambda | \Omega | }{{E}_{{\rm{F}}}\pi | \phi | }}.$$
(33)

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}_{{\rm{F}}}\xi )}^{-1}\ll \lambda\) and then go to the opposite regime \({({k}_{{\rm{F}}}\xi )}^{-1}\gg \lambda\). The operator \({\hat{\Lambda }}_{0}(\theta ,{\omega }_{n})\) is given by (14). For \({({k}_{{\rm{F}}}\xi )}^{-1}\ll \lambda\) we have

$${\hat{\Lambda }}_{0}(\theta )\Delta (\theta ,n)\approx {\left(\frac{{T}_{{\rm{c}}}(0)}{T}\right)}^{1/3}\left(1-{a}_{\xi }\frac{{({k}_{{\rm{F}}}\xi )}^{-2}}{{\lambda }^{2}}\right)\Delta (\theta ,n)\,-\,{a}_{\theta }{\theta }^{2}\Delta (\theta ,n)$$
(34)

where \({a}_{\xi }=O(1)\). In the non-local part of the equation, the finite \(\xi\) has no role since the behavior at small \(\phi\) is smooth. Consequently, the gap equation with the non-local term has the same form as Eq. (6), but with a new scaling for \(\eta\),

$$\eta {\lambda }^{2/3}\propto {\left(\frac{{T}_{{\rm{c}}}(0)}{T}\right)}^{1/3}\left(1-{a}_{\xi }\frac{{({k}_{{\rm{F}}}\xi )}^{-2}}{{\lambda }^{2}}\right)-1.$$
(35)

The new onset temperatures are,

$$\begin{array}{rcl}{T}_{{\rm{c}},n}\approx {T}_{{\rm{c}},0}{\left[\frac{1-{a}_{\xi }\frac{{({k}_{{\rm{F}}}\xi )}^{-2}}{{\lambda }^{2}}}{1+b{(\lambda n)}^{2/3}}\right]}^{3}\\ \approx {T}_{{\rm{c}}}(0)\left(1-{\left(\frac{{\bar{k}}_{F}\xi }{\lambda }\right)}^{2}\right)(1-{(\bar{\lambda }n)}^{2/3}),\end{array}$$
(36)

where \({\bar{k}}_{F} \sim {k}_{{\rm{F}}}\) and \(\bar{\lambda } \sim \lambda\). We see that, at \({({k}_{{\rm{F}}}\xi )}^{-1}\ll \lambda\), 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 \({\theta }_{{\rm{h}}}\) does not depend on \(\xi\).

In the opposite limit \({({k}_{{\rm{F}}}\xi )}^{-1}\gg \lambda\), \({\Lambda }_{0}\) has a BCS form,

$$\begin{array}{ll}&{\hat{\Lambda }}_{0}(\theta ,{\omega }_{n})\Delta (\theta ,{\omega }_{n})\approx \frac{\lambda {f}^{2}(\theta )}{2{({k}_{{\rm{F}}}\xi )}^{-1}}T\sum _{m}^{\left\lfloor \frac{{E}_{{\rm{F}}}{({k}_{{\rm{F}}}\xi )}^{-3}}{\lambda T}\right\rfloor }\frac{\Delta (\theta ,{\omega }_{m})}{| {\omega }_{m}| }\\ &\approx \frac{\lambda }{2{({k}_{{\rm{F}}}\xi )}^{-1}}\mathrm{log}\left(\frac{{E}_{{\rm{F}}}{({k}_{{\rm{F}}}\xi )}^{-3}}{\lambda T}\right)\Delta (\theta ,{\omega }_{n})-{a}_{\theta }{\theta }^{2}\Delta (\theta ,{\omega }_{n})\end{array}$$
(37)

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

$${\bar{T}}_{c}(0) \sim ({E}_{{\rm{F}}}{({k}_{{\rm{F}}}\xi )}^{-3}/\lambda )\exp \left(-{({k}_{{\rm{F}}}\xi )}^{-1}/\lambda \right).$$
(38)

This logarithm then appears in the non-local part of the gap equation as well, and affects the width of the hotspot \({\theta }_{{\rm{h}}}\). The gap equation with the non-local term becomes

$$\Delta (\theta )\frac{{\bar{T}}_{{\rm{c}}}(0)-T}{{\bar{T}}_{{\rm{c}}}(0)}\frac{\lambda }{2{({k}_{{\rm{F}}}\xi )}^{-1}}\,\propto\,-{a}_{\theta }{\theta }^{2}\Delta (\theta )+{({k}_{{\rm{F}}}\xi )}^{-1}\int d\phi \frac{\Delta (\theta +\phi )-\Delta (\theta )}{{\phi }^{2}},$$
(39)

Equation (39) can be recast into the same form as the dimensionless Eq. (6), but now

$${\theta }_{{\rm{h}}}\propto {({k}_{{\rm{F}}}\xi )}^{-1/3}$$
(40)

and

$$\eta \propto \frac{{\bar{T}}_{{\rm{c}}}-T}{{\bar{T}}_{{\rm{c}}}}\frac{\lambda }{{({k}_{{\rm{F}}}\xi )}^{-5/3}}$$
(41)

Equation (41) implies that in order for \({T}_{{\rm{c}},n}\) to be close in temperature the system must obey

$$\frac{{({k}_{{\rm{F}}}\xi )}^{-5/3}}{\lambda }\ll 1.$$
(42)

Equation (42) demonstrates that the strong evolution of \({\theta }_{{\rm{h}}}\) with temperature is a signature of QC pairing. When \({({k}_{{\rm{F}}}\xi )}^{-1} \sim {\lambda }^{3/5}\), the solutions with \(n\, {>} \,0\) do not develop at finite \(T/{T}_{{\rm{c}}}\), and the SC is of a conventional nature, although it remains strongly anisotropic till \({k}_{{\rm{F}}}\xi \sim 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 \(\lambda =0.025\), with a numerical mesh of \(512\) points in the range \(-\pi /2\ <\ \theta \ <\ \pi /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 \(\lambda =0.03,0.15,0.3\). They were obtained with a mesh of 1000 points in the range \(-\pi \ <\ \theta \ <\ \pi\) and 101 Matsubara frequencies, namely we take \({\omega }_{n}=(2n+1)\pi 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 \(\lambda =0.03\). We created Fig. 5 using the same parameters as in Fig. 6, except for \({T}_{{\rm{c}},n}\), which we estimated from solutions of the Airy equation.

Acknowledgements

We thank E. Berg, R. Fernandes, S. Kivelson, M.N. Gastiasoro, S. Lederer, and Y. Schattner for stimulating discussions. We thank the Minnesota Supercomputing Institute at the University of Minnesota for providing resources that assisted with this work. This work was supported by the NSF DMR-1523036.