Multichannel superconductivity of monolayer FeSe on SrTiO$_3$: Interplay of spin fluctuations and electron-phonon interaction

We investigate the effects of electron-phonon coupling, as well as of spin and charge fluctuations on the superconducting state in a single layer of FeSe on SrTiO$_3$ substrate. These three bosonic mediators of Cooper pairing are treated on equal footing in a multichannel, full-bandwidth, multiband, and anisotropic Eliashberg theory of the interacting state. Our self-consistent calculations show that an s-wave symmetry of the superconducting gap is compatible only with a complete absence of spin fluctuations. When spin fluctuations are present, the sign-changing nodeless d-wave pairing symmetry is always obtained, yet the essential ingredient for explaining the gap magnitude and critical temperature is still the interfacial electron-phonon interaction.

Ever since the discovery of superconductivity in monolayer FeSe on SrTiO 3 substrate (FeSe/STO) [1] no consensus has yet been reached about such fundamental questions as, for example, the Cooper pairing mechanism or the superconducting gap symmetry, despite tremendous research effort in both theory and experiment [2]. The measured superconducting critical temperature T c in FeSe/STO ranges from 50 K up to over 100 K [1,[3][4][5][6][7][8], which is an astonishing increase from the comparatively low T c (∼ 8 K) of the parent compound FeSe [9]. The bulk FeSe material has been shown to be nonmagnetic, but is poised in close vicinity to a magnetic phase transition, which is manifested in strong spin fluctuations (SFs) [10][11][12]. It is therefore commonly believed that superconductivity in bulk FeSe, as well as in other Fe-based superconductors with similar Fermi surface (FS) properties, has a magnetic origin [10][11][12][13][14].
In FeSe/STO the situation is markedly different because FS nesting conditions are changed due to doping with STO electrons at the interface. The resulting FS consists only of electron-like pockets, which makes the material inexplicable using 'standard' FS nesting arguments. There have, however, been attempts to explain superconductivity in FeSe/STO by theories for SFs that are more specialized to this particular system [15,16]. Recently, experimental [17,18] and theoretical [19] evidence for magnetic signatures have been provided, but the literature is sparse on the verification of SFs mediated superconductivity in FeSe/STO. The authors of the current work have argued in Ref. [20] that SFs possibly contribute to the pairing strength in the superconducting state, but are not the dominant 'pairing glue'.
To complicate the issue further, a sizable electronphonon interaction (EPI) between the substrate phonon and FeSe electrons has been detected in experiment [7,21], an observation reproduced by Density Functional Theory calculations [22][23][24][25]. A phonon branch of relatively large frequency Ω = 81 meV gives rise to strongly * fabian.schrodi@physics.uu.se † alex.aperis@physics.uu.se ‡ peter.oppeneer@physics.uu.se enhanced forward scattering exhibiting small momentum transfer. This characteristic feature of FeSe/STO was shown to play an important role in the superconducting state [26][27][28][29]. However, it is currently debated whether EPI and SF effects are competing or cooperative [30], and how their interplay reflects in experimentally observable quantities of the superconducting state. This is partially due to uncertainty concerning the symmetry of the superconducting order parameter. Although it has been shown that the gap function does not exhibit nodes in the Brillouin zone (BZ), an observation broadly agreed on, it is still debated whether a sign change occurs between FS pockets separated by a wave vector q = (π, π) in the unfolded BZ. In several works it has been argued that the order parameter has s-wave symmetry (no sign change), a conclusion drawn from results of impurity measurements [31] and Angular Resolved Photoemission Spectroscopy [3,7]. On the contrary, Scanning Tunneling Spectroscopy (STS) carried out in Ref. [32] revealed that a sign change of the order parameter (d-wave) is similarly possible as s-wave, while other recent STS investigations argued more determinedly in favor of a sign-changing order parameter [33][34][35], which suggested SF-mediated pairing to be most relevant [34].
In this Letter, we study the competition/cooperation between SFs and EPI in the superconducting state of FeSe/STO. This is done within a self-consistent multichannel Eliashberg formalism in which EPI, SFs and charge fluctuations (CFs) are treated on equal footing. Our results reveal that the most plausible BZ symmetry of the superconducting gap is nodeless d-wave, while an anisotropic s-wave state is also possible but only under the exclusion of any influence of spin fluctuations. Hence, our investigation shows that the EPI is responsible for the gap magnitude and high T c , but the SFs generate the unconventional pairing symmetry.
We start from the multiorbital Hubbard-Fröhlich HamiltonianĤ =Ĥ 0 +Ĥ int +Ĥ ph +Ĥ eph , expressed in an orbital basis of the five Fe-d states, witĥ HereĤ 0 andĤ ph are the electron and phonon kinetic energies, respectively withĉ † k,p,σ (b † q ) electron (phonon) creation operators, σ denotes spin, p, q, s, t orbital indices and k, k ′ , q are BZ wave vectors (we set q = k − k ′ ). We consider the 1-Fe unit cell and thus work in the unfolded BZ. Further, ξ k,p,q denotes electron energies in orbital space and Ω is a characteristic Einstein phonon frequency. The EPI is given byĤ eph with EPI scattering matrix elements g q,p,q . Electron correlations are described by the purely electronic termĤ int which carries the information about CFs and SFs mediated interactions. As usual,ˆ S i,s (n i,s ) are spin (density) operators, U , V ′ are the respective intra-and inter-orbital Hubbard interactions, J is the Hund's rule coupling, and J ′ the pair-hopping interaction with V ′ = U − 3J/4 − J ′ and J ′ = J/2 [36,37].
Our tight-binding description of FeSe/STO is adopted from Refs. [28,38], where hopping energies for bulk FeSe [39] are modified so as to account for the lattice distortion that arises when a monolayer of FeSe is deposited on the substrate. Diagonalization ofĤ 0 yields electron energies ξ k,n with band index n and matrix elements a p k,n which we utilize to derive our Eliashberg theory in band space.
Including the infinite series of Feynman diagrams for all first-order scattering processes due to EPI, SFs, and CFs, we arrive at the electron self-energy [40], where m, m ′ index Matsubara frequencies (ω m = πT (2m+1)) andρ 0(3) are Pauli matrices. The EPI kernel is derived similarly as in Refs. [20,28], where the electron-phonon scattering at the interface is modeled by the functional form g q = g 0 exp − |q|/q c with interaction strength g 0 , q c = 0.3a −1 and a the lattice constant [7]. Further, we employ an Einstein phonon with frequency Ω = 81 meV to which the FeSe electrons are coupled [7,23]. For brevity we use henceforth the notation V (eph) q,l,n,n ′ with l = m − m ′ . As is described in detail in the Supplementary Material (SM) and Ref. [20], we keep the full orbital content encoded in a p k,n when calculating band-dependent interaction kernels V q,l,n,n ′ for SFs (S) and CFs (C). Labels (+) and (−) are respectively referring to kernels as they are used in electron-energy renormalization and superconducting equations, see below.
From the Green's functionĜ −1 k,m,n = iω m Z k,m,nρ0 − ξ k,n + Γ k,m,n ρ 3 − φ k,m,nρ1 , we derive a self-consistent set of a total of 15 coupled Eliashberg equations for the mass Z k,n,m and chemical potential Γ k,n,m renormalization functions, and the superconducting order parameters φ k,n,m [41,42]. The full interaction kernels for EPI, CFs, and SFs are given by K q,l,n,n ′ , using K (+) q,l,n,n ′ in the equations for Z k,n,m and Γ k,n,m , and K (−) q,l,n,n ′ in the equation for φ k,n,m . For further details we refer to the SM and Refs. [20,28].
In the theory employed here we treat the scattering strength g 0 as a parameter to control the strength of the EPI. For CFs and SFs we are free to choose the intraorbital onsite interaction U and the Hund's rule coupling J. For convenience we set J to a fixed ratio of U , leaving us with two variational quantities g 0 and U . Hence, we have direct control over the interaction strength for all three mediators of superconductivity. Our calculations are carried out using the Uppsala Superconductivity (UppSC) code [43][44][45][46][47], in particular, by combining the advances of Refs. [28] and [20].
Motivated by our earlier work [20], we choose the temperature T = 5 K (to compare to available experiments [3,7]) and select J = U/2. For calculations of the SF and CF kernels we apply a high-energy cutoff ω cut = 0.54 eV (see [20]). We consider this pair of (J, ω cut ) as it reportedly allows for a finite superconducting gap, even in the absence of any EPI [20]. Left with two parameters U and g 0 , which respectively control the strength of CFs/SFs and EPI, we solve the full bandwidth, multiband, and anisotropic Eliashberg equations (see SM) for each pair (U, g 0 ). Computing the experimentally observable gap function as ∆ k,m,n = φ k,m,n /Z k,m,n , we plot the maximum value ∆ = max k,n |∆ k,m=0,n | in Fig. 1. The red dashed line through (U, g 0 )-space represents a gap size of ∆ ≃ 12 meV as was measured for this material [3,48].
It is apparent that the onset of superconductivity with respect to the electron-phonon scattering strength lies at g 0 ∼ 0.6 − 0.7 eV. For growing g 0 the maximum superconducting gap increases approximately linear for fixed U . On the other hand, for small g 0 the choice of U must be close to the maximally allowed value, which in term is dictated by the Stoner criterion, to obtain a finite ∆. Here we do not pay special attention to either of the limiting cases U = 0 eV or g 0 = 0 eV, because these have been analyzed in detail in previous studies [20,28,29]. When computing the results of Fig. 1 we assumed that the symmetry of the order parameter is a priori not known, which is why we performed each calculation with an initial s-wave and d-wave state. For all parameter space shown in Fig. 1 we find a d-wave symmetry as the converged solution (gap shape shown in the inset of Fig. 1), with exception of the purely phononic case (U = 0), where the symmetry is s-wave [28].
To better understand these findings let us take a closer look at the couplings in the superconducting channel, where we focus on the FS for simplicity. In Fig. 2(a) we show the FS sheets of our tight-binding model (black lines) and schematically draw all interactions included in our theory. Here we use the definitions V eph = V Choosing U = 1.07 eV as an example, we show in Fig. 2(c) and (d) the SF and CF kernel, respectively. We observe that a close-to-nesting condition between the two FS pockets leads to a leading contribution at q = (π, π) (= M ) for the spin part. Further, we find a small-q coupling for V S which has significantly lower magnitude. As indicated in panel (a), since the SF kernel enters repulsively in the equation for the superconducting order parameter (see SM), a sign change in the superconducting gap is promoted both globally between the FS sheets, and locally on each pocket. However, the repulsive small-q contribution of V S is too weak to induce such a local sign change, which has not been found here nor in Ref. [20]. As concerns CFs, we find a weak attractive interaction peaked at q = Γ, see Fig. 2(d), which has a similar order of magnitude as the repulsive SFs coupling at this wave vector.
To first order approximation we can assume that SF and CF kernels do not significantly contribute at q = Γ, due to their comparable magnitude and the fact that V S and V C are competing in the superconducting channel. Therefore, we are left with V eph peaked at Γ, and V S having a leading contribution at M . In other words, we have a locally sign-conserving EPI, and a repulsive interaction V S that induces a sign change between FS pockets. These contributions cooperatively support a signchanging d-wave solution. Note that for a global s-wave state we face a different situation. If the order parameter does not change sign between the two FS pockets, contributions V eph and V S are competing and hence reducing the size of the superconducting gap. Then, s-wave symmetry becomes energetically less favorable.
Next, we examine the temperature dependence of the superconducting gap for various pairs of (U, g 0 ). The couplings are taken from the red dashed line in Fig. 1 such that the magnitude of ∆ ≃ 12 meV corresponds to the experimental value at T ∼ 5 K. In Fig. 3(a) we show the result for the maximum superconducting gap as function of T , self-consistently computed for choices of U as written in the legend. Note that all four curves exhibit the behavior lim T →0 ∆(T ) ≃ 12 meV, while the values for T c change with U . Notably, our results for U 0.6 eV, indicating a weak coupling to SFs, fall on top of each other to a good approximation. As U increases towards the maximum value U max , which is dictated by the Stoner criterion [20], T c gradually decreases.
To elucidate this aspect further we draw the ratio ∆/k B T c as blue circles in Fig. 3(b). It is apparent that an enhancement of the SF kernel leads to a more strongly coupled system. We analyze this behavior more detailedly by fitting our results to the functional form where c 1 , c 2 , and c 3 are free parameters of the fit. The result, drawn as solid blue line in Fig. 3(b), matches our data points very accurately. Taking the limit U → 0 in Eq. (6) corresponds to the case of only considering EPI. Therefore, the value c 1 = 2.164 equals ∆/k B T c without any influence of CFs or SFs, compare Ref. [28]. Further, it was shown in Ref. [20] that the leading spin contributions to the kernel at q = M scale like U 2 /(U max − U ), which makes us associate c 3 ≡ U max = 1.183 eV. This value is remarkably close to the precise value found from the Stoner criterion [20]. With scaling constant c 2 = 0.0118 eV −1 we therefore find that enhanced contributions due to SFs drive ∆/k B T c more towards the strong coupling regime.
The red dashed lines in Fig. 3(b) serve as approximate bounding box for reasonable magnitudes of ∆/k B T c , according to existing works on FeSe/STO where the gap magnitude was found as ∼ 12 meV [3,27,48]. For our fitting curve (blue) to stay within the red shaded area, U can be chosen relatively large but not in too close vicinity of the maximally allowed value U max . This has the consequence that the influence of SFs on the superconducting gap magnitude is bounded.
Further analysis shows that the decrease of T c with growing U in Fig. 3(a) stems from a competition of SFs and EPI in mediating superconductivity. To prove this we first need to calculate the renormalized Fermi surface, defined by the condition ξ k,n + Γ k,m=0,n = 0, in the interacting state for a given T . No significant changes in comparison to the non-interacting FS are detected in the whole temperature range considered here, which goes in line with earlier predictions of a temperatureindependent FS in this system [29,49]. As the influence of SFs increases with U , we detect decreasing values of φ k,m=0,n∈FS kF , while λ (T <Tc) m = Z k,m=0,n∈FS kF − 1 grows. Therefore we observe a decrease in the superconducting gap magnitude ∆ ∼ φ/Z.
We have already seen that SFs and EPI can act cooperatively for superconductivity, a statement that holds true when considering the leading contributions at M and Γ, respectively. However, an increase in U enhances also the interaction kernel V S at Γ, which competes with the small-q EPI. Therefore, the influence of V eph gets partially suppressed, which is not compensated by an increased V S coupling at M . This is a competing aspect of these bosonic mediators of Cooper pairing in FeSe/STO. Our computed small-coupling value λ m ≃ 0.35 is furthermore consistent with experimental observations [7,50]. Importantly, our multichannel calculations show that in the presence of a sizable EPI, even modest SFs (i.e., small U ) nonetheless lead to a nodeless d-wave pairing symmetry. This unconventional symmetry is consistent with the sign-changing order parameter that was deduced in recent STS experiments [33][34][35].
In summary, we have investigated the superconducting state of FeSe/STO treating EPI, SFs and CFs on equal footing. Our self-consistent multichannel Eliashberg theory shows unambiguously that an s-wave symmetry of the order parameter is possible only with negligible magnetic contributions, i.e., in the limit U → 0. Conversely, a nodeless d-wave state is realized for any small influence of SFs. For the latter scenario the obtained maximum superconducting gap and T c are compatible with experiment, provided that U is reasonably smaller than its maximally possible value (given by the Stoner criterion). Consequently, we are led to identify the main 'pairing glue' in FeSe/STO as the EPI, but any small contribution from SFs is sufficient to cause an unconventional d-wave pairing symmetry.
This work has been supported by the Swedish Research Council (VR), the Röntgen-Ångström Cluster, the K.

SUPPLEMENTARY MATERIAL
The Hamiltonian of the system is given bŷ with kinetic termĤ 0 = k,p,q,σ ξ k,p,qĉ † k,p,σĉk,q,σ . Here we useĉ † k,p,σ andĉ k,p,σ as electron creation and annihilation operators, where σ denotes spin, k a wave vector and p, q orbital indices. Further, ξ k,p,q denotes electron energies in orbital space. The lattice vibrations are represented byĤ ph = Ω q b † qbq + 1 2 , where Ω = 81 meV is the characteristic Einstein phonon frequency of the interfacial phonon andb † q (b q ) creates (annihilates) a phonon with exchange momentum q. The electron-phonon interaction (EPI) is given byĤ eph which describes an electronĉ k,q,σ being scattered into the statê c † k ′ ,p,σ via scattering matrix elements g q,p,q and phonon displacementsû q =b † q +b −q . In this work we set q = k − k ′ .
Electron-electron interactions in the system are modeled viâ In Eq. (2) we use index i to describe the lattice site, hencê c † i,s,σ create an electron of orbital character s with spin σ at site i. We describe spin operators byˆ S i,s as defined e.g. in Ref. [37]. The occupation numbers are given byn i,s,σ =ĉ † i,s,σĉi,s,σ andn i,s = σĉ † i,s,σĉi,s,σ . The first and third term scale with intraorbital onsite coupling U and Hund's rule coupling J, respectively. For the interorbital onsite energy V ′ and the pair hopping amplitude J ′ we make the choice V ′ = U − 3J/4 − J ′ and J ′ = J/2 [36,37].
A diagonalization of the hopping energies ξ k,p,q provides us with the electronic dispersion in band space ξ k,n and matrix elements a p k,n . These are used to calculate the bare susceptibilities Im χ 0 q (ω) pq st = −π n,n ′ ,k ′ a s k,n a p, * k,n a q k+q,n ′ a t, * k+q,n ′ × n F (ξ k,n ) − n F (ξ k+q,n ′ ) δ ξ k+q,n ′ − ξ k,n + ω , as function of frequency ω, where n F (.) is the Fermi-Dirac function. Within the Random Phase Approximation (RPA) the spin and charge susceptibilities are respectively calculated as where we make use of the Stoner tensors The RPA susceptibilities can be used to find the available phase space for choosing U and J, see Ref. [20]. From here we calculate the real-frequency dependent interaction kernels for spin and charge fluctuations via The outcome of Eqs. (7) and (8) is used to compute the kernels in band space (q = k − k ′ ): (ω) n,n ′ = k,s,t,p,q a t * k,n a s * k,n V (±) q (ω) pq st a p k ′ ,n ′ a q k ′ ,n ′ .
Hereafter we obtain the Matsubara frequency-dependent interactions as q,l,n,n ′ = 1 π P ωcut −ωcut dω ω − iq m Im V (±) q (ω) n,n ′ (10) with high-energy cutoff ω cut which removes the highenergy parts of the excitation spectrum, especially the incoherent part that is irrelevant to superconductivity [20].
As stated before, the electron-phonon part of the interaction is given by an Einstein phonon spectrum, so we can write the kernel as V (eph) q,l,n,n ′ = |g q,n,n ′ | 2 2Ω Ω 2 + q 2 l , with band-dependent scattering elements g q,n,n ′ . In this work we make the approximation g q,n,n ′ = g q , with g q as given in the main text [28]. The self-consistent Eliashberg equations that we solve in this work are derived in the standard way. Starting from the non-interacting electron Green's function in Nambu space Ĝ 0 k,m,n −1 = iω mρ0 − ξ k,nρ3 , we know that bosonic interactions lead to a Green's func-tionĜ k,m,n , which obeys a Dyson equationĜ k,m,n = G 0 k,m,n +Ĝ 0 k,m,nΣk,m,nĜk,m,n . Hereρ i are Pauli matrices andΣ k,m,n represents the electron self-energy. By defininĝ G −1 k,m,n = iω m Z k,m,nρ0 − ξ k,n + Γ k,m,n ρ 3 − φ k,m,nρ1 the self-energy takes the form Σ k,m,n = iω m 1 − Z k,m,n ρ 0 + Γ k,m,nρ3 + φ k,m,nρ1 .