Introduction

It is widely believed that superconductivity in the cuprates, Fe-pnictides, heavy fermion, and other correlated electron systems is of electronic origin and at least in some portion of the phase diagram can be understood as mediated by soft fluctuations of a particle-hole order parameter, which is about to condense. The most studied scenario of this kind is pairing mediated by spin fluctuations. For the cuprates, it naturally leads to \({d}_{{{{{\rm{x}}}}}^{2}-{{{{\rm{y}}}}}^{2}}\) pairing. For Fe-pnictides, spin-mediated pairing interaction is attractive in both s-wave (s+−) and \({d}_{{{{{\rm{x}}}}}^{2}-{{{{\rm{y}}}}}^{2}}\) channels. The argument, why pairing holds despite that the electron-electron interaction is repulsive, is the same in the two cases—antiferromagnetic spin fluctuations, peaked at momentum Q, increase the magnitude of a repulsive pairing interaction at the momentum transfer Q (the pair hopping from (k, − k) to k + Q, − k − Q). A repulsive pair hopping allows for a solution for a gap function, which changes sign between Fermi points at kF and kF + Q. There is still a repulsion at small momentum transfer, which is detrimental to any superconductivity, and the bare Coulomb interaction is indeed larger at small momenta than at Q. However, when spin fluctuations are strong, a repulsion at Q gets stronger than at small momentum, and sign-changing superconducting gap does develop. This scenario has been verified by e.g., observation of a spin resonance peak below Tc1,2,3,4,5,6. Spin fluctuations were also identified as the source for spontaneous breaking of lattice rotational symmetry (nematicity) in Fe-pnictides, as nematicity there develops in the immediate vicinity of the stripe magnetic order with momenta Q = (π, 0) or (0, π). It has been argued multiple times7,8,9,10,11that spin fluctuations create an intermediate phase with a composite spin order, which breaks symmetry between (π, 0) and (0, π), but reserves O(3) spin-rotational symmetry.

Situation is different, however, in bulk Fe-chalcogenide FeSe, which has been extensively studied in the last few years using various techniques. A pure FeSe develops a nematic order at Tp ~ 85K, and becomes superconducting at Tc ~ 9K. A nematic order decreases upon isovalent doping by either S or Te (FeSe1−xSx and FeSe1−xTex) and in both cases disappears at critical xc (0.17 for S doping and 0.53 for Te doping). There is no magnetic order below Tp for any x.

The absence of magnetism lead to two conjectures: (i) that nematicity in FeSe is a d − wave Pomeranchuk order, with order parameter bilinear in fermions, rather than a composite spin order, for which an order parameter is a 4-fermion operator12, and (ii) that the origin of superconductivity may be different from the one in Fe-pnictides. On (i), there is a consistency between the Pomeranchuk scenario for nematicity and the data already in pure FeSe: a Pomeranchuk order parameter necessary changes sign between hole and electron pockets13, consistent with the data12,14,15,16, and the temperature dependence of nematic susceptibility, measured by Raman, is in line with the Pomeranchuk scenario17,18,19. On (ii), superconductivity in pure FeSe is likely still mediated by spin fluctuations20,21,22,23,24,25,26, as evidenced by the correlation between NMR 1/T1 and superconducting Tc, the consistency between ARPES data on the gap anisotropy and calculations within spin fluctuation scenario, and the fact that a magnetic order does develop under pressure27. However, near and above critical xc, magnetic fluctuations are far weaker28,29, e.g., a magnetic order does not develop until high enough pressure. This strongly reduces the strength of spin-mediated pairing as the letter requires inter-pocket interaction (pair-hopping) to be enhanced by spin fluctuations to overcome intra-pocket repulsion30. It has been argued31,32,33,34,35 based on a variety of data (see below) that superconductivity for such x is qualitatively different from the one in pure FeSe. One argument here is that the gap anisotropy changes sign, another is that Tc in FeSe1−xTex shows a clear dome-like behavior around xc.

In this communication, we address the issue whether superconductivity in doped FeSe near xc can be mediated by nematic fluctuations. It seems natural at first glance to replace spin fluctuations by soft nematic fluctuations as a pairing glue. However, there are two obstacles, both related to the fact that soft nematic fluctuations are at small momentum transfer. First, they do not affect the pair hopping term between hole and electron pockets, which is the key element for spin-mediated superconductivity. Second, the bare pairing interaction at small momentum transfer is repulsive, and dressing it by nematic fluctuations only makes the repulsion stronger.

We show that the pairing interaction Veff(k, − k; p, − p), mediated by nematic fluctuations (first two momenta are incoming, last two are outgoing), does become attractive near xc, however for a rather special reason, related to the very origin of the Pomeranchuk order. Namely, the driving force for a d − wave Pomeranchuk order is density-density interaction between hole and electron pockets. It does have a d − wave component \({U}_{{{{\rm{he}}}}}^{{{{\rm{d}}}}}\) because low-energy excitations in the band basis are constructed of dxz and dyz orbitals. A sign-changing nematic order (a spontaneous splitting of densities of dxz and dyz orbitals) develops13 when \({U}_{{{{\rm{he}}}}}^{{{{\rm{d}}}}}\) exceeds d-wave intra-pocket repulsion, much like sign-changing s+− order develops when pair hopping exceeds intra-pocket repulsion in the particle-particle channel. By itself, \({U}_{{{{\rm{he}}}}}^{{{{\rm{d}}}}}\) does not contribute to pairing, however taken at the second order, it produces an effective attractive interaction between fermion on the same pocket. We go beyond second order and collect all ladder and bubble diagrams which contain d-wave polarization bubbles at a small momentum transfer. We show that this induced attraction is proportional to the susceptibility for a d − wave Pomeranchuk order. Because a nematic susceptibility diverges at xc, the induced attraction necessary exceeds the bare intra-pocket repulsion in some range around xc, i.e., the full intra-pocket pairing interaction becomes attractive.

This attractive interaction Veff(k, − k; , p, − p) ∝ Ak,pχnem(∣k − p∣) is rather peculiar because it inherits from \({U}_{{{{\rm{he}}}}}^{{{{\rm{d}}}}}\) the d − wave form-factor \({A}_{{{{\bf{k}}}},{{{\bf{p}}}}}={\cos }^{2}({\theta }_{{{{\bf{k}}}}}+{\theta }_{{{{\bf{p}}}}})\), where θk and θp specify the location of the fermions ((in our case, this holds on the hole pocket, which is made equally out of dxz and dyz orbitals). A similar pairing interaction has been earlier suggested for one-band models on phenomenological grounds36,37,38,39 assuming that d − wave nematic coupling is attractive. We show that such an interaction emerges in the model with purely repulsive interactions, once we add the pairing component, induced by inter-pocket density-density \({U}_{{{{\rm{he}}}}}^{{{{\rm{d}}}}}\).

Because χnem(k − p) diverges at k = p, the presence of the form-factor Ak,p in Veff(k, − k; , p, − p) implies that the strength of the attraction depends on the position of a fermion on a Fermi surface. As the consequence, the gap function on the hole pocket is the largest around hot points, specified by θh = nπ/2, n = 0 − 3, and rapidly decreases in cold regions centered at θc = nπ/2 + π/4, n = 0 − 3, This has been already emphasized in the phenomenological study36. This behavior shows up most spectacularly right at a nematic QCP, where the gap emerges at Tc only at hot spots and extends at smaller T into finite size arcs. The arcs length grows as T decreases, but as long as T is finite, there exist cold regions where the gap vanishes, i.e., the system preserves pieces of the original Fermi surface. At T = 0, the gap opens everywhere except at the cold spots θc, where nematic form factor \(\cos 2\theta\) vanishes, but is still exponentially small near them, \(\Delta (\theta )\propto \exp -1/{(\theta -{\theta }_{{{{\rm{c}}}}})}^{2}\).

This, we argue, leads to highly unconventional behavior of the specific heat coefficient Cv/T, which does not display a jump at Tc and instead increases as \({({T}_{{{{\rm{c}}}}}-T)}^{1/2}\), passes through a maximum at T ~ 0.8Tc, and behaves at smaller T like there is a non-zero residual Cv/T at T → 0 (see Fig. 2). In reality, Cv/T vanishes at T = 0, but nearly discontinuously, as \(1/{(\log ({T}_{{{{\rm{c}}}}}/T))}^{1/2}\). Also, because the regions, where the gap is non-zero, are disconnected, the gap phases are uncorrelated, and s − wave, d − wave and two-component p − wave (kx + eiαky) states are degenerate.

At a finite distance from a QCP and/or in the presence of non-singular pair-hopping between hole and electron pockets, the gap function becomes continuous, but maxima at θ = nπ/2 remain. The specific heat coefficient C(T)/T acquires a finite jump at Tc, but holds the same behavior at intermediate T as in Fig. 4, within some distance to a QCP. The condensation energies for s − wave, d − wave and p − wave states split. Which order develops depends on the interplay between the attractive pairing interaction, mediated by nematic fluctuations, and non-singular repulsion. The letter is far stronger in s − wave and d − wave channels, which favors p − wave symmetry. In this case, the most likely outcome is kx ± iky state, which breaks time-reversal symmetry.

Results

Model

The electronic structure of pure/doped FeSe in the tetragonal phase consists of two non-equal hole pockets, centered at Γ, and two electron pockets centered at X = (π, 0) and Y = (0, π) in the 1FeBZ. The hole pockets are composed of dxz and dyz fermions, the X pocket is composed of dyz and dxy fermions, and the Y pocket is composed of dxz and dxy fermions. The inner hole pocket is quite small and likely does not play much role for nematic order and superconductivity. We assume that heavy dxy fermions also do not play much role and consider an effective two-orbital model with a single dxz/dyz circular hole pocket, and mono-orbital electron pockets (dyz X-pocket and dxz Y-pocket). We define fermionic operators for mono-orbital Y and X pockets as f1 and f2, respectively (f1,k,σ = dxz,k+Y,σ, f2,k,σ = dyz,k+X,σ). The band operator for the hole pocket is \({h}_{{{{\bf{k}}}},\sigma }=\cos {\theta }_{{{{\bf{k}}}}}{d}_{yz,{{{\bf{k}}}},\sigma }+\sin {\theta }_{{{{\bf{k}}}}}{d}_{xz,{{{\bf{k}}}},\sigma }\). The kinetic energy is quadratic in fermionic densities and there are 14 distinct C4-symmetric interactions40 involving low-energy fermions near the hole and the two electron pockets (see Supplementary Discussions I41 for details). We take the absence of strong magnetic fluctuations in doped FeSe as an evidence that interactions at momentum transfer between Γ and X (Y) are far smaller than the interactions at small momentum transfer and neglect them. This leaves 6 interactions with small momentum transfer: 3 within hole or electron pockets and 3 between densities of fermions near different pockets. The single interaction between hole fermions contains an angle-independent term and terms proportional to \(\cos 2{\theta }_{{{{\bf{k}}}}}\cos 2{\theta }_{{{{\bf{p}}}}}\) and \(\sin 2{\theta }_{{{{\bf{k}}}}}\sin 2{\theta }_{{{{\bf{p}}}}}\), the two interactions between hole and electron pockets contain an angle-independent and a \(\cos 2{\theta }_{{{{\bf{k}}}}}\) term, where k belongs to the hole pocket and the three interactions between fermions on electron pockets contain only angle-independent terms.

Nematic susceptibility

Like we said, we associate the nematic order with a d-wave Pomeranchuk order. In the orbital basis, this order is an orbital polarization (densities of dxz and dyz fermions split). In the band basis, we introduce two d − wave order parameters on hole and electron pockets: \({\phi }_{{{{\rm{h}}}}}={\sum }_{{{{\bf{k}}}},\sigma }\langle {h}_{{{{\bf{k}}}},\sigma }^{{\dagger} }{h}_{{{{\bf{k}}}},\sigma }\rangle \cos 2{\theta }_{{{{\bf{k}}}}}\) and \({\phi }_{{{{\rm{e}}}}}={\sum }_{{{{\bf{k}}}}}\langle {f}_{2,{{{\bf{k}}}},\sigma }^{{\dagger} }{f}_{2,{{{\bf{k}}}},\sigma }\rangle -\langle {f}_{1,{{{\bf{k}}}}\sigma }^{{\dagger} }{f}_{1,{{{\bf{k}}}},\sigma }\rangle\). The set of two coupled self-consistent equations for ϕh and ϕe is obtained by summing ladder and bubble diagrams (see Supplementary Discussions II41) and is

$$\begin{array}{rcl}{\phi }_{{{{\rm{h}}}}}&=&-{\phi }_{{{{\rm{h}}}}}\,{U}_{{{{\rm{h}}}}}^{{{{\rm{d}}}}}\,{\Pi }_{h}^{d}-{\phi }_{{{{\rm{e}}}}}\,{U}_{{{{\rm{he}}}}}^{{{{\rm{d}}}}}\,{\Pi }_{e},\\ {\phi }_{{{{\rm{e}}}}}&=&-{\phi }_{{{{\rm{e}}}}}\,{U}_{{{{\rm{e}}}}}^{{{{\rm{d}}}}}{\Pi }_{e}-2{\phi }_{{{{\rm{h}}}}}\,{U}_{{{{\rm{he}}}}}^{{{{\rm{d}}}}}{\Pi }_{h}^{d}.\end{array}$$
(1)

Here, \({\Pi }_{h}^{d}=-\int\,{G}_{p}^{h}\,{G}_{p}^{h}\cos 2{\theta }_{{{{\bf{p}}}}},\) and \({\Pi }_{e}=-(1/2){\int}_{p}\left({G}_{p}^{X}\,{G}_{p}^{X}+{G}_{p}^{Y}\,{G}_{p}^{Y}\right)\) are the polarization bubbles for the hole and the electron pockets, (\({G}_{p}^{i}={G}^{i}(p,{\omega }_{m})\) are the corresponding Green’s functions and ∫ p stands for \(T{\sum }_{{\omega }_{n}}\int\,\frac{{d}^{2}{{{\bf{p}}}}}{{(2\pi )}^{2}}\)). As defined, \({\Pi }_{h}^{d}\) and Πe are positive. The couplings \({U}_{{{{\rm{h}}}}}^{{{{\rm{d}}}}},{U}_{{{{\rm{e}}}}}^{{{{\rm{d}}}}}\) and \({U}_{{{{\rm{he}}}}}^{{{{\rm{d}}}}}\) are d − wave components of intra-pocket and inter-pocket density-density interactions. All interactions are positive (repulsive). The analysis of (1) shows that the nematic order with different signs of ϕh and ϕe develops when \({U}_{{{{\rm{he}}}}}^{{{{\rm{d}}}}}\) is strong enough with the condition \(2\,{({U}_{{{{\rm{he}}}}}^{{{{\rm{d}}}}})}^{2}\ge {U}_{{{{\rm{h}}}}}^{{{{\rm{d}}}}}\,{U}_{{{{\rm{e}}}}}^{{{{\rm{d}}}}}\).

The nematic susceptibility is inversely proportional to the determinant of (1). Evaluating it at a small but finite momentum q, we obtain χnem(q) ∝ 1/Z, where

$$Z=\left(1+{U}_{{{{\rm{h}}}}}^{{{{\rm{d}}}}}\,{\Pi }_{h}^{d}(q)\right)\left(1+{U}_{{{{\rm{e}}}}}^{{{{\rm{d}}}}}\,{\Pi }_{e}(q)\right)-2{({U}_{{{{\rm{he}}}}}^{{{{\rm{d}}}}})}^{2}{\Pi }_{h}^{d}(q){\Pi }_{e}(q).$$
(2)

Pairing interaction

Our goal is to verify whether the pairing interaction near the onset of a nematic order is (i) attractive, (ii) scales with the nematic susceptibility, and (iii) contains the d-wave form-factor \({\cos }^{2}(2{\theta }_{{{{\bf{k}}}}})\). To do this, we use the fact that χnem(q) contains polarization bubbles \({\Pi }_{h}^{d}({{{\bf{q}}}})\) and Πe(q), and obtain the fully dressed pairing interaction by collecting infinite series of renormalizations that contain \({\Pi }_{h}^{d}({{{\bf{q}}}})\) and Πe(q) with small momentum q. This can be done analytically (see refs. 41,42 for detail). Because q is small, the dressed pairing interactions are between fermions on only hole pocket or only electron pockets: \({V}_{{{{\rm{eff}}}}}^{{{{\rm{h}}}}}({{{\bf{k}}}},{{{\bf{q}}}})={V}_{{{{\rm{eff}}}}}^{{{{\rm{h}}}}}\left({{{\bf{k}}}}+{{{\bf{q}}}}/2\right.,-({{{\bf{k}}}}+{{{\bf{q}}}}/2);{{{\bf{k}}}}-{{{\bf{q}}}}/2,-({{{\bf{k}}}}-{{{\bf{q}}}}/2)\), \({V}_{{{{\rm{eff}}}}}^{{{{\rm{e}}}}}({{{\bf{k}}}},{{{\bf{q}}}})={V}_{{{{\rm{eff}}}}}^{{{{\rm{e}}}}}\left({{{\bf{k}}}}+{{{\bf{q}}}}/2\right.,-({{{\bf{k}}}}+{{{\bf{q}}}}/2);{{{\bf{k}}}}-{{{\bf{q}}}}/2,-({{{\bf{k}}}}-{{{\bf{q}}}}/2)\). We find

$${V}_{{{{\rm{eff}}}}}^{{{{\rm{h}}}}}({{{\bf{k}}}},{{{\bf{q}}}})=\frac{{U}_{h}^{d}}{1+{U}_{h}^{d}{\Pi }_{h}^{d}({{{\bf{q}}}})}-{A}_{h}{({U}_{{{{\rm{he}}}}}^{{{{\rm{d}}}}})}^{2}\,{\cos }^{2}2{\theta }_{{{{\bf{k}}}}}\,{\chi }_{{{{\rm{nem}}}}}({{{\bf{q}}}})+\ldots$$
(3)
$${V}_{{{{\rm{eff}}}}}^{e}({{{\bf{k}}}},{{{\bf{q}}}})=\frac{{U}_{e}^{d}}{1+{U}_{e}^{d}{\Pi }_{e}({{{\bf{q}}}})}-{A}_{e}{({U}_{{{{\rm{he}}}}}^{{{{\rm{d}}}}})}^{2}\,{\chi }_{{{{\rm{nem}}}}}({{{\bf{q}}}})+\ldots ,$$
(4)

where \({A}_{h}=\frac{{\Pi }_{e}}{1+{U}_{{{{\rm{h}}}}}^{{{{\rm{d}}}}}\,{\Pi }_{h}^{d}({{{\bf{q}}}})}\) and \({A}_{e}=\frac{1}{2}\frac{{\Pi }_{h}^{d}}{1+{U}_{{{{\rm{e}}}}}^{{{{\rm{d}}}}}\,{\Pi }_{e}({{{\bf{q}}}})}\). The dots stand for other terms which do not contain \({\Pi }_{h}^{d}({{{\bf{q}}}})\) and \({\Pi }_{e}^{d}({{{\bf{q}}}})\) and are therefore not sensitive to the nematic instability.

We see that each interaction contains two terms. The first is the dressed intra-pocket pairing interaction. It does get renormalized, but remains repulsive and non-singular at the nematic instability. The second term is the distinct interaction, induced by \({U}_{{{{\rm{he}}}}}^{{{{\rm{d}}}}}\). It is (i) attractive, (ii) scales with the nematic susceptibility, and (iii) contains the d − wave nematic form-factor \({\cos }^{2}2{\theta }_{{{{\bf{k}}}}}\). We emphasize that the attraction is induced by inter-pocket density-density interaction, despite that relevant nematic fluctuations are with small momenta and the pairing interactions involve fermions from the same pocket.

Gap equation

Near a nematic QCP, χnem(q) is enhanced and the interaction, induced by \({U}_{{{{\rm{he}}}}}^{{{{\rm{d}}}}}\), is the dominant one. In the absence of pair-hopping, the gap equation decouples between hole and electron pockets. The most interesting case is when the gap develops first on the hole pocket (the case Ah > Ae). We use Ornstein-Zernike form χnem(q) = χ0/(δ2 + q2), where δ is the distance to a nematic QCP in units of momentum. At small δ, relevant q are of order δ. To first approximation, the non-linear equation for Δh(k) then becomes local, with angle-dependent coupling:

$$1=g\,{\cos }^{2}2{\theta }_{{{{\bf{k}}}}}\mathop{\int}\nolimits_{0}^{\Lambda }dx\frac{\tanh \left(\frac{\sqrt{{x}^{2}+| {\Delta }_{{{{\rm{h}}}}}({{{\bf{k}}}}){| }^{2}}}{2T}\right)}{\sqrt{{x}^{2}+{\Delta }_{{{{\rm{h}}}}}{({{{\bf{k}}}})}^{2}}}.$$
(5)

where \(g=m{({U}_{{{{\rm{he}}}}}^{{{{\rm{d}}}}})}^{2}{\chi }_{0}/(4\pi {k}_{{{{\rm{F}}}}}\delta )\). Because the coupling is larger at θk = nπ/2, n = 0 − 3, the gap appears at \({T}_{{{{\rm{c}}}}}=1.13\Lambda \exp (-1/g)\) only at these points. As T decreases, the range, where the gap is non-zero, extends to four finite arcs with the width \({\theta }_{0}(T)=0.5\arctan \sqrt{g\log {T}_{{{{\rm{c}}}}}/T}\) (see Fig. 1c). In the areas between the arcs, the original Fermi surface survives. We emphasize that this is the original Fermi surface, not the Bogoliubov one, which could potentially develop inside the superconducting state43,44. We plot ∣Δh(k)∣ along the Fermi surface at T = 0 and a finite T in Fig. 1a, b, and plot θ0(T) as a function of T/Tc in Fig. 1c. The phases of the gap function in the four arcs are not correlated, hence s − wave, d − wave (\({d}_{{x}^{2}-{y}^{2}}\)) and two-component p − wave (kx + eiαky with arbitrary α) are all degenerate. At T = 0, the arcs ends merge at θk = nπ/2 + π/4, n = 0 − 3 and the gap becomes non-zero everywhere except these cold spots(red dots in Fig. 1a). In explicit form, \(| {\Delta }_{{{{\rm{h}}}}}(k)| =1.76{T}_{{{{\rm{c}}}}}\exp \{-\,{\tan }^{2}2{\theta }_{{{{\bf{k}}}}}/g\}\). The gap near cold spots becomes a bit smoother if we keep the Landau damping in χnem and solve the dynamical pairing problem, but still Δh(k) remains highly anisotropic.

Fig. 1: Superconductivity at the nematic QCP(δ = 0).
figure 1

a We plot the absolute value of the SC gap scaled by the transition temperature Tc, Δh/Tc(blue) on the hole pocket at zero temperature as a function of the angle on the Fermi surface(green disk) with radius 3 unit. The gap is exponentially small near the cold spots(red dots) while maximum along the kx and ky axis. In (b), we plot the angular variation of the gap, ∣Δh∣/Tc on the hole pocket for three different reduced temperatures, t = 0.9, 0.7 and 0 below the transition point Tc. At finite temperature, the gap vanishes on the four patches of the Fermi surface arcs(flat segments of the gap). In (c), we plot the width of the gap, θ0(T) as a function of the reduced temperature t = T/Tc. We keep the coupling strength g = 1 here.

Full size image

Specific heat

We split the specific heat coefficient γc = Cv(T)/T into contributions from the gapped and ungapped regions of the Fermi surface: \({\gamma }_{{{{\rm{c}}}}}(T)={\gamma }_{{{{\rm{c}}}}}^{n}(T)+{\gamma }_{{{{\rm{c}}}}}^{s}(T)\). The first term, \({\gamma }_{{{{\rm{c}}}}}^{n}(T)=\frac{8{N}_{0}\pi }{3}\left[\frac{\pi }{4}-{\theta }_{0}(T)\right]\) which at small T becomes: \({\gamma }_{c}^{n}(T)\approx \frac{4\,{N}_{0}\,\pi }{3\,\sqrt{g\,\log {T}_{{{{\rm{c}}}}}/T}}\). It evolves almost discontinuously: vanishes at T = 0, but reaches 1/3 of the normal state value already at T = 0.01Tc. We obtained \({\gamma }_{c}^{s}(T)\) numerically and show the result for the full γc(T) in Fig. 2. We see that γc(T) does not jump at Tc. Instead, it increases from its normal state value as \(\sqrt{{T}_{{{{\rm{c}}}}}-T}\), passes through maximum at T ≈ 0.8Tc and nearly linearly decreases at smaller T, apparently with a finite offset at T = 0. It eventually drops to zero at T = 0, but only at extremely small T, as \(1/{(\log {T}_{{{{\rm{c}}}}}/T)}^{1/2}\). We emphasize that γc(T) is a function of a single parameter T/Tc, i.e., the smallness of the range, where γc(T) drops, is purely numerical.

Fig. 2: Specific heat at the nematic QCP(δ = 0).
figure 2

We plot the specific heat coefficient scaled by the density of state on Fermi surface, N0 as a function of the reduced temperature, T/Tc for the coupling strength g = 1.0 in the main frame and for g = 0.6, 0.4 and 0.2 in the inset. The black dashed line represents the normal state contribution and is equal to 2π2/3.

Full size image

Away from a nematic QCP

At a finite δ, s − wave, d − wave, and p − wave solutions for the gap function are no longer degenerate. If we keep only the interaction induced by \({U}_{{{{\rm{he}}}}}^{{{{\rm{d}}}}}\) (the second term in (4)), we find that s-wave solution has the lowest condensation energy. We show the eigenvalues λs,p,d and the gap functions in Fig. 3a, b. The gap function is smooth and finite for all angles, but remains strongly anisotropic up to sizable δ. We define the gap anisotropy α as the ratio of the gap function on the hole fermi surface at θ = π/4(kx − ky axis) to θ = 0(kx axis): α = Δh(π/4)/Δh(0) and show its variation with the nematic mass parameter δ in Fig. 3c. The specific heat coefficient γc(T) has a finite jump at Tc, whose magnitude increases with δ, yet the low temperature behavior remains nearly the same as at a QCP up to sizable δ (Fig. 4). If we consider the full pairing interaction in (3), situation may change as the first term in (3) has comparable repulsive s − wave and d − wave harmonics, but a much smaller p − wave harmonic. As the consequence, p − wave may become the leading instability. The condensation energy for a p-wave state is the lowest for kx + iky and kx − iky gap functions. A selection of one of these states breaks time-reversal symmetry.

Fig. 3: Superconductivity away from the nematic QCP(δ ≠ 0).
figure 3

In (a), we plot the largest eigenvalue λ of the linearized gap equation defined in Supplementary Discussions VIII41 in different angular momentum channels labeled as λs, λd and λp for the s, d and p − wave respectively as a function of the nematic mass parameter δ. We show how the ratio of these eigenvalues vary with δ in the inset and find λp/λs < λd/λs < 1 which indicates that s − wave is the leading instability. With δ going to zero, the ratios become closer to each other, indicating possible degeneracy among different pairing channels. In (b) we plot the angular variation(θ) of the gap function, Δh(θ)/Tc on the hole pocket for a set of reduced temperatures, t = T/Tc below the transition point Tc for δ = 0.01. In (c) we plot the gap anisotropy α = Δh(θ = π/4)/Δh(θ = 0) as a function of the nematic mass δ and fit (blue dashed line) our result upto second order in δ with the fitting parameters α(δ) = 2.12 δ2 + 0.44 δ.

Full size image
Fig. 4: Specific heat away from the nematic QCP(δ ≠ 0).
figure 4

We plot the variation of the specific heat coefficient scaled by the density of state on Fermi surface, N0 with the reduced temperature T/Tc for a set of values of the nematic mass δ = 0.01 and 0.1. Here, Tc is the superconducting transition temperature. The black dashed line represents the normal state contribution and is equal to 2π2/3. There is a finite specific heat jump at the transition point represented by the vertical line for both values of δ.

Full size image

Comparison with experiments

We argued in this work is that pairing in doped FeSe near a nematic QCP is mediated by nematic fluctuations rather than by spin fluctuations. This is generally consistent with the observations in refs. 31,32,34 of two distinct pairing states in pure FeSe and in doped FeSe1−xSx and FeSe1−xTex at x ≥ xc. More specifically, one can distinguish between magnetic and nematic pairing scenarios by measuring the angular dependence of the gap along the hole dxz/dyz pocket. We argued that a nematic-mediated pairing gives rise to an anisotropic gap, with maxima along kx and ky directions. Within spin-fluctuation scenario, the gap \({\Delta }_{{{{\rm{h}}}}}(k)=a+b\cos 4\theta\) is the largest along the diagonal directions kx ± ky (b < 0, see e.g., ref. 45). The angular dependence of the gap in pure and doped FeSe has been extracted from ARPES and STM data in ref. 23,46,47,48,49,50. For pure and weakly doped FeSe, an extraction of \(\cos 4\theta\) dependence is complicated because superconductivity co-exists with long-range nematic order, in which case the gap additionally has \(\cos 2\theta\) term due to nematicity-induced mixing of s − wave and d − wave components24,51. Still, the fits of the ARPES data in refs. 23,46 yielded a negative b, consistent with spin-fluctuation scenario. A negative b is also consistent with the flattening of the gap on the hole pocket near θ = π, observed in the STM study47. A negative prefactor for \(\cos 4\theta\) term was also reported for Fe-pnictides, e.g., Ba0. 24K0.76Fe2As2, ref. 52. In contrast, gap maximum along ky has been reported in a recent laser ARPES study of FeSe0.78S0.22 (ref. 48). Further, recent STM data for FeSe0.81S0.19 (ref. 49,50) detected a clear gap maxima along kx and ky. For Te-doped case, STM data for tetragonal FeSe0.45Te0.55 (ref. 53) also found the maximal gap along kx and ky directions, consistent with the pairing by nematic fluctuations. This STM data is in contradiction to earlier ARPES data for the same material, which reported a near-isotropic gap on the hole pocket54. However, the gap magnitude is only 2 meV, which calls for laser ARPES analysis. The gap anisotropy has also been analyzed in ref. 55, using angle-resolved specific heat data, but the authors of that work focused on the gap on the electron pockets. Taken together, these data strongly support the idea about different pairing mechanisms in pure FeSe and in doped ones at x≥xc, and are consistent with the change of the pairing glue from spin fluctuations at x < xc to nematic fluctuations at x ≥ xc.

Next, we argued that right at a nematic QCP, the gap vanishes in the cold regions on the Fermi surface, and this leads to highly unconventional behavior of the specific heat coefficient γc(T). This holds when we neglect pair hopping between hole and electron pockets. In the presence of pair hopping, the gap becomes non-zero everywhere except, possibly, special symmetry-related points. Still, in the absence of magnetism nearby, pair-hopping is a weak perturbation, and the gap in cold regions is small. The specific heat of FeSe1−xSx has been measured in refs. 33,56. The data clearly indicate that the jump of γc(T) at Tc decreases with increasing x and vanishes at around xc. At smaller T, γc(T) passes through a maximum at around 0.8Tc and then decreases nearly linearly towards apparently a finite value at T = 0. The authors of ref. 31 argued that this behavior is not caused by fluctuations because residual resistivity does not exhibit a noticeable increase around xc (ref. 57). Other experiments58 also indicated that fluctuation effects get weaker with increasing x.

The behavior of γc(T) around xc was first interpreted as potential BCS-BEC crossover32 and later as a potential evidence of an exotic pairing that creates a Bogolubov Fermi surface in the superconducting state43,48,59. We argue that the specific heat data are consistent with the nematic-mediated pairing, in which near xc the gap develops in the arcs near kx and ky and nearly vanishes in between the arcs. This explanation is also consistent with recent observation60 that superfluid density in FeSe1−xSx drops at x≥xc, indicating that some fermions remain unpaired.

Finally, recent μSR experiments60,61 presented evidence for time-reversal symmetry breaking in FeSe. The μSR signal is present below Tc for all x, however in FeSe1−xTex it clearly increases above xc. This raises a possibility that the superconducting state at x > xc breaks time-reversal symmetry, at least in FeSe1−xTex. Within our nematic scenario, this would indicate a p − wave pairing with kx ± iky gap structure. We argued that p − wave pairing, mediated by nematic fluctuations, is a strong competitor to s+− pairing.

There is one recent data set, which we cannot explain at the moment. Laser ARPES study of FeSe0.78S0.22 (ref. 48) detected superconducting gap in the polarizarion of light, which covers momenta near the X direction, but no gap in polarization selecting momenta near Y. Taken at a face value, this data implies that superconducting order strongly breaks C4 symmetry. In our nematic scenario, pure kx (or ky) order is possible, but has smaller condensation energy than kx ± iky. More analysis is needed to resolve this issue.

Discussion

In this paper, we derived an effective pairing interaction near the onset of a nematic order in a 2D two-orbital/three band system of fermions and applied the results to doped FeSe. The model consists of a hole band, centered at Γ and made equally of dxz and dyz fermions, and two electron bands, centered at X and Y and made out of dyz and dxz fermions, respectively. The nematic order is a spontaneous polarization between dxz and dyz orbitals, which changes sign between hole and electron pockets. We found the pairing interaction as the sum of two terms: a dressed bare interaction, which remains non-singular and repulsive, and the term, induced by inter-pocket density-density interaction \({U}_{{{{\rm{he}}}}}^{{{{\rm{d}}}}}\). This last term contains the square of the nematic form-factor and scales with the nematic susceptibility, and is the dominant pairing interaction near the onset of a nematic order. We obtained the gap function and found that it is highly anisotropic with gap maxima along kx and ky directions. This is in variance with pairing by spin fluctuations, for which the gap has maxima along diagonal directions kx ± ky. Right at the nematic QCP, the gap develops in four finite arcs around kx and ky, while in between the arcs the original Fermi surface survives. Such a gap function, degenerate between s-wave, d − wave, and p − wave, gives rise to highly unconventional behavior of the specific heat coefficient with no jump at Tc and seemingly finite value at T = 0 (the actual Cv(T)/T vanishes at T = 0, but drops only at extremely low T ~ 10−2Tc). In the tetragonal phase away from a QCP, the degeneracy is lifted, and there is a competition between s − wave and kx ± iky, the latter breaks time-reversal symmetry. In both cases, the gap remains strongly anisotropic, with maxima along X and Y directions. We compared our theory with existing experiments in some details.

Methods

Analytical calculations

Analytic calculations have been performed using diagrammatic theory. We obtained the susceptibility of the sign-changing d-wave nematic order by summing bubble and ladder series of diagrams for the renormalization of the nematic vertex, and obtained the pairing vertex by summing up ladder, bubble, and maximally crossed diagrams for the bare intra-pocket pairing interaction and for the pairing interaction induced by inter-pocket density-density interaction. We present the details of the calculations in Supplementary Discussion I–VIII.

Numerical calculations

We numerically solved the non-linear integral equation for the superconducting gap and used the results to compute the specific heat. The details are presented in the Supplementary Discussion VIII.