Origin of pressure induced second superconducting dome in A_y Fe_2−xSe₂ [A = K, (Tl,Rb)]

High-T_c superconductivity often occurs when the system is driven from its pristine phase to the verge of a magnetic quantum critical point via external parameters, such as chemical doping, magnetic field or pressure (P) in most of the cuprates, heavy fermions, pnictides and organic superconductors. However, several recent breakthrough discoveries of a second superconducting (SC) dome—completely isolated or slightly connected to the first SC dome—without the intervention of any competing order as a function of P in several families of high-T_c superconductors [1–3], and/or extreme chemical doping in KFe₂As₂ [4], LaFeAsO₁₋₂H_x [5] or strain [6] have questioned this expectation. A much higher optimal T_c value than that of the first dome, as well as a substantial increase in electronic mass with P, as obtained in iron-based compounds, suggest an interesting and exotic phenomena of superconductivity along this tuning axis. It is notable that in stochiometric SrFe₂As₂, a SC dome appears both as a function of P [7], as well as crystallographic strain [6] near the quantum critical point of the spin-density wave, as often observed in other high-T_c superconductors, further supporting the notion that the P induced superconductivity is unconventional. In this paper, we present a model that allows us to capture the onset of the second SC phase with different pairing symmetry as a function of P.

Our approach is based on the weak or intermediate coupling scenario in which the shape of the Fermi surface (FS) topology plays a key role in creating pairing instability at the 'hot-spot' [8]. For such a case, Cooper pairing arises from repulsive interaction with sign-reversal pairing symmetry constrained by the FS topology and crystal symmetries. This theoretical framework consistently describes d-wave pairing in cuprates [8], and Ce-based heavy fermions [9], s^±-pairing in iron–pnictides and chalcogenides [10] and nodeless d-wave in A_yFe_2−xSe₂ [A = K, Cs, Rb, (Tl,Rb), (Tl,K)] families [11–13].

Here, we focus on the latter family in which a second isolated dome is observed, as shown in figure 1(a) (reproduced from [14]). In an analogous 122 compound BaFe₂As₂, it is established that the effective mass, m*, gradually increases as a function of P, see the plot of 1/m* in figure 1(b). Constrained by these experimental facts we postulate that uniform P renormalizes the bands in a way that goes beyond standard Fermi-liquid (FL) behavior, as also shown in first-principles calculations [15], and thereby affects the FS topology. In A_yFe_2−xSe₂ systems, both first principles calculations and angle-resolved photoemission spectroscopy (ARPES) have demonstrated that a hole-pocket lies slightly below the Fermi level (E_F) at ambient P. With band renormalization this hole pocket appears on the FS above a critical P, and the overall FS topology changes from only electron-pockets at P = 0 to the coexisting electron and hole pocket as in pnictides. This topological FS transition induces a crossover from FS nesting along the diagonal Q₁ → (π,π) direction to along the zone direction Q₂ → (π,0)/(0,π), which makes a pairing symmetry transition from nodeless d-wave to s^±-pairing.

Zoom In Zoom Out Reset image size

We explore these postulates via first-principles band structure and pairing eigenvalue calculations within a random-phase approximation (RPA). The first-principles calculation is performed for KFe₂Se₂ within the atomic sphere approximation using the TB-LMTO code [16, 17]. For a given P, the uniform volume contraction is evaluated using the Birch–Murnaghan equation of state [18] formalism (given below). The key here is to optimize the Se atomic position with respect to the Fe plane by minimizing total energy at each P. We find z_Se = 0.3452c at P = 0 GPa and z_Se = 0.351c at P = 12 GPa. The obtained band structure in figure 1(c) shows that at P = 12 GPa, the bands are renormalized by Z ∼ 0.9 (see inset), which allows the hole-pockets to appear on the FS. The change of FS via band renormalization can be understood this way. As deduced below, the pressure-induced renormalization affects the hopping energies, not the chemical potential. Thus at each pressure a new Fermi energy appears which allows the shifting of the band along the energy direction. This result is unexpected from the standard FL picture. Therefore the dominant FS instability arises along Q₂ at a critical pressure, and thus s^± pairing with large coupling constant commences. We also semi-quantitatively reproduce the phase diagram of the d and s^±-pairings within the conditions for best nestings along each channels.

To grasp further insight into how P modifies the electronic structure, and also to enable adding correlations for pairing symmetry calculations, we use a low-energy five-bands tight-binding (TB) formalism from [12] at zero pressure, and include pressure effects by band renormalizations. The TB hopping integrals t_i in a given crystal is defined as t = 〈Ψ|V_c|Ψ〉, V_c is the Coulomb interaction between lattice and electron. Since the TB hopping involves integration over the unit-cell volume, the simplest approximation to account for the change in TB parameters due to the change in unit cell volume is t_i(P)/t_i(0)∝V (P)/V₀, where t_i(P) and V (P) are the TB parameters and lattice volume at any given P P, and t_i(0) and V₀ are their corresponding values at ambient P. This effective theory, which is reasonably justified by the first-principles calculations, is valid when changes in the wavefunction and Coulomb potential V_c as a function of pressure are negligible. For simplicity, we take t_i(P)/t_i(0) ∼ V (P)/V₀ ≈ m/m* = Z, for all bands. From the value of Z, or more strictly from the volume ratio, we can obtain the value of P using the Birch–Murnaghan equation of state [18]

$$\begin{equation} P=\frac{3B}{2}\left( x^{7/3}-x^{5/3}\right)\left[ 1-\frac{3}{4}(4-B^{\prime})(x^{2/3}-1) \right], \end{equation} \tag{ 1 }$$

where x = V₀/V =1/Z, and B, B' are the bulk modulus and its derivative with respect to pressure. By fitting the calculated optimal T_c of the second dome to the experimental value of  ∼ 48 K, we find B ≈ B'P = 67 GPa, which is close to the first-principle value of B = 64–70 GPa for this sample, and also to the available experimental data of B = 62 GPa for other iron–pnictides [19]⁵.

In what follows, the renormalization, applied to the TB parameters, reflects in the band structure as Eⁱ_k = Zξⁱ_k − E_F, where ξⁱ_k is the ith TB band taken from [12]. Unlike in FL-theory, Z does not renormalize the spectral weight, and thus at each pressure a new Fermi level occurs to keep the number of electron unchanged. This constraint allows the change in FS upon renormalization of the bands. We compute the Fermi level E_F self-consistently, beyond a simple rigid band shift approximation, by integrating the density of state up to E_F. This electron number can be related to the FS volume in a slightly revised Luttinger theory as

$$\begin{equation} n=\frac{1}{N}\sum_{{\boldsymbol{k}},i=1-5}\int \frac{{\rm d}\omega}{\pi} \delta(\omega-ZE_{\boldsymbol{k}}^i) = \frac{1}{N}\sum_{{\boldsymbol{k}},i=1-5}\frac{\delta(E_{\boldsymbol{k}}^i)}{Z}=\frac{V_L(0)}{Z}, \end{equation} \tag{ 2 }$$

where V_L(0) is the unitless FS volume at P = 0. This formula can be contrasted with the conventional Luttinger theorem, which says that the FS volume, after applying a FL-like renormalization, does not change, because $n=\frac {1}{N}\sum _{{\boldsymbol{k}},i=1-5}\int \frac {{\mathrm { d}}\omega }{\pi } Z\delta (\omega -ZE_{\boldsymbol{k}}^i) = \frac {1}{N}\sum _{{\boldsymbol{k}},i=1-5}Z\frac {\delta (E_{\boldsymbol{k}}^i)}{Z}=V_{\mathrm { L}}(0)$, where V_L should be read as bare FS volume. The important difference between the present formalism and FL thus comes from the fact that while in FL Z renormalizes both spectral weight and dispersion, in the present case it only renormalizes the band, not the spectral weight. This allows the change in FS topology and volume as a function of pressure, while V_L/Z remains constant, constrained by the number of electron.

Figure 2 reveals the evolution of electronic states at ambient P and at P = 12 GPa. At P = 0, the hole pocket lies at 60 meV below E_F at Γ-point, as is also seen in ARPES data [20]. In this case, the presence of electron pockets at X-points (in one Fe per unit cell notation) is well established [20–22]. With the uniform band renormalization by Z = 0.9 at P = 12 GPa, according to equation (1), we find that two concentric hole pockets are fully formed on the FS, as seen from figures 2(b1) and (b2).

Zoom In Zoom Out Reset image size

To determine the dominant nesting 'hot-spot', we calculate multiband susceptibility using a full orbital overlap matrix-element. The many-body correction is incorporated within RPA by explicitly including intraorbital interaction U, interorbital interaction V =U − 2J, Hund's coupling J = U/4 and the pair hopping energy J' = J as defined in footnote⁶ , and the details of the calculation can, for example, be found in [11, 13, 23, 24]. We present results as a function of U as given in figure 4(a), where the other parameters change accordingly following these equations. We plot the static spin RPA susceptibility (trace of the susceptibility tensor) at the two representative P in figure 2(a3) for P = 0 and in figure 2(b3) for P = 12 GPa. As expected, the dominant nesting in the former case is aligned along the inter-electron-pocket direction (consistent with earlier calculation in [11–13]). On the other hand, at P = 12 GPa, the dominant nesting changes to Q₂ → (π,0) between electron and hole pockets, as obtained for many iron–pnictide superconductors [10, 23].

Next we evaluate the effective pairing vertex in the singlet channel for scattering between two FSs i,j within the spin and charge fluctuations exchange approximation

$$\begin{equation}\fl {\hat \Gamma}({\boldsymbol{k}},{\boldsymbol{k}}^{\prime})=\frac{1}{2}\left[3{\hat U}_{\rm s}{\hat \chi}_{\rm s}({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime}){\hat U}_{\rm s} - {\hat U}_{\rm c}{\hat \chi}_{\rm c}({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime}){\hat U}_{\rm c}+ \frac{1}{2}({\hat U}_{\rm s}+{\hat U}_{\rm c}){\hat \chi}_0({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime})({\hat U}_{\rm s}+{\hat U}_{\rm c}) \right]. \end{equation} \tag{ 3 }$$

It should be noted that the above equation is similar to the one used in [23–27], with differences only being in the bare bubble and onsite Coulomb repulsion terms. However, since the dominant contributions to the pairing channel come from the RPA term, the absence or presence of the other comparatively weaker terms does not alter the resulting pairing symmetry, except giving an overall shift in the value of the pairing strength. We now study the evolution of the pairing strength λ for a given gap function g(k) from the following equation [23]:

$$\begin{equation} \lambda_{ij} [g] = - \frac{\oint_{c_i} \frac{{\rm d} k}{v_{\rm F} (k)} \oint_{c_j} \frac{{\rm d} k^{\prime}}{v_{\rm F} (k^{\prime})} g({\boldsymbol{k}}){\rm Re}\big[\Gamma_{ij}({\boldsymbol{k}},{\boldsymbol{k}}^{\prime})g({\boldsymbol{k}}^{\prime})\big]}{(2\pi)^2\oint_{c_j} \frac{{\rm d} k}{v_{\rm F} (k)} [g({\boldsymbol{k}})]^2}. \end{equation} \tag{ 4 }$$

Here the line integrals over C_i are performed on each FS loops, and v_F is the Fermi velocity. The key point of equation (4) is if the gap function g(k) possesses an opposite sign at k and k', mediated by a large peak in Γ at the 'hot-spot' q = k − k'. Based on this framework, we now study the evolution of λ for g(k) = cos k_x − cos k_y for d_x²−y² and g(k) = cos k_x cos k_y for s^±-pairing channels (in the one Fe per unit cell notation) as a function of P in figures 3 and 4 (in two Fe per unit cell, the form of the pairing structure g(k) transforms by the same unitary transformation as the FS so that the macroscopic properties such as nodeless and isotropic gap structure remains the same in any notation [13, 28]). To understand in detail the origin of the pairing symmetry transition, we study each intra- and inter-band component of λ_ij in figure 3 at two representative P. These results confirm our initial assumption.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For d_x²−y²-wave pairing at P = 0 in figure 3(a) and at P = 12 GPa in figure 3(b), we obtain λ > 0 and <0, respectively. This can be understood from the corresponding FS topological changes. At P = 0, the two electron pockets at (π,0) and (0,π) posses an opposite sign of SC gap, and this phase is supported by nesting along Q₁ → (π,π), as deduced in figure 2(a3) [11–13]. So we get λ_γ−γ > 0. But as the dominant 'hot-spot' changes to Q₂ → (0,π) at P = 12 GPa, the positive λ_γ−γ component is overturned by large negative λ_α−γ and λ_β−γ components being supported by Q₂ 'hot-spot' without sign-reversal of g(k). Therefore, the total λ becomes negative, making the d_x²−y² pairing unstable at this P. For s^±-pairing, the situation is reversed in that λ_γ−γ < 0 at all values of P, but λ_α/β−γ > 0. It is important to note that due to the presence of the two hole pockets α and β, the total value of λ becomes positive and large at some critical value of P, otherwise, s^±-pairing would have been favorable.

In figure 4, we show the full P and interaction U dependence of the total pairing strength λ, and the corresponding calculated T_c. The d-wave pairing survives up to P ∼ 5 GPa, slightly less than the termination point of the first SC shown in figure 1(a). Of course, these quantitative consistencies rely strongly on the exact shape of FS topology and the value of bulk modulus B used in equation (1). We focus on fitting B so that the optimum T_c for the second dome matchs with the experimental value of P [3, 29]. With a separation of about 1–2 GPa, we find that the s^±-pairing channel appears abruptly for a large range of U. Although we obtain an optimum λ as a function of P at which the FS nesting between α/β to γ is strong for all values of U considered, it survives to a larger P range than the experimental data of K_0.8Fe_1.7Se₂ [3]. However, for other samples within the same family, new data shows that T_c in the second dome is very much P independent and survives up to P as large as 40 GPa measured so far [29].

For a purely electronic mechanism, we compute T_c using the spin-fluctuation exchange formula within the weak-coupling limit [30–34]

$$\begin{equation} T_{\rm c}=\frac{\omega_{\rm sf}}{1.2}\exp{(-1.04/\lambda)} \end{equation} \tag{ 5 }$$

with ω_sf = 55 K, we obtain optimum T_c ∼ 38 and 45 K for the first and second dome, respectively for U = 1.2 eV, which are close to the experimental values of 37 and 48 K. It is obvious that the absolute value of T_c depends on the interaction parameter, however, the ratio between the optimum value at two domes is always maintained. Interestingly, we find that the computed T_c for the second dome is even flatter than that for λ and thus agrees well with the new data [29].

Relating the glue function to the spin-resonance mode that appears in the SC state, we obtain a resonance ω_res ∼ 15 meV at the optimum P. The phenomena of spin-resonance in the SC state is well known [8–11, 23] and is essentially similar to the condition for obtaining a positive pairing eigenvalue λ as discussed above: given that the sign reversal of the SC gap is connected by momentum transfer q, a spin-resonance appears at an energy ω(q) ∼ |Δ(k)| + |Δ(k + q)| (RPA correction shifts the mode to a slightly lower energy). We compute the spin-resonance spectra at the same representative P values as in figures 2 and 3 for d-wave and s^±-pairing and the results are shown in figure 5. We immediately see that for d-wave pairing we get a nearly commensurate spin-resonance mode at ω/2Δ ∼ 0.8, while it becomes incommensurate and shifts to higher energy for the s^±-pairing at higher P. Although, ARPES and INS measurements are difficult to perform under P, a neutron diffraction experiment can be done here to test the changes of Q vector as a function of P to verify our results.

Zoom In Zoom Out Reset image size

The FS topology modeled here is constrained by the consistency between band-structure calculations and ARPES measurements, and thus we expect that correlation effects driving either coexistence with vacancy order and/or magnetic phase [28, 35], or phase separation between them [36] will not dramatically change our results, but more studies are needed to address this question. Furthermore, we also note that a recent ARPES measurement [37] has shown that a tiny hole pocket develops around a Z-point at ambient P with an isotropic SC gap. In our present scenario of d-wave pairing, one would expect a node on this FS. However, for such tiny FS elevated along the k_z direction and small nodal quasiparticle weight, one can expect the node to become filled due to various extraneous effects such as disorder, final state scattering of ARPES measurement. A recent experimental study on LaFeAsO₁₋₂H_x found two SC dome with a dip in T_c as a function of doping x [5]. Based on first-principle band-structure calculations, it has been argued that the pairing changes from s^± to s⁺⁺ in that case due to orbital degeneracy.

In conclusion, we present an analysis of pressure evolution of pairing interaction in a A_yFe_2−xSe₂ family of superconductors. We argue that observed two SC domes can be naturally explained by changes in FS topology driven by mass enhancements. At lower P, we expect no Γ-point hole pocket on the FS, and therefore d-wave SC state [37]. Upon increase in P, we expect two hole pockets to develop at the Γ-point, leading to pairing symmetry transition to s^±-pairing. To test these predictions, we suggest the use of magnetic field dependent tunneling spectroscopies to investigate gap changes [38]. Neutron scattering and neutron diffraction measurements would reveal the distinct pattern of spin resonance and nesting properties in these phases.

This work was supported, in part, by UCOP-TR, Nordita and Los Alamos National Laboratory, of the US Department of Energy under contract DE-AC52-06NA25396, and benefited from the allocation of supercomputer time at NERSC.

The spin-fluctuation mediated pairing interaction is well studied in the literature [8, 23–27], which sometimes include the bare bubble term in equation (3) and/or the onsite interaction term (not included here). The spin-fluctuation spectrum in the unconventional superconductors (SCs) obtains a sharp peak at a resonance energy, ω_res, and it then falls off sharply on the energy scale (see figure 5 of main text). Therefore, although the SC gap equation depends on the full vertex, Im  $[{\hat \Gamma }({\boldsymbol{k}},{\boldsymbol{k}}^{\prime },\omega )]$, the relevant k and k' values are restricted by this energy cutoff to remain in the vicinity of the FSs. In this paper, just as for the electron–phonon case, the strength of the pairing interaction is characterized by an energy integral over ${\hat \Gamma }$ weighted by ω⁻¹ following the Kramers–Kronig relation as

$$\begin{equation} \int_0^{\infty}\frac{{\rm Im}[{\hat \Gamma}({\boldsymbol{k}},{\boldsymbol{k}}^{\prime},\omega)]}{\pi\omega}{\rm d}\omega={\rm Re}[{\hat \Gamma}({\boldsymbol{k}},{\boldsymbol{k}}^{\prime},\omega=0)]. \end{equation} \tag{ A.1 }$$

This allows us to evaluate the pairing strength by considering only the real part of the static pairing interaction. This leads to solving the pairing eigenvalue problem by integrating over a closed FS for a given pairing symmetry, say g(k), as

$$\begin{equation} -\sum_i\oint_{C_i}\frac{{\rm d}k}{2\pi}\frac{{\rm Re}[\Gamma_{ij}({\boldsymbol{k}},{\boldsymbol{k}}^{\prime})]}{2\pi v_{\rm F}}g({\boldsymbol{k}})=\lambda g({\boldsymbol{k}}^{\prime}), \end{equation} \tag{ A.2 }$$

where v_F is the Fermi velocity and C_i gives the closed FS for the ith-band. If the gap function can be decomposed by its amplitude Δ₀ and structure factor g(k), then the above eigenvalue problem can be reduced to a dimensionless pairing strength functional [8] given in equation (4). The total pairing strength is then obtained by summing over all bands.

Combining equations (4) and (A.1), we can recast the eigenvalue problem in the typical Eliasberg framework as

$$\begin{equation} \lambda[g] = \frac{2}{\pi}\int_0^{\infty}\frac{\alpha^2F(\omega)[g]}{\omega}{\rm d}\omega, \end{equation} \tag{ A.3 }$$

where the electron–boson spectral function α²F is nothing but a momentum average over the dynamical pairing interaction weighted by the gap function g as

$$\begin{equation} \alpha^2F(\omega)[g] = - \frac{\oint_{c_i} \frac{{\rm d} k}{v_{\rm F} (k)} \oint_{c_j} \frac{{\rm d} k^{\prime}}{v_{\rm F} (k^{\prime})} g({\boldsymbol{k}}){\rm Im}\big[\Gamma_{ij}({\boldsymbol{k}},{\boldsymbol{k}}^{\prime},\omega)\big]g({\boldsymbol{k}}^{\prime})}{(2\pi)^2\oint_{c_j} \frac{{\rm d} k}{v_{\rm F} (k)} [g({\boldsymbol{k}})]^2}. \end{equation} \tag{ A.4 }$$

For any general electron–boson interaction (including phonon [31] and spin-fluctuation [30, 32–34]), SC transition temperature, T_c can be calculated from the pairing strength λ in the weak-coupling scenario as

$$\begin{equation} T_{\rm c}=\frac{\omega_{\rm sf}}{1.2}\exp{(-1.04/\lambda)}. \end{equation} \tag{ A.5 }$$

Here the spin-fluctuation cutoff frequency ω_sf is given by [31]

$$\begin{equation} \omega_{\rm sf} = \exp{\left(\frac{2}{\lambda}\int_0^{\infty}\log{\omega}\frac{\alpha^2F(\omega)}{\omega}{\rm d}\omega\right)}. \end{equation} \tag{ A.6 }$$

Finally, as discussed before, we assume that the α²F has a strong peak at ω_res, and falls off rapidly away from this energy. Using equation (A.3), we obtain ω_sf ≈ ω_res. In figure 4 of the main text, we evaluate T_c by using equations (A.5) and (A.6) with the coupling constant evaluated from equation (4).