Vertex dominated superconductivity in intercalated FeSe

Abstract

Bulk FeSe becomes superconducting below 9 K, but the critical temperature (T_c) is enhanced almost universally by a factor of ~4–5 when it is intercalated with alkali elements. How intercalation modifies the structure is known from in-situ X-ray and neutron scattering techniques, but why T_c changes so dramatically is not known. Here we show that there is one-to-one correspondence between the enhancement in magnetic instabilities at certain q vectors and superconducting pairing vertex, even while the nuclear spin relaxation rate 1/(T₁T) may not reflect this enhancement. Intercalation modifies electronic screening both in the plane and also between layers. We disentangle quantitatively how superconducting pairing vertex gains from each such changes in electronic screening. Intercalated FeSe provides an archetypal example of superconductivity where information derived from the single-particle electronic structure appears to be insufficient to account for the origins of superconductivity, even when they are computed including correlation effects. We show that the five-fold enhancement in T_c on intercalation is not sensitive to the exact position of the d_xy at Γ point, as long as it stays close to E_F. Finally, we show that intercalation also significantly softens the collective charge excitations, suggesting the electron-phonon interaction could play some role in intercalated FeSe.

Introduction

With the rise of layered materials^{1,2,3,4,5,6,7}, intercalation^8,9,10 and exfoliation^{11,12,13,14,15} became two of the most commonly used methods to alter their structural and electronic properties. Such intercalation with elements or molecular moieties changes the separation between layers, screening environment and electronic structure, including Fermi surface properties. Indeed, shortly after their first realisation^5,16 in bulk crystalline form, intercalation has become one of the most popular methods to modify magnetic and superconducting characteristics of the iron-based superconductors (IBS).

Bulk FeSe superconducts up to 9 K, deep inside an orthorhombic phase that sets in at a much higher temperature, 90 K¹⁷. The tetragonal phase, can be made to superconduct in various ways, e.g. through doping^{18,19,20,21,22}, or pressure^23,24,25, as a monolayer^26,27,28, or when intercalated with alkali elements^{29,30,31,32,33,34} (Li, Na, K, Cs), under surface doping³⁵ and ionic liquid gating^36,37,38,39. Remarkably, the critical temperature T_c of intercalated FeSe is enhanced by roughly 4–5 times over the bulk, for many kinds of intercalated variants. What leads to such dramatic enhancement in T_c remains unanswered even after a decade of investigation.

Theoretical attempts to answer this question face some challenges as well. One of the primary challenges for ab-initio theoretical calculations in such situation relates to the lack of proper information of the crystalline structure post intercalation. In recent years, significant progress has been made on the front of determination of the crystal structure using in situ X-ray and neutron powder diffraction techniques^40,41,42. Further, magnetometry and muon-spin rotation techniques are used to determine the superconducting properties of the same sample, leading to unambiguous determination of both its structural and superconducting properties post intercalation. Reliable information of the crystal structure helps in performing ab-initio theoretical calculations for these materials.

FeSe is a strongly correlated material, which is reflected in its spectral properties and spin fluctuations^43,44,45,46. Its particularly large Hund’s coupling^45,47 drives orbital differentiation⁴⁸ and large electronic vertex corrections to its two particle instabilities^45,49. Previous theoretical works⁵⁰ on intercalated FeSe connect the enhancement in T_c to the enhancement in Fermi surface nesting and subsequent enhancement in density of states at Fermi energy ρ(E_F). This is the usual line of argument in superconductors where attempts have been made to connect the enhancement in T_c to increase in ρ(E_F) in the spirit of BCS theory⁵¹. In BCS theory T_c has an exponential dependence on ρ(E_F). The positive correlation of ρ(E_F) and T_c was broadly discussed for both conventional superconductors such as e.g. A15 or C15 families⁵² and for high-temperature cuprates^53,54.

On the other hand, DFT calculations of intercalated FeSe finds rather weak electron doping compared to the parent compound⁵⁰ and thus no significant change of ρ(E_F) is expected. We observe the same in both DFT and many body-perturbative approaches (quasi-particle self-consistent GW, QSGW⁵⁵, Fig. 2). While QSGW predicts ρ(E_F) to be weakly suppressed in the intercalated material, it is often noted that FeSe is not non-magnetic, rather paramagnetic. We take this into account by augmenting QSGW with Dynamical Mean Field Theory (DMFT), QSGW+DMFT⁵⁶. Paramagnetic DMFT simulates site-local magnetic fluctuations that are crucial for FeSe^{45,49,57,58,59,60,61}. When QSGW is augmented with DMFT⁵⁶, the primary conclusion remains, namely intercalation weakly suppresses ρ(E_F). In the absence of two-particle vertex corrections, T_c can only drop in such situation. Thus, the dramatic enhancement in T_c can not be explained based on any theory that is reliant on electronic density of states alone, and a primary aim of this work is to show the quintessential driving force for enhancement of T_c originates from the two-particle sector. This is revealed through careful examination of the orbital, frequency and momentum dependence of two-particle vertex functions, in magnetic and superconducting channels, and how they quantitatively determine the nature of T_c enhancement.

It is not totally unexpected since even in the BCS theory the constant λ determining T_c is the product of ρ(E_F) and the effective inter-electron interaction, and the primary reason why common attention is focused on ρ(E_F) is that it is much easier to calculate. This simplification can be sometimes dangerous, and as we will show here, intercalated FeSe is a very clear example where all the essential changes happen in the effective interaction only. In a more formal language, this effective interaction is nothing but the vertex, that is, an essentially two-particle characteristic of correlated systems.

When materials are intercalated, the separation between active layers for superconductivity (in this case the Fe–Fe square planes) enhances and most experimental studies attempt to establish a relation between T_c and the separation between layers. Such phenomenological attempts are crucial to advance the technology to the next level, where we will be able to intercalate samples in a controlled manner to tune superconducting properties as desired. However, most such attempts turn out to be failures. In intercalated FeSe, there are several samples with very large inter-layer separation where T_c nevertheless remains invariant. That being said, when the inter-layer separation enhances it reduces significantly the electronic screening and that leads to larger Hund’s coupling. This aspect of the correlated Fe-3d Hamiltonian is often entirely discarded from theoretical analysis, but it turns out to be crucial. Within our ab-initio constrained QSGW-RPA theory, we compute the changes to the bare two-particle parameters entering into the Anderson impurity model, namely the Hubbard U and Hund’s J, in an unbiased fashion, and show how they affect T_c.

We quantitatively analyse the consequences of the following: changes to the Fe–Se–Fe bond angles, changes to the Hund’s J driven by enhancement in layer separations in building the correlated Hamiltonian for intercalated FeSe (see Fig. 1). We solve the correlated two-particle Hamiltonian with Bethe–Salpeter equations in the magnetic and superconducting channels to show in detail how the five-fold T_c enhancement in Li/NH₃ intercalated FeSe originates. We analyse the parent non-intercalated bulk tetragonal FeSe (n − i), a bulk FeSe structure simulated only using the planar distortion (p − d) parameters (Fe–Se bond lengths and Fe–Se–Fe bond angles reported in the supplemental materials Table 2 of ref. ⁴¹); and also intercalated FeSe with the full structural (f − d) distortion as reported in the same work. The primary difference between the planar-distortion (p − d) and full-distortion (f − d) is the enlargement of the c-axis, that reduces the electronic screening perpendicular to the Fe–Fe square plane (see Fig. 1), leading to about 10% enhancement to the Hubbard parameters for the Fe-3d orbitals (in n − i and p − d, U = 3.5 eV and J = 0.6 eV and in f − d, U = 3.76 eV and J = 0.68 eV). While several prior studies have explored the relation between T_c and structural parameters in iron-based superconductors^62,63,64,65, here we are able to disentangle the effects coming from planar distortions and inter-layer separations, and show quantitatively how each of these structural changes modifies the strength of the pairing instability. Changes in crystal structure leads to changes in electronic screening which affects both the single-particle and two-particle sectors, and the enhancement in T_c can only be explained when both are treated in the presence of higher-order vertex corrections. Previous works have analysed the screened Coulomb vertex corrected pairing instabilities. in a one-band model for cuprates⁶⁶ and a two-band model for nickelates⁶⁷. However, our solutions for Bethe–Salpeter equations in pairing channel keeps full orbital-, momentum- and two-frequency-dependence (fermionic Matsubara frequencies) of all five Fe-3d correlated orbitals. What is key for our Cooper pairing are the reducible vertex functions⁶⁸ in the particle-hole magnetic and charge channels that contribute to superconducting pairing equation. These reducible quantities pick up explicit momentum dependence (and the sign for the pairing interaction) when the outgoing propagator lines are flipped. The only approximation we make is to put the centre-of-mass co-ordinates to zero (centre-mass momentum q = 0 and the bosonic frequency Ω = 0) since the linearised Eliashberg equation simplifies to an eigenvalue problem that can be solved then at different temperatures in the normal phase looking for all possible pairing instabilities and their competitions^69,70 in a fully unbiased fashion.

Fig. 1: Separability of conditions for dramatic enhancement in T_c on intercalation.

A two-step increment in T_c on intercalation is explored; the first is due to smaller Fe–Se–Fe bond angle that reduces the Fe–Se–Fe hopping and makes Fe 3d states more correlated, and the second is due to enhanced c-axis length that reduces electronic screening and enhances correlations. The rest of the letter discusses how these two mechanisms strongly modify the pairing vertex while the one-particle density of states remain nearly unchanged.

Results and discussion

Effect of intercalation on electronic density of states

We start by discussing the one-particle properties: the electronic band structure, density of states (DOS) and the Fermi surfaces computed using LDA and QSGW (see Fig. 2). ρ(E_F) enhances slightly going from n − i to p − d, since the Fe–Se–Fe bond angle α reduces slightly, leading to reduced Fe–Se–Fe hopping and narrower Fe bands. However, the effect is nullified in f − d owing to the enlarged c-axis (see Fig. 2d–f). In this case ρ(E_F) decreases. The electron pockets at M becomes slightly larger, while the hole pockets at Γ become smaller, leading to net electron doping of the system (see Fig. 2a–c). Nevertheless, we find that the electron doping within the QSGW approximation remains as small as ~1%. The situation is similar in QSGW+DMFT where ρ(E_F) in f − d also drops relative to n − i and p − d. Thus, the effective electron doping of the system remains rather weak. If the vertex were to remain constant, a five-fold enhancement in T_c would require about 80% enhancement in electronic density of states in a BCS picture (${T}_{c,{{{\rm{BCS}}}}} \sim \exp [-1/V\rho ({E}_{F})]$), (assuming the correlation parameters V remain unchanged).

Fig. 2: Failure of one-particle properties from LDA, QSGW and QSGW+DMFT in explaining the enhancement in T_c.

The orbital projected electronic QSGW Fermi surfaces for (a) n − i FeSe, (b) p − d FeSe and (c) f − d FeSe are shown. The Fe-3d${}_{{x}^{2}-{y}^{2}}$+d${}_{{z}^{2}}$ orbitals are shown in blue, d_xz,yz in green and d_xy in red. Finally, the total density of states (d–f) from different levels of the theory are shown. In QSGW+DMFT the total density of states decrease in intercalated sample, and yet T_c enhances. k_x,k_y are in units of $\frac{2\pi}{a}$.

Magnetic susceptibilities with and without vertex corrections

We now turn our attention to two-particle instabilities, particularly the frequency and momentum resolved magnetic susceptibility Im χ^m(ω, q) (Fig. 3). It is computed by solving a non-local Bethe–Salpeter equation (BSE) that dresses the non-local polarisation bubble in the particle-hole magnetic channel by the local irreducible dynamic vertex^49,69. Both the inputs for the BSE, the single-particle vertex which enters the self-energy and the two-particle vertex entering into response functions, are computed using DMFT. Our computed Im χ^m(ω, q) has been rigorously benchmarked against Inelastic Neutron Scattering measurements³⁰ in a prior study⁴⁹ over all relevant energies and momenta. We analyse Im χ^m(ω, q) along the reciprocal lattice vectors q = (q_x, q_y, q_z), in units of 2π/a. We fix q_z at 0 and explore the susceptibilities in the q = (q_x, q_y) plane. In the vicinity of q = (1/2,1/2), intercalation (p − d and f − d) significantly enhances the strength of Imχ^m(ω, q) at lower energies, making spin-fluctuation mediated superconductivity more favourable. We also observe that this tendency is independent of temperature; nevertheless, with lower temperatures the strength of spin-glue at low-energies for superconductivity becomes more prominent. To understand the essential ingredient of such spin susceptibilities in the static limit, we recompute χ^m(q) = χ^m(ω = 0, q) with a vertex and without (RPA).

Fig. 3: Enhancement in low energy magnetic glue for superconductivity on intercalation.

The vertex corrected dynamic and momentum resolved magnetic susceptibility Imχ^m(ω, q) is shown for (a) n − i, (b) p − d and (c) f − d at 100 K. On intercalation, Im χ^m(ω, q) becomes more intense at low energies, particularly at q = (1/2, 1/2) which corresponds to the anti-ferromagnetic instability vector in 2-Fe atom unit cell of FeSe. Also with lowering of temperature the low energy structure of Imχ^m(ω, q) becomes more prominent.

Superconducting susceptibilities with and without vertex corrections

Within the RPA, χ^m(q) is resolved in different intra-orbital channels (Fig. 4a, b) and it remains the largest in n − i at all q. Intercalation weakly suppresses χ^m(q) in different inter-orbital channels suggesting a weak suppression in superconductivity in a purely RPA picture of spin-fluctuation mediated superconductivity. Our computed χ^m(q) has similarities with the susceptibility computed in ref. ⁵⁰. Particularly striking is the very large intra-orbital elements of χ^m(q) at q = 0. This is primarily a shortcoming of the RPA approximation: without vertex corrections the momentum dependence of χ^m(q) is rather weak and there is no clear separation between χ^m at q = 0 and q = (1/2,1/2). However, experimentally, in bulk FeSe magnetic fluctuations at q = 0 are rather weak compared to the main instability at q = (1/2,1/2)³⁰ and it is the vertex that brings in the needed momentum-dependent variation in χ^m, as we have shown previously⁴⁹.

Fig. 4: Shortcoming of the RPA theory for magnetic and superconducting susceptibilities.

χ^m and χ^pp are resolved in different intra-orbital channels. RPA theory predicts magnetic (a, b) and superconducting (c, d) susceptibilities to be maximal in the non-intercalated sample, while vertex corrections predict a nearly five-fold enhancement in magnetic instability (e, f) in the intercalated sample.

We compute the RPA particle-particle superconducting susceptibility χ^pp(q), and resolve it by intra-orbital components (Fig. 4c, d). Indeed χ^pp(q) gets weakly reduced on intercalation. This conclusion appears qualitatively slightly different from earlier work⁵⁰, and the main reason for that, we think, is in how the RPA pairing susceptibility was computed. In Ref. ⁵⁰ it was computed from a tight-binding model derived from DFT, while in the present case the eigenfunctions are computed directly from QSGW+DMFT. From this earlier work⁵⁰ it can also be seen that the superconducting eigenvalue λ gets weakly reduced for most of the doping range, and gets weakly enhanced at the maximal electron doping. However, the true intercalated sample in the experiment does not correspond to the maximal doping and it is rather difficult to understand from their figure whether λ for the experimentally intercalated sample should increase or decrease. In any case, such minor differences between these two different RPA calculations are not relevant for the essential fact of a five-fold enhancement in T_c on intercalation.

When χ^m(ω, q) is calculated with the vertex, it shows a systematic enhancement on intercalation (Fig. 4e, f). The static χ^m(1/2, 1/2) in the d_xy channel gets enhanced by a factor of nearly 5, and ~ 3 in the d_yz,xz channels. However, ${\chi }_{xy}^{m}(1/2,1/2)$ always remains at least a factor of two larger than ${\chi }_{yz,xz}^{m}(1/2,1/2)$, suggesting that the low energy spin fluctuations originate primarily in the d_xy channel. Orbital-differentiation is a signature of Hund’s correlations; mass enhancement factors can be quite different for different orbitals^{45,49,57,58,59,71}. The effect on χ is a two-particle analogue of the orbital-differential for the self-energy. The fact that d_xy is the primary source of magnetic fluctuations in a range of strongly correlated chalcogenide and pnictide superconductors was discussed in previous works^45,49,70. To understand the role of the vertex functions in this remarkable enhancement in χ^m(q) we analyse the magnetic vertex functions Γ^ph,m (Fig. 5a, b) and their energy, momentum and orbital dependence. Γ^ph,m(ω₁, ω₂, Ω) is a dynamic quantity that depends on two Matsubara frequency indices (ω_1,2) and one bosonic frequency (Ω). We observe that the five-fold enhancement in χ^m is directly related to the five-fold enhancement in Γ^ph,m in the d_xy channel in the static limit (Ω = 0 and ω₁ = ω₂). Γ^ph,m also gets enhanced on intercalation in the d_yz,xz channels but only by a moderate amount. Also, the magnetic vertex corrections always remain about factor of two larger in the d_xy channel than the d_xz,yz channels, consistent with the magnetic susceptibilities. This shows that the enhancement in magnetic fluctuations at low energies on intercalation is purely a phenomenon emerging from the two-particle electronic vertex, and it is not contained in the bare RPA polarizability, even when computed using Green’s functions G_k,ω dressed with the DMFT self-energy Σ(ω).

Fig. 5: Vertex mediated five-fold enhancement in T_c on intercalation.

The orbital projected magnetic Γ^ph,m (a, b) and pairing Γ^pp (c, d) vertex functions show a factor of ~5 enhancement in the d_xy channel, while a factor of ~3 enhancement in the d_yz,xz channels on intercalation. The enhancement in T_c directly correlates with enhancement in pairing vertex strength in the d_xy channel.

We further analyse the momentum and orbital structure of the pairing vertex Γ^pp at Ω = 0 (after all internal frequencies are integrated). In complete consistency with χ^m(1/2, 1/2) we observe a similar nearly five-fold enhancement in Γ^pp(q = 1/2, 1/2)_xy (see Fig. 5c, d). Weaker enhancements in the d_yz,xz pairing vertex can be observed as well. Intriguingly enough, Γ^pp(q = 1/2, 1/2)_xy gets enhanced in p − d compared to n − i by a factor of ~ 3, and this is purely due to the changes in Fe–Se–Fe bond angle. The reduction in screening from the increased c-axis length that enhances U,J by about 10% accounts for about a third of the total enhancement in T_c from ~ 9 K to ~ 45 K. We also observe that a proper treatment of the Γ^pp reduces the contribution to superconducting instability originating from q = 0 and makes it mostly dominated by q = (1/2,1/2). Further, we compute the superconducting order parameters and find that while the leading instability is of extended s-wave in nature^45,49 in n − i and f − d, it has d${}_{{x}^{2}-{y}^{2}}$ symmetry in p − d. This is another testimony to the importance of reliable computation of vertex and Hubbard correlation parameters in such strongly correlated systems, where moderate changes to these quantities can lead to significant qualitative changes to collective instabilities.

Impact of missing hole pocket on superconductivity

We next address the major limitation of the present ab initio theory, namely its prediction of a hole pocket in the d_xy channel at the Γ point, which is missing in photoemission experiments^72,73, this band is found to fall slightly below E_F. As we have discussed elsewhere⁴⁵, the Fermi surfaces computed from our QSGW+DMFT approach, significantly reduces the hole pocket compared to DFT+DMFT approaches but it does not drive the d_xy state below E_F. The origin of this discrepancy originates from a non-local self-energy, likely magnetic fluctuations or the electron-phonon interaction. Neither of these are yet built into the present theory, but whatever the cause it is important to assess its effect on the conclusions of this paper. The absolute position of the d_xy states at Γ point can be sensitive to doping^74,75, intercalation and other structural changes.

Here we model how the role of proximity of d_xy state to E_F affects superconductivity by adding an external potential to shift to the d_xy state only, leaving the remainder of the system intact. It is similar in spirit to QSGW+U+DMFT where U is applied only to the d_xy orbital. We use U as a free parameter to create a potential shifts δ_xy to d_xy and explore its consequence in superconducting pairing instability by solving the BSE in the pairing channel, exactly as before. Here we show results for four cases, with δ_xy of 0, − 40, − 100 and − 130 meV (see Fig. 6a–d). In the first three cases, (a)–(c), the d_xy hole pocket survives, although it keeps getting smaller with larger δ_xy. We find that the leading eigenvalue λ for superconductivity in s^± channel remains nearly invariant for δ_xy= − 40 and −100 meV where the d_xy hole pocket still survives. When finally δ_xy= − 130 meV, the d_xy hole pocket is pushed below E_F by ~ 25 meV, close to the position observed in photoemission. In that case we find λ is reduced only slightly compared to δ_xy=0. However, note that in all cases in Fig. 6a–d, the d_xy character survives at the electron pockets. Superconductivity is a low-energy phenomena and it is known that for multi-orbital superconductivity, it is important that both the narrow and dispersive orbitals are present close to the Fermi energy and it is not necessary that all have them have to be right at E_F^76,77,78. In some of our older works, we discussed the impact of anomalous screening on properties like, Kondo physics, Mott transitions and superconductivity, even when the van Hove singularities in electronic density of states sit beyond the thermal broadening energy scales from the Fermi surface^79,80. Our fully ab-initio framework establishes that the proximity of the most correlated d_xy state to the Fermi energy, is key to superconductivity⁴⁵. If it gets pushed too far below the Fermi energy, the superconducting instability is suppressed. This is understandable since the paramagnon dispersion in bulk FeSe is ~100 meV and the superconducting gap energy scale is only about 1 meV. Bands that are pushed far beyond these typical relevant energy scales, would have less impact on the pairing mechanism.

Fig. 6: d_xy hole pocket and its impact on superconductivity.

For both non-intercalated (a–d) and intercalated (e–h) the electronic band structure is shown along M–Γ–A path for different values of shift δ_xy in potential. Red represents the d_xy orbital character and green represents d_xz,yz orbital characters. The leading superconducting eigenvalues λ remains roughly 4–5 times larger in intercalated samples for all δ_xy (the reported values for λ are at 300 K, which is the lowest temperature where we could compute λ reliably in all cases). The pairing vertex Γ^pp uniformly enhances by similar factors in intercalated variants for all such δ_xy, as the electron pockets around M and A, always contain significant d_xy component.

We now explore the consequence of similar shifts δ_xy in intercalated (f–d) FeSe. Note that in Fig. 4 and Fig. 5, we have shown that for δ_xy=0, intercalation (f–d) produces a ~ 4-5 factor enhancement is pairing instability compared to the non-intercalated (n–i) bulk FeSe. As we add, δ_xy= − 40, − 100 and − 130 meV (see Fig. 6f–h)), we find that in all cases the relative pairing strengths still enhances ~4–5 on intercalation. This intriguing result supports the essential claim of our work that the nearly five-fold T_c enhancement on intercalation is a robust fact and is not sensitive to the presence or absence of the d_xy hole pocket. This is possible because d_xy orbital character is still present in the electron pockets and the pairing instability mediated by d_xy still gets enhanced by a similar factor on intercalation. A significant challenge for future ab-initio studies will be to explain the physical mechanism that slightly suppresses the d_xy hole pocket. In this context, it is also relevant to note that a self-energy or vertex function simulated within a beyond DMFT mechanism can have impact on both the band energies and the susceptibilities. For the moment, beyond DMFT approaches like DΓA^81,82 are limited by the maximum number of orbitals that can included in the correlated Hamiltonian, nevertheless, such approaches have potential to qualitatively better describe the electronic properties of materials like FeSe.

ARPES, NMR and superconductivity in intercalated FeSe

Experiments offer valuable insights into the origins of superconductivity, and ideally experiments such as ARPES, NMR and Knight shift could provide hard tests that either lend support to the conclusions we have drawn, or be at variance with them. We present here some key experimental findings that broadly support the theoretical findings; however there are enough gaps in both experiment and theory that some caution is needed: we cannot exclude the possibility that a boson we have not considered, (e.g. the electron-phonon interaction) also makes some contribution to superconductivity in intercalated FeSe. For example, an ARPES study on Li-intercalated FeSe⁸³ at inside the superconducting phase at 20K suggests that all the hole pockets are pushed slightly below E_F. This study, as well as valence analysis of Fe^40,41, indicate that the intercalated samples are slightly electron-doped. Even in intrinsic FeSe d_xy state is below E_F^84,85 and d_xz,yz only leads to a very small hole pocket. Some experiments suggest suppression of the hole pocket suppresses superconductivity, in accord with a traditional nesting picture. Two instances of heavily electron-doped systems that do not superconduct are Fe_1.03Se⁸⁶ and Fe_0.9Co_0.1Se. Li-intercalated FeSe⁸³, on the other hand are the counterexamples (alkali doped FeSe^87,88,89 and intercalated FeSe) that establishes a hole pocket crossing E_F at Γ, and therefore a static one-particle nesting picture, is not essential. Our theoretical treatment indicate that provided those hole pockets lie close enough to E_F they can mediate strong pairing even if they do not cross it. A similar argument was made based on a two-band model Hamiltonian⁹⁰, where the authors showed that spin fluctuations can mediate pairing in the strong coupling limit when the electron-like pocket remains at E_F while the hole-like pocket becomes ‘incipient’. This suggests that the traditional static concept of nesting needs to be generalised to a dynamical one, where the frequency-dependence of χ plays a key role.

Other key experiments are NMR and the Knight shift. We compute $\sum {{{\rm{Im}}}}\,{\chi }^{m}({{{\bf{q}}}},\omega )/\omega$, which is the main factor in determining 1/(T₁T) measured by NMR; and also χ^m(q = 0, ω = 0), which is the main factor controlling the Knight shift K_S. Figure 7a shows that $\sum {{{\rm{Im}}}}\,{\chi }^{m}({{{\bf{q}}}},\omega )/\omega$ for f–d and n–i differ widely at room temperature, but the difference shrinks with temperature, to about 20% at 77 K. This is consistent with our interpretation that while at a particular q vector χ^m can get enhanced by a factor of 4–5 in the f–d, the local quantity may not show similar enhancement. Turning to the Knight shift, the theory predicts (see Fig. 7b) almost no difference between f–d and n–i in χ^m(q = 0, ω = 0) at any temperature. Both of these observations are found in experiments as well^86,91. Together, they suggest that our computed magnetic susceptibilities and their momentum and energy structures are of good fidelity and reasonably consistent with NMR data. However, we can not fully interpret the absence of the build up of 1/(T₁T) right above T_c as observed in experiments. We believe, to an extent we understand why that is the case, but a complete understanding is lacking. We discuss this in detail in the next paragraph. As an additional caveat, the QMC solver available to us limits the temperatures we can reach. The lowest temperature for which we could compute vertex corrected susceptibilities was 77 K, somewhat above the critical region around 45 K.

Fig. 7: q-integrated and q = 0 spin susceptibilities for comparison against the NMR $\frac{1}{{T}_{1}T}$ and Knight Shift data.

DMFT vertex corrected $\sum \frac{{{{\rm{Im}}}}{\chi }^{m}({{{\bf{q}}}},\omega )}{\omega }$ and χ(q = 0, ω = 0) are computed in the temperature range 300 K to 77 K. a While $\sum \frac{{{{\rm{Im}}}}{\chi }^{m}({{{\bf{q}}}},\omega )}{\omega }$ remains almost twice large in the intercalated phase at high temperatures compared to the non-intercalated phase, at lower temperatures it is only about 20% larger. b χ^m(q = 0, ω = 0) remains almost invariant over all temperatures between the non-intercalated and intercalated phases.

In bulk FeSe, the primary nesting vector is (1/2,1/2) (in the unit cell that contains two Fe atoms). The primary nesting vectors are different in FeSe and FeSe_1−xTe_x⁹². With Te doping more than one nesting vector emerges. A similar situation occurs in alkali doped FeSe⁹³ where the primary nesting vectors are entirely different^94,95 from bulk undoped FeSe. A new primary nesting vector, often, does co-exist with the old one and also with other additional nesting vectors (several vectors with similar strengths of spin fluctuations). What it primarily suggests is that spectral weights get redistributed in the material over different q and energies, on doping. One can imagine that what the material chooses as the primary nesting vector from many, in doped systems, can depend on multitude of factors and the balance can shift easily (depending on parameters like atomic co-ordinates and changes in hopping parameters). Intriguingly enough, in almost all alkali doped FeSe compounds, the hole pockets appear to be pushed below E_F^87,88,89, much like its intercalated counterpart. This is of extreme relevance to our case, because 1/(T₁T) is a local quantity that sums over dynamical spin susceptibility over all q (in the ω → 0 limit). In complicated cases like these, where there is spectral weight redistribution over various q points, the summed over quantity over a certain temperature window may still appear very similar from different materials (for example, in bulk FeSe, alkali doped FeSe⁹⁶ and intercalated FeSe⁹¹). Also, when we check the magnetic susceptibility χ^m at the vector q = (1/2,1/4,1/2) (unfolded) we do find a strong enhancement in spin fluctuations compared to non-intercalated case. This is the same vector where the primary nesting appears when FeSe is doped with Rb or K. Additionally, from our 1/(T₁T) calculations we also observe that as we include more q-points to compute $\sum {{{\rm{Im}}}}\,{\chi }^{m}({{{\bf{q}}}},\omega )/\omega$, the curves for the n–i and f–d phases come closer at lower temperatures. This is a clear indication for the fact that when spin fluctuations are distributed over several q vectors we need to be more careful about obtaining this quantity. A counter-argument to this could be, APRES suggests that in (alkali doped and) intercalated samples nesting from (1/2,1/2) is removed and data from Neutron scattering suggests that a new nesting vector appears at (1/2,1/4,1/2), then would not it be reasonable to argue that magnetic fluctuations and superconductivity chooses this new vector (1/2,1/4,1/2), instead of (1/2,1/2)? The problem with this argument is that if this new vector nests the Fermi surface and mediate spin fluctuations it still does not answer why the build up in 1/(T₁T) remains absent right above T_c. This supports our observation for the need of more careful interpretation of 1/(T₁T) in cases where multiple (finely balanced) nesting vectors can emerge under doping or intercalation. Taken together with our another observation that superconductivity does not necessarily gain from a static nesting picture, but rather from ‘incipient’ bands, they stress the need of more careful computation and interpretation of 1/(T₁T), since the relevant ‘incipient’ bands for superconductivity are sitting few meV below E_F. (Note that a similar situation occurs in uni-axially strained Sr₂RuO₄, where a new spin susceptibility peak emerges at q = (0.5,0.25,0)⁹⁷ while the peak at incommensurate q = (0.3,0.3,0) from the unstrained material also survives with almost equal intensity and the T_c enhances on strain.)

Soft collective charge excitations and electron-phonon vertex

While our findings are largely consistent with key experimental data, the theory does not take into account the electron-phonon interaction, and we cannot exclude the possibility that it can also play some role in the enhancement of T_c on intercalation. As a hint towards addressing this question, we compute the charge susceptibilities in both QSGW (within RPA approximation) and also with local vertex corrections from DMFT. We find that intercalation causes the real part of the charge susceptibility χ^c to drop significantly within either approximation: Fig. 8a–c shows the vertex corrected χ^c calculated from DMFT. (Note that χ^c in the f–d phase is multiplied a factor of three to bring them to the same scale.) A strong suppression of χ^c in the f–d phase suggests that collective charge excitations themselves can not drive pairing, nevertheless, if χ^c becomes small enough (~0), then χ⁻¹ can diverge. The electron–phonon vertex is linear in χ⁻¹ and can diverge too. This suggests that while the electron-phonon plays little role in intrinsic FeSe, it may play some role in the intercalated case. A recent study explores the role of charge criticality in FeSe_1−xTe_x⁹⁸. It is likely that a charge mechanism, as elucidated above, is at play in doped, monolayer, ultra-thin films⁹⁹ and intercalated variants. A more definitive answer is beyond the scope of our present study.

Fig. 8: Softening of charge susceptibility χ^c on intercalation.

DMFT vertex corrected χ^c is plotted in different intra-orbital channels (a) d_xy (b) d_yz and (c) d_xz. For all the channels real part of χ^c reduces on intercalation. We plot three times the χ^c for the intercalated phase to bring them to the same scale with non-intercalated phase.

To summarise, we perform ab-initio all Green’s function calculation for non-intercalated and intercalated FeSe in the presence of vertex corrections. QSGW supplies a good reference one-body hamiltonian; this with the dynamical local self-energy and vertex from DMFT yields a very good ab-initio description of the spin susceptibility⁴⁹ and charge- and spin-fluctuation- mediated superconductivity. The vertex functions, with their orbital, momentum and energy dependence are directly computed out of the theory and no form factor is assumed. We rigorously establish that the essential component of the superconducting pairing vertex is the magnetic vertex in FeSe. Changes in the electronic density of states cannot explain the enhancements to T_c in these systems. In the absence of the vertex, superconductivity in intercalated materials either would not increase or get weakly suppressed.

Intercalation enhances the superconducting instability by a factor of five, primarily because the magnetic vertex gets enhanced by a similar factor at some particular q vectors in the Fe-d_xy channel. The d_yz,xz channels are also enhanced by a factor of two, but they always remain the secondary source of pairing glue. Such clear orbital differential in two-particle channels is a hallmark of large Hund’s coupling.

Further, we show that incorporating the effects of reduced electronic screening due to enhanced layer-separation post-intercalation is crucial when constructing a realistic many-body Hamiltonian for these intercalated materials and it is the enhancement in pairing vertex driven by such reduced electronic screening that can account for about a third of the total enhancement in T_c. We believe, our work establishes the foundation for tetragonal FeSe where similarly enhanced superconducting T_c is realised on intercalation, alkali-doping, under-pressure, on ionic gating and surface doping. Finally, we address the outstanding problem of missing d_xy hole pocket from Fermi surfaces in bulk FeSe. We show by creating an artificial potential shift to the d_xy state, that even in the extreme case when the d_xy band energy is pushed nearly 100 meV below the Fermi energy, on intercalation, the pairing instability still enhances by nearly a factor of 4. We also show that on intercalation, electronic spectral weight gets redistributed over various q vectors and in such cases, 1/(T₁T) may not enhance significantly above T_c compared to the non-intercalated variant, in reasonable agreement with NMR measurements⁹¹. Finally, we show that while q-selective enhancements in the pairing vertex are closely connected to enhancements in the magnetic susceptibility and concomitant superconducting instability, intercalation also induces a significant softening in collective charge excitations. This raises the possibility that electron-phonon coupling may also contribute to superconductivity in intercalated FeSe.

Methods

One-particle calculations using LDA and QSGW

Single particle calculations (LDA, and energy band calculations with the static quasiparticlized QSGW self-energy Σ⁰(k)) were performed on a 16 × 16 × 16 k-mesh while the (relatively smooth) dynamical self-energy Σ(k) was constructed using a 8 × 8 × 8 k-mesh and Σ⁰(k) is extracted from it. The charge density was made self-consistent through iteration in the QSGW self-consistency cycle: it was iterated until the root mean square change in Σ⁰ reached 10⁻⁵ Ry. Thus the calculation was self-consistent in both Σ⁰(k) and the density. At the end of QSGW cycles, we use the quasi-particlised electronic band structures as the starting point of our DMFT calculations.

One-particle calculations using DMFT

The impurity Hamiltonian is solved with continuous time Quantum Monte Carlo solver^100,101. For projectors onto the Fe d subspace, we used projectors onto augmentation spheres, following the method described in this reference¹⁰². This approach is sightly different from the approach used to compute U,J parameters in a previous work by Miyake et al.¹⁰³. Further, those numbers¹⁰³ are computed while building the Hubbard Hamiltonian on top a DFT bath, while ours is a QSGW bath. QSGW already takes into long-range charge correlations missing from DFT, so it is only natural that the correlations (mostly of spin fluctuations origin) that our QSGW calculations miss out would be lesser compared to DFT, leading to smaller U,J estimations. The double counting correlations are implemented using fully localised limit approximation. The DMFT for the dynamical self energy is iterated, and converges in 30 iterations. Calculations for the single particle response functions are performed with 10⁹ QMC steps per core and the statistics is averaged over 128 cores.

Two-particle calculations using DMFT

The two particle Green’s functions are sampled over a larger number of cores (40,000–50,000) to improve the statistical error bars. The local effective interactions for the correlated impurity Hamiltonian are given by U and J. These are calculated within the constrained RPA¹⁰⁴ from the QSGW Hamiltonian using an approach⁶⁹ similar to that of ref. ¹⁰⁵, using projectors from ref. ¹⁰².