Electronic Origin of T_c in Bulk and Monolayer FeSe

Abstract

FeSe is classed as a Hund’s metal, with a multiplicity of d bands near the Fermi level. Correlations in Hund’s metals mostly originate from the exchange parameter J, which can drive a strong orbital selectivity in the correlations. The Fe-chalcogens are the most strongly correlated of the Fe-based superconductors, with $d_{xy}$ the most correlated orbital. Yet little is understood whether and how such correlations directly affect the superconducting instability in Hund’s systems. By applying a recently developed ab initio theory, we show explicitly the connections between correlations in $d_{xy}$ and the superconducting critical temperature $T_c$. Starting from the ab initio results as a reference, we consider various kinds of excursions in parameter space around the reference to determine what controls $T_c$. We show small excursions in J can cause colossal changes in $T_c$. Additionally we consider changes in hopping by varying the Fe-Se bond length in bulk, in the free standing monolayer M-FeSe, and M-FeSe on a SrTiO${}_3$ substrate (M-FeSe/STO). The twin conditions of proximity of the $d_{xy}$ state to the Fermi energy, and the strength of J emerge as the primary criteria for incoherent spectral response and enhanced single- and two-particle scattering that in turn controls $T_c$. Using c-RPA, we show further that FeSe in monolayer form (M-FeSe) provides a natural mechanism to enhance J. We explain why M-FeSe/STO has a high $T_c$, whereas M-FeSe in isolation should not. Our study opens a paradigm for a unified understanding what controls $T_c$ in bulk, layers, and interfaces of Hund’s metals by hole pocket and electron screening cloud engineering.

1. Introduction

Iron-pnictogen and iron-chalcogen based superconductors (IBS) are classed as Hund’s metals, meaning correlations mostly originate from the Hund’s exchange parameter J. In recent years a consensus has evolved that strong Hund’s correlations drive the ubiquitous bad metallicity observed in their phase diagrams [1,2]. Such metals have a multiplicity of bands near the Fermi level $E_F$; in particular FeSe [3,4,5] has all five Fe d states active there. Correlations are observed to be highly orbital-selective [5,6] (a signature of “Hundness”) with $d_{xy}$ the most strongly correlated orbital. In recent years, the role of Hund’s exchange J is explored in determining the orbital-selectivity [7] and triplet pairing in Uranium based superconductors [8]. However, very little is known whether Hund’s correlation can generate glue for superconducting pairing and control $T_c$ in singlet-pairing channel.

$T_c$ is a mere 9 K in bulk, but it has been observed to increase to ∼75 K when grown as a monolayer on SrTiO${}_3$[9] (M-FeSe/STO), and 109 K on doped SrTiO${}_3$ [10]. Thus while “Hundness” has been found to be important in controlling the single- and two-particle spectral properties of bulk FeSe, the multiplicity of factors (orbital character, spin-orbit coupling, shape of Fermi surface and dispersion of states around it, differences in susceptibilities, nematicity, electron-phonon interaction) obfuscate to what extent Hundness, or other factors, drive superconductivity, and whether ‘Hundness’ can at all explain the jump in $T_c$ going from bulk to M-FeSe/STO.

In this work we calculate the superconducting instability using a new high fidelity, ab initio approach [11,12]. For the one-particle Green’s function it combines the quasiparticle self consistent GW (QSGW) approximation [13] with CTQMC solver [14,15] based dynamical mean field theory (DMFT) [16]. This framework [17,18] is extended by computing the local vertex from the two-particle Green’s function by DMFT [19,20], which is combined with nonlocal bubble diagrams to construct a Bethe-Salpeter equation [21,22,23]. The latter is solved to yield the essential two-particle spin and charge susceptibilities ${\chi }^d$ and ${\chi }^m$ — physical observables which provide an important benchmark. Moreover they supply ingredients needed for the particle–particle Bethe–Salpeter equations, which yields eigenvalues and eigenfunctions that describe instabilities to superconductivity. We will denote QSGW⁺⁺ as a shorthand for the four-tier QSGW+DMFT+BSE.

The numerical implementation is discussed in Pashov et al. [12].

QSGW⁺⁺ has high fidelity because QSGW captures non-local dynamic correlation particularly well in the charge channel [12,24], but cannot adequately capture effects of spin fluctuations. DMFT does an excellent job at the latter, which are strong but mostly controlled by a local effective interaction given by Hubbard parameters U and J. These are calculated within the c-RPA [25] from the QSGW Hamiltonian using an approach [12] similar to that of Ref. [26]. That it can well describe superconductivity in a parameter-free manner has now been established in several Hund’s materials [21,22]. In FeSe, we have also shown [27] that it reproduces the main features of neutron structure factor [28,29].

Sr${}_2$RuO${}_4$ was the first system where we implemented our QSGW⁺⁺ ability. We found a dense spectrum of triplet and singlet superconducting instabilities therein [21]. Lacking at that time the ability to calculate U, J, we borrowed values from an earlier LDA+DMFT study [30]. Later we were able to build our own c-RPA calculated from QSGW (QSGW+cRPA), and found U, J to be approximately 2/3 of the original values. Our initial study concluded that the triplet instability was slightly stronger than the B${}_{1g}$ − d${}_{x^2-y^2}$ singlet in unstrained Sr${}_2$RuO${}_4$. Redoing the procedure with U, J calculated from QSGW+cRPA [31] the conclusions did not qualitatively change, except the singlet eigenvalue became larger than the triplet. Thus the original study [21] did not predict the correct ground state. From a prior study we discovered that the superconducting instability can be sensitive to the exact choices of U, J. Here we show small changes in J have a dramatic effect on spin susceptibility and superconductivity in FeSe. For that reason we recalculate QSGW+cRPA for for all cases we report here, on bulk and layered variants of FeSe. This is a step forward in making the calculations as ab-initio as possible, and as we discuss later in the paper, this provides us with a remarkable unified understanding of both bulk and layered variants of FeSe, the correlations and superconducting instabilities therein.

To isolate the effect of Hundness we make excursions about the ab initio reference point, by treating J as a free parameter. One of our primary conclusions is that intense and broad low energy spin fluctuations in the vicinity of an antiferromagnetic ordering vector is the primary glue for pairing and the controlling element for $T_c$ in several variants of FeSe (both bulk and layered); and that this in turn is controlled by J. We further show, using c-RPA, that J can be tuned by varying the screening through, e.g., changes in the geometry such as the change from bulk to monolayer. We also consider excursions in Fe-Se bond length, and can find correlation can be enhanced or suppressed by changes in it. For Hund’s coupling to be effective in driving $T_c$ it needs a certain ‘universal’ band feature: the Fe $d_{xy}$ state must be in close proximity to the Fermi energy.

The paper organised as follows. We:

• Present the mathematical formulation and implementation of different susceptibilities in spin, charge and superconducting channels;
• List the structural and Hubbard parameters used for different materials, and present ab initio results bulk FeSe;
• By treating J and the Se height as free parameters, show how they affect correlations and superconductivity in bulk FeSe;
• Show results for a standing monolayer of FeSe (M-FeSe);
• Show results for a monolayer of FeSe on SrTiO${}_3$ (M-FeSe/STO);
• Interpret these results to explain what controls T${}_c$.

2. Methodology

The QSGW+DMFT implementation is discussed in our method paper [12]. The one-body part of QSGW is performed on a $12\times 12\times 12$ k-mesh and charge has been converged up to ${10}^{-6}$ accuracy, while the (relatively smooth) many-body static self-energy ${\Sigma }^0\left(\mathbf{k}\right)$ is constructed on a $8\times 8\times 8$ k-mesh from the dynamical $GW$$\Sigma (\mathbf{k},\omega )$. ${\Sigma }^0\left(\mathbf{k}\right)$ is iterated until convergence (RMS change in ${\Sigma }^0<{10}^{-5}$ Ry). DMFT is solved for all five Fe-$3d$ orbitals using a Continuous time Quantum Monte Carlo technique (CTQMC) [14] on a rotationally invariant Coulomb interaction. From CTQMC-DMFT we sample the local full and connected two-particle Green’s functions. The DMFT for the dynamical self energy is iterated, and converges in ∼20 iterations. Calculations for the single particle response functions are performed with ${10}^9$ QMC steps per core and the statistics is averaged over 64 cores. The two particle Green’s functions are sampled over a larger number of cores (30,000-50,000) to improve the statistical error bars. We compute the local polarisation bubble from the local single-particle Green’s function. In order to extract ${\Gamma }_{loc}^{irr}$, we employ the local Betlpeter equation which relates the local two-particle Green’s function (${\chi }_{loc}$) sampled by CTQMC, with both the local polarisation function (${\chi }_{loc}^0$) and ${\Gamma }_{loc}^{irr}$.

$${\Gamma }_{loc\genfrac{}{}{0pt}{}{{\alpha }_1,{\alpha }_2}{{\alpha }_3,{\alpha }_4}}^{irr,m\left(d\right)}{(i\nu ,i{\nu }^{\prime })}_{i\omega }={[{\left({\chi }_{loc}^0\right)}_{i\omega }^{-1}-{\chi }_{loc}^{m\left(d\right)-1}]}_{\genfrac{}{}{0pt}{}{{\alpha }_1,{\alpha }_2}{{\alpha }_3,{\alpha }_4}}{(i\nu ,i{\nu }^{\prime })}_{i\omega }.$$

(1)

$\Gamma$ is the local irreducible two-particle vertex functions computed in magnetic (m) and density (d) channels. $\Gamma$ is a function of two fermionic frequencies $\nu$ and ${\nu }^{\prime }$ and the bosonic frequency $\omega$.

The non-local polarisation bubble in the p-h channel is computed from single-particle DMFT Green’s functions embedded into the QSGW bath.

$${\chi }_{\genfrac{}{}{0pt}{}{{\alpha }_1{\sigma }_1,{\alpha }_2{\sigma }_2}{{\alpha }_3{\sigma }_3,{\alpha }_4{\sigma }_4}}^0{(i\nu ,i{\nu }^{\prime })}_{\mathbf{q},i\omega }=-\frac{T}{N_k}\sum _{\mathbf{k}}{G}_{{\alpha }_2{\alpha }_1,{\sigma }_1}(\mathbf{k},i\nu )·G_{{\alpha }_3{\alpha }_4,{\sigma }_3}(\mathbf{k}+\mathbf{q},i\nu +i\omega )·{\delta }_{i\nu ,i{\nu }^{\prime }}·{\delta }_{{\sigma }_1{\sigma }_2}·{\delta }_{{\sigma }_3{\sigma }_4}$$

(2)

Spin (${\chi }^m$) and charge (${\chi }^d$) susceptibilities are computed from momentum dependent Bethe-Salpeter in magnetic (spin) and density (charge) channels.

$${\chi }_{\genfrac{}{}{0pt}{}{{\alpha }_1,{\alpha }_2}{{\alpha }_3,{\alpha }_4}}^{m\left(d\right)}{(i\nu ,i{\nu }^{\prime })}_{\mathbf{q},i\omega }={[{\left({\chi }^0\right)}_{\mathbf{q},i\omega }^{-1}-{\Gamma }_{loc}^{irr,m\left(d\right)}]}_{\genfrac{}{}{0pt}{}{{\alpha }_1,{\alpha }_2}{{\alpha }_3,{\alpha }_4}}^{-1}{(i\nu ,i{\nu }^{\prime })}_{\mathbf{q},i\omega }.$$

(3)

The susceptibilities ${\chi }^{m\left(d\right)}(\mathbf{q},i\omega )$ are computed by closing ${\chi }_{\genfrac{}{}{0pt}{}{{\alpha }_1,{\alpha }_2}{{\alpha }_1,{\alpha }_2}}^{m\left(d\right)}{(i\nu ,i{\nu }^{\prime })}_{\mathbf{q},i\omega }$ with spin or charge bare vertex $\gamma$ ($\gamma$=1/2 for spin and $\gamma$=1 for charge) and summing over frequencies ($i\nu$,$i{\nu }^{\prime }$) and orbitals (${\alpha }_{1,2}$).

$${\chi }^{m\left(d\right)}(\mathbf{q},i\omega )=2{\gamma }^2\sum _{i\nu ,i{\nu }^{\prime }}\sum _{{\alpha }_1{\alpha }_2}{\chi }_{\genfrac{}{}{0pt}{}{{\alpha }_1,{\alpha }_2}{{\alpha }_1,{\alpha }_2}}^{m\left(d\right)}{(i\nu ,i{\nu }^{\prime })}_{\mathbf{q},i\omega }.$$

(4)

The superconducting pairing susceptibility ${\chi }^{p-p}$ is computed by dressing the non-local pairing polarisation bubble ${\chi }^{0,p-p}(\mathbf{k},i\nu )$ with the pairing vertex ${\Gamma }^{irr,p-p}$ using the Bethe-Salpeter equation in the particle–particle channel.

$$\begin{array}{c}\hfill {\chi }^{p-p}={\chi }^{0,p-p}·{[\mathbf{1}+{\Gamma }^{irr,p-p}·{\chi }^{0,p-p}]}^{-1}\end{array}$$

(5)

The particle–particle vertex in the singlet channel has odd-symmetry under exchange of two external spins (and even symmetry in the triplet channel),

$$\begin{array}{c}\hfill {\Gamma }^{p-p,s}=\frac{1}{2}[{\Gamma }_{\genfrac{}{}{0pt}{}{\uparrow \downarrow }{\downarrow \uparrow }}^{p-p}-{\Gamma }_{\genfrac{}{}{0pt}{}{\downarrow \uparrow }{\downarrow \uparrow }}^{p-p}]\end{array}$$

(6)

$$\begin{array}{c}\hfill {\Gamma }^{p-p,t}=\frac{1}{2}[{\Gamma }_{\genfrac{}{}{0pt}{}{\uparrow \downarrow }{\downarrow \uparrow }}^{p-p}+{\Gamma }_{\genfrac{}{}{0pt}{}{\downarrow \uparrow }{\downarrow \uparrow }}^{p-p}]\end{array}$$

(7)

The irreducible particle–particle vertex function channel ${\Gamma }^{p-p,irr}$ which provides the pairing glue to form Cooper pairs, consists of the fully irreducible vertex function ${\Gamma }^{f,irr}$ and the reducible vertex functions computed in the particle-hole channels

$$\begin{array}{c}\hfill {\tilde{\Gamma }}^{p-h}={\Gamma }^{full,p-h}-{\Gamma }^{irr,p-h}\end{array}$$

(8)

The reducible magnetic/charge vertex ${\tilde{\Gamma }}^{irr,p-h}$ is obtained from the non-local magnetic/charge susceptibilities and magnetic/charge irreducible vertex functions by

$$\begin{array}{cc}\hfill & {\tilde{\Gamma }}_{\genfrac{}{}{0pt}{}{{\alpha }_2,{\alpha }_4}{{\alpha }_1,{\alpha }_3}}^{p-h,m/d}{(i{\nu }^{\prime },i\nu )}_{\mathbf{q},i\omega }={\sum }_{i\nu 1,i\nu 2}{\sum }_{\genfrac{}{}{0pt}{}{{\alpha }_2^{\prime },{\alpha }_4^{\prime }}{{\alpha }_3^{\prime },{\alpha }_1^{\prime }}}{{\Gamma }_{loc}^{irr,p-h,m/d}}_{\genfrac{}{}{0pt}{}{{\alpha }_2,{\alpha }_2^{\prime }}{{\alpha }_1,{\alpha }_1^{\prime }}}{(i\nu ,i\nu 1)}_{i\omega }{{\chi }_{\mathbf{q}}^{p-h,m/d}}_{\genfrac{}{}{0pt}{}{{\alpha }_2^{\prime },{\alpha }_4^{\prime }}{{\alpha }_1^{\prime },{\alpha }_3^{\prime }}}{(i\nu 1,i\nu 2)}_{i\omega }\\ \hfill & {{\Gamma }_{loc}^{irr,p-h,m/d}}_{\genfrac{}{}{0pt}{}{{\alpha }_4^{\prime },{\alpha }_4}{{\alpha }_3{,}^{\prime }{\alpha }_3}}{(i\nu 2,i{\nu }^{\prime })}_{i\omega }\end{array}$$

(9)

The irreducible particle–particle vertex function ${\Gamma }^{irr,p-p}$ is finally written in terms of the reducible magnetic/charge vertex ${\tilde{\Gamma }}^{m/d}$ functions.

$$\begin{array}{cc}& {\Gamma }_{\genfrac{}{}{0pt}{}{{\alpha }_2\uparrow ,{\alpha }_4\downarrow }{{\alpha }_1\uparrow ,{\alpha }_3\downarrow }}^{irr,p-p}(\mathbf{k},i\nu ,{\mathbf{k}}^{\prime },i{\nu }^{\prime })={\Gamma }_{\genfrac{}{}{0pt}{}{{\alpha }_2\uparrow ,{\alpha }_4\downarrow }{{\alpha }_1\uparrow ,{\alpha }_3\downarrow }}^{f-irr}(i\nu ,i{\nu }^{\prime })\\ & -\frac{1}{2}[{\tilde{\Gamma }}^{p-h,\left(d\right)}\\ & -{\tilde{\Gamma }}^{p-h,\left(m\right)}{]}_{\genfrac{}{}{0pt}{}{{\alpha }_2,{\alpha }_3}{{\alpha }_1,{\alpha }_4}}\left({\mathbf{k}}^{\prime }-\mathbf{k},i{\nu }^{\prime }-i\nu \right)\\ & +{\tilde{\Gamma }}_{\genfrac{}{}{0pt}{}{{\alpha }_4,{\alpha }_3}{{\alpha }_1,{\alpha }_2}}^{p-h,\left(m\right)}\left(-{\mathbf{k}}^{\prime }-\mathbf{k},-i{\nu }^{\prime }-i\nu \right)\end{array}$$

(10)

$$\begin{array}{cc}& {\Gamma }_{\genfrac{}{}{0pt}{}{{\alpha }_2\downarrow ,{\alpha }_4\uparrow }{{\alpha }_1\uparrow ,{\alpha }_3\downarrow }}^{irr,p-p}(\mathbf{k},i\nu ,{\mathbf{k}}^{\prime },i{\nu }^{\prime })={\Gamma }_{\genfrac{}{}{0pt}{}{{\alpha }_2\downarrow ,{\alpha }_4\uparrow }{{\alpha }_1\downarrow ,{\alpha }_3\uparrow }}^{f-irr}(i\nu ,i{\nu }^{\prime })\\ & -\frac{1}{2}[{\tilde{\Gamma }}^{p-h,\left(d\right)}\\ & -{\tilde{\Gamma }}^{p-h,\left(m\right)}{]}_{\genfrac{}{}{0pt}{}{{\alpha }_4,{\alpha }_3}{{\alpha }_1,{\alpha }_2}}\left(-{\mathbf{k}}^{\prime }-\mathbf{k},-i{\nu }^{\prime }-i\nu \right)\\ & -{\tilde{\Gamma }}_{\genfrac{}{}{0pt}{}{{\alpha }_2,{\alpha }_3}{{\alpha }_1,{\alpha }_4}}^{p-h,\left(m\right)}\left({\mathbf{k}}^{\prime }-\mathbf{k},i{\nu }^{\prime }-i\nu \right)\end{array}$$

(11)

Finally, exploiting the Equations (6) and (7) and Equations (10) and (11) we obtain the ${\Gamma }^{irr,p-p}$ in the singlet (s) and triplet (t) channels from the magnetic and density particle-hole reducible vertices,

$$\begin{array}{cc}& {\Gamma }_{\genfrac{}{}{0pt}{}{{\alpha }_2,{\alpha }_4}{{\alpha }_1,{\alpha }_3}}^{irr,p-p,s}(\mathbf{k},i\nu ,{\mathbf{k}}^{\prime },i{\nu }^{\prime })={\Gamma }_{\genfrac{}{}{0pt}{}{{\alpha }_2,{\alpha }_4}{{\alpha }_1,{\alpha }_3}}^{f-irr}(i\nu ,i{\nu }^{\prime })\\ & +\frac{1}{2}[\frac{3}{2}{\tilde{\Gamma }}^{p-h,\left(m\right)}\\ & -\frac{1}{2}{\tilde{\Gamma }}^{p-h,\left(d\right)}{]}_{\genfrac{}{}{0pt}{}{{\alpha }_2,{\alpha }_3}{{\alpha }_1,{\alpha }_4}}{(i\nu ,-i{\nu }^{\prime })}_{{\mathbf{k}}^{\prime }-\mathbf{k},i{\nu }^{\prime }-i\nu }\\ & +\frac{1}{2}[\frac{3}{2}{\tilde{\Gamma }}^{p-h,\left(m\right)}\\ & -\frac{1}{2}{\tilde{\Gamma }}^{p-h,\left(d\right)}{]}_{\genfrac{}{}{0pt}{}{{\alpha }_4,{\alpha }_3}{{\alpha }_1,{\alpha }_2}}{(i\nu ,i{\nu }^{\prime })}_{-{\mathbf{k}}^{\prime }-\mathbf{k},-i{\nu }^{\prime }-i\nu }\end{array}$$

(12)

$$\begin{array}{cc}& {\Gamma }_{\genfrac{}{}{0pt}{}{{\alpha }_2,{\alpha }_4}{{\alpha }_1,{\alpha }_3}}^{irr,p-p,t}(\mathbf{k},i\nu ,{\mathbf{k}}^{\prime },i{\nu }^{\prime })={\Gamma }_{\genfrac{}{}{0pt}{}{{\alpha }_2,{\alpha }_4}{{\alpha }_1,{\alpha }_3}}^{f-irr}(i\nu ,i{\nu }^{\prime })\\ & -\frac{1}{2}[\frac{1}{2}{\tilde{\Gamma }}^{p-h,\left(m\right)}\\ & +\frac{1}{2}{\tilde{\Gamma }}^{p-h,\left(d\right)}{]}_{\genfrac{}{}{0pt}{}{{\alpha }_2,{\alpha }_3}{{\alpha }_1,{\alpha }_4}}{(i\nu ,-i{\nu }^{\prime })}_{{\mathbf{k}}^{\prime }-\mathbf{k},i{\nu }^{\prime }-i\nu }\\ & +\frac{1}{2}[\frac{1}{2}{\tilde{\Gamma }}^{p-h,\left(m\right)}\\ & +\frac{1}{2}{\tilde{\Gamma }}^{p-h,\left(d\right)}{]}_{\genfrac{}{}{0pt}{}{{\alpha }_4,{\alpha }_3}{{\alpha }_1,{\alpha }_2}}{(i\nu ,i{\nu }^{\prime })}_{-{\mathbf{k}}^{\prime }-\mathbf{k},-i{\nu }^{\prime }-i\nu }\end{array}$$

(13)

With ${\Gamma }^{irr,p-p}$ in hand we can solve the p-p BSE to compute the p-p susceptibility ${\chi }^{p-p}$.

$$\begin{array}{c}\hfill {\chi }^{p-p}={\chi }^{0,p-p}·{[\mathbf{1}+{\Gamma }^{irr,p-p}·{\chi }^{0,p-p}]}^{-1}\end{array}$$

(14)

The critical temperature $T_c$ is determined by the temperature where ${\chi }^{p-p}$ diverges. The pairing susceptibility diverges when the leading eigenvalue approaches unity. The corresponding eigenfunction represents the momentum structure of ${\chi }^{p-p}$. For such divergence the sufficient condition is that at least one eigenvalue of the pairing matrix ${\Gamma }^{irr,p-p}·{\chi }^{0,p-p}$ approaches unity. Hence $T_c$, eigenvalues $\lambda$ and eigenfunctions ${\varphi }^{\lambda }$ associated with different superconducting gap symmetries (in the singlet channel) can all be computed by solving the eigenvalue equation,

$$\frac{T}{N_k}\sum _{k^{\prime },i{\nu }^{\prime }}\sum _{\genfrac{}{}{0pt}{}{{\alpha }_2{\alpha }_4}{{\alpha }_5,{\alpha }_6}}{\Gamma }_{\genfrac{}{}{0pt}{}{{\alpha }_2,{\alpha }_4}{{\alpha }_1,{\alpha }_3}}^{irr,p-p,s}(k,i\nu ,k^{\prime },i{\nu }^{\prime })·{\chi }_{\genfrac{}{}{0pt}{}{{\alpha }_5,{\alpha }_6}{{\alpha }_2,{\alpha }_4}}^{0,p-p}(k,i{\nu }^{\prime }){\varphi }_{{\alpha }_5{\alpha }_6}^{\lambda }=\lambda ·{\varphi }_{{\alpha }_5{\alpha }_6}^{\lambda }$$

(15)

The gap function can be written in a symmetric and Hermitian form by

$$\frac{T}{N_k}\sum _{k^{\prime },i{\nu }^{\prime }}\sum _{\genfrac{}{}{0pt}{}{{\alpha }_2{\alpha }_4{\alpha }_5}{{\alpha }_6{\alpha }_7{\alpha }_8}}{\left({\chi }_{\genfrac{}{}{0pt}{}{{\alpha }_2,{\alpha }_4}{{\alpha }_1,{\alpha }_3}}^{0,p-p}(k,i\nu )\right)}^{1/2}·{\Gamma }_{\genfrac{}{}{0pt}{}{{\alpha }_5,{\alpha }_7}{{\alpha }_2,{\alpha }_4}}^{irr,p-p,s}(k,i\nu ,k^{\prime },i{\nu }^{\prime })·{\left({\chi }_{\genfrac{}{}{0pt}{}{{\alpha }_6,{\alpha }_8}{{\alpha }_5,{\alpha }_7}}^{0,p-p}(k^{\prime },i{\nu }^{\prime })\right)}^{1/2}·{\varphi }_{{\alpha }_6{\alpha }_8}^{\lambda }(k^{\prime },i{\nu }^{\prime })=\lambda ·{\varphi }_{{\alpha }_1{\alpha }_3}^{\lambda }(k,i\nu )$$

(16)

However, it can also be explicitly shown that the eigenvalues of the non-Hermitian gap equation are the same as eigenvalues of the Hermitian gap equation.

Finally, ${\chi }^{p-p}$ can be represented in terms of eigenvalues $\lambda$ and eigenfunctions ${\varphi }^{\lambda }$ of the Hermitian particle–particle pairing matrix.

$$\begin{array}{ccc}\hfill {\chi }^{p-p}(k,i\nu ,k^{\prime },i{\nu }^{\prime })& =& \sum _{\lambda }\frac{1}{1-\lambda }·(\sqrt{{\chi }^{0,p-p}(k,i\nu )}·{\varphi }^{\lambda }(k,i\nu ))\hfill \\ & ·(& \sqrt{{\chi }^{0,p-p}(k^{\prime },i{\nu }^{\prime })}·{\varphi }^{\lambda }(k^{\prime },i{\nu }^{\prime }))\hfill \end{array}$$

(17)

To solve this eigenvalue equation, the most important approximation we make is to take the static limit of ${\Gamma }^{irr,p-p}$ in the bosonic frequency i$\omega$ = 0 (real frequency axis). The explicit dependence on the fermion frequencies are kept, as are all the orbital and momentum indices.

As is apparent from Equations (13) and (14) at what wave vector spin and charge fluctuations are strong is of central importance to the kind of superconducting pairing symmetry they can form. The entire momentum, orbital and frequency dependence of the vertex functions are computed explicitly and the BSE equations are solved with them. The crucial point is that the vertex structure has no predefined form-factor, so the emergent superconducting gap symmetry is calculated in an unbiased manner. This provides an unbiased insight into the superconducting gap symmetries, the strength of the leading eigenvalues in different systems and, most importantly, allows for a fair comparison of the relative strength of the leading superconducting instabilities in bulk FeSe and the monolayer of FeSe/STO. Thus, our ability to predict these properties is limited mostly by the fidelity of the Green’s functions that determine the vertices and $\chi$.

There is a practical limitation, however. Since we compute the vertex functions from CTQMC, which limits the temperatures down to which the vertex can be computed. We have observed in different materials that the leading eigenvalue $\lambda$ does not have a simple, analytic dependence on temperature [22], and hence $\lambda$ can not be reliably extrapolated to very low temperatures. For that reason, we avoid estimating $T_c$ (the temperature at which $\lambda$ reaches 1) for different systems from our method, rather, we compare the strength of $\lambda$ for a given temperature in different materials, which is free from any ambiguities.

In some recent studies [32,33,34,35,36], the shortcomings of the standard Eliashberg formalism in explaining some of the fundamental observation related to unconventional superconductivity has been highlighted and alternative formulation for a better theory is proposed. While our particle–particle ladder approximation does not use the standard Eliashberg formalism, it is an important task for future efforts, also for our own approach, to see how these proposed formulations can be adapted for better description of unconventional superconductivity.

3. Results and Discussion

3.1. Structural Parameters and c-RPA Estimates for U and J

Table 1. Structural parameters, chalcogen height $h_{Se}$ and computed U and J for the correlated many body Hamiltonian from our QSGW+c-RPA implementation. We also list the bare electronic bandwidths from QSGW which show how electronic correlations enhance in M-FeSe/STO in comparison to the bulk as the bands get narrower. References indicate where the structural inputs were taken.

Variants	a (Å)	c (Å)	${\mathit{h}}_{\mathit{se}}$(A${}^0$)	U (eV)	J (eV)	d${}_{\mathit{xy}}$ (eV)	d${}_{\mathit{yz},\mathit{xz}}$ (eV)
B-FeSe [37]	3.779	5.5111	1.463	3.4	0.60	2.5	3.3
B-FeSe (reduced h${}_{se}$)	3.779	5.511	1.27	3.9	0.69	3.1	4.3
M-FeSe [38]	3.905		1.39	4.3	0.71	2.4	3.2
M-FeSe/STO [38]	3.905		1.40	3.8	0.67	2.35	2.95

provides structural parameters used for simulating different systems, and the correlation parameters calculated using QSGW+constrained-RPA. In the bulk crystal, Se sits above and below the Fe plane with height $h_{Se}=1.463$ Å and bond length $l_{Fe-Se}=2.39$ Å. cRPA yields U = 3.4 eV and J = 0.6 eV, calculated from the QSGW band structure.

Over the large variety on iron based superconductors it is observed [6] that larger chalcogen/pnictogen height ($h_{se}$ in our case) reduces the Fe-chalcogen/pnictogen-Fe hopping and hence makes the systems more strongly correlated. We observe that $h_{se}$ is reduced in both M-FeSe/STO and free-standing monolayer M-FeSe relative to the bulk FeSe (B-FeSe), however, the Fe-Se bond lengths remain almost invariant in both monolayer variants in comparison to the bulk (as a increases in layer). It suggests that the monolayer variants have similar or less electronic correlations (larger electronic hopping). However, as we will show, the free standing monolayer is non-superconducting while both the bulk FeSe and M-FeSe/STO are superconducting. Further, the c-RPA calculations for each materials indicate that both U and J increase in the monolayer variants relative to the bulk.

As we will show below, the degree of correlation and resulting superconductivity are highly sensitive to small changes in both quasiparticle levels and the Hubbard parameters driving the correlations. Thus, the electronic origins of unconventional superconductivity in FeSe can only be understood when the interplay between crystal structure, one-particle properties, and Hubbard parameters is taken into account.

3.2. Fermi Surfaces and Spectral Functions

QSGW Fermi surfaces are generally smaller than DFT Fermi surfaces in the Fe superconductors [39]. However, in FeSe, both DFT and QSGW predict the $d_{xy}$ band to cross the Fermi level forming a pocket around $\Gamma$ (see Figure 1), apparently in contradiction to ARPES measurements which place it at approximately $-17$ meV [40,41,42]. When $\Sigma \left(\omega \right)$ is further dressed by DMFT this pocket shrinks further, but it still remains, leaving a discrepancy with ARPES of ∼20 meV. The $d_{xy}$ orbital character emerges as the most incoherent (see Figure 2a) but its position is slightly higher than where ARPES places it.

3.3. Bulk FeSe

The most direct way to isolate the contributors to superconductivity is to make parametric adjustments to the ab initio results. Thus to assess the role of ‘Hundness’ we consider, in addition to the ab initio QSGW⁺⁺ calculations, how excursions in J between 0 and 1 eV affect the single-particle scattering rate $\Gamma$, and the inverse Z factor, a measure of mass or bandwidth normalisation.

The single-site DMFT Im$\Sigma \left(i\omega \right)$ is fit to a fourth order polynomial in $i\omega$ for low energies (first 6 matsubara points at $1/\beta$ = 1/40 eV = 290 K). The mass enhancement, related to the coefficient ($s_1$) of the linear term in the expansion $m_{DMFT}/m_{QSGW}=1+\left|s_1\right|$ [43], and the intercept $|s_0|$ = $\Gamma m_{DMFT}/m_{QSGW}$. m${}_{DMFT}$/m${}_{QSGW}$ = $Z^{-1}$ is resolved in different intra-orbital channels. $\Gamma$ is found to be insensitive to J for $J<0.4$ eV, while for $0.4<J<0.7$ eV it increases monotonically for all orbitals (see Figure 3). The size of the local moment also increases steadily and reaches a value of 2.2 ${\mu }_B$ for $J=0.6$ eV which is the atomic moment experimentally observed in Fe. This observed moment in bulk FeSe is also very close to the spin local moment estimated from NMR, INS and x-ray emission spectroscopy studies [44]. Thus there is smooth transition from coherence to incoherence with increasing J. $1/Z$ increases from 1.33 at $J=0$, reaching a maximum of 4.5 in the $d_{xy}$ channel at $J\sim 0.68$ eV. Correlation increases for all states, but $d_{xy}$ is the heaviest and most incoherent, followed by $d_{xz}+d_{yz}$ (see Figure 3). For still larger J, both ${\Gamma }_{xy}$ and $1/Z_{xy}$ begin to slowly decrease (see Figure 3). A similar non-monotonic behaviour was observed in a previous study of Hund’s metals [45] Similar conclusions were drawn in a recent orbitally-resolved quasi-particle scattering interference measurement by Kostin et al. [5] in the low temperature orthorhombic phase of FeSe, and the orbital selectivity was emphasised in Ref. [46].

Small increases in J induce remarkable changes in the transverse spin susceptibility $\chi (q,\omega )$ (see Figure 4a,b). The peak near the antiferromagnetic nesting vector $(1/2,1/2)$ (in 2-Fe unit cell) and $q_z$ = 0 becomes markedly more intense as shown in Figure 2e and Figure 4b. In Figure 4, Im $\chi (q,\omega )$ is plotted along the $(0,0)\to (1/2,0)\to (1/2,1/2)\to (0,0)$ lines in the Fe${}_2$Se${}_2$ unit cell. The energy dispersion in Im$\chi$ also becomes strongly compressed. Elsewhere [27] we perform a rigorous benchmarking of $\chi (q,\omega )$ against inelastic neutron scattering measurements [28,29,47] where we show that our calculations reproduce all intricate structures of $\chi (q,\omega )$ for all energies and momenta. Resolving Im $\chi$ into orbital channels, $d_{xy}$ is seen to be the leading component. Along $(1/2,0)\to (1/2,1/2)$ it contributes about 50% of the total with $d_{z^2}$, $d_{x^2-y^2}$ and $d_{xz,yz}$ combining to contribute the rest.

How do these striking changes in Im $\chi$ correlate with superconductivity? We compute the full two particle scattering amplitude in the particle–particle channel within our DMFT framework, and solve Equation (17) in the BCS low energy approximation [19,20,21,22]. Resolving the eigenfunctions of the gap equation into inter- and intra-orbital channels, two dominant eigenvalues $\lambda$ are found. Both of them increase with J up to the point of maximum intensity in $\chi$ (J = 0.68 eV) and then begin to decrease, as shown in Figure 2f. The corresponding eigenfunctions have extended s-wave- cos($k_x$)+cos($k_y$) (leading eigenvalue) and $d_{x^2-y^2}$− cos($k_x$) − cos($k_y$) structures (second eigenvalue) [48]. Calculations show that these instabilities reside primarily in the intra-orbital $d_{xy}$ − $d_{xy}$ channel and the inter-orbital components are negligible. In the bulk crystal, varying J from the ab initio value ($\lambda$ = 0.067 at J = 0.60 eV), we find $\lambda$ reaches its maximum 0.9 at the point where the spin susceptibility is most intense, J=0.68 eV (see Figure 2f). We attribute the decrease for $J>0.68$ eV to the softening of electron masses and loss of spin fluctuations at $q=(1/2,1/2)$ as can be seen in local moment ${〈M^2〉}^{1/2}$ plotted in Figure 3.

3.4. Excursion in Fe-Se bond length in Bulk FeSe

We next consider how parametric changes in the one-body Hamiltonian alter the spectral function, $\chi$ and $T_c$. We first vary the Fe-Se bond length $l_{Fe-Se}$. When it is reduced from its experimental value of 2.39 Å, the $d_{xy}$ band initially near $E_F$ at $\Gamma$, gets pushed down well below $E_F$. It is particularly easy to see at the QSGW level (blue band in Figure 2b), reaching about $E_F-160$ meV when h = 1.27 or $l_{Fe-Se}$ = 2.275 Å. From c-RPA we compute the corresponding U and J as 3.9 eV and 0.69 eV respectively. A similar shift is found in the spectral function calculated by QSGW+DMFT (Figure 2b). Further, quasi-particles become more coherent: orbital-selective $1/Z$ values range between 1.6 and 1.3 and scattering rates become small (○, Figure 5). The system behaves as an itinerant metal, and the peak in Im$\chi (q,\omega )$ at $(1/2,1/2)$ shifts to higher $\omega$ and becomes gapped (see Figure 4c). It also becomes very weak and broad (Figure 2e). The leading eigenvalue of Equation (17) become negligibly small (see Figure 2f and Figure 5e), suggesting extremely weak or no superconducting instability. $\lambda$ also becomes completely insensitive to J. A similar observation was made in our recent work on LaFe${}_2$As${}_2$ [22], where the collapsed tetragonal phase with lesser $l_{Fe-As}$ in comparison to its uncollapsed phase, loses superconductivity [49] as bands become itinerant due to increased Fe-Fe hopping mediated via As. It suggests that a small pnictogen/chalcogen height over the Fe-Fe square plane is not conducive for superconductivity as the system gains significant amount of kinetic energy and electronic scatterings are remarkably reduced. Similar observations are made for Fe based pnictides containing Phosphorus, where the system has a Fermi liquid normal phase with dispersive electrons and either they do not superconduct or have extremely small $T_c$’s [20].

3.5. Free Standing Monolayer of FeSe

We next turn to a free standing monolayer of FeSe, M-FeSe. For this study we take the structural inputs from recent work by Mondal et al. [38] which finds the minimum-energy value for h to be 1.39 Å within a combined DFT+DMFT framework, while the lattice parameter a is somewhat larger, close to that of SrTiO${}_3$. As a benchmark, the same work found the equilibrium h to be close to the measured value in bulk crystalline FeSe [50]. DFT has a tendency to underestimate h; the error is not easily fixed by other kinds of density-functionals. However, DFT does predict a similar change in $\delta h$ between bulk and monolayer which gives us some confidence in the value of h.

At the QSGW level, the $d_{xy}$ band is pushed to $E_F-300$ meV on the $\Gamma$–M line (band structure in Figure 2c and SM, bottom panel). c-RPA calculations yield U = 4.3 eV and J = 0.71 eV, the increase arising from reduced screening. $1/\mathrm{Z}$ is found from QSGW+DMFT to be 1.5, 1.35, 1.3, 1.25 on $xy$, $z^2$, $yz+xz$, $x^2-y^2$ respectively (Figure 5b); also a negligible scattering rate is found (Figure 5a). As a further indicator of a good metal, Im $\chi (q,\omega )$ shows negligible spin fluctuations in the $d_{xy}$ channel; and at $q=(1/2,1/2)$ spin excitations are gapped and vanishingly small (see Figure 2e and Figure 4e). The superconducting instability is almost entirely suppressed (see Figure 2f and Figure 5f). It is noteworthy that the reduction in electronic screening reflects in a marked increment in J. Unfortunately, this beneficial effect is more than counterbalanced by the fact that $d_{xy}$ is pushed far below $E_F$ on the scale of magnetic excitation energies. Suppression of low energy one- and two-particle scattering is not conducive for superconductivity.

3.6. Monolayer of FeSe/SrTiO${}_3$

How does the effect of a SrTiO${}_3$ substrate modify superconductivity in M-FeSe? M-FeSe/STO is a subject of intense debate, since as noted above, $T_c$ of M-FeSe/STO has been measured to be an order of magnitude higher than bulk FeSe (In the Supplementary Materials we include a small section on uniformly electron doped bulk FeSe). Many explanations have been put forward, e.g., that superconductivity is boosted by large electron-phonon coupling [51,52,53] as SrTiO${}_3$ is close to a ferroelectric instability (for a contrary view see [54]), but the simplest explanation is that SrTiO${}_3$ modifies M-FeSe to restore $d_{xy}$ to be proximate to $E_F$. M-FeSe/STO is a partially formed Schottky barrier; we can expect the Fermi level to sit in the SrTiO${}_3$ bandgap. SrTiO${}_3$ modifies M-FeSe in two important ways: the interfacial dipole controls the Schottky barrier height and changes the electron count in M-FeSe; also the STO (especially the O-p-derived bonding states) couple to the Fe-d in an orbital-selective manner. Both effects are accurately incorporated by a direct QSGW calculation of M-FeSe/STO. We consider 5-ML slab of SrTiO${}_3$ terminated on the Sr side by M-FeSe. The structure is relaxed with DFT, subject to the constraint that h is fixed to 1.40 Å, as predicted in Ref. [38]. This value is close to h = 1.39 Å ($l_{Fe-Se}$ = 2.403 Å) found for free-standing M-FeSe. Its value is critical, as we have seen in the bulk case, and we cannot rely on DFT for it.

A c-RPA calculation for M-FeSe/STO yields U = 3.8 eV and J = 0.67 eV. Using QSGW + DMFT, we extract the Fermi surfaces spectral functions (see Figure 1d and Figure 2d). Remarkably, while $d_{xy}$ is pushed far below $E_F$ in M-FeSe, its position returns near to $E_F$ in the M-FeSe/STO case ( to $E_F-50$ meV in QSGW and $E_F-100$ meV in QSGW + DMFT as shown in Figure 2d). This is fully consistent with ARPES studies [55]. Furthermore, its bandwidth is slightly reduced relative to bulk, in keeping with a stretched a.

In the QSGW⁺⁺ calculation, the hole pockets of $d_{xz,yz}$ character survive although are significantly narrowed. This is a discrepancy with ARPES, which finds no such pockets. However, their small size suggests that they should be sensitive to electron-phonon coupling that can further renormalise them and push them down.

Using QSGW+DMFT we compute the orbital dependent $1/Z$ (4, 3.3, 3.2 2.7 on $xy$, $yz+xz$, $z^2$, $x^2-y^2$) and $\Gamma$, which show significantly enhanced incoherence in the quasiparticle spectrum relative to M-FeSe (see Figure 5a,b). Further, the dispersion in Im $\chi (q,\omega )$ is significantly narrower than the bulk, and the intensity is spread over momenta in the region around $(1/2,1/2)$ (see Figure 2e and Figure 4f). These signatures suggest that M-FeSe/STO is more correlated than either the bulk or free standing monolayer. The leading extended-s wave instability from the particle–particle ladder BSE equation survives but the second ($d_{x^2-y^2}$) instability gets suppressed. Whenever the $d_{xy}$ state is pushed below $E_F$ at $\Gamma$ and yet continues to be present at the electron pockets, the only surviving superconducting instability is the $s\pm$ state, consistent with earlier predictions [56,57,58].

Experimentally, the superconducting gap of FeSe/STO show maxima on the bands with $d_{xy}$ character [59]. Furthermore, in surface doped bulk FeSe where superconductivity enhances, the appearance of the $d_{xy}$ electron pocket on the Fermi surface coincides with the beginning of the second superconducting dome with higher $T_c$ [60]. The leading $\lambda$ for superconducting instability is 0.34, which is five times larger that the estimate for bulk FeSe. These ab initio results do not support a recently published results from Eliashberg theory where it is claimed that spin fluctuations can, at most, account for only two fold increment in $T_c$. [61] We also find that suppressing the $d_{xz,yz}$ contributions in Equation (17) lead to only 6% reduction in $\lambda$. Thus, the $d_{xy}$ orbital is the determinative one.

Our estimate of five-fold increment is less than the eight-fold increment in $T_c$ = 75 K observed for M-FeSe/STO, but there can be many reasons for this. The comparison to experiment is not direct: experimentally M-FeSe/STO has a buffer layer between M-FeSe and STO, which our calculation omits. Furthermore, the calculated $T_c$ is extremely sensitive to J and h. Both are theoretically calculated; moreover h is assumed to be the same for the Se planes above and below Fe, while there should be small differences. Finally the interface can have several other effects we omitted such as an enhanced phonon contribution to $T_c$, as others have suggested [51,52,53].

We close with a note on the origins of the discrepancies in QSGW⁺⁺ spectral function with ARPES. We have recently developed the ability to incorporate the electron-phonon self-energy into QSGW⁺⁺ via a field-theoretic technique. While initial results are very preliminary, they suggest that much of the discrepancy we see with ARPES in bulk FeSe, originates from this interaction. It is perhaps not surprising, given that the electron-phonon interaction is expected to be stronger when pockets are small. These new findings, however, are beyond the scope of this work.

4. Conclusions

To summarise, a unified picture of the origins of superconductivity in FeSe emerges from evidence drawn from several parametric studies of FeSe around a high-fidelity ab initio theory. Superconducting glue mainly originates from low-energy spin fluctuations concentrated in a region near the antiferromagnetic ordering vector. The instability and the single- and two-particle correlations characterised by the band renormalisations $1/Z$, scattering rate $\Gamma$ and Im$\chi (q,\omega )$ are all closely linked as summarised in Figure 2 and Figure 5. The Fe $d_{xy}$ orbital is the most strongly correlated and contributes maximally to the pairing glue as long as it is auspiciously near the Fermi level and further, that the Hund’s J is sufficiently large to induce a high degree of ‘bad metallic’ behaviour. Further, the superconducting instability is found to lie predominately in the $d_{xy}$ channel.

We show that in the bulk, $T_c$ is extremely sensitive to J over a certain energy window. At the optimum J $\lambda$ increases by nearly 15 times its value at the ab initio J, even though the increment in J is only 10–15%. This indicates that a possible way to enhance $T_c$ is to reducing the electron screening. We presented two scenarios where a structural modification induces a change in J (M-FeSe and M-FeSe/STO), and further that this enhancement drives a five fold enhancement in $\lambda$, in M-FeSe/STO relative to the bulk. M-FeSe and M-FeSe/STO differ in the latter preserves the proximity of the $d_{xy}$ channel to $E_F$ while M-FeSe does not. This provides a natural explanation for the enchantment in superconducting M-FeSe/STO as a purely electronic one: the two dominant factors are J and the proximity of the $d_{xy}$ band to $E_F$.

Controlling number of layers, applying pressure to tune Fe-chalcogenide bond length, doping and intercalation are some other possibilities. At the same time it is important to realise it is necessary to control both the screening and specific features of the single-particle spectrum. Our calculations show that one promising directions for reaching an optimised $T_c$ appears to be controlling number of layers and interfaces to simultaneously satisfy both conditions: lesser electron screening leading to a larger Hund’s correlation and larger $d_{xy}$ contribution to the Fermiology.