Theory of potential impurity scattering in pressurized superconducting La3Ni2O7

Disordered bilayer model

We consider a bilayer Hamiltonian with point-like non-magnetic impurities \({{{\mathcal{H}}}}={{{{\mathcal{H}}}}}_{0}+{{{{\mathcal{H}}}}}_{{{{\rm{int}}}}}+{{{{\mathcal{H}}}}}_{{{{\rm{dis}}}}}\), where the interaction term gives rise to an instability towards unconventional s_±– or d-wave superconductivity, respectively. The non-interacting part is given by

where \({c}_{l,o}^{{{\dagger}} }({{{\bf{k}}}})\) creates an electron in layer l and orbital o with momentum k. In the clean case, there is reflection symmetry between the two layers of a bilayer sandwich, such that the non-interacting Hamiltonian can be written in terms of intra-layer \({\hat{H}}_{0}^{\parallel }\) and interlayer blocks \({\hat{H}}_{0}^{\perp }\)

where the hats refer to the matrices in orbital space. In the following, we drop explicit momentum and frequency dependence and exchange subscripts and superscripts whenever it is convenient for readability. The bilayer phase factors \({e}^{\pm i{k}_{z}d}\) arise from the Fourier transform for a bilayer sandwich with thickness d. The presence of reflection symmetry between the two layers allows to make a transformation into bonding-antibonding space (ba-space).In particular, assuming negligible interbilayer hopping, the non-interacting Hamiltonian is diagonalized by the simple unitary transformations

Here, the bilayer phase factors have been factored out and are no longer present in the ba-space. Within this framework, we thus have a two-dimensional non-interacting Hamiltonian in the ba-space. One should mention, however, that the k_z dependence is implicitly included as the physical observables like charge and spin response functions can be classified into even and odd ones with respect to the k_z dependence, see for example, ref. ⁴⁶ for further details.

We next introduce superconductivity on a mean-field level. After applying mean-field approximation, the interaction part reads

and we can again introduce the ba-block structure

which is again diagonalized using the ba-transformation \(\hat{V}\). Therefore, the superconducting gaps are also two-dimensional in ba-space and given by \({\hat{\Delta }}^{{{{\rm{b/a}}}}}={\hat{\Delta }}^{\parallel }\pm {\hat{\Delta }}^{\perp }\).

To simulate the effect of the electron irradiation, we consider randomly distributed point-like impurities at Ni positions in both NiO₂ layers. Such impurities locally break the reflection symmetry and consequently mix bonding and antibonding blocks. Impurity averaging will re-introduce the ba-block structure, but one has to be careful to perform the averaging in both layers separately⁵⁴. The impurity matrices for the upper and lower layers read

where \({\hat{\tau }}_{i}\) denotes the ith Pauli matrix in the Gor’kov-Nambu space. We consider the self-energy arising due to impurity scattering in the non-crossing approximation and assume the same impurity concentration in both layers. For simplicity, we focus on the intraorbital scattering, which preserves C₄ symmetry, i.e. \({\hat{W}}_{{{{\bf{k}}}}o,{{{{\bf{k}}}}}^{{\prime} }{o}^{{\prime} }}=W\hat{{\mathbb{1}}}\). After impurity averaging, we obtain the self-energy in the form

Importantly, the above expression corresponds to an averaging over bonding and antibonding blocks. To see this, the relation \(2{\hat{G}}_{{{{\bf{k}}}}}^{\parallel }={\hat{G}}_{{{{\bf{k}}}}}^{b}+{\hat{G}}_{{{{\bf{k}}}}}^{a}\) for the Green’s functions can be inserted, which can again be shown by applying the ba-transformation. This ba-averaging is the direct consequence of the local breaking of reflection symmetry. Note that an impurity matrix of the form \({\hat{W}}_{1}+{\hat{W}}_{2}\) is not creating such an averaging⁵⁴. The Green’s functions are self-consistently calculated via

Note, the self-energy is unity in the layer space and, therefore, transforms trivially to ba-space, which results from impurity averaging restoring the global reflection symmetry.

Single-orbital model

To study the effect of the averaging over ba-space on the T_c-suppression, it is instructive to consider first a simple bilayer model with only one orbital, considered e.g. in ref. ⁵⁶. To be specific, we employ a t − J– like bilayer model with in-plane nearest neighbor hopping t and interlayer hopping t_⊥

Similarly, we consider the superexchange in-plane interaction J and inter-layer interaction J_⊥

The summation brackets 〈i, j〉 indicate summation over nearest neighbors only and r is the layer index. Unless stated otherwise, we assume \({J}_{\perp }={({t}_{\perp }/t)}^{2}J\), which is the expected relation once the t − J model is derived from a Hubbard-type model in the strong-coupling limit. The non-interacting parts yields ba-dispersions \({\epsilon }_{{{{\bf{k}}}}}^{{{{\rm{b/a}}}}}=-2t(\cos ({k}_{x})+\cos ({k}_{y}))\pm {t}_{\perp }-\mu\). We perform a mean-field decoupling from which we obtain the following gap equations (neglecting possible in-plane triplet part)⁵⁷:

with \({\gamma }_{{{{\rm{d/s}}}}}=(\cos ({k}_{x})\mp \cos ({k}_{y}))/2\) and Δ_b/a = γ_sΔ_s + γ_dΔ_d ± Δ_⊥ and V = − 3J and V_⊥ = − (3/8)J_⊥. Here, Δ_s/d denotes the in-plane extended s-wave and \({d}_{{x}^{2}-{y}^{2}}\)-wave gap functions, respectively, whereas Δ_⊥ is the interlayer s-wave gap and in the case of opposite signs between bonding and antibonding bands refers to the bonding-antibonding s_±-wave solution. The d_xy-wave solution would have a \(\sin {k}_{x}\sin {k}_{y}\) functional form in the plane. We do not consider it specifically as it behaves similarly to the \({d}_{{x}^{2}-{y}^{2}}-\) wave solution in the presence of non-magnetic point-like impurities. In addition, the \({d}_{{x}^{2}-{y}^{2}}-\) wave solution can also have an interlayer component \(\pm {\Delta }_{\perp }^{{{{\rm{d}}}}}{\gamma }_{{{{\rm{d}}}}}\) and yield either in-phase or out-of-phase locking between two neighboring planes. At the same time, its magnitude is expected to be much smaller than that in the s-wave case, since the d-wave solution is only dominant for the smaller values of the interlayer hopping. For simplicity, we assume it to be small and consider the in-phase locking of the d-wave solution between the layers.

In the clean case, the \({d}_{{x}^{2}-{y}^{2}}\)-wave solution for the given model is found at half-filling for t_⊥ ≲ 1.1t. For larger t_⊥, the interlayer ba-s_± solution becomes dominant at half-filling. Importantly, as Δ_s and Δ_⊥ correspond to gap structures of A_1g symmetry, they always occur simultaneously. In the following, we will focus on the pure interlayer ba-s_± solution, and the mixed solution is discussed in Supplementary Note 1. In this one-orbital model, the bonding and antibonding Green’s functions in Eq. (8) are given by

and the renormalized quantities \(i{\tilde{\omega }}_{n}=i{\omega }_{n}-{\Sigma }_{0}\), \({\tilde{\epsilon }}_{{{{\bf{k}}}}}^{{{{\rm{b/a}}}}}={\epsilon }_{{{{\bf{k}}}}}^{{{{\rm{b/a}}}}}+{\Sigma }_{3}\) and \({\tilde{\Delta }}_{{{{\bf{k}}}}}^{{{{\rm{b/a}}}}}={\Delta }_{{{{\bf{k}}}}}^{{{{\rm{b/a}}}}}+{\Sigma }_{1}\) are defined as

Recall, that for the \({d}_{{x}^{2}-{y}^{2}}\)-wave gap function, Σ₁ vanishes by symmetry in momentum space as can be readily seen from inserting \({\Delta }_{{{{\bf{k}}}}}^{{{{\rm{b/a}}}}}={\gamma }_{{{{\rm{d}}}}}{\Delta }_{{{{\rm{d}}}}}\) in Eq. (16). The interlayer ba-s_±-wave solution with Δ_b/a = ± Δ_⊥ has no dependence on momentum k. Eq. (15) describes the renormalization of the chemical potential. It is included for completeness, however, it has almost negligible effect on our numerical calculations. This is expected for the case of a superconducting BCS-like state. Note also that Σ₃ vanishes exactly in case of particle-hole symmetry. As we are interested in the evolution of the superconducting transition temperature, T_c, when the superconducting order parameter approaches the limit Δ → 0, we rewrite Eq. (16) as

and obtain two equations \({\tilde{\Delta }}_{{{{\rm{b/a}}}}}={\Delta }_{{{{\rm{b/a}}}}}+{\tilde{\Delta }}_{{{{\rm{b}}}}}{\Omega }_{{{{\rm{b}}}}}+{\tilde{\Delta }}_{{{{\rm{a}}}}}{\Omega }_{{{{\rm{a}}}}}\). A similar set of equations was discussed in the context of iron-based superconductors^49,51 with one important difference. In particular, here the inter-band coefficients are not symmetric because in general Ω_a ≠ Ω_b. Solving for \({\tilde{\Delta }}_{{{{\rm{b/a}}}}}\), we find

and after inserting Δ_b/a = ± Δ_⊥ explicitly, we obtain

As follows from Eq. (19), the T_c suppression for the ba-s_±-wave depends on the difference between Ω_b and Ω_a. For the specific case Ω_b = Ω_a, we find \({\tilde{\Delta }}_{{{{\rm{b/a}}}}}=\pm {\Delta }_{\perp }\) and this version of the ba-s_±-wave solution will be suppressed by potential point-like impurities in the same way as the d-wave superconducting state with Σ₁ = 0. However, if there is an imbalance (in this particular case it is particle-hole asymmetry) between bonding and antibonding sub-spaces, there will be  deviation from the AG behavior. Moreover, finite Ω_b ≠ Ω_a induces different magnitudes for the bonding and superconducting gap functions \({\tilde{\Delta }}_{{{{\rm{b/a}}}}}\) as the second term on the r.h.s. of Eq. (19) is of different sign. In particular, for increasing impurity density, Ω_b + Ω_a approaches unity and \({\tilde{\Delta }}_{{{{\rm{b/a}}}}}\) must change its sign, if Ω_a/b > Ω_b/a, respectively. Therefore, a s_±-wave to s₊₊-wave crossover happens, which already was discussed previously in the context of iron-based superconductors^49,51. As Ω_b + Ω_a approaches unity, one might be worried about a divergence due to the denominator. However, this is of course not the case because of the concurrent suppression of Δ_⊥. Note that this can be understood with the renormalized Matsubara frequency \({\tilde{\omega }}_{n}={\omega }_{n}/(1-{\Omega }_{{{{\rm{b}}}}}-{\Omega }_{{{{\rm{a}}}}})\) entering the gap equation Eq. (12). The renormalization of the quasiparticle energies is typically negligible. Note that the numerator in the correction term in Eq. (19) depends linearly on the impurity density. Consequently, the deviation from AG behavior is small for low impurity densities but sizable for higher impurity densities.

For our simple model, we employ \({\epsilon }_{{{{\bf{k}}}}}^{{{{\rm{b/a}}}}}=-2t(\cos ({k}_{x})+ \cos ({k}_{y}))\pm {t}_{\perp }-\mu\). At half-filling, μ = 0, and the relation \({\epsilon }_{{{{\bf{k}}}}}^{{{{\rm{b}}}}}=-{\epsilon }_{{{{\bf{k}}}}+(\pi ,\pi )}^{{{{\rm{a}}}}}\) guarantees Ω_b = Ω_a (note that for this case particle-hole symmetry is present). Therefore, moving away from half-filling allows us to tune the particle-hole asymmetry between bonding and antibonding bands and, consequently, between Ω_b and Ω_a. In Fig. 1 the T_c suppression for exemplary band fillings n = 1, n = 0.75, and n = 0.5 are compared for the pure interlayer ba-s_±-wave and \({d}_{{x}^{2}-{y}^{2}}\)-wave Cooper-pairing scenarios. Clearly, the qualitative behavior is changing as the imbalance between bonding and antibonding bands increases. At half-filling, both scenarios give the same qualitative AG shape (Fig. 1a). Upon lowering the filling to n = 0.75, the T_c curve for the interlayer s_±-wave starts to deviate from AG form close to the critical impurity density at which superconductivity is completely suppressed (Fig. 1(b)). At quarter-filling, n = 0.5, the initial slope is still comparable for both scenarios but the T_c curve for the interlayer ba-s_±-wave state changes from convex to concave shape and persists up to more than three times larger impurity densities. We adopt J and t_⊥ such that we have roughly the same T_c in the clean state. The results are not sensitive to the variation of t_⊥. Fixing t_⊥ and varying J and J_⊥ independently gives qualitatively identical results as shown in Supplementary Fig. S4. More details on the numerical calculation can be found in the Supplementary Note 4.

The superconducting transition temperature T_c is calculated for different impurity concentrations for bonding-antibonding s_±-wave (with bonding/antibonding superconducting order parameter given by Δ_b/a = ± Δ_⊥, where Δ_⊥ is the momentum independent interlayer order parameter) and d-wave (\({\Delta }_{{{{\rm{b/a}}}}}={\Delta }_{{{{\rm{d}}}}}(\cos ({k}_{x})- \cos ({k}_{y}))/2\) with Δ_d denoting the intralayer d-wave order parameter) Cooper-pairings for a bilayer model with single orbital for various band fillings n = 1 (a), n = 0.75 (b), and n = 0.5 (c). The temperature is normalized to the calculated transition temperature of the clean system T_c0. The impurity concentrations are normalized using also the impurity strength W and density of states at the Fermi level N(0). The insets show the corresponding Fermi surface with bonding and antibonding bands as solid blue and dashed red lines, respectively. The imbalance (here it is controlled by the particle-hole asymmetry) between bonding and antibonding bands increases from left to right. The intralayer interaction J and interlayer hopping t_⊥ are modified such that T_c0 is roughly the same for all cases.