Pairing Mechanism of the Heavily Electron Doped FeSe Systems: Dynamical Tuning of the Pairing Cutoff Energy

We studied pairing mechanism of the heavily electron doped FeSe (HEDIS) systems, which commonly have one incipient hole band -- a band top below the Fermi level by a finite energy distance $\epsilon_b$ -- at $\Gamma$ point and ordinary electron bands at $M$ points in Brillouin zone (BZ). We found that the system allows two degenerate superconducting solutions with the exactly same $T_c$ in clean limit: the incipient $s^{\pm}_{he}$-gap ($\Delta_h^{-} \neq 0$, $\Delta_e^{+} \neq 0$) and $s_{ee}^{++}$-gap ($\Delta_h =0$, $\Delta_e^{+} \neq 0$) solutions with different pairing cutoffs, $\Lambda_{sf}$ (spin fluctuation energy) and $\epsilon_b$, respectively. The $s_{ee}^{++}$-gap solution, in which the system dynamically renormalizes the original pairing cutoff $\Lambda_{sf}$ to $\Lambda_{phys}=\epsilon_b$ ($<\Lambda_{sf}$), therefore actively eliminates the incipient hole band from forming Cooper pairs, but without loss of $T_c$, becomes immune to the impurity pair-breaking. As a result, the HEDIS systems, by dynamically tuning the pairing cutoff and selecting the $s_{ee}^{++}$-pairing state, can always achieve the maximum $T_c$ -- the $T_c$ of the degenerate $s^{\pm}_{he}$ solution in the ideal clean limit -- latent in the original pairing interactions, even in dirty limit.

We studied pairing mechanism of the heavily electron doped FeSe (HEDIS) systems, which commonly have one incipient hole band -a band top below the Fermi level by a finite energy distance ε b -at Γ point and ordinary electron bands at M points in Brillouin zone (BZ). We found that the system allows two degenerate superconducting solutions with the exactly same T c in clean limit: the incipient s ± he -gap (∆ − h = 0, ∆ + e = 0) and s ++ ee -gap (∆ h = 0, ∆ + e = 0) solutions with different pairing cutoffs, Λ s f (spin fluctuation energy) and ε b , respectively. The s ++ ee -gap solution, in which the system dynamically renormalizes the original pairing cutoff Λ s f to Λ phys = ε b (< Λ s f ), therefore actively eliminates the incipient hole band from forming Cooper pairs, but without loss of T c , becomes immune to the impurity pair-breaking. As a result, the HEDIS systems, by dynamically tuning the pairing cutoff and therefore selecting the s ++ ee -pairing state, can always achieve the maximum T c -the T c of the degenerate s ± he solution in the ideal clean limit -latent in the original pairing interactions, even in dirty limit. The discovery of the FeSe/SrTiO 3 monolayer system (T c ≈ 60 − 100K) [1][2][3] and other heavily electron-doped iron selenide (HEDIS) compounds such as A x Fe 2−y Se 2 (A=K, Rb, Cs, Tl, etc.) (T c ≈ 30 − 40K) [4][5][6], (Li 1−x Fe x OH)FeSe (T c ≈ 40K) [7], and pressurized bulk FeSe (T c ≈ 37K) [8] are posing a serious challenge to our understanding of the Iron-based superconductors (IBS). The main puzzles are two: (1) why is T c so high, up 100K ? (2) what is the pairing mechanism and pairing solution with only electron pockets at M point?
Among the HEDIS systems, we think, FeSe/SrTiO 3 monolayer system [1][2][3] has one extra mechanism, which is the long-sought small-angle scattering phonon boost effect [9,10]. Lee et al. [11] have measured the presence of the ferroelectric polar phonon in the SrTiO 3 (STO) substrate and its strong coupling with the conduction band of the FeSe monolayer. Subsequently, theoretical works [12][13][14][15] have elaborated this phonon boost effect specifically to the FeSemonolayer system. This phonon boost effect is theoretically trivial to understand. When there exist a large momentum exchange repulsive interaction V Q -provided by antiferromagnetic(AFM) spin fluctuations -and a small momentum exchange attractive phonon interaction V ph , two pairing potentials do not interfere each other living in the different sectors of momentum space but work together to boost T c of the s ±or d-wave gap solutions as following way [9,10]: 14 Λ˜λ s f s f · Λ˜λ ph ph e −1/λ tot (1) whereλ s f = λ s f /λ tot ,λ ph = λ ph /λ tot , and λ tot = (λ s f + λ ph ). A small momenta q exchanging phonon attractive interaction λ ph and the large momenta Q exchanging AFM repulsive interaction λ AFM , living in different momentum space, do not interfere but cooperate to enhance the total pairing interaction λ tot = λ AFM + λ ph for the s ± -and d-wave pairing channels.
phonon, respectively. The Eq. (1) shows that the phonon coupling λ ph -even if weak strength -entering into the exponent of the exponential as λ tot = (λ s f + λ ph ), its boosting effect of T c can be far more efficient than a simple algebraic addition. Apart from this phonon boost effect, which exists -or particularly strong -only in the FeSe/STO monolayer system, theoretically more challenging question of the HEDIS superconductors is: even without the phonon boosting effect, how these systems can achieve such high T c of 30-40K only with the electron pockets. To clarify this question, we would like to understand the following more specific questions: (1) without the hole pocket at Γ point, what is the pairing mechanism and the pairing state ? (2) does the incipient (sunken) hole band have any specific role or mechanism for such high T c ? (3) why do all HEDIS compounds seem to have the optimal incipiency distance ε optimal b of 60-80 meV for the maximum T c arXiv:1605.01509v3 [cond-mat.supr-con] 3 Jun 2016 ?
In this paper, we provide a noble mechanism as an unified answer to all these questions, that is dynamical tuning of the pairing cutoff energy. We studied a simple model for the HEDIS systems which consists of one incipient hole band and ordinary electron band(s) with a dominant repulsive interband interaction V inter between them. The original pairing interaction cutoff energy Λ s f is assumed to be larger than the incipient energy ε b (< Λ s f ). This model system has been recently studied, and was shown to have the s ± he -gap solution [16,17]. This gap solution forms non-zero gaps both on the incipient hole band and electron band (∆ h = 0, ∆ e = 0), hence it is specifically named as "incipient" s ± he -gap solution. In this paper, we showed that this incipient band model allows another degenerate pairing solution with the same T c as the incipient s ± he -gap solution, but only with the electron bands (∆ h = 0, ∆ e = 0), hence named as s ++ ee -gap solution. For this pairing solution, the incipient hole band is dynamically eliminated by the renormalization process, and the bare pairing interactions V 0 ab (Λ s f ) > 0(a, b = h, e) get renormalized and in particular the bare repulsive intraband interaction (V 0 ee (Λ s f ) > 0) for the electron band evolves into an attractive interaction (V ren ee (Λ = ε b ) < 0) with a reduced pairing cutoff energy Λ = ε b , with which the system can form a s-wave pairing with the electron bands only, namely, the s ++ ee -gap. The idea of the s-wave pairing state (s ++ ee ) in the electron bands only was also proposed in the previous works [18,19] using the functional renormalization group (FRG) technique. Although the FRG approach correctly captured the leading instability of the s ++ ee -pairing state, this method could not distinguish the s ++ ee -state and the incipient-s ± he -state, nor did it clarify the subtle relation between these two pairing states. The FRG techniques traces the RG flow of the leading pairing form factor φ α (k) by reducing the scaling energy Λ until its eigenvalue diverges. However, the form factor φ α (k), being a lattice harmonics in the C 4 symmetric lattice, is the same both for the incipient s ± he -state and the s ++ ee -state. Therefore when the eigenvalue of this form factor diverges at Λ * , the FRG result just indicates T c ∼ Λ * in the pairing channel φ α , so that it cannot determine which of two pairing states -the incipient s ± he -state and the s ++ ee -state -is the real ground state because the FRG technique itself does not determine the physical pairing cutoff Λ phys .
Our important finding is that in the clean limit of the genuine incipient band model these two pairing states -the incipient s ± he -state (∆ h = 0, ∆ e = 0) and s ++ ee -state (∆ h = 0, ∆ e = 0) -are degenerate two physical solutions with the exact same T c but with different physical pairing cutoff energies, Λ phys = Λ s f and Λ phys = ε b , respectively. These two states are both physical and distinct states. On the one hand, our finding is in accord with the important principle of the physical invariance of the renormalization group (RG) transformation [20]. But on the other hand, our finding reveals a new potential of the RG transformation: the RG transformation can go beyond a mathematical technique to conveniently study the low energy physics [21] and the system can actively utilize its own RG flow to determine the physical cutoff energy scale Λ phys and optimize its best ground state. This active mechanism of RG is not new but already known with the renormalization of the Coulomb pseudopotential [22]. The degeneracy of these two pairing solutions of the HEDIS system becomes broken by impurities and the system choose the s ++ ee -gap solution and extracts the maximum T c , potentially stored in the system by avoiding the impurity pair-breaking.
In section II, we illustrated the concept of this dynamical tuning of the cutoff energy scale with the well known example of the Coulomb pseudopotential and the phononmediated BCS superconductor [22]. We then showed that even the phonon-mediated BCS superconductor with a phonon energy ω D needs not to have the physical cutoff Λ phys to be ω D , as commonly believed, but can take arbitrary scale Λ without affecting T c , ∆ sc , and the condensation energy, within the RG scheme. And higher order correction is necessary to determine the physical cutoff as Λ phys = ω D .
In section III, equipped with this new concept of dynamical tuning of cutoff energy by RG, we studied the incipient band model with one incipient hole band with ε b and electron band(s) mediated by dominant interband repulsive potential In subsection III.A, first we studied the minimal two band model and demonstrated that the T c is invariant with scaling Λ < Λ s f . Therefore when the scaling energy Λ crossovers from above to below ε b , the pairing solution continuously changes from the incipient s ± he -state to the s ++ ee -state keeping the same symmetry and same T c . In section III.B, we showed that the non-magnetic impurity scattering severely weakens the incipient s ± he -state, but would not affect the T c of the s ++ eestate with the physical cutoff Λ phys = ε b .
We propose that this is the key mechanism why the HEDIS systems can achieve reasonably high T c of 30-40K with the sunken (incipient) hole band. They can avoid the impurity pair-breaking by dynamical tuning the pairing cutoff energy to ε b . The standard IBS systems with both hole and electron bands crossing the Fermi level cannot have the choice of the s ++ ee -solution but only the s ± he -solution, therefore the standard IBS systems would suffer severe reduction of T c from the inevitable impurities introduced by dopings. Otherwise all standard IBS systems could have achieved much higher T c . On the other hand, the s ++ ee -solution with the sunken hole band has a different drawback; i.e. increasing the incipient distance ε b weakens the pair susceptibility, hence reduces T c . Therefore, we could expect that the maximum T c of the s ++ ee -solution should occur with ε b → 0, but it is not the case in experiments. We found that there exists a mechanism to determine an optimal value of ε optimal b . The above mentioned RG scaling breaks down when ε b becomes too small because the pair susceptibility for the s ++ ee -state, n +ε 2 (k) , (ω n = πT c (2n + 1)), becomes saturated as ε b → T c , which sets the stable minimum cutoff energy scale Λ phys = ε b ≈ T c . In real system, the optimal Λ phys should increase further by other broadening processes as Λ phys = ε optimal b ≈ (πT c + Γ imp + Γ inela ) (where Γ imp is the impurity scattering rate and Γ inelas is the inelastic scattering rate). This is the reason why optimal incipient energy ε optimal b is about 60-80 meV in all HEDIS systems.
In section III.C and D, we studied a more realistic three band model with one incipient hole band and two electron bands e1 and e2. With this model, we could consider another possible pairing solution: the s +− e1e2 -state, also called, "nodeless" d-wave state, in which ∆ e1 = −∆ e2 [23][24][25]. We showed in general that the s +− e1e2 -state is favored when the incipient hole band is deep (larger ε b ) but the s ++ e1e2 -state becomes winning when the incipient hole band becomes intermediate to shallow (optimal ε b ). With non-magnetic impurity scattering, we showed that the s +− e1e2 -state (nodeless d-wave) is most rapidly destroyed, but the s ++ e1e2 -state is immune to it and can survive with high T c in the region of optimal values of ε optimal b . On the other hand, the "incipient" s −++ he1e2 -state -which had the same T c as the s ++ e1e2 -state in clean limit -can survive with much reduced T c in the region of small values of ε b if the impurity scattering is not strong enough to completely kill this pairing state. This double-dome structure of the phase diagram (see Fig.10(B)) is quite similar with recent experiments of electron doped FeSe systems [26][27][28][29][30].

II. RENORMALIZATION OF PAIRING INTERACTIONS
The bare Coulomb repulsion µ 0 , operating up to the plasma frequency ω pl , is renormalized to become pseudopotential µ * , operating up to the phonon frequency ω D (< ω pl ), as follows [22], Using Cooperon propagator ω D ], we obtain the well known result, where N 0 is the density of states (DOS) at Fermi level. Since ω pl ω D , the strong repulsive Coulomb potential µ becomes much weakened Coulomb pseudo-potential µ * µ. Therefore the total BCS pair potential V pair (ω D ) = [V ph + µ * ](< 0) can now be attractive with the common pairing cutoff energy Λ phys = ω D . This is the well known mechanism of how weak phonon attraction V ph (< 0) can overcome the stronger Coulomb repulsion µ 0 (> 0) by retardation (ω D ω pl ). But the important message for us is that the physical pairing cutoff is not necessarily fixed by the boson energy scale of the corresponding interaction (in this example, ω pl ).
A standard theory of the renormalization of the pairing interactions stops here. But now we would like to perform a thought-experiment, namely, we continue to scale down the effective cutoff energy Λ below ω D to see what happens. It is straightforward to continue the RG scaling as For an attractive interaction (V pair < 0), Eq.(4) indicates that the strength of |V pair (Λ)| increases as the cutoff energy Λ decreases (Λ < ω D ). But it is straightforward to show that the T c is invariant with RG flow as T c = 1.14ω D e −1/N 0 |V pair (ω D )| = 1.14Λe −1/N 0 |V pair (Λ)| . Furthermore, as far as the BCS limit (∆ sc /Λ 1) holds, the superconducting (SC) gap size ∆ sc is invariant as ∆ sc (ω D ) = ∆ sc (Λ) and the total condensation energy ∆E = 1 2 N 0 ∆ sc (Λ) is also invariant with respect to the RG scaling. Therefore all physical quantities are not affected by this level of RG scaling.
Therefore, at this level, the system has no reason to choose ω D , the Debye frequency, as a physical pairing cutoff energy and higher order corrections in O(∆ sc /Λ) are necessary to determine the true physical cutoff energy Λ phys . Exact calculations of all higher order corrections are difficult but a simple hint comes from the total condensation energy expression, more precise form of which is given as [31]. The second term is the next order correction to the BCS approximation ∆E BCS = − 1 2 N 0 ∆ 2 sc and tells us that more condensation energy is gained with larger cutoff energy Λ when the gap size ∆ sc is the same. Therefore, among the degenerate T c (Λ) solutions, the system would choose the largest cutoff solution that is Λ phys = ω D in this particular case.
With the above exercises we would like to propose the key concept of this paper, i.e., the physical pairing cutoff energy Λ phys of the SC transition is not automatically determined by the physical boson energy scales of the pairing interactions such as ω pl , ω D , or ω s f . But it can be dynamically tuned by the system to maximize the T c and the condensation energy.

III. INCIPIENT BAND MODEL
In this section, we would like to apply the concept of "dynamical tuning of cutoff energy" discussed in the previous section to the incipient band model for the HEDIS systems. At the moment, the most accepted theory of the superconductivity for the Fe-pnictide superconductors is the sign-changing s-wave pairing state (s ± he ) between hole band(s) around Γ and electron band(s) around M points in BZ, mediated by the antiferromagnetic (AFM) fluctuations with the wave vector Q connecting Γ and M points (C-type) [32,33]. However, the HEDIS systems commonly have only electron pockets(s) at M points and the hole pocket is missing at Γ point, which exists only as an incipient hole-band (see Fig.2). However, even without the hole pocket, experimental evidences are that the AFM spin correlation is dominantly the C-type with the characteristic wave vector Q connecting the incipient hole band and the electron band [34]. Therefore, it is a natural attempt to extend the standard paradigm of the s ± he pairing mechanism to the HEDIS systems: a dominant pairing interaction, V he > 0, between the incipient hole band and electron band(s). There is possibly some fraction of deviation from C-type AFM correlation toward G-type AFM correlation with the wave vector Q connecting two electron bands at Q x = (π, 0) and Q y = (0, π) [35][36][37]. Phenomenologically, we will consider this deviation by introducing another weaker repulsive interaction between two electron bands V e1e2 > 0.
The same model but with the pairing cutoff scaled down as Λ < Λ s f . For each value of Λ, H ren (Λ) is defined and the best SC gap solutions are calculated as: incipient s −+ he -gap for ε b < Λ < Λ s f , and s ++ ee -gap for Λ < ε b , with the same T c for all Λ.

A. Incipient Two Band Model
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 Fig.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 Fig.2, only s −+ he -gap ( Fig.2(a)) and s ++ ee -gap ( Fig.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.
In this paper, we assume the incipient energy ε b to be smaller than the spin-fluctuation pairing interaction cutoff Λ s f as illustrated in Fig.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 ( Fig.2(a)) was investigated. Here we will consider the incipient s −+ he -gap solution and s ++ ee -gap solution ( Fig.2(b)) on an equal footing in the two band model. The s +− e1e2 -gap solution (Fig.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.
where the pair susceptibilities are defined as where N h,e are the density of states (DOS) of the hole band and electron band, respectively. Assuming all repulsive pair , the coupled gap equations Eq.(5) with the susceptibilities Eq.(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 |∆ h |/|∆ e | is insensitive to the value of ε b and the gap size of the incipient hole band |∆ h | is comparable to the size of |∆ 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 Λ s f to produce the effective pair potentials V 0 ab (Λ s f ) with cutoff energy Λ s f . In particular, the characteristic AFM spin fluctuation energy scale Λ s f 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. Eq.(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 Λ < Λ s f and define the renormalized model H ren (Λ), as depicted in Fig.3(b). The renormalized pair potential V ab (Λ) is defined by a standard RG process as follows (see Fig.4 and the formal solution is defined bŷ whereV (Λ) is the renormalized 2x2 matrix pair potential with the new cutoff energy Λ, andV 0 is the bare pair potential with the cutoff energy Λ s f defined aŝ Accordingly, the Cooper susceptibilityχ is also 2x2 matrix defined asχ with and Notice that the factor 2 is missing in χ h since the hole band exists only below the Fermi level. Also χ h (Λ s f ; Λ) is defined only up to Λ → ε b , meaning that when the scaling cutoff Λ runs below ε b as Λ < ε b , only the electron band will contribute to the RG flow.
In Fig.5(a), we show the results of the renormalizedV ab (Λ) for the whole range of Λ < Λ s f . In this representative case, we chose ε b = 0.5Λ s f and the bare potentials: all four components ofV ab (Λ) 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 0 hh = N e V 0 ee = 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 .
In Fig.5(b), we showed the calculated T c of the renormalized model H ren (Λ) for all values of Λ < Λ s f depicted in Fig.3(b). Namely, we solved the gap equations Eq.(5) in the limit of ∆ h,e → 0 with the renormalized pair potentials V ab (Λ) and the reduced cutoff Λ. For ε b < Λ < Λ s f , the pairing gap solution is the incipient s −+ he -solution, forming SC OPs both on the hole-(∆ − ) and electron band (∆ + ) (see Fig.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 (Λ), i.e. s ++ ee -solution (see Fig.2(b)). Figure 5(b) shows that the T c (Λ) remain the same all the time with the scaling for Λ < Λ s f . This is a consistent result with the invariance principle of RG [20] but still surprising to be confirmed with the multiband SC model with a complicated pair potentialV ab , in particular, despite the change of the pairing solutions from s −+ he to s ++ ee [21]. 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 cequation Eq.(5) together with the renormalized pair potentialŝ V ab (Λ; T = T c ) Eq.(8) at same temperature, which needs selfconsistent 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 α (Λ)φ * α (k)φ α (k ) 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 ++ eechannel and the s −+ he -channel are symmetry-wise identical and have the same form factor φ α (k) ∼ (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 = Λ s f (the spinfluctuation energy scale). And when Λ s f > ε b (Fig.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 Λ s f < ε b and V 0 ee (Λ s f ) < 0, the pairing solution can be the s ++ ee -state, forming gaps only on the electron bands. However, we have shown in Fig.5(b) that even when Λ s f > ε b , there is no reason to fix the physical pairing cutoff as Λ phys = Λ s f because the whole range of Λ < Λ s f produces the pairing states with the same T c and 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 (Λ)state for ε b < Λ < Λ s f , the system should choose the incipient s ± he (Λ s f ) state with cutoff energy, Λ phys = Λ s f , as a physical ground state. By the same reasoning, among the degenerate solutions of the s ++ ee (Λ)-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 Fig.5(b) and specifically illustrated in Fig.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 Λ s f that is larger than ε b . B. Incipient s ± he -state, s ++ ee -states, and impurity scattering Now we compare two physical ground state solutions of a given incipient two band model with Λ s f > ε b : the incipient s ± he -state with Λ phys = Λ s f and the s ++ ee -state with Λ phys = ε b (see Fig.6). We numerically calculate the T c of these two gap solutions with varying ε b for 0 < ε b < Λ s f . 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 ( T c Λ phys 1). Therefore, we numerically calculate  Figure.7 shows that the T c s of the incipient s ± he -state with Λ phys = Λ s f (black stars) and the T c s 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. One important new feature is that while the T c of the incipient s ± he -state continuously increases as ε b → 0, the T c of the s ++ ee -state becomes ill-defined when ε b → 0. This is because the Cooper susceptibility χ e (Λ = ε b ; T ) gets saturated  below which the s ++ ee -pairing state is not stabilized regardless of the strength of the pair potentialsV ab (Λ). 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 Λ mim ∼ (πT c + Γ imp + Γ inela ), where Γ imp is static impurity scattering rate, and Γ inelas is inelastic scattering rate that can be provided by spin fluctuations as Γ inelas = ImΣ s f (T ).
We now consider the impurity pair-breaking effect on T c of both pairing states and in this paper we considered nonmagnetic 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 scal-ing because the Cooperon propagators (pair susceptibilities) χ h,e (Λ s f ; Λ) defined in Eq.(10) are invariant with the nonmagnetic impurity scattering. The fundamental reason for this invariance has the same origin as the Anderson's theorem [38]: 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 (Λ s f ; Λ) 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 [39] is given in Appendix A. Therefore, the renormalized pair potentialsV (Λ) are the same with and without the non-magnetic impurities.
Then given the same pairing potentialsV (Λ), the nonmagnetic impurity scattering would not change T c of the swave pairing according to Anderson's theorem [38], so that the T c of the s ++ ee -state in Fig.7 will not be affected with nonmagnetic impurities. On the other hand, it is well known that non-magnetic impurities will suppress the T c of the "standard" s ± he -state almost as strong as in the d-wave [40]. We might expect that the T c suppression with non-magnetic impurities on the "incipient" s ± he -state will be weakened compared to the case of the standard s ± he -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 -state is almost as strong as the case of the "standard" s ± hestate. 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 -state is described in Appendix B.
In Fig.7, we plotted the results of T c suppression of the incipient s ± he -state with the different impurity pair-breaking rate Γ/Λ s f = 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 imp he,eh = V imp hh,ee (see Appendix A). As expected, the T c suppression rate is strong and comparable to the case of a standard s ± he -state. The key messages of Fig.7 is: (1) The T c s of two degenerate pairing solutions, the incipient s ± he -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 -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 ε b < (πT c + Γ imp + Γ inela ). This implies that there exists an optimal incipient energy distance ε optimal b for the s ++ ee -pairing state. c , respectively, but crossings are not possible.

C. Incipient three band model
Now we consider more realistic three band model: one incipient hole band (h), plus two electron bands (e1, e2) (see Fig.8). In particular, we can study the effect of the interelectron band interaction V e1,e2 which might have a nonnegligible strength when the magnetic correlation has a deviation from the standard C-type (Q = (π, 0), (0, π)) toward G-type (Q = (π, π)) [35][36][37]. If the G-type magnetic correlation is dominant, i.e., when V e1,e2 > V he1 ,V he2 , the leading pairing solution should always be the one in which ∆ e1 and ∆ e2 have the opposite signs but the same sizes (Fig.2(c), Fig.8(b), Fig.8(d)) [19,[23][24][25]. 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 [34][35][36][37]. Therefore, in the three band model studied in this paper, we assumed V e1e2 < V he (= V he1 = V he2 ) along with the assumption of all repulsive (V ab > 0) inter-and intraband 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.
In Fig.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 Λ s f > ε b . In Fig.8(a) and (b), without RG scaling, two possible pairing solutions, s −++ he1e2 and s 0+− he1e2 , are illustrated. The s −++ he1e2 -state is the same state as the incipient s −+ he -state of the two band model in the previous section. The s 0+− he1e2 -state in Fig.8(b) is subtle. First, we found that this pairing state with OPs ∆ + e1 and ∆ − e2 , with opposite signs each other on the electron bands e1 and e2, is possible  even with V e1,e2 < V he1 ,V he2 , as will be shown with numerical calculations. Second, having the OPs ∆ + e1 and ∆ − 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 Fig.8(d) because they have different physical pairing cutoffs, Λ phys = Λ s f and Λ phys = ε b , respectively. In Fig.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 flow between them are not possible. Therefore the T c s of (a) and (c) are equal and the T c s of (b) and (d) are equal, respectively.
he2 . In Fig.9(a), the renormalized potentials N ab V ab (Λ = ε b ) are plotted as functions of ε b . These potentials are used to calculate T c s of the s ++ e1e2 and s +− e1e2 states illustrated in Fig.8(c) and (d) with the physical cutoff Λ phys = ε b . The results of T c versus ε b are plotted in Fig.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 (ε b ) turns into negative as indicated by black arrow in Fig.9(a).
Interestingly, Fig.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 Fig.8(b) and the T c of Fig.8(d) should be the same due to the RG invariance. In the pairing state of Fig.8(b), i.e. s 0+− he1e2 -state, we argued ∆ h = 0. Then it is easy to note that with the given bare values of N ab V 0 ab and Λ s f , varying ε b has no effect on the gap equation and T c . On the other hand, in the cases of s −++ he1e2 (Fig.8(a)) and s ++ e1e2 (Fig.8(c)) -they have the same T c in clean limit -varying ε b should strongly affect T c because ∆ h = 0 in the s −++ he1e2 -state. Fig.9(b) shows this behavior of T c versus ε b for the s −++ he1e2 (black stars) and s ++ e1e2 (solid red pentagons) states. As in the case of two band model, the T c s of the incipient s −++ he1e2 -state with the fixed pairing cutoff energy Λ phys = Λ s f continuously increases as ε b → 0. However, the T c s of the s ++ e1e2 -state with the pairing cutoff energy Λ phys = ε b becomes ill-defined when ε b < T c .
D. Impurity effect on three band model Now we would like to consider the non-magnetic impurity effect on the T c of s −++ he1e2 -, 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 −++ he1e2 -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). Figure10(a) is the same plot as Fig.9(b) -T c versus ε b for three pairing states: s −++ he1e2 , s ++ e1e2 , and s +− e1e2 -with the same parameters as in Fig.9  (black stars) and s ++ e1e2 -state (solid symbols) decreases with increasing N e1e2 V 0 e1e2 , as expected. Secondly, as already shown in Fig.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 0 e1e2 . Finally, the T c s of the s ++ e1e2 -state is ill-defined when shows the same calculations of Figure10(a) but including non-magnetic impurity scattering. We assumed unitary limit scatterers (c=0; see Appendix A) for all interand intra-band scatterings and the impurity scattering rate Γ = 0.3Λ s f . 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 c s of the s ++ e1e2 -state remain the same as in Fig.10(a), the T c s of the incipient s −++ he1e2 -state are strongly suppressed. As a result, the incipient s −++ he1e2 -state (open symbols) survives only in the region of the left corner of small values of ε b in Fig.10(b). The results of Fig.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 [26][27][28][29][30] in a typical FeSe system, the shallow incipient hole band system (small ε b ) can support the incipient s −++ he1e2 -state with low T c , while the s ++ e1e2 -state can appear with much higher T c with increasing ε b to the optimal value ε optimal b . Further increasing impurity scattering rate, all incipient s −++ he1e2 -state disappears and only the s ++ e1e2 -state will survive with high T c in the region of ε optimal b . Incidently, these results of Fig.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 [26][27][28].

IV. SUMMARY AND CONCLUSIONS
We believe that the FeSe/STO monolayer system is an exception among all HEDIS system, in that it has an extra phonon boost effect to achieve the exceptionally high T c up to ∼100K [3]. In this paper, we focused on the common electronic pairing mechanism with the electron pockets only, the characteristics shared by all HEDIS systems including the FeSe/STO monolayer system. We studied the pairing mechanism of the phenomenological incipient band models: one incipient hole band (h) plus one electron band (e) or two electron bands (e1, e2) with the pairing interactions V 0 ab > 0 -possibly provided by the AFM spin-fluctuations -with the original pairing cutoff energy Λ s f (> ε b ).
We introduced the concept of dynamical tuning of physical cutoff by RG. Using this concept and direct numerical calculations, we found that the incipient band model allows two degenerate SC solutions with the exactly same T c in clean limit: the s ± he -gap (∆ − h = 0, ∆ + e = 0) and s ++ ee -gap (∆ h = 0, ∆ + e = 0) solutions with different pairing cutoffs, Λ phys = Λ s f and Λ phys = ε b , respectively. The s ++ ee -gap solution, with Λ phys = ε b , actively eliminates the incipient hole band from forming Cooper pairs, and becomes immune to the impurity pair-breaking. As a result, the HEDIS systems, by dynamically tuning the pairing cutoff and by selecting the s ++ eepairing state, can always achieve the maximum T c -the T c of the degenerate s ± he solution in the ideal clean limit -latent in the original pairing interactions, even in dirty limit. We also found that there exist an optimal incipient energy ε optimal b of the hole band, below this value the s ++ ee -pairing state cannot be stabilized. We estimated ε optimal b ≈ (πT c + Γ imp + Γ inelas ). With more realistic three band model with one incipient hole band (h) and two electron bands (e1, e2), we also considered additional pairing state: s +− e1e2 , also called as "nodeless d-wave" state. We showed in general that the s +− e1e2 -state is favored when the incipient hole band is deep (large ε b ) but the s ++ e1e2 -state becomes favored when the incipient hole band becomes intermediate to shallow depth (ε b ∼ ε optimal b ). Including non-magnetic impurity scattering, the s +− e1e2 -state (nodeless d-wave) becomes most rapidly destroyed and the "incipient" s −++ he1e2 -state might barely survive with very low T c in the region of small values of ε b . However, the s ++ e1e2 -state, being immune to the non-magnetic impurity pair-breaking, can exist with much higher T c in the region of optimal values of ε optimal b . This double-dome structure of the phase diagram shown in Fig.10(b) and the general trend of the transition of the pairing states: "incipient" s −++ he1e2 -state → s ++ e1e2 -state, with electron doping, accompanied with increasing T c , looks very much similar to the recent experiments of electron doped FeSe systems [26][27][28][29][30].
In conclusion, we showed: (1) The standard paradigm of the IBS superconductivity [32,33] -the s ± he -pairing mediated by a dominant interband repulsion between the hole band(s) around Γ and the electron band(s) around M -continues to operate in the HEDIS systems. (2) The new ingredient of the pairing mechanism in the HEDIS systems is the dynamical tuning of the pairing cutoff Λ phys from Λ s f to ε b by RG scaling, which allows the s ++ e1e2 -pairing state as the ground state of the HEDIS systems. By this, the HEDIS systems can avoid the impurity pair-breaking effect. But the drawback is that the hole band has to sink below Fermi level by ε b , which has to reduce T c . In view of this picture and relatively high T c ∼ 30K − 40K of the HEDIS, all pnictide and chalcogenide IBS seem to suffer severe T c -suppression due to intrinsic impurities, inevitably introduced with dopings; otherwise, the IBS systems in general could have had much higher T c . (3) The FeSe/STO monolayer system is special among other HEDIS systems, which has the additional phonon boost effect [9][10][11][12][13][14][15] on top of the above mentioned common pairing mechanism of the HEDIS systems.
, where ω n = T π(2n + 1) is the Matsubara frequency, n imp the impurity concentration, and N tot = N h + N e is the total DOS.
The T -matrices T 0,1 are the Pauli matrices τ 0,1 components in the Nambu space, and are written for two band superconductivity as where c = cot δ 0 = 1/[πN tot V 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 imp → ∞) 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 (∆ h,e → 0 limit), all differences come from the restricted integration of the incipient hole band below Fermi level as Being N e f f h < N h , it succinctly captures all the effects arising from the sunken hole band by ε b .
To determine T c , we take T → T c limit and linearize the gap equations (5)  It is immediately clear that if η ω = δ a , as in a single band s-wave state, there is no renormalization of the pair susceptibilityχ 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 η ω = 0, hence results in the maximum T c -suppression. In the case of a standard s +− hewave, it was shown that |δ a | ≈ 0 because of the inverse relation of N h N e ≈ |∆ e | |∆ h | and Eq.(B.15), in the limit of the dominant interband pairing (V he,eh >> V hh,ee ) [40].
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 → N e f f h , 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 N e f f h N e ≈ |∆ e | |∆ 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 Eq.(B.16) together with the gap equations (5) in the main text in the limit ∆ a → 0, we have calculated T c in Fig.7 with non-magnetic impurity scattering. For convenience of parametrization, we used the impurity scattering rate parameter Γ = n imp /(πN tot ) with N tot = (N h + N e ) in Fig.7 instead of using the physically more relevant parameter Γ = n imp /(πN tot ) with N tot = (N e f f h + N e ), which is a complicate function of ε b .

The three band incipient s −++
he1e2 -state Noticing that N e1 = N e2 and ∆ e1 = ∆ e2 (see Fig.2(a)), the above formulas for the two band incipient s −+ he -case can be used only with the replacement of N e = 2N e1,e2 . And the result of T c -suppression is qualitatively the same as the incipient two band s +− he -state.