Eliashberg theory for spin-fluctuations mediated superconductivity -- Application to bulk and monolayer FeSe

We present a novel method for embedding spin and charge fluctuations in an anisotropic, multi-band and full-bandwidth Eliashberg treatment of superconductivity. Our analytical framework, based on the random phase approximation, allows for a selfconsistent calculation of material specific characteristics in the interacting, and more specifically, the superconducting state. We apply this approach to bulk FeSe as representative for the iron-based superconductors and successfully solve for the superconducting transition temperature $T_c$, the gap symmetry and the gap magnitude. We obtain $T_c \approx 6$ K, consistent with experiment ($T_c \approx 8$ K), as well as other quantities in good agreement with experimental observations, thus supporting spin fluctuations mediated pairing in bulk FeSe. On the contrary, applying our approach to monolayer FeSe on SrTiO$_3$ we find that spin fluctuations within the full Eliashberg framework give a $d$-wave gap with $T_c\le 11$ K and therefore cannot provide an explanation for a critical temperature as high as observed experimentally ($T_c \approx 70$ K). Our results hence point towards interfacial electron-phonon coupling as the dominant Cooper pairing mediator in this system.


I. INTRODUCTION
Ever since the discovery of superconductivity in ironbased compounds [1][2][3][4] an enormous effort has been made to understand the prevailing mechanism responsible for Cooper pairing in these materials, both experimentally and theoretically (see Refs. [5][6][7][8][9][10] and references therein). For most members of this family the superconducting transition temperature is rather small compared to the high-T c cuprates, with a few exceptions such as monolayer FeSe on SrTiO 3 (FeSe/STO). There, the onset of superconductivity has been reported at temperatures as large as 60 − 100 K [11][12][13][14][15][16][17], which is an order of magnitude higher than in bulk FeSe (∼ 8 K) [2,18]. For many of the iron-based compounds there is consensus about an underlying unconventional mechanism responsible for superconductivity, although many theoretical investigations are so far based on the linearized Bardeen-Cooper-Schrieffer (BCS) equations at the Fermi level [5,[19][20][21][22][23]. In this way the superconducting gap symmetry can presumably be obtained correctly, but other important experimental aspects, such as T c or the gap magnitude, remain partially elusive depending on the level of approximation, hampering hence also the unambiguous identification of the pairing mechanism. Thus there is an extensive need for a microscopic theory directly applicable to unconventional pairing, which can naturally provide the experimentally accessible characteristics of the system. Superconductivity in bulk FeSe has been intensively studied since its discovery [2]. As one of the intriguing properties of this material, a small nematic distortion of the tetragonal unit cell at low temperatures has been * fabian.schrodi@physics.uu.se † alex.aperis@physics.uu.se ‡ peter.oppeneer@physics.uu.se observed [24] and its role to the complex superconductivity has been much discussed [25][26][27][28][29][30][31]. A second important feature of FeSe is that there is no long range magnetic order at ambient pressure [6,7] yet there are enhanced spin fluctuations observed that signal a proximity to a magnetic phase transition [32][33][34][35]. A similar behavior is observed in several related Fe-based superconductors [9,[36][37][38] but the spin fluctuations without local ordered moment are particularly strong in FeSe [35]. The spin fluctuations in Fe-based bulk superconductors have consequently been investigated theoretically with several approaches [39][40][41][42]. A further hint for an unconventional Cooper pairing mechanism comes from ab initio calculations that predict a small effect of electronphonon coupling in bulk iron pnictides [43,44]. Moreover, the appearance of superconductivity with an unusual s ± gap symmetry in the Fe-based superconductors has drawn much theoretical attention (see, e.g. [45][46][47][48]) and has in particular strengthened the picture of a close connection between spin fluctuations and superconductivity [5,[49][50][51][52][53]. The most plausible scenario for bulk FeSe is dominant spin-fluctuations mediated pairing while nematicity in itself only modifies the magnetic fluctuations and thereby modifies the superconductivity indirectly [30,31].
The case of monolayer FeSe on STO is markedly different from bulk FeSe. Due to substrate doping the Fermi surface of FeSe/STO does not exhibit the same nesting properties as FeSe [13,14,16,54] and the gap symmetry is plain s-wave instead of s ± [55]. As a consequence, there has been an intensive discussion for FeSe/STO recently about the role of spin fluctuations and substrate phonons for the superconducting state [56,57]. On the one hand, it has been argued that a high-energy interfacial phonon mode can give rise to an enhanced coupling to FeSe electrons [16,58,59]. By imposing this assumption in a multiband full-bandwidth Eliashberg formalism many experimentally measured quantities can arXiv:2004.01539v1 [cond-mat.supr-con] 3 Apr 2020 indeed be explained [60,61]. On the other hand, arguments have been presented in favor of an unconventional, spin-fluctuations mechanism [62][63][64], e.g. via an incipient band scenario [65] or orbital selective modifications in quasiparticle weights [66]. However, these previous theories have so far been formulated on the basis of approximations that are tailored to address one specific aspect of the problem. For example, predictions within the incipient band scenario were based on solving isotropic two-band Eliashberg equations where only interband coupling was assumed and the s ± symmetry of the gap was imposed as the only selfconsistent solution, which resulted in a very high T c . In the case of the orbital selective scenario, the focus was on explaining the momentum anisotropy of the superconducting gap on the Fermi surface. This was achieved through a combination of static random phase approximation (RPA) and linearized BCS theory calculations [66]. Evidently, to tackle the multiorbital spin-fluctuation problem at its full capacity, there is a need for a more generally applicable theoretical framework.
We develop here a full Eliashberg theory generalization of the multiorbital Hubbard-type model and show that it provides the amplitude, symmetry and momentum dependence of the superconducting gap over multiple bands, the renormalization of the electron mass and energy for all Brillouin zone (BZ) momenta, electronic energies and temperatures and therefore also the superconducting T c ; all calculated on the same footing. Our theory hence opens the door for treating both phononic and electronic pairing mechanisms on the same footing to settle the question about the dominant pairing mediator.
In the following we introduce a generic way of selfconsistently solving the anisotropic, full-bandwidth and multi-band Eliashberg equations for spin and charge fluctuations on an RPA level. Subsequently, we apply our microscopic theory to bulk FeSe as a representative example for the iron-based superconductors. The only ingredient that is needed to selfconsistently solve for the critical temperature, the gap magnitude and its associated symmetry, is a tight-binding model reliably reproducing density functional theory (DFT) calculations for the electronic dispersions [67]. In addition we examine the case of monolayer FeSe using modified electronic energies [68], while neglecting any possible influence of the substrate phonon. Our results for this material reveal a strong mismatch to experimental findings, in particular, a computed T c similar only to that of bulk FeSe. This leads us to the conclusion that spin fluctuations play a minor role in FeSe/STO only for temperatures characteristic for superconductivity in the parent compound. We hence attribute the high T c to the characteristic interfacial electron-phonon coupling as extensively studied in previous works [16,58,60,61,69].
In the following we introduce in Sec. II the methodology used to compute spin-and charge-fluctuations mediated superconductivity in a full Eliashberg framework. We then directly apply the method to bulk FeSe and compute several properties representative of its superconducting state. In Sec. III we apply the same methodology to monolayer FeSe on STO. Finding that spinfluctuations mediated pairing can explain superconductivity in bulk FeSe but not in FeSe/STO, we analyze deeper the origin of this result and compare to simplified approaches within BCS theory. Our conclusions on the plausible mechanisms for superconductivity in bulk FeSe and FeSe/STO are given in Sec. IV.

II. METHODOLOGY
In this section we present a recipe of how to embed spin and charge fluctuations in a full-bandwidth, multiband and anisotropic Eliashberg theory. As a representative example for the iron-based superconductors we investigate bulk FeSe and solve for the main characteristics of this system in the superconducting state.

A. Electronic energies of bulk FeSe
The full system is modeled asĤ =Ĥ 0 +Ĥ int , where the interacting partĤ int is explained below in Section II B. For the noninteracting part we consider a tightbinding model as introduced in Ref. [67] that describes the electron band energies for the five iron d-orbitalŝ where we use k, p and σ as labels for momentum, orbital character and spin, respectively.ĉ † kpσ (ĉ kpσ ) are electronic creation (annihilation) operators and ξ kp describes the dispersion in orbital space. From the diagonalization ofĤ 0 we find the band dependent energies ξ kn as shown in Fig. 1(a), and retrieve the matrix elements a p kn , which serve as connection between band and orbital space. The tight-binding model used here has tetragonal symmetry, hence we assume that the influence of nematicity on the superconducting state is to first order negligible (cf. Ref. [31]). This assumption is to some extend justified because superconductivity does not compete with orthorhombicity in FeSe, as has been shown via hydrostatic pressure measurements [26]. Here we work in the unfolded Brillouin zone corresponding to the one-FeSe unit cell in real space. The mapping to the folded BZ is indicated in Fig. 1(b) by cyan dashed lines and the corresponding high-symmetry points. As explicitly shown in Fig. 1(c) the hole band at Γ has pure d xy orbital character, while d xz and d yz play important roles on the remaining FS sheets.
With the just discussed dispersions ξ kn and orbital weights a p kn we have the necessary tools at hand to calculate, within linear-response theory, the bare susceptibility of the system which is the same for both spin and charge channels. This in turn will be used for the calculation of the RPA interacting susceptibilities as we show in the next section. Note that the basis vectors fulfill the orthonormality condition p a p * kn a p kn = δ n,n .

B. Linear response in the RPA approximation
The interaction part of our Hamiltonian is given by the intrasite Hubbard-type terms [21,70,71] whereˆ S is is the spin operator for orbital index s at site i, compare Ref. [22,50,72]. The occupation at site i with electrons of orbital index s, for spin σ, isn isσ =ĉ † isσĉisσ , which can be used to getn is = σn isσ . Above we use U (V ) as intraorbital (interorbital) onsite interaction, J is the Hund's rule coupling and J the pair hopping energy. The interactions are related via J = J/2 and V = U − 3J/4 − J , a choice consistent with related works [21][22][23]. We start by calculating the imaginary part of the bare susceptibilities for real frequency ω where we set T = 5 K to evaluate the Fermi-Dirac functions n F (·), and keep this temperature from here on unless noted otherwise. As briefly discussed in Appendix A it is a reasonable assumption to keep T fixed with respect to bare susceptibilities. The delta function in Eq. (3) is treated by using an adaptive smearing method [73], see Appendix A for details. The intraorbital sum is shown in Fig. 2(a) as function of momenta and frequencies. It is apparent that the elementary excitations of our system range to rather large frequencies, that introduce a scale inappropriate to our low-energy theory. As we show below and in Section II D, in the interacting case the high-energy region gives rise to the Stoner continuum, which generally suppresses Cooper pair formation. Next, we want to use the above results to obtain the real part of the bare susceptibility, which can be done by using the Kramers-Kronig relation for spectral functions, where P denotes the principal value. From here we define the static bare susceptibility as which, as in Eq. (4), is calculated from the intra-orbital components only. In Fig. 2(b) we observe dominant peaks at q = X (=(0, π)) when plotting χ 0,stat q in the first quadrant of the BZ. This is easily explained by the relatively enhanced nesting between electron and hole FS pockets, compare Fig. 1. Similarly, hole-hole and electron-electron nesting features at momenta slightly smaller than (π, π) give rise to the two rings around M (=(π, π)) in panel 2(b). We note that such a bare susceptibility as shown in Fig. 2 is rather generic for the family of Fe-based superconductors, see for example Refs. [22,66], and that qualitatively comparable results have been obtained by DFT calculations, too [74].
Within the RPA the spin and charge susceptibilities are defined via Dyson equations, with Stoner tensors U S for spin and U C for charge. The nonzero elements are given by We solve Eqs. (7) and (8) by mapping all four-rank tensors involved to usual (two-rank) matrices, for example U S pq st → U S ij ≡Û S . This leads to matrix equationŝ which allow us to directly identify the Stoner instabilities. More explicitly, we define static susceptibilities, in analogy to Eq. (6), as where the aforementioned mapping is inverted to retrieve back the four-rank tensors corresponding to outcomes of Eqs. (10) and (11). The requirement χ S,stat q > 0 and χ C,stat q > 0 ∀q, i.e. staying below the first Stoner instability, defines a finite region of allowed values in (U, J)-space. In Fig. 3(a) we plot the resulting phase diagram, where the allowed region is drawn in blue. For all remaining parts the Stoner criterion is violated due to spin (cyan), charge (green) or spin and charge (yellow). We test three different ratios for bulk FeSe, each indicated by a solid line in Fig. 3(a): J = U/10 (gray line) describes strongly localized electrons, J = U/2 (purple) resembles a Hund-metal situation, and J = U/6 (red) is a reasonable intermediate choice.
The RPA susceptibilities do not change appreciably with (U, J). This can be explicitly seen in Fig. 3(b) where we show spin results for all three ratios and varying distance from the border of the allowed region in panel (a), using similar color code for U/J. In all three cases we observe two peaks near the X point. These have their origin in the bare susceptibility (compare Fig. 2(b)) and are enhanced when approaching the Stoner instability. From this behavior one can directly conclude that momenta around q = X give rise to the leading instability and will approximately become delta-peaks in the vicinity of the spin-border in (U, J)-space. Further we find increased susceptibilities with growing U , while a change in J leads barely to noticeable modifications. Since we are not interested in cases where J > U there is no need of explicitly discussing the RPA charge susceptibility. The changes are minor in this quantity because we stay always well separated from χ C,stat q < 0. In the following we want to look at the frequency dependence of Eqs. (10) and (11). Let us therefore define the dynamic spin and charge susceptibilities as Note, that these dynamical susceptibilities always pertain to the imaginary parts in our notation. For J = U/2 and a rather critical value of U = 0.827 eV we plot in Fig. 3(c) the spin structure factor calculated from Eq. (14).
We can compare the result of our calculations with the outcome of DFT-dynamical mean field theory (DFT-DMFT) calculations carried out in Ref. [41]. Although the two approaches are rather different, the main characteristics of S q (ω) are actually similar. Starting with the maximum value of S q (ω), which is 22 eV −1 in Fig. 3(c) and 16 eV −1 as obtained in Ref. [41]. Note that our results are scalable with respect to criticality, hence we could fine-tune U to achieve the same maximum structure factor. Both our calculation and that of Ref. [41] reveal an enhanced contribution at X for small frequencies, as well as substantial values for S q (ω) at (π, π) and (π/2, π/2) for larger ω. There are differences along the frequency axis that can be attributed partially to deviating choices for T and U/J, but mainly stems from the different ways how S q (ω) is calculated [41,75]. From the above comparison we can conclude that our results for the spin and charge susceptibilities show the correct main features.

C. Coupling via the spin and charge sectors
Before turning to the full Eliashberg problem we need to calculate band dependent interaction kernels in Matsubara space. To achieve this we first define where we distinguish between a kernel for electron mass and energy dispersion renormalization (+) and the superconducting pairing (−), corresponding to diagonal and off-diagonal elements of the Green's function, respectively. Setting q = k − k and averaging over k we transform these kernels from orbital into band space: Since we are interested in imaginary frequencies on the Matsubara axis, we transform the result of Eq. (19) via the Kramers-Kronig relation where q m = 2πT m is a bosonic frequency. In analogy to Fig. 3(b) we want to understand the influence of U and J on the interaction kernels. To this end we define static intra-and inter-band contributions as where the focus lies on the kernel of the superconducting channel. For similar choices of (U, J) as in Fig. 3(b) (same color code) we show the outcome of Eqs. (21) and (22) in Fig. 4, panels (a) and (b), respectively.
As an overall trend one observes increasing kernel values with growing U . In contrast to the spin susceptibilities plotted in Section II B, this increase applies to values throughout the whole BZ. While peaks at X are still dominant in all panels of both figures, V s,intra q grows large also at Γ with increasing J and U . Furthermore, a nonnegligible background develops as we slowly approach the Stoner instability (the values shown correspond to noncritical regions). This is caused by taking the full Stoner continuum into account when transforming from real to Matsubara frequencies and causes potential problems for Cooper pair formation due to e.g. frustration effects, see Section II D for details.
Next we focus on the Matsubara axis and define frequency dependent counterparts to Eqs. (21) and (22) via The results are drawn in Fig. 4, where panel (c) and (d) corresponds to the intra-and inter-band terms, respectively. In the first, second and third row of these two panels we show respectively the pairs (q, J) = (X, U/2), (Γ, U/6), and (M, U/10). The aforementioned influence of the Stoner continuum is reflected along q m in all panels. With growing U the kernels are enhanced along the full frequency axis resulting in slowly decaying tails, which in turn increases the computational load of Eliashberg calculations significantly. We further observe that a moderate Hund's rule coupling of J = U/6 can lead to attractive intra-band kernels at q = Γ, as seen in the central panel of Fig. 4(c). This sign change in V d,intra q=Γ (iq m ) is lost as soon as U approaches the Stoner instability. Such properties, i.e. the interplay of repulsive and attractive couplings can have an important impact on Cooper pair formation and the superconducting gap symmetry.
To resolve the aforementioned difficulty concerning the Stoner continuum we introduce a frequency cutoff ω cut > 0 to truncate the integration in Eq. (20). The Matsubara frequency dependent kernels we treat from here on are therefore At this stage ω cut can be considered a variational parameter. Its main effect is to controllably remove high energy parts of the magnetic excitation spectrum, especially the incoherent part which should be irrelevant to superconductivity. The need for this cutoff will become clear later below. To make contact between ω cut and the real-frequency dependence of kernels as obtained from Eq. (19), we define Note that V d,· q is used for kernels as function of ω, while V d,· q in Eqs. (23) and (24) are Matsubara frequency dependent.
To show possible consequences arising from the cutoff in Eq. (25) we choose (U, J) = (0.8 eV, U/6) and plot Eq. (26) in Fig. 5(a). At momenta/frequencies where V d,intra q (ω) < 0 the charge contributions dominate, since they enter with negative sign in Eq. (18), and hence the kernel becomes attractive in the superconducting channel. Wherever the spin dominates over the charge content, the coupling is repulsive. For ω ≤ 0.5 eV the spectrum is rather discrete, making it possible to identify clear features at Γ (attractive) and X (repulsive). At larger frequencies substantial increases of V d,intra q (ω) are observed throughout the full BZ, a similar feature as discussed in connection with Fig. 4. Going to ω ∼ 3 eV, we see an enhanced influence of charge fluctuations which make the kernel attractive. A qualitatively similar picture is found in Fig. 5(b) when considering the dynamic inter-band kernel. These graphs indicate that taking the full Stoner continuum into account is not favorable for unconventional superconductivity. Intuitively this becomes clear in the light of FS nesting, which becomes combined with an incoherent background at all q.
In Fig. 5(c) we use different cutoffs for Eq. (25) and show in panel i (ii) the resulting static intra-(inter-) band kernels. From the lower graph we learn that the incoherent background can be directly controlled via decreasing ω cut , which makes sense having in mind the discrete nature at low frequencies of V d,intra q (ω) and V d,inter q (ω) in panels (a) and (b). Similarly interesting, the kernels in i show a sign change at Γ with increasing cutoff. Concerning the superconducting gap, see Section II D, this can lead to different tendencies concerning the favored gap symmetry since small ω cut leads to attractive intra-band coupling on the FS pockets; the limit ω cut → ∞ on the other hand induces a sign change on individual FS sheets, which overall can lead to a different momentum structure of the order parameter. We consider the former situation as more physical, since we employ a low-energy theory, hence any large frequency effects should not drastically change the qualitative picture of our results.
As another direct consequence of ω cut we note the changes along Matsubara frequencies in Fig. 5(d), where we show the dynamic intra-(i) and inter-band (ii) kernels found from Eqs. (23) and (24). Plotting our results at q = Γ we observe in panel (c)i again the sign change with increasing cutoff, as already discussed. Additionally the tails for small ω cut are decaying much faster with q m , which makes the computations more efficient.

D. The superconducting state
We are now in the position to address the selfconsistent Eliashberg problem for spin-fluctuation mediated pairing. The interaction kernels introduced in the previous section are used to solve the following set of coupled and selfconsistent equations: Here we use Z kn (iω m ) as the electronic mass renormalization with fermionic frequencies ω m = πT (2m + 1), Γ kn (iω m ) is the chemical potential renormalization and φ kn (iω m ) the superconducting order parameter. The gap function is found via ∆ kn (iω m ) = φ kn (iω m )/Z kn (iω m ), the zero-frequency component of which is accessible in experiment. Note that all functions in Eqs. (28)-(31) are explicitly momentum, Matsubara frequency, and band dependent. Within our five-band model, this gives rise to 15 coupled selfconsistent Eliashberg equations in total. The above-presented Eliashberg equations are solved selfconsistently without any further approximation, see Appendix A for details. The full mathematical modeling presented in this work has been implemented in the Uppsala Superconductivity code (UppSC) [60,61,[76][77][78][79].
We now consider all three U/J ratios highlighted in Fig. 3(a) and perform a variation in U and ω cut . For each configuration we test several symmetries for initializing the order parameters to make sure that we capture all possible solutions. Interestingly, our calculations show that a sufficiently large Hund's rule coupling is needed for finding φ = 0, i.e. not a single nonzero gap is found for J = U/10 and J = U/6. Contrarily, we find that selfconsistent solutions are possible when choosing J = U/2. In Fig. 6 we plot the maximum zero-frequency gap in (U, ω cut ) space, keeping J = U/2 fixed. The vertical border drawn in white represents the first Stoner instability. Self-consistent solutions are found only in a very confined region of the phase space, pointing towards two characteristic cutoffs that we identify as 0.42 eV and 0.66 eV. This corresponds to an energy range where nothing (0.42 eV) or only a very small fraction (at 0.66 eV) of the Stoner continuum is included when calculating the Matsubara frequency dependent kernels. We therefore see here explicitly what is already discussed in Sec. II C, namely, that including the Stoner continuum does not allow for a selfconsistently obtained superconducting state. This is due to frustration effects [23], caused by an incoherent distribution of the kernel among nearly all possible exchange wave vectors in the BZ. It is also worth noting that lowering ω cut too much again leads to disappearance of superconducting solutions.
Next we perform a temperature variation for the aforementioned two cutoffs to obtain the corresponding transition temperatures. Following the evolution of maximal zero frequency gaps in Fig. 7(a) we find T c 5.6 K for the cutoff ω cut = 0.66 eV (drawn in red), that is closer to the onset of the Stoner continuum. A slightly larger value of T c 6 K is possible for ω cut = 0.42 eV, represented by the blue solid curve in the same graph. As guide for the eye the onset of superconductivity is marked in both cases by dashed lines with respective color code. These results resemble the experimental value of ∼ 8 K remarkably well [2].
We conclude the discussion of bulk FeSe by turning to the gap symmetry, choosing (U, J) as before, ω cut = 0.42 eV, and T = 5 K. Having access to the fully momentum dependent zero-frequency component ∆ kn (0), we project our results on the FS, drawn in Fig. 7(b). As directly evidenced, there is a sign change between electron and hole pockets without any FS nodes. Since the sign on each individual pocket is constant we obtain a global s ± symmetry. Further we note that max kF ∆ kF (0) = |min kF ∆ kF (0)|, which suggests an additional pure s-wave component. Hence, our result for the gap symmetry of superconducting bulk FeSe is s ± + s, matching experimental findings [18]. Concerning the magnitude of our calculated superconducting gap we deviate only very slightly from measured values of ∆ ∼ 1.67 meV [18,80]. The difference to our result of ∆ ∼ 1.4 meV directly explains the small mismatch in the calculated critical temperature. From our selfconsistent results we furthermore observe that superconducting gap values as found experimentally are not primarily related to nematicity [26,31]. Although our FS obeys C 4 symmetry, the main features measured in the nematic (orthorhombic) state are reproduced reliably.
We conclude this part by referring again to Appendix A, where several further aspects of our calculations are discussed. These details concern the mathematical and numerical steps in all subsections presented so far.

III. MONOLAYER FeSe ON STO
Having introduced our method in Sec. II we now want to apply it to FeSe/STO, imposing that spin fluctuations are the only relevant ingredient for the superconducting state in this system. Any influence of the interfacial phonon that presumably plays an important role for superconductivity [16,58,60,61] is hence neglected here. After describing some characteristic properties of FeSe/STO within our framework in Sec. III A, we directly go to the discussion of our selfconsistent results for the superconducting state, in Sec. III B. We continue in Sec. III C by examining the effect of changes in our tight-binding model on the superconducting properties. We compare our results to predictions of BCS theory in Sec. III D and conclude by commenting on how our solutions scale with respect to the proximity to antiferromagnetic criticality in Sec. III E.

A. Basic properties of FeSe/STO
The electronic energies are given by a tight-binding model that we take from Ref. [68]. We show the corresponding energy bands along high-symmetry lines in Fig. 8(a). The lattice distortion arising from the deposition of monolayer FeSe on the substrate is explicitly taken into account [68]. In addition, the hopping parameters are modified in order to move the hole bands present in bulk FeSe to below the Fermi level [60]. The FS, plotted with its orbital content in Fig. 8(b), consists of two electron pockets, which are dominated by d xy character. Compared to experiment the FS sheets are slightly smaller [81], this aspect is further addressed in Sec. III C.
With these energy dispersions we calculate the real and imaginary parts of the bare susceptibility, Eqs. (5) and (3), which serve as input for obtaining the static and dynamic bare susceptibilities of Eqs. (6) and (4), respectively. From Fig. 8(d) we observe that χ 0,stat q is peaked near q = M , which is the wave vector connecting the FS pockets. Compared to Fig. 2(b) we no longer have pronounced contributions at X since the hole bands can no longer be statically connected to the electron sheets at the FS. A small ring around Γ is found due to small wavevectors connecting states within the electron pockets. Turning to the dynamic susceptibility, shown in Fig. 8(c), confirms the aforementioned picture clearly, namely, that the leading excitations are located at M . Small and distinct branches can be seen at small frequencies as more explicitly shown in Fig. 18(a) in Appendix B; for ω > 0.6 eV a continuum of nonnegligible contributions occurs throughout the BZ. The only exception to this is q ∼ Γ, where a minimal χ 0,dyn q is found for all frequencies.
Inserting the bare susceptibilities into Eqs. (10) and (11), we perform a variation in (U, J) to find the onset of charge and magnetic order. The blue region in Fig. 8(e) represents allowed choices for the Hubbard U and the Hund's rule J coupling. Contrarily, the green, yellow and cyan parts of the diagram correspond to a violation of the Stoner criterion due to χ C,stat q < 0, χ C,stat q < 0 and χ S,stat q < 0 or χ S,stat q < 0, respectively. We show in the same panel three ratios considered in the following: J = U/2 (purple line), J = U/6 (red line) and J = U/10 (gray line). When comparing this phase diagram with bulk FeSe, Fig. 3(a), remarkable similarities are observed.
Since the bandwidths of both tight-binding models are comparable, the scales of U and J do not differ much. Further, the overall shapes are very alike, except for a small piece of allowed phase space missing in Fig. 8(e) in the limit J U . From discussions in the previous section we know that a nonzero solution for the superconducting gap is to be expected close to magnetic order. We hence show the dynamic spin susceptibility Eq. (14) in Fig. 8(f) for a rather critical pair (U, J) = (1.5802 eV, U/10). A zoom into the low-energy region is presented in Fig. 18(b) in Appendix B. At frequencies of approximately 0.6 eV the Stoner continuum begins, introducing contributions for all momenta sufficiently away from Γ. For smaller ω we find enhanced susceptibilities at X, which partially represent the connection of FS electron pockets with hole bands below the Fermi level, made possible by a relatively small energy exchange. This aspect is treated in more detail in Sec. III C. Comparing Fig. 8(f) with the dynamic bare susceptibility in panel (c) we identify the leading instabilities at ω ∼ 0.8 eV and ω ∼ 0.1 eV, both around q = M .
The appearance of a Stoner continuum is reflected in interaction kernels similarly as discussed in Sec. II C. We calculate full frequency, momentum and orbitaldependent couplings from Eqs. (17) and (18), which are then used as input for Eq. (25), additionally as function of ω cut . Keeping parameters (U, J) = (1.5802 eV, U/10) as in Fig. 8(f), we show the outcomes at momenta Γ, M , and X in Fig. 9(a), (b) and (c) for several cutoffs. From panel (a) we observe a similar behavior as in bulk FeSe: for ω cut sufficiently small, i.e. below the onset of the Stoner continuum, small-q couplings are attractive in the superconducting channel for Matsubara frequencies close to zero. As soon as the cutoff gets larger than a threshold of roughly 0.8 eV, V d,intra q=Γ becomes repulsive for all q m , thus favoring a sign change on the FS pockets. Combined with dominant contributions at q = M , see Fig. 9(b), such a configuration could still lead to a nonvanishing gap, possibly with unconventional symmetry, but we did not find it in any of our calculations. From these considerations one can, even without solving the Eliashberg equations, qualitatively predict that a nonvanishing order parameter is likely to be found only when one excludes the Stoner continuum; as shown in the next subsection this is indeed what we observe. In the above discussion we omit showing explicitly the interband kernels, since these do not provide further insights.

B. Spin-fluctuations mediated pairing
We solve the coupled set of Eliashberg equations (28)- (30) for the three ratios U/J as indicated in Fig. 8(e). Varying U and ω cut , compare Sec. II D, we are able to find the available phase space for a nonvanishing order parameter. The selfconsistently calculated maximum superconducting gap is shown in Fig. 10(a), (b), and (c) for As function of ω cut there are three dome-like regions allowing for a finite gap. Such behavior points towards distinct branches in V (±) q (ω) nn , which can be constructive or destructive for Cooper pairing. Without the need of plotting these kernels we can understand the underlying mechanism already on the level of the dynamic spin susceptibility, see Fig. 8(f), using U = 1.5802 eV and J = U/10 as example. Note that a direct comparison of frequency values is not appropriate, since the ω-scale in general is shifted to smaller frequencies when going from susceptibilities to interaction kernels. Clearly, the smallω contribution at M in Fig. 8(f) is responsible for the first dome of panel 10(a). For growing cutoff the pairing is suppressed due to substantial contributions that extends to large parts of the BZ, see Fig. 8(f) at ω ∼ 0.6 eV. At the high-frequency end we find that the leading instability at M , in competition with the Stoner continuum, is responsible for the last region (at largest ω cut ) of nonzero gap. The intermediate frequency regime is less easily understood, since it lies directly within the onset of the Stoner continuum.
When we compare the three panels of Fig. 10 we dis- cern that an increase in Hund's rule coupling leads to a smaller phase space, as well as enhanced gap values. For J = U/2 we find that the phase diagram is similar to that observed in bulk FeSe, see Fig. 6. Although the gap sizes are different, the characteristic frequencies at which a nonzero solution is possible, seem to be almost the same. Besides a small region with ∆ = 0 near ω cut ∼ 0.6 eV, we find the most relevant frequency cutoff at around 0.45 eV. This agreement might be explained by similar choices of (U, J) and the fact that both tight-binding models are derived from Ref. [67].
As we show in the inset of Fig. 10(a) the gap size increases approximately linear with U going towards its critical value. We observe similar trends for the other choices of Hund's rule coupling (not shown). The behavior of ∆ close to criticality is further discussed in Sec. III D. Focusing now on panels (a) and (c) of Fig. 10, i.e. on choices J = U/10 and J = U/2, we choose three representative cutoffs for both. Choosing the respective value of U very close to the Stoner instability, we calculate the temperature dependence of the maximum gap and show the results in Fig. 11(a) and (b).
The largest critical temperature for J = U/10, panel (a), is found for ω cut = 0.45 eV as T c = 10.6 K (∆/T c ≈ 4). Both other cutoffs lead to smaller values of T c = 9 K (∆/T c ≈ 3.6) and 6.2 K (∆/T c ≈ 6), as seen from the blue and yellow curves. Results for J = U/2 do not change much in this respect: The maximal critical temperature is found for the intermediate ω cut = 0.54 eV, and is again T c = 10.6 K (∆/T c ≈ 5.6). From Fig. 11 we also learn that a change in cutoff can lead to changes in the ratio of ∆ = lim T →0 max ∆ kn (0) over T c which is usually taken as a measure of how strongly coupled superconductivity is in a system. This can for example be seen in panel (b) when comparing the yellow and the red lines. At T = 5 K the former shows a larger gap size, although the T c is higher for the latter, hence the gap over T c ratio at T → 0 K can not be expected to be the same. We note that a precise limit of zero temperature can not be calculated here due to the associated computational costs.
Next we examine the FS momentum dependence and hence the symmetry of ∆ kn (0), which is a direct result from our selfconsistent calculations. To achieve this we calculate the renormalized FS as which does not noticeably deviate from the bare electron pockets. The band-dependent superconducting gap is projected ontoξ kn in Fig. 11(c), where we choose U = 1.5802 eV and J = U/10 at T = 5 K. As is easily observed the gap follows d-wave symmetry, which is similarly true for all solutions presented for FeSe/STO in this work. Although this particular symmetry has been proposed for monolayer FeSe on STO [82], the magnitude of ∆ kn at the Fermi level as found here is not sufficiently large to account for experimental findings [81,83]. From the specifics of the interaction kernels discussed in connection to Fig. 9 it is worth doing a simplified treatment to explain our calculated gap symmetry. For simplicity we might consider Z kn (iω m ) = 1 and Γ kn (iω m ) = 0 and focus on the order parameter only. φ kn and V (−) q nn are both largest at the zeroth frequency, hence we omit the dependence on the Matsubara axis. In this sense we can write Eq. (30) as Since we are interested in FS properties, and only a single band in our dispersion crosses the Fermi level, we can remove the band index and the associated sum on the right-hand side. As a drastic simplification we may write the potential as sum of delta peaks at the high-symmetry points, such that where V (Γ,M,X) is the approximate kernel size at the respective momentum. We directly can neglect the contribution at X, since no FS points can be connected by this reciprocal space vector. Evaluating the sum over k leads to where we use the sign change of V d,intra q=Γ , as discussed in Sec. III A, and set V (Γ) = −|V (Γ) |. From here it is easy to see that φ kF is maximal if φ kF−M = −φ kF , which is just the symmetry observed in Fig. 11(c).
From the results presented in this section one can conclude that spin fluctuation mediated pairing cannot be the prevailing mechanism in monolayer FeSe on STO. The onset of superconductivity as found here does not compare well with experimentally obtained values of T c ∼ 60 − 100 K [11,12,17], and the maximum calculated superconducting gap is also too small [81,83]. A possible counter argument could be, however, that the specifics of our tight-binding model are not accurately resembling the experimentally observed situation. We therefore check below the influence of the electronic dispersions on our results.

C. Influence of the tight-binding model and role of incipient pairing
Now we want to examine the influence of two specific aspects within the tight-binding model used for the calculation. On the one hand we dope the system by means of a global chemical potential, ξ kn → ξ kn − µ with µ = 140 meV. This significantly increases the size of our electronic FS pockets, while placing the hole bands at rather deep energies, which are not confirmed by experiment, see Ref. [16] and Appendix C. In addition, we want to explore the role that the distance between electron bands at X and hole bands at M plays in the superconducting state. This is achieved by an artificial nonrigid shift of only the hole bands by δµ = −48 meV. For further details on the modified energies see Appendix C. In Fig. 12 we show the maximum superconducting gap as function of U and ω cut . The upper row panels (a), (b), and (c) are obtained by using ξ kn − µ, while the lower ones (d), (e), and (f) represent ξ kn − δµ. The columns correspond to different choices of the Hund's rule coupling: (a) and (d) are for J = U/10, (b) and (e) for J = U/6, and (c) and (f) for J = U/2. Note that we do not show the full frequency range of Fig. 10, since in the current simulations we always find a vanishing gap for ω cut > 1 eV.
Starting from results for ξ kn − µ, we see a slight enhancement of allowed phase space for a selfconsistent solution in the superconducting state, compare Fig. 10. The scales along U are slightly different now due to a change in boundaries regarding the first Stoner instability. This aspect is explicitly discussed in Appendix C. As in the case of our original dispersion, an increase in Hund's rule coupling leads to larger superconducting gaps near the onset of magnetic order, see especially Fig. 12(c). The largest values possible for a rigidly shifted dispersion with increased FS pockets is about 7 meV. This particular configuration (U, J, ω cut ) = (1.16 eV, U/2, 0.42 eV) leads to a critical temperature of T c ∼ 8.3 K, which is the maximal value found in all the parameter space explored. Compared to Fig. 10 these observations point clearly towards an enhanced ratio ∆/T c ≈ 9.8, which can not be expected to be compatible with experiment. Generally we observe that the enhancement in the FS size leads to increased superconducting gaps, which can be explained by the larger DOS at the Fermi level, but not to higher critical temperatures. The FS pocket sizes are very well comparable to the experimentally observed situation, although the resulting upper (lower) position of the hole (electron) bands are not [16]. However, we still observe that spin fluctuations are not sufficient to give the correct characteristics of FeSe/STO.
When examining a decrease in distance between the bottom of the electron band at X and the top of hole bands at M we do not find any superconductivity for J = U/10 in Fig. 12(d). Further, the available phase space for J = U/6 (panel (e)) and J = U/2 (panel (f)) is significantly smaller than in corresponding graphs of Fig. 10. The maximum gap possible for using a nonrigid shift of δµ occurs at J = U/6 and is around 8 meV. We can explain these results by a qualitative comparison to bulk FeSe. Although the hole bands at M do not cross the Fermi level for ξ kn − δµ, the associated coupling to electron pockets is significantly enhanced through an exchange momentum q = X. For the bulk material we find a hole band also at Γ, which in the current case is still too far away from the Fermi level to contribute significantly, compare Fig. 8(a). From above arguments it follows that the leading instability for the nonrigidly shifted dispersion is still at q = M , favoring a sign change between FS pockets, and hence the d-wave solution as observed in Fig. 11(c). Combining this with an enhanced coupling at X there might be a need of changing sign on the incipient hole band pockets, introducing a node. As apparent, such a solution is hard to accomplish, which means we need to choose U very close to its critical value, such that the relative significance of couplings at X is suppressed. This qualitative argument explains why selfconsistent solutions for the superconducting state are found only very close to magnetic order in panels Fig. 12(e) and (f), while the gap size for these confined regions is enhanced, in comparison to Fig. 10, to a maximum of almost 8 meV. For dispersion ξ kn − δµ and J = U/6 we calculate the largest critical temperature as T c ∼ 11.4 K by choosing (U, ω cut ) = (2.48, 0.18) eV, which is a slight enhancement when compared to results obtained for our actual dispersion, though by far not large enough to account for experimental values. We note that this particular value is rather artificial in any case, since corresponding positions of the hole bands are far from angular resolved photoemission spectroscopy (ARPES) measurements [16], see also Appendix C.
It has been proposed that the just-discussed incipient band coupling can provide an explanation for the high critical temperature in FeSe/STO, imposing an uncon-ventional pairing mechanism [65]. Our calculations presented in this work point, however, towards the opposite direction, i.e. that incipient interaction leads to a decrease in the phase space available for superconductivity, while mainly increasing the gap size and hence the ∆/T c ratio. This discrepancy can be attributed to neglected intra-band interactions in Ref. [65], which in our framework provide the dominant contributions, compare Fig. 9, since FS points can be connected mainly via q ∼ Γ and q ∼ M . Moreover, here we do not constrain the symmetry of the superconducting order parameter [65], instead we obtain it from the selfconsistent calculation.
To further examine the relevance of incipient pairing within our theory, we now remove all hole bands in the electron dispersions and couplings artificially on the level of Matsubara kernels from Eq. (25). With only the electron bands (for simplicity referred to as ξ kn ) left we perform our selfconsistent Eliashberg calculations and follow the maximum superconducting gap with temperature, see Fig. 13. We use (U, J) = (1.5802 eV, U/10) and run the simulations for our dispersion ξ kn (dashed) and the reduced ξ kn (solid). Two characteristic cutoff frequencies ω cut = 0.21 eV and ω cut = 0.45 eV are shown in blue and red colors, respectively. The onset of superconductivity is marked via dotted (dashed-dotted) straight lines for ξ kn (ξ kn ).
For cutoff ω cut = 0.21 eV we find a decrease of T c from 9.0 K to 7.6 K due to neglecting the hole bands, while the critical temperature goes from 10.6 K to 8.8 K for frequency ω cut = 0.45 eV. Hence we find that in both cases almost 85% of the pairing stems from the electron bands only. This in turn reveals that the role of hole bands is rather minor in monolayer FeSe, compared to its bulk parent compound. From the results in Fig. 13 and the discussion about changing the dispersion via a nonrigid δµ we can safely conclude that the incipient band scenario put forward in Ref. [65] does not lead to a satisfac- tory enhancement of T c to values near the experimentally observed ones in a full bandwidth Eliashberg theory for spin fluctuations.

D. Simplified calculations
We now want to compare our results to computationally less demanding and more effective theories, as they are rather commonly employed [19,22,23]. For this purpose we use the zero-frequency kernel V nn q = V (−) q (0) nn and start with a linearized BCS equation at the Fermi surface. By decoupling ∆ kn = ∆g k , with g k a global symmetry form factor, we solve for the coupling via [19,22] (36) In Eq. (36) we integrate over FS sheets C n and test form factors g k ∈ {1, cos(k x ) + cos(k y ), cos(k x ) − cos(k y ), cos(k x ) cos(k y ), sin(k x ) sin(k y )}. Among the g k , the prevailing symmetry is determined by the largest coupling. We test this setup for all cutoffs shown in Fig. 10 with an enlarged range for U towards smaller values, and with all three ratios of U/J as in panels (a), (b) and (c) of this figure. Our results reveal that the leading λ has exclusively d-wave symmetry, which implies that our selfconsistent result in Fig. 11(c) is governed by FS contributions. Since on this level of approximation any value λ > 0 leads to a finite T c [20], solutions are found for all U and ω cut , which generally overestimates largely the available phase space allowing for superconductivity. In addition, we find hardly any dependence on ω cut , which is easily explained by the fact that this approximation neglects the frequencies and all momenta away from the Fermi level. A less drastic approximation can be made by neglecting only the Matsubara frequency components but keeping the full momentum dependence. The selfconsistent BCS equation for the superconducting gap is then found as with E kn = ∆ 2 kn + ξ 2 kn [20]. We solve the above equation at T = 5 K as function of U and ω cut . Results for the maximum superconducting gap in case of J = U/10, J = U/6, and J = U/2 are shown in panels (a), (b) and (c) of Fig. 14. The inset of (a) shows the dependence of max ∆ kn on U = 10J for two frequency cutoffs ω cut = 0.5 eV in blue and ω cut = 2.0 eV in red. Comparing these results to Fig. 10 immediately reveals a large increase of the superconducting gap sizes in all three panels. Close to the Stoner instability ∆ reaches values as large as few hundred meV's, which is by no means compatible with experiment. Such a large enhancement is not particularly surprising given that BCS theory has no quantitative use beyond the very weak-coupling limit where the mass renormalization is negligible. In addition the onset of superconductivity occurs at smaller U , while the dependence on ω cut , as observed in Fig. 10, is almost completely lost. As the inset of Fig. 14(a) clearly shows, the gap diverges for U close to the Stoner instability. This behavior translates into almost arbitrary gap magnitudes and corresponding critical temperatures, which is clearly a deficiency of BCS theory. In the following Sec. III E we discuss the aspects of criticality in U within both, BCS and Eliashberg treatments, and explicitly show how the aforementioned problem is cured when solving the full set of Eqs. (28)- (30).
Our calculations in the current section reveal that spinfluctuation mediated pairing is drastically overestimated within theories as BCS that neglect any frequency dependence [62,65,66]. This holds true for the gap magnitude (and hence T c ) and the available phase space with respect to the choice of (U, J).

E. An upper bound for the gap magnitude
The inset of Fig. 10 shows an approximately linear decrease of the superconducting gap as function of U , obtained from our Eliashberg formalism. Contrarily we find approximately a 1/U trend when solving the BCS equation in Sec. III D, compare with the inset of Fig. 14(a). Below we use some qualitative arguments to reproduce these trends, which in case of the Eliashberg theory allows us to give a realistic estimation of the maximally possible gap size in FeSe/STO within our unconventional theory presented here.
Let us start from the BCS Eq. (37). If we confine ourselves to the FS we can set ξ kn = 0, hence E kn = ∆ kn .
Since in the unfolded BZ only a single band is crossing the Fermi level we omit the band index n and get The leading instability occurs at q = M , so we write V q = V δ(q − M ) with V a function of the intraorbital coupling only. Using the sign change among FS pockets which can be solved graphically as function of V = V (U ), see Fig. 15 As an alternative one can start from an estimate of the scaling for V (U ) by considering the spin susceptibility in Eq. (10). Close to the Stoner instabilityχ S q (0) grows as 1/(U crit − U ), with leading contributions at q = M . Inserting this in Eq. (18) and neglecting all but the first contribution, which represents the coupling via the spin degree of freedom, gives a scaling V ∝ U 2 /(U crit − U ).
We fit the zero-frequency maximum kernel for the superconducting channel by and find U crit = 1.594 eV. If ∆ kF is sufficiently large we can neglect the hyperbolic tangent in Eq. (39), and hence the gap grows as ∆ kF ∝ V ∝ U 2 /(U crit − U ). This approach is depicted by blue dashed lines in Fig. 15(a) and (c). In panel (a) we compare our actual data (red, solid) for the maximum static Matsubara kernels with the fitted form of Eq. (40) in blue and dashed, leading to an accurate agreement. This particular functional form serves as input for the comparison in panel (c), using similar color code. As expected the scaling is very accurate for critical U → U crit , while less precise for smaller U . From the above discussions we learn that the scaling of the interaction kernel is directly translated to the order parameter, leading to a 1/U divergence of the superconducting gap as U → U crit in the BCS approximation. Such a behavior does not correspond to the actual physical situation and is an artifact of a too simplistic modeling. Below we show that such a divergence does not occur in our Eliashberg formalism due to the explicit inclusion of the electronic mass renormalization.
We use similar assumptions as before, i.e. we confine Matsubara frequency indices to m = m = 0 and k at the Fermi level. The latter condition within Eliashberg theory translates as ξ kn +Γ kn (0) = 0, hence we only consider Eqs. (28) and (30). Let us assume for the moment that T > T c , then φ kF = 0 and Taking the kernel again as V q = V δ(q − M ), together with Z kF−M = Z kF leads to Z kF ∼ 1 + V /π 2 T Z kF . The solution to this second order polynomial is given by Next, we assume that the mass renormalization does not change significantly when going to the superconducting state. The order parameter can then be simplified as φ kF = T V − π 2 T 2 Z 2 kF , where we use similar arguments that led to Eqs. (39) and (42). Finally, the superconducting gap function then reads In Fig. 16 we present the solutions to fitting our actual data to the simplified dependencies of Eqs. (42) and (43). Panel (a) shows the maximum mass renormalization on the FS as function of intraorbital coupling, where selfconsistently obtained results from the Eliashberg equations which is close to results already found in Fig. 11(a) due to a rather critical choice of U therein.
From these results we learn two important aspects: First, the Eliashberg formalism employed in this work removes the unphysical divergence of the superconducting gap (and hence T c ) as function of U . This is due to a mild divergence of the mass renormalization, which counteracts the scaling of the interaction kernel, such that ∆ grows only linear with U in the instability region close to magnetic order. On the other hand, we can consider our results for the superconducting gap, and hence the transition temperatures, of Fig. 11 as upper bounds already, since all values of U employed for associated calculations have been chosen very close to the Stoner instability.

IV. DISCUSSION AND CONCLUSIONS
Before formulating our conclusions it is appropriate to discuss possible limitations of the here-developed formalism.
First, to start with, there is the influence that the employed tight-binding energy bands could have. For example, the tight-binding models used here do not account for spin-orbit coupling. Nonetheless, we have taken care to simulate the influences due to changes in the near-Fermi energy bands of FeSe/STO in Sec. III C and did not find an appreciable change in the superconductivity characteristics. We have furthermore performed our simulations in the unfolded, instead of the folded BZ due to computational performance, but changes due to this aspect are presumably small as well.
Second, in the theory presented in Sec. II we introduce a frequency cutoff to truncate the real-frequency dependent kernels when transforming into Matsubara space. This particular cutoff is one of the key ingredients to find selfconsistent solutions in the superconducting state. With this procedure, we are able to controllably remove the Stoner continuum-like incoherent end of the magnetic spectrum and the concomitant high energy diverging tendencies which are not relevant to the low-energy superconducting phenomena. Ideally, these high-energy degrees of freedom should be integrated out and used to renormalize the remaining interactions, in a manner similar to the well-known treatment of the Coulomb interaction in the electron-phonon problem that leads to the low-energy Coulomb pseudopotential µ * [20]. Such a procedure would however increase considerably the already high complexity and computational cost of the here-proposed framework and is therefore out of the scope of the present work. The influence of such corrections to our method is expected to be rather minor since our results for bulk FeSe are very accurate already. By enabling the systematic numerical solution of the Eliashberg equations over a broad energy range, our method provides a way of mapping out the relevant parts of the spin-fluctuation spectrum that are important to the pairing and therefore provides new insights to electronic mechanisms of superconductivity.
Third, as a further step closer to the definite answer to the superconductivity mediating mechanisms in bulk FeSe and FeSe/STO, Eliashberg theory calculations, wherein both spin fluctuations and electron-phonon coupling are treated on the same footing are required. Nonetheless, a clear hint towards the outcome of such an investigation is provided in this work, together with studies [16,58,60,61] of phonon-mediated superconductivity in FeSe/STO. The evidence collected in the latter studies speaks clearly in favor of interfacial electron-phonon coupling as dominant contribution to the superconducting gap, and hence its high T c .
Summarizing, we have presented a full Eliashberg treatment of spin-fluctuations mediated superconductivity with broad applicability. While our results for bulk and monolayer FeSe are representative for a wide spectrum in the class of iron-based superconductors, the herepresented Eliashberg formalism provides a way of examining the influence of magnetic fluctuations in arbitrary materials, with and without superconductivity. With a faithful tight-binding model for bulk FeSe as the only input needed for our treatment, we achieve for FeSe good agreement with known experimental quantities in the superconducting state: our maximum gap is 1.4 meV, compared to measurements that yield 1.67 meV [18]. The selfconsistently computed gap symmetry of s ± +s type is also obtained correctly. Further, we find a critical temperature of 6 K, which compares well to the measured T c ∼ 8 K [2]. As the outcomes presented are consistent with the main characteristics of superconductivity in this material, we are confident that our microscopic treatment captures the important physics and thus provides strong support for spin-fluctuations mediated superconductivity in bulk FeSe.
Applying the same methodology to monolayer FeSe on STO we find a clear discrepancy between computed results for the superconducting state and well established experimental facts as e.g. the critical temperature. Although the tight-binding model we use for our main findings might deviate to some degree from reality, we have explicitly tested the influence of changing the FS size and distance between electron and hole bands. For the former situation, i.e. a rigid shift in the electronic dispersions, the maximum gap size is increased to ∼ 7 meV. Changing the electron-hole distance on the other hand leads to even ∼ 8 meV, which is not far from experimental values. The critical temperature, however, does not noticeably increase and stays at the order of 10 K. Our results further indicate that band incipiency has little effect on T c and therefore, solely spin fluctuation in combination with band incipiency cannot explain superconductivity in FeSe/STO. As we showed, the T c value computed with full-bandwidth Eliashberg theory is notably lower than that predicted by BCS theories based on spin fluctuations, due to approximations involved in the latter. We have explicitly shown that within Eliashberg theory there exists an upper limit to the obtained T c as U → U crit , in contrast to BCS theory calculations where the gap values and critical temperatures can be almost arbitrarily increased by choosing U close to U crit of the antiferromagnetic transition. This emphasizes an important quantitative limitation of the RPA-BCS approach one should be aware of when comparing with experiments. Overall, our results lead us to surmise that superconductivity due to spin and charge fluctuations could be possible in FeSe/STO, but only at temperatures slightly larger than the T c ∼ 8 K of bulk FeSe. Consequently, we conclude that a spin-fluctuations mechanism alone cannot explain the observed high T c and that another, dominant pairing mechanism such as interfacial electron-phonon coupling must be responsible.
where we used k = k − q in the first step and replaced k → k in the second step. This expression simply means that the system's linear response respects causality. We hence calculate the imaginary susceptibilities only for ω < 0. The delta function in Eq. (3) is approximated by a Gaussian, where we make use of an adaptive smearing method to obtain reasonable broadenings [73]. Written explicitly, we approximate with the broadening matrix that adapts to the electronic velocities. The parameter α can be chosen close to unity and ∆k is the spacing of the momentum grid [73]. As stated in the main text, the real bare susceptibility is found by using a Kramers-Kronig relation. The integration bounds of ±∞, see Eq. (5), are truncated as follows: All momentum and orbital symmetries valid for the imaginary part translate directly to the real part. The cutoff C > 0 must be chosen such that no high-frequency information is lost. Since the delta function in Eq. (3) peaks at frequencies ω = ξ kn − ξ k+qn we keep all contributions by choosing C > 2 max k,n |ξ kn |. Note that C is not in any way related to the truncation parameter ω cut , which we extensively use in the main text. Let us now turn to characteristic properties of the bare susceptibility. The real part in Eq. (5) can be rewritten analytically by inserting Eq. (3): Such a form is used e.g. in Ref. [22]. It is well known that for q → 0 and vanishing frequencies the susceptibility is equal to the density of states at the Fermi level, provided that also T → 0. It can be shown from Eq. (A7) that we recover this limit. Considering χ 0,stat q from Eq. (6) we find At zero temperature one gets χ 0 ≡ 2χ 0,stat q → n N n (0), where N n (0) is the band-resolved density of states at the Fermi level [22]. Note that we do not use a p kn ∈ R in the above calculation, i.e. the result holds for general tightbinding models.
We plot χ 0 as function of temperature in Fig. 17 for both bulk FeSe (red line) and the monolayer case (blue line). No significant changes for temperatures T ∈ [5, 100] K are observed for both electronic dispersions. We hence assume that thermal broadening effects do not noticeably alter the bare susceptibility found from Eqs. (3) and (5). Therefore, taking T = 5 K, we calculate the susceptibility only once for given electron energies and use the result also for larger temperatures. A distinction with respect to T enters therefore when transforming the interaction kernels from real to Matsubara frequencies, Eq. (25).
For RPA susceptibilities we keep only the imaginary parts from solutions to Eqs. (7) and (8), since the real parts are not needed for subsequent calculations. The momentum and frequency symmetries are similar to the imaginary bare susceptibility. Mapping the four-rank tensors to simple matrices, as mentioned in the main text, is a two-step procedure, which we explain here briefly by using U S pq st as an example. First, the orbital indices within the tensors must be rearranged according to U S pq st → U S qp ts . Note, that this transformation is nontrivial since it does not fall into any orbital symmetry of the bare susceptibility or the Stoner tensors, compare Eqs. (9) and (A2). The second step is the actual mapping, which is not uniquely defined. In our implementation we use where L o is the number of orbitals and Note that this mapping is unidirectional, i.e. from knowing p, q, s, t one cannot solve for the associated i, j analytically. One can, however, at any time convert the matrices back to four-rank tensors numerically, so the mapping can be considered invertible. For implementing Eqs. (17) and (18) we employ again the tensor-matrix mapping as discussed above. The Coulomb terms are present only for V (−) q (ω) pq st , since this is the kernel used for off-diagonal Green's function elements. As is common practice we assume the Coulomb interaction to be taken into account in the tight-binding model, which is why no additional terms appear in the coupling used for diagonal elements of the Green's function. In real-frequency space the terms U S /2 and U C /2 in Eq. (18) are needed for double-counting. As we show below, no such contributions enter the Matsubara space calculation. Inserting the Coulomb terms into Eq. (19) and using a p kn ∈ R gives In Eq. (20) we need to explicitly take the imaginary part of the above, which is identically zero; hence the Coulomb terms do not enter into the kernels used for the Eliashberg equations.
The integral in Eq. (20) is treated similarly as in Eq. (A6). We take only the negative frequency axis into account and set the lower integration bound to −C < −2 max k,n |ξ kn |: When calculating Eq. (25) we replace C by the truncation cutoff ω cut in the above. The kernels in Matsubara space are even in q m and, as all quantities considered in this work except the matrix elements, tetragonal in q. Additionally, Re V (±) q (iq m ) nn is invariant under exchanging the band indices, which is a symmetry directly translated from the real-frequency kernel. Below we show that the latter is invariant under exchanging n ↔ n : where we used a p kn ∈ R and k = k + q. Due to inversion symmetry of V where we also rename k to k. Reshuffling the orbital dummy indices then gives V (±) q (ω) n n = k stpq a t * kn a s * kn V (±) q (ω) st pq a p k−qn a q k−qn The Eliashberg Eqs. (28)- (30) are solved iteratively with a convergence criterion of 10 −9 as threshold for the maximal absolute change in all three functions. Always taking more than 2000 points on the Matsubara axis we are confident to be converged in the number of frequencies. Note that a much larger number is partially employed (and needed), depending on the cutoff used for calculating the interaction kernels, see Secs. II C and III A. In the calculations presented we adjust the number of frequencies as required depending on ω cut . For increasing the numerical performance we employ a Fourier convolution scheme in momenta and frequencies. Since the BZ symmetry of the order parameter is a priori not known one needs to be careful with the initialization before starting the iterative loop. If the initial guess does not obey the favored symmetry the algorithm might not converge or give a zero-solution. For all results presented in the main text we tested several alternative form factors as initial guess and did not find any nonzero solutions different from the ones presented. This shows that the symmetries discussed in the main text are the only possible scenarios for the gap within the associated setups.
Appendix B: Susceptibilities of FeSe/STO In Fig. 8, panels (c) and (f), of the main text we show the dynamic bare and spin susceptibilities respectively. As is apparent from this graph, dominant contributions of both quantities appear at wave vector q = M . Due to the magnitudes of χ 0,dyn q and χ S,dyn q in close vicinity of this momentum, the low frequency region is not resolved very clearly. For this purpose we re-plot both functions in Fig. 18, zoomed into the low energy regime and with modified contrast. In this way the distinct branches below the onset of the Stoner continuum are better visible. In Section III C we show effects on our results for FeSe/STO when changing either the size of the FS pockets or the distance between hole and electron bands. The former is achieved by introducing a rigid chemical potential µ, such that ξ kn → ξ kn −µ. This operation commutes with all constituents of the Hamiltonian, hence the matrix elements are not affected. On the other hand, when shifting only the hole bands by a nonrigid δµ we need to recalculate a p kn . Consider the initial kinetic termĤ 0 , which is diagonalized asĤ 0 =âξâ † , orâ †Ĥ 0â =ξ. For simplicity we omit the momentum dependence here and write the eigenvalues and eigenvectors in matrix notation. Now we add a selective shift on both sides, which only affects the hole bands: In Eq. (C1) we take the band ordering such that first all hole bands, and then all electron bands are listed. This is indicated by subscripts h and e, denoting the respective subspaces. On the right hand side we get a modified dispersionξ , which are the eigenvalues of a Hamiltonian H 0 =Ĥ 0 . To find the corresponding eigenvectors we calculate the new Hamiltonian aŝ which leads to the modified eigenvectorsâ .
To summarize, we realize a relative shift of the hole bands by diagonalizingĤ 0 , which gives the matrix elementsâ. These are used to calculateĤ 0 from Eq. (C2), which again needs to be made diagonal to have access to the desired dispersionξ and the associatedâ .
The discussions in Sec. III C are made for energies as plotted in Fig. 19(a). Starting from the initial ξ kn shown in solid black, we modify in two different ways. The influence of the FS pocket size is studied by rigidly shifting the energies as ξ kn − µ, shown in dotted red, with µ = 140 meV. To study changes with the distance between electron and hole bands we use δµ = −48 meV, shown in dashed purple. We choose these two modifications of the tight-binding model, since the effects we want to study are perfectly decoupled in this way.
Both rigid and nonrigid shifts produce changes in the bare susceptibilities due to the evaluation of Fermi-Dirac functions in Eq. (3). This in turn leads to altered boundaries in the (U, J) phase diagram, as we sketch in Fig. 19(b). As seen from the red and dotted curves, a rigid shift by µ introduces slight changes in the boundary for magnetic order, but the effect is rather minor. Contrarily, bringing electron and hole bands closer together allows for significantly larger values of U and J as can be seen from the purple dashed lines. When performing our selfconsistent calculations in (U, ω cut ) space, Fig. 12, we adjust our parameter choices according to the changed boundaries of Fig. 19(b).