Enhanced superconductivity due to forward scattering in FeSe thin films on SrTiO3 substrates

We study the consequences of an electron-phonon ($e$-$ph$) interaction that is strongly peaked in the forward scattering ($\mathbf{q} = 0$) direction in a two-dimensional superconductor using Migdal-Eliashberg theory. We find that strong forward scattering results in an enhanced $T_c$ that is linearly proportional to the strength of the dimensionless $e$-$ph$ coupling constant $\lambda_m$ in the weak coupling limit. This interaction also produces distinct replica bands in the single-particle spectral function, similar to those observed in recent angle-resolved photoemission experiments on FeSe monolayers on SrTiO$_3$ and BaTiO$_3$ substrates. By comparing our model to photoemission experiments, we infer an $e$-$ph$ coupling strength that can provide a significant portion of the observed high $T_c$ in these systems.

Determining the origin of the T c enhancement in these interface systems is critical. At the moment, proposals include charge transfer between the substrate and FeSe [2][3][4]20], electric field [6] and strain effects due to the substrate [5,8], and phononic related effects such as enhanced e-ph coupling in the FeSe layer [1,13,16] or across the interface [7,19]. Strong evidence for the latter has been provided by a recent angle-resolved photoemission spectroscopy (ARPES) study [7], which observed replica bands in the single-particle spectral function of the FeSe monolayer. These replicas are interpreted as being produced by coupling between the FeSe 3d electrons and an optical oxygen phonon branch in the STO substrate. Moreover, the replica bands are complete copies of the corresponding main bands, which implies that the responsible e-ph interaction is strongly peaked in the forward scattering direction (small momentum transfers). Such momentum dependence is notable because it can enhance superconductivity in most pairing channels [23][24][25][26][27][28][29][30][31]. As such, this cross-interface coupling provides at the same time a suitable mechanism for the T c enhancement in the FeSe/STO and FeSe/BTO systems [7,8].
In this Letter we explore this possibility and examine the consequences of strong forward scattering in the e-ph interaction for superconductivity and the spectral properties of a two-dimensional system. By solving the momentum dependent Eliashberg equations, we show that a pronounced forward scattering results in a T c that scales linearly with the dimensionless e-ph coupling constant λ m (see below) in the weak coupling limit. This is in stark contrast to the usual exponential dependence predicted by BCS theory. Furthermore, this coupling produces distinct replica structures in the spectral function similar to those observed experimentally. By comparing our model to experiments [7], we infer a significant e-ph contribution to the total T c observed in the FeSe/STO system with a modest value of λ m . Formalism -To model the FeSe monolayer we consider a single-band model for the FeSe electron pockets, which includes coupling to an oxygen phonon branch in the STO substrate. The Hamiltonian is given by where c † k,σ (c k,σ ) and b † q (b q ) are electron and phonon creation (annihilation) operators, respectively, ξ k is the band dispersion, Ω q is the phonon dispersion, and g(k, q) is the momentum-dependent e-ph coupling constant.
We calculate the single-particle self-energy due to the e-ph interaction using Migdal-Eliashberg theory. Using the Nambu notation with fermionic Matsubara frequencies ω n = (2n + 1)π/β, where β = 1/T is the inverse temperature, the self-energy isΣ(k, iω n ) = iω n [1 − Z(k, iω n )]τ 0 + χ(k, iω n )τ 3 + φ(k, iω n )τ 1 , whereτ i are the Pauli matrices, Z(k, iω n ) and χ(k, iω n ) renormalize the single-particle mass and band dispersion, respectively, and φ(k, iω n ) is the anomalous self-energy, which is zero in the normal state. In Migdal-Eliashberg theory, the self-energy is computed by self-consistently evaluating the one-loop diagram and is given bŷ In what follows we parameterize the electronic dispersion as ξ k = −2t[cos(k x a) + cos(k y a)] − µ with t = 75 meV and µ = −235 meV. This choice in parameters produces at Γ an electron-like Fermi pocket with k F = 0.97/a and a Fermi velocity v F = 0.12 eV·a/h along the k y = 0 line, where a is the in-plane lattice constant. This closely resembles the electron pocket at M point measured by ARPES experiment. Since first principles calculations indicate that the relevant oxygen phonon branch in STO is relatively dispersionless near the Γ-point [14,32,33], we approximate the phonon with a flat Einstein mode Ω q = Ω = 100 meV (h = 1), which is consistent with the observed energy separation of the replica bands [7]. Furthermore, as we are interested in the case of forward scattering, we neglect any potential fermion momentum dependence in the e-ph interaction and set g(q) = g 0 exp(−|q|/q 0 ), as microscopically derived before [7,19]. Here, q 0 sets the range of the coupling in momentum space. For different values of q 0 we adjust g 0 to obtain the desired value of the dimensionless e-ph coupling constant λ m , which is computed from the Fermi surface averaged mass enhancement in the normal state λ m = − ∂ReΣ(k,ω) ∂ω ω=0 [38]. Throughout we assume an s-wave symmetry for the gap function, consistent with the observations of a fully gapped state on the Fermi level [1,7,10,18]. Finally, we neglect the Coulomb pseudopotential µ * . One can therefore regard our T c values as upper bounds for the e-ph contribution to the FeSe/STO system.
Analytical Results -Before proceeding to full numerical solutions, we can gain some insight by first considering the case of perfect forward scattering, where the e-ph matrix element is a delta function |g(q)| 2 = g 2 0 δ q N with g 2 0 = λ m Ω 2 [34]. In the weak coupling limit, we further set Z(k, iω n ) = 1, χ(k, iω n ) = 0, and therefore φ(k, iω n ) = ∆(k, iω n ). With these approximations, the gap function on the Fermi surface is given by To determine T c we take the ansatz ∆(iω n ) = ∆ 0 /[1 + (ω n /Ω) 2 ] and follow the usual steps [36]: the gap equa- tion is linearlized by setting ∆ 2 0 = 0 for T ∼ T c and we set ω n=1 /Ω = 0. This results in the condition for T c .
The Matsubara sum can be performed exactly, yielding our final expression For the case of FeSe, T c Ω, and the hyperbolic functions dominate. To the leading order, the critical temperature is quasi-linear in the coupling strength in the weak coupling limit, T c = λm 2+3λm Ω. (A similar result was obtained in Ref. 23 in the context of the cuprates using square-well models.) For λ m = 0. 16 and Ω = 100 meV one obtains T c = 75 K, which is a remarkably high temperature for such a modest value of λ m .
The increased T c should be compared to the standard BCS value obtained for a momentum-independent coupling. In this case, the linearized gap equation simplifies to [34][35][36] where we have expanded at large Ω D /T c and ψ(z) is the digamma function [34]. This form produces the usual exponential behavior for the critical temperature, T c = 1.13Ω D exp(−1/λ m ), which predicts a T c = 2.5 K for λ m = 0. 16 and Ω D = 100 meV.  Comparing these two results, one sees that the origin of the enhanced T c lies in the momentum decoupling [24] that occurs in the Eliashberg equations when the interaction is strongly peaked at q = 0. In the BCS case, the integration over the Fermi surface is equally weighted at all momenta, leading to a n 1 |ωn| term in the BCS gap equation and subsequently a leading logarithmic behavior. In the forward scattering case, there is no integration over momentum so the ω −2 m term remains, resulting in a leading behavior that scales like 1/T c [34]. Thus, strong forward scattering serves as an ideal mechanism for producing high-T c superconductivity [30]. Furthermore, a strong forward scattering peak in the coupling constant means that this interaction will contribute in most pairing channels [7,[23][24][25][26][27][28][29][30][31]. It can therefore act in conjunction with other active unconventional channels, providing another means to increase T c further.
Numerical Results for T c and the superconducting gap -In real materials the e-ph interaction is expected to have a finite range q 0 in momentum space [7,33]. Therefore we now consider an interaction with a finite width by numerically solving the full Eliashberg equations for an e-ph coupling constant g(q) = g 0 exp(−|q|/q 0 ). Fig. 1 shows the superconducting gap at the lowest Matsubara frequency ∆(k F , iπ/β) as a function of temperature for several values of λ m and q 0 = 0.1/a. We find that the superconducting T c is already large for a modest value of λ m and increases approximately linearly with λ m in the weak coupling limit; however, the finite range of the coupling in momentum space reduces the total T c slightly with respect to the perfect forward scattering limit (see the inset of Fig. 1). The linear dependence of T c with respect to λ m may account for the wide variation of reported T c values in the literature, as differences in sample preparation or doping are likely to result in differences in the screening of the e-ph coupling and subsequently T c .
Replica Bands -The above results show that, in principle, a modest coupling to a phonon with a forward scattering peak is capable of accounting for the large T c enhancement observed in FeSe on STO and BTO. The natural question is then how much of the experimental T c is accounted for by this coupling? The observed shape and intensity of the replica bands [7,8] provide us with a direct means to estimate this by comparing our model to experiment. To do so, we calculate the single particle spectral function A(k, ω) = −ImG 11 (k, ω)/π, which requires the analytic continuation of the self-energy to the real frequency axis using the method of Ref. 39 (see also [34]). Fig. 2 plots the temperature evolution of the spectral function obtained from a full numerical solution to our model for several values of λ m , as indicated on the left, and q 0 = 0.1/a. In all cases clear replica bands are produced by the coupling, offset in energy from the main band by a fixed energy, which is Ω for small values of λ m . The separation, however, grows for increasing λ m . This is due to χ(k, ω), which shifts the main band upward in energy. This is most clearly seen in the λ m = 0.33 results, where the value of k F has visibly shrunk in the main band. In addition, for stronger values of λ m we begin to see the formation of a second replica band located at ∼ 2Ω below the main band. Thus the observation of only a single replica band in the bandstructure of FeSe/STO is consistent with a small λ m .
An intuitive picture for the intensity and energy splitting of the replica band can again be obtained in the limit of perfect forward scattering. On the real axis, the zerotemperature self-consistent equation for the self-energy in the normal state can be written as Σ(ω) = g 2 0 G(ω + Ω).  The spectral weight ratio and energy splitting between the main and replica bands can be extracted from our numerical simulations for finite values of q 0 . Fig. 3 shows A(k, ω) for k = (0, 0) as a function of λ m with q 0 = 0.1/a. The behavior matches our expectations gained from the perfect forward scattering limit: both the distance between the bands and the relative spectral weight grow with increasing λ m , though the rate of increase is slower than for the case of perfect forward scattering. ARPES experiments on the FeSe/STO system [7] observe a spectral weight ratio of ∼ 0.15 − 0.2 [34]. Comparing to our model calculations, we extract a value of λ m ∼ 0.15−0.2. This corresponds to a T c ∼ 60 − 70 K and a gap magnitude of ∆ ∼ 10 − 15 meV, which are consistent with measurements [1,7,10,18].
In Fig. 4 we present the evolution of the spectral function for increasing values of q 0 where λ m is fixed to give the same value of Z − /Z + . As expected, the replica bands are observed to smear both in energy and momentum as the value of q 0 is increased. This shows that a weakly momentum-dependent coupling (large q 0 ) to an optical mode does not reproduce the observation of a perfect replica band, with the same effective mass and termination points in the Brillouin zone. Consequently, strong forward scattering is a necessary ingredient to understand the experimental observations [7]. Summary and Conclusions -We have examined the consequences of e-ph coupling that is strongly peaked in the forward scattering direction on the spectral properties and superconducting transition of a two-dimensional electronic system. We demonstrated that such a coupling produces distinct replica bands in the electronic bandstructure consistent with recent ARPES measurements on FeSe/STO and FeSe/BTO interface systems. In order to reproduce the experimentally observed spectral function, we find that relatively modest values of the e-ph coupling are needed with λ m ∼ 0.15 − 0.2. Strong forward scattering results in a momentum decoupling of the Eliashberg equations, which subsequently produces a larger superconducting T c in comparison to the predictions of conventional BCS theory. As a result, the inferred values of λ m predict T c values on the order of 60 − 70 K due to e-ph coupling alone.
We stress that our results do not exclude the presence of another unconventional pairing channel such as spin fluctuations. The predicted values of T c and ∆ will be reduced somewhat by the inclusion of the Coulomb pseudopotential µ * . This reduction, however, can be overcome by the combination of the e-ph and unconventional interactions, since forward scattering will contribute to Cooper pairing in most channels [7]. An obvious way to distinguish between these possible scenarios is to measure the oxygen isotope effect. If a purely phononic mechanism is present then T c should have an isotope coefficient α = −∂ log(T c )/∂ log(M ) = 1/2, while the energy separation between the replica bands should decrease by ∼ 0.5(18 − 16)/16 ∼ 6% for 18 O rich substrates. Alternatively, in a multi-channel scenario, the isotope coefficient α will be reduced from 1/2 when the unconventional channel is significant in comparison to the e-ph interac-tion [37]. This provides a clear means to distinguish between these scenarios.
Finally, we note that e-ph coupling with a pronounced forward scattering peak has been studied in several contexts related to of unconventional superconductivity in the cuprates [24][25][26][27][28][29][30] and pnictides [31]. Moreover, it is also now being addressed in the context of nematic fluctuations [40,41]. This suggests forward scattering has a broader applicability in enhancing superconducting beyond interface systems.
Acknowledgments -We thank E. Dagotto, T. P. Devereaux, D.-H. Lee, R. G. Moore, D. Scalapino, and J. Zaanen for useful discussions. L. R. acknowledges funding from Rubicon (Dutch Science Foundation). T. B. is supported as a Wigner Fellow at the Oak Ridge National Laboratory. A portion of this research was conducted at the Center for Nanophase Materials Sciences, which is a DOE Office of Science User Facility. CPU time was provided in part by resources supported by the University of Tennessee and Oak Ridge National Laboratory Joint Institute for Computational Sciences (http://www.jics.utk.edu).
In this work we have solved the momentum-dependent Eliashberg equations and obtained the electron Green's function on the real axis, following the method in Ref. 1. Here, we briefly outline the main aspects of this approach.
In the main text we introduced the electron self-energy in the Eliashberg equation on the imaginary axis with the Nambu notationΣ whereτ i are the usual Pauli matrices, Z(k, iω n ) and χ(k, iω n ) renormalize the single-particle mass and band dispersion, respectively, and φ(k, iω n ) is the anomalous self-energy, which is zero in the normal state. In Migdal-Eliashberg theory, the self-energy is computed by self-consistently evaluating the one-loop diagram as shown in Fig. 1 and is given bŷ where is the dressed electron propagator, N is number of momentum grid points, and β = 1/(k B T ) is the inverse temperature. For the dispersionless optical phonon mode we study here, Ω q = Ω. For forward scattering the e-ph matrix element is independent of the fermion momentum k and is given by where f (|q|) is a function peaked at |q| = 0. For the perfect forward scattering, we take f (|q|) = N δ q where δ q is the Kronecker delta function. For a general forward scattering we take f (|q|) ∝ exp(−2|q|/q 0 ) with a proper normalization factor so that the result in the q 0 → 0 limit can be compared with the result obtained using the Kronecker delta function.
To obtain the electron self-energy for real frequencies, one needs either to perform an analytical continuation by Padé approximants, or to solve the Eliashberg equations directly on the real axis. The latter method is cumbersome while the former is unreliable; however, Ref. 1 introduced a method whereby the real-frequency self-energy can bê 1: The electron self-energy within Midgal-Eliashberg theory is given by the one-loop diagram above, where the Green's function in the diagram is not the usual bare one but the fully dressed Green's function. By solving the self-energy selfconsistently, we implicitly include higher order terms in the one-loop diagram, but the vertex corrections are neglected according to Migdal's theorem.
arXiv:1507.03967v2 [cond-mat.supr-con] 25 Aug 2015 found iteratively by using the results on the imaginary axis. Given the imaginary axis Green's functionĜ(k, iω n ), the real frequency self-energy is found by solving the iterative equation where B(z) = δ(z − Ω) − δ(z + Ω) is the spectral function of the phonon, and η > 0 is a small number. For a specific momentum dependence f (|q|) in the coupling function |g(k, q)| 2 = g 2 0 f (|q|), we set the coefficient g 0 to fix the total dimensionless coupling constant λ m , formally defined as (see, for example, Ref. 2), For the relevant low temperature range T c T Ω, the temperature dependencies of the self-energy in the normal state and λ m are negligible, so we calculate λ m at T = 100 K, if not otherwise specified. Here, . . . denotes the average over the Fermi surface, and δ(ξ q ) is approximated by a Gaussian λ m is a measure of the mass enhancement due to electron-phonon coupling, m * /m = 1 + λ m , where m * is the effective electron mass. For a normal metal (that is, with a smooth coupling function g(q) and Ω E F ), Eq. (4)- (6) gives the textbook result for the total electron-phonon coupling constant λ where N F = 1 N k δ(ξ k ) is the density of states per spin. For the forward scattering case, however, the mass enhancement λ m following Eq. (6) differs significantly from the λ defined in Eq. (8) for a normal metal. For the forward scattering with a width q 0 in momentum space, λ m ∼ q 0 λ as q 0 k F . For one set of parameters used in our numerical simulation, η = 15 meV, N = 64×64 discrete momentum points and q 0 = 0.1/a, we obtain λ m ≈ 0.24λ from numerical calculation. The Eq. (8) is therefore unsuitable to characterize the strength of electron-phonon coupling for forward scattering and to predict the corresponding T c due to electron-phonon interaction. Consequentially, we use λ m as the strength of electron-phonon coupling in the main text.
Numerically, we solve the Eliashberg equations with a momentum grid of N = 64 × 64 grid points, a cut-off energy of 600 meV for the Matsubara and real frequencies with the real frequency discretized into 1 meV intervals, and a small-scale parameter η = 15 meV. For a given set of parameters, we start the computation with initial self-energies Z = 1, χ = 0 and a small anomalous self-energy φ = 0.7 meV. The Eliashberg equations are then iterated, where upon each iteration the chemical potential is adjusted as to fix the total number of electrons. When the difference between two consecutive self-energies is less than 10 −3 meV, we have reached the convergence.

BCS GAP EQUATION
For comparison, we derive the standard BCS gap equation from the Eliashberg theory. An extended review of Eliashberg theory and the BCS limit can be found, for example, in Ref. 3. In the weak coupling limit, we can neglect the normal electron self-energy by setting Z(k, iω n ) = 1 and χ(k, iω n ) = 0. The anomalous self-energy is now equal to the gap function, φ(k, iω n ) = ∆(k, iω n ). For an isotropic s-wave gap, the gap function is momentum-independent. Under these assumptions the Eliashberg equation reduces to Unlike for the forward scattering, we now assume that |g(k, q)| 2 D (0) (q, iω n − iω m ) is momentum independent. Now the only momentum-dependent term in the Eliashberg equation is the electron dispersion ξ k+q . In this case, the sum over momenta can be replaced by an integral over the free electron density of states N ( ), Since most of the relevant physics happens close to the Fermi level, we replace the N ( ) by its approximate value N (0). The resulting integration over electron energies changes the 1/ω 2 m behavior into a 1/|ω m |, The coupling strength as a function of energy can be expressed as λ(iω ν ) = −N (0)|g| 2 D (0) (iω ν ) so that we recover the BCS gap equation, The critical temperature can be found by taking the limit ω n = 0 and assuming that λ(iω n − iω m ) and ∆(iω m ) are constants λ and ∆ 0 with a hard cut-off at the Debye frequency Ω D . The linearized gap equation with ∆ 2 0 = 0 now takes the form where we have expanded at large Ω D /T c and ψ(z) is the digamma function. The critical temperature is given by the exponential function of the total coupling constant λ This exponential suppression can be traced back to the integration over the electron energies. This, in turn, is possible because the vertex and phonon part of the gap equation are momentum independent. When the vertex has strong momentum-dependence, the integration over electron energies cannot be performed and the leading behavior is still 1/ω 2 m instead of 1/|ω m |, as will be shown when we consider perfect forward scattering.

PERFECT FORWARD SCATTERING
For the forward scattering with g(q) ∝ g 0 exp(−|q|/q 0 ), the Eliashberg equations can be solved analytically in the limit of perfect forward scattering, q 0 → 0. The results of this limit are discussed in the main text; here we provide a more detailed derivation of the proper coupling strength λ m and the corresponding critical temperature T c .
The coupling for the perfect forward scattering is on a discrete lattice with total N grid points, where N δ q is the Kronecker delta function normalized to unity in the sum 1 N q N δ q = 1. (In the large N and continuum limit, N δ q → (2π) 2 δ(q), and the Dirac delta function δ(q) can be simulated with a peak function with finite width q 0 , as for the general forward scattering case.) The one-loop electron self-energy now reads where B(z) = δ(z − Ω) − δ(z + Ω). The subscript "(1)" of the self-energy means we have performed a single iteration of the solving the Eliashberg equation, where we substitute in the bare electron Green's function instead of the dressed Green's function. At the Fermi surface, the electron energy vanishes (ξ k = 0) so the self-energy simplifies to which is momentum independent. In the low temperature limit, coth(Ω/2T ) ≈ 1, and by taking the derivative of the self-energy at ω = 0 we obtain the total coupling constant as defined in Eq. (6) The numerical value of λ m from the self-consistently calculated self-energy is very close to the value given by the above equation. For a nearly perfect forward scattering with a peak width q 0 = 0.1/a, even after considering the small temperature dependence, the numerical value differs from the above result by a few percent. We now continue to derive the linear dependence of transition temperature T c in terms of the coupling strength λ m in the weak coupling limit for the perfect forward scattering. In this limit, we can neglect the normal electron self-energies, that is Z(k, iω n ) = 1 and χ(k, iω n ) = 0. We assume an isotropic s-wave gap. After summing over the momentum q, we obtain the self-consistent equation for the gap function φ(k, iω n ) = ∆(k, iω n ), To find the critical temperature, we take an ansatz for the gap function with the form ∆(iω n ) = ∆ 0 1 + (ω n /Ω) 2 .
Our numerical results verify that this shape properly represents the Matsubara frequency dependence of the gap.
Setting ω n=1 /Ω = 0 and linearizing the gap equation with ∆ 2 0 = 0, we find the equation for the critical temperature, where the sum over Matsubara frequencies can be performed exactly by the method of contour integration. We finally obtain 1 = λ m 2T c 2Ω + Ω cosh Ω Tc − 3T c sinh Ω Tc 1 + cosh Ω Tc .
It is reasonable to assume that the critical temperature T c is much less than the optical phonon frequency, as Ω is typically of the order of 100 meV ≈ 1000 K. In this limit, the hyperbolic functions dominate and cancel each other in the numerator and denominator of the above equation. We find a critical temperature of This gives a remarkably high critical temperature, even for a modest coupling constant λ m because of the quasi-linear dependence. For example, λ m = 0.16 with Ω = 100 meV yields a critical temperature of T c = 75 K.

SPECTRAL WEIGHT RATIO AND λm FROM ARPES DATA
By fitting ARPES data from Lee et al. [4], we found that the ratio of spectral weights at the peak of the main band and the mirror band is approximately Z − /Z + ∼ 0.15-0.2, as shown in Fig. 2