Pairing mechanism of heavily electron doped FeSe systems: dynamical tuning of the pairing cutoff energy

In order to illustrate the essence of the dynamical tuning mechanism and RG scaling of the cutoff energy, we first study a minimal two band model with one incipient hole band and one electron band as depicted in figure 3(a). Here we ignore the interaction between electron bands e1 and e2, therefore two electron bands, e1 and e2, can be considered as identical one e-band as far as the SC pairing mechanism is concerned. As a result, among the three possible pairing states depicted in figure 2, only ${s}_{{he}}^{-+}$-gap (figure 2(a)) and ${s}_{{ee}}^{++}$-gap (figure 2(b)) are possible with the two band model. This simplification is only for the proof of concept and we will consider the full three band model later.

Zoom In Zoom Out Reset image size

In this paper, we assume the incipient energy ε_b to be smaller than the spin-fluctuation pairing interaction cutoff Λ_sf as illustrated in figure 3. This incipient band model [16] and its extended models [17] were recently studied as a possible candidate model for the HEDIS SC systems but only the incipient ${s}_{{he}}^{-+}$-type solution (figure 2(a)) was investigated. Here we will consider the incipient ${s}_{{he}}^{-+}$-gap solution and ${s}_{{ee}}^{++}$-gap solution (figure 2(b)) on an equal footing in the two band model. The ${s}_{e1e2}^{+-}$-gap solution (figure 2(c)) will be considered with three band model later.

The coupled gap equations are the same as the ordinary two band SC model as follows.

$$\begin{eqnarray}{{\rm{\Delta }}}_{h} & = & [{V}_{{hh}}{\chi }_{h}]{{\rm{\Delta }}}_{h}+[{V}_{{he}}{\chi }_{e}]{{\rm{\Delta }}}_{e},\\ {{\rm{\Delta }}}_{e} & = & [{V}_{{ee}}{\chi }_{e}]{{\rm{\Delta }}}_{e}+[{V}_{{eh}}{\chi }_{h}]{{\rm{\Delta }}}_{h}\end{eqnarray} \tag{ 5 }$$

where the pair susceptibilities are defined as

$$\begin{eqnarray}{\chi }_{h}(T) & = & -\displaystyle \,\frac{{N}_{h}}{2}{\displaystyle \int }_{-{{\rm{\Lambda }}}_{{sf}}}^{-{\varepsilon }_{b}}\displaystyle \frac{{\rm{d}}\xi }{\xi }\tanh \left(\displaystyle \frac{\xi }{2T}\right)\approx -\displaystyle \frac{{N}_{h}}{2}\mathrm{ln}\left[\displaystyle \frac{1.14{{\rm{\Lambda }}}_{{sf}}}{{\varepsilon }_{b}}\right],\\ {\chi }_{e}(T) & = & -\,{N}_{e}{\displaystyle \int }_{-{{\rm{\Lambda }}}_{{sf}}}^{{{\rm{\Lambda }}}_{{sf}}}\displaystyle \frac{{\rm{d}}\xi }{\xi }\tanh \left(\displaystyle \frac{\xi }{2T}\right)\approx -{N}_{e}\mathrm{ln}\left[\displaystyle \frac{1.14{{\rm{\Lambda }}}_{{sf}}}{T}\right]\end{eqnarray} \tag{ 6 }$$

where N_h,e are the density of states (DOS) of the hole band and electron band, respectively. Assuming all repulsive pair potentials V_ab, but with ${V}_{{\rm{inter}}-{\rm{band}}}(={V}_{{he}},{V}_{{eh}})\gt {V}_{{\rm{intra}}-{\rm{band}}}(={V}_{{hh}},{V}_{{ee}})$, the coupled gap equations equation (5) with the susceptibilities equation (6) produce the incipient ${s}_{{he}}^{+-}$-gap solution with the OPs Δ_h and Δ_e with opposite signs [16, 17]. T_c decreases with increasing ε_b; however the relative size of $| {{\rm{\Delta }}}_{h}| /| {{\rm{\Delta }}}_{e}|$ is insensitive to the value of ε_b and the gap size of the incipient hole band $| {{\rm{\Delta }}}_{h}|$ is comparable to the size of $| {{\rm{\Delta }}}_{e}|$.

This model is already a low energy effective model in which all the high energy interaction processes originating from U (Hubbard on-site interaction), J (Hund coupling), etc are renormalized down to Λ_sf to produce the effective pair potentials ${V}_{{ab}}^{0}({{\rm{\Lambda }}}_{{sf}})$ with cutoff energy Λ_sf. In particular, the characteristic AFM spin fluctuation energy scale Λ_sf is an experimentally measurable—by neutron experiments—physical excitations just like a phonon energy ω_D in the BCS theory. In a standard theory of superconductivity, the renormalization stops here and remained is to solve the gap equation(s) (e.g. Equation (5)) by a mean field method (BCS theory) or by its dynamical extension (Eliashberg theory).

In this paper, we continue the RG scaling down to arbitrary low energy scale Λ < Λ_sf and define the renormalized model H_ren(Λ), as depicted in figure 3(b). The renormalized pair potential V_ab(Λ) is defined by a standard RG process as follows (see figure 4),

$$\begin{eqnarray}&&\hat{V}({\rm{\Lambda }})={\hat{V}}^{0}+{\hat{V}}^{0}\cdot \hat{\chi }({{\rm{\Lambda }}}_{{sf}};{\rm{\Lambda }})\cdot \hat{V}({\rm{\Lambda }})\end{eqnarray} \tag{ 7 }$$

and the formal solution is defined by

$$\begin{eqnarray}&&\hat{V}({\rm{\Lambda }})={[1-{\hat{V}}^{0}\cdot \hat{\chi }({{\rm{\Lambda }}}_{{sf}};{\rm{\Lambda }})]}^{-1}\cdot {\hat{V}}^{0},\end{eqnarray} \tag{ 8 }$$

where $\hat{V}({\rm{\Lambda }})$ is the renormalized 2 × 2 matrix pair potential with the new cutoff energy Λ, and ${\hat{V}}^{0}$ is the bare pair potential with the cutoff energy Λ_sf defined as

$$\begin{eqnarray}{\hat{V}}^{0}=\left(\begin{array}{cc}{V}_{{hh}}^{0} & {V}_{{he}}^{0}\\ {V}_{{eh}}^{0} & {V}_{{ee}}^{0}\end{array}\right).\end{eqnarray} \tag{ 9 }$$

Accordingly, the Cooper susceptibility $\hat{\chi }$ is also 2 × 2 matrix defined as

$$\begin{eqnarray}\hat{\chi }({{\rm{\Lambda }}}_{{sf}};{\rm{\Lambda }})=\left(\begin{array}{cc}{\chi }_{h} & 0\\ 0 & {\chi }_{e}\end{array}\right)\end{eqnarray} \tag{ 10 }$$

with ${\chi }_{e}({{\rm{\Lambda }}}_{{sf}};{\rm{\Lambda }})=-{T}_{c}{\sum }_{n}2{N}_{e}{\int }_{{\rm{\Lambda }}}^{{{\rm{\Lambda }}}_{{sf}}}{\rm{d}}\xi \tfrac{1}{{\omega }_{n}^{2}+{\xi }^{2}(k)}$ and ${\chi }_{h}({{\rm{\Lambda }}}_{{sf}};{\rm{\Lambda }})=-{T}_{c}{\sum }_{n}{N}_{h}{\int }_{-{{\rm{\Lambda }}}_{{sf}}}^{-{\rm{\Lambda }}}{\rm{d}}\xi \tfrac{1}{{\omega }_{n}^{2}+{\xi }^{2}(k)}$. Notice that the factor 2 is missing in χ_h since the hole band exists only below the Fermi level. Also χ_h(Λ_sf; Λ) is defined only up to ${\rm{\Lambda }}\to {\varepsilon }_{b}$, meaning that when the scaling cutoff Λ runs below ε_b as Λ < ε_b, only the electron band will contribute to the RG flow.

Zoom In Zoom Out Reset image size

In figure 5(a), we show the results of the renormalized ${\hat{V}}_{{ab}}({\rm{\Lambda }})$ for the whole range of Λ < Λ_sf. In this representative case, we chose ε_b = 0.5 Λ_sf and the bare potentials: ${N}_{h}{V}_{{hh}}^{0}={N}_{e}{V}_{{ee}}^{0}=0.5$ and $\sqrt{{N}_{h}{N}_{e}}{V}_{{he}}^{0}=\sqrt{{N}_{h}{N}_{e}}{V}_{{eh}}^{0}=2.0$. For ε_b < Λ < Λ_sf, all four components of ${\hat{V}}_{{ab}}({\rm{\Lambda }})$ get renormalized. The key features are: (1) the repulsive interband pair potentials V_{he, eh}(Λ) slightly decrease at the beginning but eventually become more repulsive; this is very different behavior compared to the standard single band repulsive potential under RG such as the Coulomb potential. (2) The weaker repulsive intraband pair potentials V_{hh, ee}(Λ) turn quickly into attractive ones; the different flow tracks of V_hh(Λ) and V_ee(Λ), despite the same starting bare potential ${N}_{h}{V}_{{hh}}^{0}={N}_{e}{V}_{{ee}}^{0}=0.5$, is due to the difference of χ_h and χ_e. For Λ < ε_b, only V_ee(Λ) continues to scale and other pair potentials stop scaling because χ_h stops scaling at Λ = ε_b.

Zoom In Zoom Out Reset image size

In figure 5(b), we showed the calculated T_c of the renormalized model H_ren(Λ) for all values of Λ < Λ_sf depicted in figure 3(b). Namely, we solved the gap equations equation (5) in the limit of ${{\rm{\Delta }}}_{h,e}\to 0$ with the renormalized pair potentials V_ab(Λ) and the reduced cutoff Λ. For ε_b < Λ < Λ_sf, the pairing gap solution is the incipient ${s}_{{he}}^{-+}$-solution, forming SC OPs both on the hole- (Δ⁻) and electron band (Δ⁺) (see figure 2(a)). When Λ < ε_b, the hole band is outside the pairing cutoff Λ, hence cannot participate forming Cooper pairs. Therefore the gap solution consists of the electron band only with attractive ${V}_{{ee}}({\rm{\Lambda }})$, i.e. ${s}_{{ee}}^{++}$-solution (see figure 2(b)).

Figure 5(b) shows that the T_c (Λ) remain the same all the time with the scaling for Λ < Λ_sf. This is a consistent result with the invariance principle of RG (see footnote 1) but still surprising to be confirmed with the multiband SC model with a complicated pair potential ${\hat{V}}_{{ab}}$, in particular, despite the change of the pairing solutions from ${s}_{{he}}^{-+}$ to ${s}_{{ee}}^{++}$ (see footnote 2). We would like to emphasize that, in order to achieve this T_c invariance with RG transformation, it is crucially important to calculate the T_c-equation equation (5) together with the renormalized pair potentials ${\hat{V}}_{{ab}}({\rm{\Lambda }};T={T}_{c})$ equation (8) at same temperature, which needs self-consistent iterative calculations of two equations, because the RG flow itself depends on temperature.

In fact, this electron band only pairing state, ${s}_{{ee}}^{++}$, has been suggested by previous works [18, 19] of FRG studies, applied to the tight binding model with local interactions U, U', and J_H, designed for the FeSe systems. The FRG technique traces the RG flow of the effective pairing channels ${V}_{\alpha }({\rm{\Lambda }}){\phi }_{\alpha }^{* }(k){\phi }_{\alpha }(k^{\prime} )$ and identifies the most diverging channel 'α' as the winning instability of the system as the RG scale Λ runs to arbitrarily low energy. However, because the ${s}_{{ee}}^{++}$-channel and the ${s}_{{he}}^{-+}$-channel are symmetry-wise identical and have the same form factor ${\phi }_{\alpha }(k)\sim (\cos {k}_{x}+\cos {k}_{y})$ (in two Fe/cell BZ), these two channels are running on the same track of the FRG flow. Therefore, identifying the most diverging form factor ϕ_α(k) itself does not distinguish which of these two pairing states is a true ground state. What distinguishes these two pairing states is neither 'channel' nor 'symmetry' (they are always the same), but the physical pairing cutoff Λ_phys, however the FRG scheme itself does not determine the physical pairing cutoff Λ_phys.

This is exactly the point we are addressing in this paper: how to determine the physical pairing cutoff Λ_phys. A common sense would be to identify as Λ_phys = Λ_sf (the spin-fluctuation energy scale). And when Λ_sf > ε_b (figure 3(a)), the pairing solution should be the 'incipient' ${s}_{{he}}^{-+}$-state forming gaps on the sunken hole band as well as on the electron bands [16, 17]. And only if Λ_sf < ε_b and ${V}_{{ee}}^{0}({{\rm{\Lambda }}}_{{sf}})\lt 0$, the pairing solution can be the ${s}_{{ee}}^{++}$-state, forming gaps only on the electron bands. However, we have shown in figure 5(b) that even when Λ_sf > ε_b, there is no reason to fix the physical pairing cutoff as Λ_phys = Λ_sf because the whole range of Λ < Λ_sf produces the pairing states with the same T_c and the same symmetry but with different pairing cutoffs.

In the previous section, we have argued that among the degenerate pairing solutions with different pairing cutoffs Λ, the physical ground state should be chosen as the one with the largest Λ because the condensation energy (CE) gain is larger with a larger cutoff energy Λ when considering the higher order corrections ∼O(T_c/Λ, Δ_h,e/Λ) to the CE. Therefore, among the continuously degenerate solutions of the ${s}_{{he}}^{\pm }({\rm{\Lambda }})$-state for ε_b < Λ < Λ_sf, the system should choose the incipient ${s}_{{he}}^{\pm }({{\rm{\Lambda }}}_{{sf}})$ state with cutoff energy, Λ_phys = Λ_sf, as a physical ground state. By the same reasoning, among the degenerate solutions of the ${s}_{{ee}}^{++}({\rm{\Lambda }})$-state for Λ < ε_b, the system will choose the ${s}_{{ee}}^{++}$ state with Λ_phys = ε_b as a physical ground state. These two physical ground states are marked as black star symbols in figure 5(b) and specifically illustrated in figure 6. Therefore, our work provides the rationale for justifying the electron band only pairing state, ${s}_{{ee}}^{++}$, with the physical cutoff Λ_phys = ε_b even when the original pairing cutoff is Λ_sf that is larger than ε_b.

Zoom In Zoom Out Reset image size

Now we compare two physical ground state solutions of a given incipient two band model with ε_b < Λ_sf: the incipient ${s}_{{he}}^{\pm }$-state solution with Λ_phys = Λ_sf and the ${s}_{{ee}}^{++}$-state solution with Λ_phys = ε_b (see figure 6). We numerically calculate the T_c of these two gap solutions with varying ε_b for 0 < ε_b < Λ_sf. For all these calculations, it is important to calculate the Cooper susceptibilities without approximation because the cutoff energy Λ_phys = ε_b of the ${s}_{{ee}}^{++}$-state can be very low to violate the BCS limit $\left(\tfrac{{T}_{c}}{{{\rm{\Lambda }}}_{{\rm{phys}}}}\ll 1\right)$. Therefore, we numerically calculate ${\chi }_{e}({\rm{\Lambda }};0)=-2T{\sum }_{n}{N}_{e}{\int }_{0}^{{\rm{\Lambda }}}{\rm{d}}\xi \tfrac{1}{{\omega }_{n}^{2}+{\xi }^{2}(k)}$, which can be very different from the BCS approximation ${\chi }_{e}({\rm{\Lambda }};0)\sim -{N}_{e}\mathrm{ln}\tfrac{1.14{\rm{\Lambda }}}{T}$ when Λ < T. We chose the same interaction parameters, ${N}_{h}{V}_{{hh}}^{0}={N}_{e}{V}_{{ee}}^{0}=0.5$ and $\sqrt{{N}_{h}{N}_{e}}{V}_{{he}}^{0}=\sqrt{{N}_{h}{N}_{e}}{V}_{{eh}}^{0}=2.0$ as in the case of figure 5.

Figure 7 shows that the T_cs of the incipient ${s}_{{he}}^{\pm }$-state with Λ_phys = Λ_sf (black plus symbols) and the T_cs of the ${s}_{{ee}}^{++}$-state with Λ_phys = ε_b (red solid circles) are exactly overlayed on top of each other, in clean limit (Γ = 0), as explained before. T_c decreases with increasing ε_b as the hole band sinks deeper below Fermi level. In the simple case when ${V}_{{hh}}^{0}={V}_{{ee}}^{0}=0.0$, we can derive an analytic formula for T_c without impurities as [16]

$$\begin{eqnarray}&&{T}_{c}({\varepsilon }_{b})\approx 1.14{{\rm{\Lambda }}}_{{sf}}\exp [-1/{\lambda }_{{\rm{eff}}}({\varepsilon }_{b})]\end{eqnarray} \tag{ 11 }$$

with ${\lambda }_{{\rm{eff}}}({\varepsilon }_{b})=[{V}_{{he}}^{0}{N}_{e}{V}_{{eh}}^{0}{N}_{h}]\cdot \mathrm{ln}\left[\tfrac{1.14{{\rm{\Lambda }}}_{{sf}}}{{\varepsilon }_{b}}\right]$. This formula shows that T_c (ε_b)/Λ_sf is a function of ε_b/Λ_sf.

Zoom In Zoom Out Reset image size

One important new feature is that while the T_c of the incipient ${s}_{{he}}^{\pm }$-state continuously increases as ${\varepsilon }_{b}\to 0$, the T_c of the ${s}_{{ee}}^{++}$-state becomes ill-defined when ${\varepsilon }_{b}\to 0$. This is because the Cooper susceptibility ${\chi }_{e}({\rm{\Lambda }}={\varepsilon }_{b};T)$ gets saturated and loses its logarithmic divergence as ε_b ∼ O(T_c); in fact χ_e(Λ = ε_b; T) decreases, instead of increasing, with lowering temperature when ε_b ∼ O(T_c) so that the T_c-equation fails to find a solution. This behavior provides an important clue why there exists a minimum incipient distance energy ε_b (=Λ_mim) below which the ${s}_{{ee}}^{++}$-pairing state is not stabilized regardless of the strength of the pair potentials ${\hat{V}}_{{ab}}({\rm{\Lambda }})$. In clean limit, this minimum cutoff energy scale is Λ_mim ∼ T_c, but it will increase with additional relaxation processes, thermal or dynamical origins, such as ${{\rm{\Lambda }}}_{{\rm{mim}}}\sim (\pi {T}_{c}+{{\rm{\Gamma }}}_{{\rm{imp}}}+{{\rm{\Gamma }}}_{{\rm{inela}}})$, where Γ_imp is static impurity scattering rate, and Γ_inelas is inelastic scattering rate that can be provided by spin fluctuations as Γ_inelas = Im Σ_sf(T).

We now consider the impurity pair-breaking effect on T_c of both pairing states and in this paper we considered non-magnetic impurities only. First, we have to investigate the impurity effect on the RG scaling itself. The answer is that the non-magnetic impurities has no effect on the RG scaling because the Cooperon propagators (pair susceptibilities) ${\chi }_{h,e}({{\rm{\Lambda }}}_{{sf}};{\rm{\Lambda }})$ defined in equation (10) are invariant with the non-magnetic impurity scattering. The fundamental reason for this invariance has the same origin as the Anderson's theorem [36]: the s-wave pairing is not affected by the non-magnetic impurity scattering. This is easy to understand by noting that the Cooperon propagators χ_h,e(Λ_sf; Λ) entering the RG equation are nothing but the s-wave pair susceptibility. Diagrammatic derivation of this proof following the formalism of Abrikosov and Gor'kov [37] is given in appendix A. Therefore, the renormalized pair potentials $\hat{V}({\rm{\Lambda }})$ are the same with and without the non-magnetic impurities.

Then given the same pairing potentials $\hat{V}({\rm{\Lambda }})$, the non-magnetic impurity scattering would not change T_c of the s-wave pairing according to Anderson's theorem [36], so that the T_c of the ${s}_{{ee}}^{++}$-state in figure 7 will not be affected with non-magnetic impurities. On the other hand, it is well known that non-magnetic impurities will suppress the T_c of the 'standard' ${s}_{{he}}^{\pm }$-state almost as strong as in the d-wave [38]. We might expect that the T_c suppression with non-magnetic impurities on the 'incipient' ${s}_{{he}}^{\pm }$-state will be weakened compared to the case of the standard ${s}_{{he}}^{\pm }$-state because the incipient band, being sunken below Fermi level, should be less effective for any scattering process. However, we have found that this weakening effect due to the incipiency is only marginal and the non-magnetic impurity T_c-suppression rate of the 'incipient' ${s}_{{he}}^{\pm }$-state is almost as strong as the case of the 'standard' ${s}_{{he}}^{\pm }$-state. The physical reason is because the relative size of OPs, Δ_h and Δ_e are not affected by the incipient distance energy ε_b as far as the pairing interactions is dominated by the interband potentials as V_he,eh > V_hh,ee [16, 17]. The detailed formalism of the impurity effect on the incipient ${s}_{{he}}^{\pm }$-state is described in appendix B.

In figure 7, we plotted the results of T_c suppression of the incipient ${s}_{{he}}^{\pm }$-state with the different impurity pair-breaking rate Γ/Λ_sf = 0.1, 0.2, 0.3, 0.4, and 0.5, respectively. In this paper, we used the equal strength unitary scatterers (c = 0) both for inter- and intra-band impurity scatterings as ${V}_{{he},{eh}}^{{\rm{imp}}}={V}_{{hh},{ee}}^{{\rm{imp}}}$ (see appendix A). As expected, the T_c suppression rate is strong and comparable to the case of a standard ${s}_{{he}}^{\pm }$-state. The key messages of figure 7 is:

• (1)  
  The T_cs of two degenerate pairing solutions, the incipient ${s}_{{he}}^{\pm }$-state and the ${s}_{{ee}}^{++}$-state, of the incipient two band model tract each other in clean limit.
• (2)  
  When non-magnetic impurities exist, this degeneracy breaks down and the incipient ${s}_{{he}}^{\pm }$-pairing state quickly becomes suppressed with impurities, while the ${s}_{{ee}}^{++}$-state is robust against the non-magnetic impurity pair-breaking.
• (3)  
  Most interestingly, however, the ${s}_{{ee}}^{++}$-state cannot be stabilized when the incipient hole band approaches too close to Fermi level such as ${\varepsilon }_{b}\lt (\pi {T}_{c}+{{\rm{\Gamma }}}_{{\rm{imp}}}+{{\rm{\Gamma }}}_{{\rm{inela}}})$. This implies that there exists an optimal incipient energy distance ${\varepsilon }_{b}^{{\rm{optimal}}}$ for the ${s}_{{ee}}^{++}$-pairing state.

Now we consider more realistic three band model: one incipient hole band (h), plus two electron bands ($e1,e2$) (see figure 8). In particular, we can study the effect of the inter-electron band interaction ${V}_{e1,e2}$ which might have a non-negligible strength when the magnetic correlation has a deviation from the standard C-type (${\bf{Q}}=(\pi ,0),(0,\pi )$) toward G-type (${\bf{Q}}=(\pi ,\pi )$) [33–35]. If the G-type magnetic correlation is dominant, i.e., when ${V}_{e1,e2}\gt {V}_{{he}1},{V}_{{he}2}$, the leading pairing solution should always be the one in which ${{\rm{\Delta }}}_{e1}$ and ${{\rm{\Delta }}}_{e2}$ have the opposite signs but the same sizes (figures 2(c), (b), (d)) [19, 21–23]. This gap solution is also called as 'nodeless d-wave' or '${d}_{{x}^{2}-{y}^{2}}$', etc in the literature, but in this paper we denote it as '${s}_{e1e2}^{+-}$' to be contrasted to '${s}_{e1e2}^{++}$-state'. In real FeSe systems, it is most likely that the C-type correlation is dominant, but mixed with some fraction of the G-type correlation [32–35]. Therefore, in the three band model studied in this paper, we assumed ${V}_{e1e2}\lt {V}_{{he}}(={V}_{{he}1}={V}_{{he}2})$ along with the assumption of all repulsive (V_ab > 0) inter- and intra-band pair potentials. In this case, we found that there is an interesting transition from the ${s}_{e1e2}^{+-}$-state to the ${s}_{e1e2}^{++}$-state as the incipient energy ε_b varies.

Zoom In Zoom Out Reset image size

In figure 8, we have sketched the typical band structure of the three band model with its possible pairing solutions. As in two band model, we assumed Λ_sf > ε_b. In figures 8(a) and (b), without RG scaling, two possible pairing solutions, ${s}_{{he}1e2}^{-++}$ and ${s}_{{he}1e2}^{0+-}$, are illustrated. The ${s}_{{he}1e2}^{-++}$-state is the same state as the incipient ${s}_{{he}}^{-+}$-state of the two band model in the previous section. The ${s}_{{he}1e2}^{0+-}$-state in figure 8(b) is subtle. First, we found that this pairing state with OPs ${{\rm{\Delta }}}_{e1}^{+}$ and ${{\rm{\Delta }}}_{e2}^{-}$, with opposite signs each other on the electron bands e1 and e2, is possible even with ${V}_{e1,e2}\lt {V}_{{he}1},{V}_{{he}2}$, as will be shown with numerical calculations. Second, having the OPs ${{\rm{\Delta }}}_{e1}^{+}$ and ${{\rm{\Delta }}}_{e2}^{-}$ with opposite signs, the OP on incipient hole band Δ_h can be anything: positive, negative, or zero, without altering the pairing energy with the pair potentials V_he1 and V_he2. But when considering V_hh > 0, Δ_h = 0 is the best solution. However, this state is still not exactly the same state as the ${s}_{e1e2}^{+-}$-state in figure 8(d) because they have different physical pairing cutoffs, Λ_phys = Λ_sf and Λ_phys = ε_b, respectively. In figures 8(c) and (d), two obviously possible solutions, ${s}_{e1e2}^{++}$ and ${s}_{e1e2}^{+-}$, with electron bands only are illustrated. As indicated by the vertical red arrows, the RG flows are: (a) to (c), and (b) to (d), but crossing the flow between them are not possible. Therefore the T_cs of (a) and (c) are equal and the T_cs of (b) and (d) are equal, respectively.

Figure 9 shows the numerical results of a representative case of the three band model with the dimensionless bare pairing potentials ${N}_{{ab}}{V}_{{ab}}^{0}=\sqrt{{N}_{a}{N}_{b}}{V}_{{ab}}^{0}$: $\sqrt{{N}_{h}{N}_{e}}{V}_{{he}}^{0}=1.0$, $\sqrt{{N}_{e1}{N}_{e2}}{V}_{e1e2}^{0}=0.6$, ${N}_{h}{V}_{{hh}}^{0}=0.25$ and ${N}_{e}{V}_{{ee}}^{0}=0.25$, where we assumed ${N}_{e1}={N}_{e2}$, ${V}_{{ee}}^{0}={V}_{e1e1}^{0}={V}_{e2e2}^{0}$, and ${V}_{{he}}^{0}={V}_{{he}1}^{0}={V}_{{he}2}^{0}$. In figure 9(a), the renormalized potentials N_abV_ab(Λ = ε_b) are plotted as functions of ε_b. These potentials are used to calculate T_cs of the ${s}_{e1e2}^{++}$ and ${s}_{e1e2}^{+-}$ states illustrated in figures 8(c) and (d) with the physical cutoff Λ_phys = ε_b. The results of T_c versus ε_b are plotted in figure 9(b): ${s}_{e1e2}^{++}$ (solid red pentagons) and ${s}_{e1e2}^{+-}$ (open red pentagons), respectively. when the incipient hole band is deep (for large ε_b), ${s}_{e1e2}^{+-}$-state (nodeless d-wave) has the highest T_c, and as the incipient hole band becomes shallow (for smaller ε_b), ${s}_{e1e2}^{++}$-state becomes winning. This transition happens when the renormalized inter-electron band potential ${V}_{e1e2}({\varepsilon }_{b})$ turns into negative as indicated by black arrow in figure 9(a).

Zoom In Zoom Out Reset image size

Interestingly, figure 9(b) shows that the T_c of the ${s}_{e1e2}^{+-}$-state (nodeless d-wave; open red pentagons) doesn't change as ε_b varies. This behavior can be understood, however, if we remember that the T_c of figure 8(b) and the T_c of figure 8(d) should be the same due to the RG invariance. In the pairing state of figure 8(b), i.e. ${s}_{{he}1e2}^{0+-}$-state, we argued Δ_h = 0. Then it is easy to note that with the given bare values of ${N}_{{ab}}{V}_{{ab}}^{0}$ and Λ_sf, varying ε_b has no effect on the gap equation and T_c. On the other hand, in the cases of ${s}_{{he}1e2}^{-++}$ (figure 8(a)) and ${s}_{e1e2}^{++}$ (figure 8(c))—they have the same T_c in clean limit—varying ε_b should strongly affect T_c because ${{\rm{\Delta }}}_{h}\ne 0$ in the ${s}_{{he}1e2}^{-++}$-state. Figure 9(b) shows this behavior of T_c versus ε_b for the ${s}_{{he}1e2}^{-++}$ (black stars) and ${s}_{e1e2}^{++}$ (solid red pentagons) states. As in the case of two band model, the T_cs of the incipient ${s}_{{he}1e2}^{-++}$-state with the fixed pairing cutoff energy Λ_phys = Λ_sf continuously increases as ${\varepsilon }_{b}\to 0$. However, the T_cs of the ${s}_{e1e2}^{++}$-state with the pairing cutoff energy Λ_phys = ε_b becomes ill-defined when ε_b < T_c.

Now we would like to consider the non-magnetic impurity effect on the T_c of ${s}_{{he}1e2}^{-++}$-, ${s}_{e1e2}^{++}$-, and ${s}_{e1e2}^{+-}$-states. As we have argued in the two band model, the T_c of the ${s}_{e1e2}^{++}$-state is immune to the non-magnetic impurity scattering. The impurity effect on the ${s}_{e1e2}^{+-}$-state (nodeless d-wave) is mathematically equivalent to the case of the d-wave state, hence we expect the strongest T_c suppression. Finally, the impurity effect on the incipient ${s}_{{he}1e2}^{-++}$-state is the same to the case of the incipient ${s}_{{he}}^{-+}$-state in two band model. The impurity theory for the two band model can be straight forwardly used for the three band model with a replacement of ${N}_{e}={N}_{e1}+{N}_{e2}$ (see appendix A).

Figure 10(a) is the same plot as figure 9(b)—T_c versus ε_b for three pairing states: ${s}_{{he}1e2}^{-++}$, ${s}_{e1e2}^{++}$, and ${s}_{e1e2}^{+-}$—with the same parameters as in figure 9 except for varying values of ${N}_{e1e2}{V}_{e1e2}^{0}=0.5,0.6,0.7,$ and 0.8, respectively. It is busy plot but easy to understand the general trend. First, the T_c of the ${s}_{e1e2}^{+-}$-state (nodeless d-wave; open symbols) increases with increasing ${N}_{e1e2}{V}_{e1e2}^{0}$. On the contrary, the T_cs of the ${s}_{{he}1e2}^{-++}$ (black stars) and ${s}_{e1e2}^{++}$-state (solid symbols) decreases with increasing ${N}_{e1e2}{V}_{e1e2}^{0}$, as expected. Secondly, as already shown in figure 9(b), there are crossovers of the higher T_c pairing state as ε_b decreases from the ${s}_{e1e2}^{+-}$-state (nodeless d-wave; open symbols) to the ${s}_{e1e2}^{++}$-state (solid symbols) for each value of ${V}_{e1e2}^{0}$. Finally, the T_cs of the ${s}_{e1e2}^{++}$-state is ill-defined when ε_b < T_c.

Zoom In Zoom Out Reset image size

Figure 10(b) shows the same calculations of figure 10(a) but including non-magnetic impurity scattering. We assumed unitary limit scatterers (c = 0; see appendix A) for all inter- and intra-band scatterings and the impurity scattering rate Γ = 0.3Λ_sf. This is quite a large scattering rate and this value was chosen to kill all T_c of the ${s}_{e1e2}^{+-}$-state (nodeless d-wave) just for illustration. While the T_cs of the ${s}_{e1e2}^{++}$-state remain the same as in figure 10(a), the T_cs of the incipient ${s}_{{he}1e2}^{-++}$-state are strongly suppressed. As a result, the incipient ${s}_{{he}1e2}^{-++}$-state (open symbols) survives only in the region of the left corner of small values of ε_b in figure 10(b).

The results of figure 10(b) imply the following simple picture. If we can change the depth of the incipient hole band, ε_b, by electron doping, pressure, or dosing [24–28] in a typical FeSe system, the shallow incipient hole band system (small ε_b) can support the incipient ${s}_{{he}1e2}^{-++}$-state with low T_c, while the ${s}_{e1e2}^{++}$-state can appear with much higher T_c with increasing ε_b to the optimal value ${\varepsilon }_{b}^{{\rm{optimal}}}$. Further increasing impurity scattering rate, all incipient ${s}_{{he}1e2}^{-++}$-state disappears and only the ${s}_{e1e2}^{++}$-state will survive with high T_c in the region of ${\varepsilon }_{b}^{{\rm{optimal}}}$. Incidently, these results of figure 10(b) looks very similar to the recent experimental observations of the phase diagram of electron doped FeSe systems, which shows the curious double dome and single dome structure of the T_c versus electron doping phase diagram [24–26].

The impurity effect enters the pair susceptibilities $\tilde{\chi }{(T)}_{h,e}$ for T < T_c in the gap equations equation (5) in the main text as follows,

$$\begin{eqnarray}&&{\tilde{\chi }}_{h}(T)=-T\sum _{n}{N}_{h}{\int }_{-{{\rm{\Lambda }}}_{{sf}}}^{-{\varepsilon }_{b}}{\rm{d}}\xi \displaystyle \frac{{\tilde{{\rm{\Delta }}}}_{h}(k)}{{\tilde{\omega }}_{n}^{2}+{\xi }^{2}+{\tilde{{\rm{\Delta }}}}_{h}^{2}(k)},\end{eqnarray} \tag{ B.1 }$$

$$\begin{eqnarray}&&{\tilde{\chi }}_{e}(T)=-T\sum _{n}{N}_{e}{\int }_{-{{\rm{\Lambda }}}_{{sf}}}^{{{\rm{\Lambda }}}_{{sf}}}{\rm{d}}\xi \displaystyle \frac{{\tilde{{\rm{\Delta }}}}_{e}(k)}{{\tilde{\omega }}_{n}^{2}+{\xi }^{2}+{\tilde{{\rm{\Delta }}}}_{e}^{2}(k)},\end{eqnarray} \tag{ B.2 }$$

where ${\tilde{\omega }}_{n}$ and ${\tilde{{\rm{\Delta }}}}_{h,e}$ contain all selfenergy corrections due to the impurity scattering. Notice the difference of the quasiparticle energy integration domains, ${\int }_{-{{\rm{\Lambda }}}_{{sf}}}^{-{\varepsilon }_{b}}$ and ${\int }_{-{{\rm{\Lambda }}}_{{sf}}}^{{{\rm{\Lambda }}}_{{sf}}}$, between the hole and electron band, respectively. For a standard ${s}_{{he}}^{-+}$-state with ordinary hole- and electron-band, where both bands has the integration domains of ${\int }_{-{{\rm{\Lambda }}}_{{sf}}}^{{{\rm{\Lambda }}}_{{sf}}}$, ${\tilde{\omega }}_{n}$ and ${\tilde{{\rm{\Delta }}}}_{h,e}$ can be calculated using the ${ \mathcal T }$-matrix method as [38–41]

$$\begin{eqnarray}&&{\tilde{\omega }}_{n}={\omega }_{n}+{{\rm{\Sigma }}}_{h}^{0}({\omega }_{n})+{{\rm{\Sigma }}}_{e}^{0}({\omega }_{n}),\end{eqnarray} \tag{ B.3 }$$

$$\begin{eqnarray}&&{\tilde{{\rm{\Delta }}}}_{h,e}={{\rm{\Delta }}}_{h,e}+{{\rm{\Sigma }}}_{h}^{1}({\omega }_{n})+{{\rm{\Sigma }}}_{e}^{1}({\omega }_{n}),\end{eqnarray} \tag{ B.4 }$$

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{h,e}^{0,1}({\omega }_{n})={\rm{\Gamma }}\cdot {{ \mathcal T }}_{h,e}^{0,1}({\omega }_{n}),\,\,{\rm{\Gamma }}=\displaystyle \frac{{n}_{{\rm{imp}}}}{\pi {N}_{{\rm{tot}}}},\end{eqnarray} \tag{ B.5 }$$

where ${\omega }_{n}=T\pi (2n+1)$ is the Matsubara frequency, n_imp the impurity concentration, and ${N}_{{\rm{tot}}}={N}_{h}+{N}_{e}$ is the total DOS. The ${ \mathcal T }$-matrices ${{ \mathcal T }}^{\mathrm{0,1}}$ are the Pauli matrices τ^0,1 components in the Nambu space, and are written for two band superconductivity as

$$\begin{eqnarray}&&{{ \mathcal T }}_{a}^{i}({\omega }_{n})=\displaystyle \frac{{G}_{a}^{i}({\omega }_{n})}{D}\quad (i=0,1;\,\,a=h,e),\end{eqnarray} \tag{ B.6 }$$

$$\begin{eqnarray}&&D={c}^{2}+{[{G}_{h}^{0}+{G}_{e}^{0}]}^{2}+{[{G}_{h}^{1}+{G}_{e}^{1}]}^{2},\end{eqnarray} \tag{ B.7 }$$

$$\begin{eqnarray}&&{G}_{a}^{0}({\omega }_{n})=\displaystyle \frac{{N}_{a}}{{N}_{{\rm{tot}}}}\displaystyle \frac{{\tilde{\omega }}_{n}}{\sqrt{{\tilde{\omega }}_{n}^{2}+{\tilde{{\rm{\Delta }}}}_{a}^{2}(k)}},\end{eqnarray} \tag{ B.8 }$$

$$\begin{eqnarray}&&{G}_{a}^{1}({\omega }_{n})=\displaystyle \frac{{N}_{a}}{{N}_{{tot}}}\displaystyle \frac{{\tilde{{\rm{\Delta }}}}_{a}}{\sqrt{{\tilde{\omega }}_{n}^{2}+{\tilde{{\rm{\Delta }}}}_{a}^{2}(k)}},\end{eqnarray} \tag{ B.9 }$$

where $c=\cot {\delta }_{0}=1/[\pi {N}_{{\rm{tot}}}{V}_{{\rm{imp}}}]$ is a convenient measure of scattering strength and we assumed the equal strength for both inter- and intra-band impurity scatterings. c = 0 means the unitary limit (${V}_{{\rm{imp}}}\to \infty$) and c > 1 is the Born limit (V_imp ≪ 1) scattering.

Now we need a modification of the above formulas for the incipient hole band model. For T_c-equation (${{\rm{\Delta }}}_{h,e}\to 0$ limit), all differences come from the restricted integration of the incipient hole band below Fermi level as

$$\begin{eqnarray}&&{N}_{h}{\int }_{-{{\rm{\Lambda }}}_{{sf}}}^{-{\varepsilon }_{b}}{\rm{d}}\xi \displaystyle \frac{1}{{\omega }_{n}^{2}+{\xi }^{2}}\approx {N}_{h}\displaystyle \frac{1}{| {\omega }_{n}| }{\tan }^{-1}\left(\displaystyle \frac{{{\rm{\Lambda }}}_{{sf}}}{{\varepsilon }_{b}}\right),\end{eqnarray} \tag{ B.10 }$$

compared to the standard two band case with the ordinary hole band,

$$\begin{eqnarray}&&{N}_{h}{\int }_{-{{\rm{\Lambda }}}_{{sf}}}^{{{\rm{\Lambda }}}_{{sf}}}{\rm{d}}\xi \displaystyle \frac{1}{{\omega }_{n}^{2}+{\xi }^{2}}\approx {N}_{h}\displaystyle \frac{\pi }{| {\omega }_{n}| }.\end{eqnarray} \tag{ B.11 }$$

Therefore, by replacing N_h by ${N}_{h}^{{\rm{eff}}}={N}_{h}{\tan }^{-1}\left(\tfrac{{{\rm{\Lambda }}}_{{sf}}}{{\varepsilon }_{b}}\right)/\pi$, all the above formulas for the standard two band model can be continuously used without further changes. We only need to redefine: ${N}_{{\rm{tot}}}=({N}_{h}^{{\rm{eff}}}+{N}_{e})$, and ${\tilde{N}}_{h}={N}_{h}^{{\rm{eff}}}/{N}_{{\rm{tot}}}$, ${\tilde{N}}_{e}={N}_{e}/{N}_{{\rm{tot}}}$, and the strength of the impurity scattering rate ${\rm{\Gamma }}=\tfrac{{n}_{{\rm{imp}}}}{\pi {N}_{{\rm{tot}}}}$ should be understood with newly defined ${N}_{{\rm{tot}}}=({N}_{h}^{{\rm{eff}}}+{N}_{e})$. Being ${N}_{h}^{{\rm{eff}}}\lt {N}_{h}$, it succinctly captures all the effects arising from the sunken hole band by ε_b.

To determine T_c, we take $T\to {T}_{c}$ limit and linearize the gap equation (5) in the main text with respect to the OPs ${{\rm{\Delta }}}_{h,e}$. First, the impurity renormalized Matsubara frequency equation (B.3) and the OPs equation (B.4) are written as

$$\begin{eqnarray}&&{\tilde{\omega }}_{n}={\omega }_{n}(1+{\eta }_{\omega }),\end{eqnarray} \tag{ B.12 }$$

$$\begin{eqnarray}&&{\tilde{{\rm{\Delta }}}}_{h,e}={{\rm{\Delta }}}_{h,e}(1+{\delta }_{h,e}),\end{eqnarray} \tag{ B.13 }$$

with

$$\begin{eqnarray}&&{\eta }_{\omega }=\displaystyle \frac{{\rm{\Gamma }}}{1+{c}^{2}}\displaystyle \frac{1}{| {\omega }_{n}| },\end{eqnarray} \tag{ B.14 }$$

$$\begin{eqnarray}&&{\delta }_{h,e}=\displaystyle \frac{{\rm{\Gamma }}}{1+{c}^{2}}\displaystyle \frac{1}{| {\omega }_{n}| }\displaystyle \frac{[\tilde{{N}_{h}}{{\rm{\Delta }}}_{h}+\tilde{{N}_{e}}{{\rm{\Delta }}}_{e}]}{{{\rm{\Delta }}}_{h,e}}.\end{eqnarray} \tag{ B.15 }$$

And the pair susceptibilities equations (B.1) and (B.2) are now simplified as

$$\begin{eqnarray}&&{\tilde{\chi }}_{h,e}(T)=-\pi {{TN}}_{{\rm{tot}}}{\tilde{N}}_{h,e}\sum _{n}\displaystyle \frac{{{\rm{\Delta }}}_{h,e}(k)(1+{\delta }_{h,e})}{| {\omega }_{n}(1+{\eta }_{\omega })| }.\end{eqnarray} \tag{ B.16 }$$

It is immediately clear that if η_ω = δ_a, as in a single band s-wave state, there is no renormalization of the pair susceptibility ${\tilde{\chi }}_{a}(k)$ with the impurity scattering. This is just the Anderson theorem of T_c-invariance of the s-wave SC. For a d-wave case, obviously δ_a = 0 and ${\eta }_{\omega }\ne 0$, hence results in the maximum T_c-suppression. In the case of a standard ${s}_{{he}}^{+-}$-wave, it was shown that $| {\delta }_{a}| \approx 0$ because of the inverse relation of $\tfrac{{N}_{h}}{{N}_{e}}\approx \tfrac{| {{\rm{\Delta }}}_{e}| }{| {{\rm{\Delta }}}_{h}| }$ and equation (B.15), in the limit of the dominant interband pairing (${V}_{{he},{eh}}\gg {V}_{{hh},{ee}}$) [38].

In our incipient two band case, it is more complicated to draw a simple conclusion. However, as we have shown above, after replacing ${N}_{h}\to {N}_{h}^{{eff}}$, all the formulas are the same as in the case of the standard ${s}_{{he}}^{+-}$-wave state. The remaining question is whether the inverse relation $\tfrac{{N}_{h}^{{eff}}}{{N}_{e}}\approx \tfrac{| {{\rm{\Delta }}}_{e}| }{| {{\rm{\Delta }}}_{h}| }$ is still hold or not, and we found that this relation is still hold in the incipient ${s}_{{he}}^{+-}$-state with numerical calculations [16, 17]. With the susceptibilities equation (B.16) together with the gap equations (5) in the main text in the limit ${{\rm{\Delta }}}_{a}\to 0$, we have calculated T_c in figure 7 with non-magnetic impurity scattering. For convenience of parametrization, we used the impurity scattering rate parameter ${\rm{\Gamma }}={n}_{{\rm{imp}}}/(\pi {N}_{{\rm{tot}}})$ with ${N}_{{\rm{tot}}}=({N}_{h}+{N}_{e})$ in figure 7 instead of using the physically more relevant parameter ${\rm{\Gamma }}={n}_{{\rm{imp}}}/(\pi {N}_{{\rm{tot}}})$ with ${N}_{{\rm{tot}}}=({N}_{h}^{{\rm{eff}}}+{N}_{e})$, which is a complicate function of ε_b.

Noticing that ${N}_{e1}={N}_{e2}$ and ${{\rm{\Delta }}}_{e1}={{\rm{\Delta }}}_{e2}$ (see figure 2(a)), the above formulas for the two band incipient ${s}_{{he}}^{-+}$-case can be used only with the replacement of ${N}_{e}=2{N}_{e1,e2}$. And the result of T_c-suppression is qualitatively the same as the incipient two band ${s}_{{he}}^{+-}$-state.

Noticing that this is a d-wave (nodeless) state (see figure 2(c)), it is always ${\delta }_{a=e1,e2}=0$, hence results in the maximum T_c-suppression as in a standard d-wave case.

Finally, the most focused pairing state of this paper, the three band ${s}_{e1e2}^{++}$-state, having ${{\rm{\Delta }}}_{e1}={{\rm{\Delta }}}_{e2}$ with the same sign and Δ_h = 0, this state is the same as the single band s-wave state, hence should satisfy the Anderson's theorem [36] of the T_c-invariance with the time-reversal invariant disorders such as the non-magnetic impurities.