Supporting data and preprint "Superconductivity in the spin-state crossover materials: Nickelates with planar-coordinated low-spin Ni^2+ ions"

Superconductivity in the spin-state crossover materials: Nickelates with planar-coordinated low-spin Ni2+ ions Jiˇr´ı Chaloupka1and Giniyat Khaliullin2 1Department of Condensed Matter Physics, Faculty of Science, Masaryk University, Kotl´aˇrsk´a 2, 61137 Brno, Czech Republic 2Max Planck Institute for Solid State Research, Heisenbergstrasse 1, D-70569 Stuttgart, Germany (Dated: April 26, 2025) We theoretically study quasi-two-dimensional nickel compounds, where the nickel ions assume Ni2+ d8valence state and feature a low-spin S= 0 ground state quasidegenerate with S= 1 ionic excitations. Such a level structure is supported by square-planar coordination of nickel ions or a suitable substitution of apical oxygens. We construct the corresponding singlet-triplet exchange model and explore its phase diagram and excitation spectrum. By hole doping, we further introduce mobile Ni3+ d7ionic configurations with effective spins S= 1/2, and analyze their interactions with the d8singlet-triplet background. The interplay with the triplet excitations in the d8sector is found to have a deep impact on the propagation of the doped hole-like charge carriers and is identified as a powerful source of Cooper pairing among them. I. INTRODUCTION After decades of efforts to find cuprate-like superconductivity (SC) in nickel-based compounds, unconventional SC has finally been discovered in a number of nickelate families in recent years [1–8]. It is truly remarkable that the Ni valence states, orbital degeneracy, and fermiology vary broadly among the nickelate superconductors: While NdNiO2discovered in 2019 [1] apparently falls well into the cuprate’s “single-layer, single-orbital, spin one-half” paradigm, the bilayer La3Ni2O7[2] and trilayer La4Ni3O10 [6,7] nickelate SCs seem greatly differ from cuprates in all respects. Indeed, (i) in contrast to only weakly coupled CuO2planes in cuprates, the NiO2 planes are tightly bound into bilayers or trilayers in them; (ii) distinct from cuprates, both egorbitals x2−y2and 3z2−r2are essential; and, finally, (iii) Ni-valence is far from the cuprate-like spin one-half d9configuration, but rather in the mixed-valence regime with ∼7.5 electrons in the 3dshell. Such a diversity of Ni-based SCs should be related to the rich spin-and-orbital physics typical for nickelates, and suggests that more superconductors may be found in the future among the nickel-based materials with different lattice and orbital structures. One well-known consequence of the orbital degeneracy in compounds of transition metal ions such as Fe, Co, and Ni are the famous spin-state crossover phenomena [9]. In spin-crossover materials, a close competition between the lattice crystalline fields and intraionic Hund’s coupling results in the quasi-degeneracy of the ionic spin states. For example, the Ni2+ ion with d8 configuration may adopt either S= 0 singlet or S= 1 triplet ground state, as illustrated in Fig. 1(a). Importantly, the ground state spin value depends sensitively on the lattice structure details, covalency, etc, and can thus be tuned by external factors such as pressure or strain [10]. The spin-state quasi-degeneracy is quite common phenomena among the Ni2+ compounds with various lattice structures, e.g., Sr2NiO2Cu2(S,Se)2[11–13], (Sr,Ba)2NiO2Ag2Se2[14], BaNiO2[15], LaNiO2.5[16]. In this paper we perform a model exploration of singlelayer nickelates in the spin-state crossover regime. Specifically, we consider low-spin Ni2+ d8Mott insulators doped by mobile hole-like charge carriers, and address the question of pairing mechanisms and possible superconductivity in them. We find that the doped holes strongly interact with the spin-state fluctuations, namely, while moving they emit and absorb singlet-triplet excitations. This leads to polaron physics similar to that in the electron-phonon problem, and gives rise to an effective attraction between the charge carriers. We analyze and evaluate various pairing processes, and conclude that the nickel-based spin-crossover materials are promising candidates to host high-temperature superconductivity. The paper is organized as follows: Section II derives the exchange model for the undoped Mott insulator with low-spin Ni2+ d8ions (Sec. II A), and presents its phase diagram along with a cross-check against exact diagonalization of the underlying Hubbard model (Sec. II B). The mobile d7hole-like carriers are considered in Sec. III: we first introduce the interactions between the holes and the d8background (Sec. III A), then analyze the singleparticle propagation of the doped holes at low densities (Sec. III B), and finally move on to the central topic of the paper—a study of the mechanisms of Cooper pairing among the holes (Sec. III C). Section IV concludes the paper. II. LOW-SPIN Ni2+d8INSULATOR We construct our model as a low-energy effective model of the Hubbard model based on egorbitals. We assume that crystal field splitting ∆ between x2−y2and 3z2−r2 orbitals [Fig. 1(b)] is sufficient to compensate the Hund’s coupling JHand stabilize spin-singlet ground state. It is important to note at this point that JHin a solid is well reduced from its free-ion value due to covalency effects. 2 S=1 S=1 S=0 S=0 ∆>JH ∆<JH Ni2+ Ni2+ JH s = = T+1 f= 2 −r3z 2 −y 2 x2 (a) high−spinlow−spin eg ∆ (b) (c) spin-state crossover FIG. 1. (a) Ground-state spin configuration of a Ni2+ ion with d8(t6 2ge2 g) electronic filling, tuned by the competition of Hund’s coupling JHand the crystal field splitting ∆. If the Ni2+ ion is planar-coordinated or placed into a highly elongated oxygen octahedra, it would assume a singlet ground state with doubly occupied 3z2−r2orbital. (b) Splitting of egorbital levels in a crystal field promoting the above low-spin configuration. (c) Low-energy states selected as a basis of our effective model: d8singlet swith doubly occupied 3z2−r2 orbital, d8triplet T, where the CF splitting ∆ is nearly compensated by Hund’s coupling JH, and d7state corresponding to the hole-like particle fwith spin one-half. We further assume that the balance between ∆ and JHis such that while the ground state of Ni2+ is nonmagnetic, the excited spin-triplet states are still nearby in energy. In this regime, the relevant basis states are those depicted in Fig. 1(c). As a first step, we consider the non-doped case with d8(t6 2ge2 g) electronic configuration, and study the corresponding exchange model based on the ionic singlet sand triplet Tstates. A. Exchange model Built on top of a t6 2gconfiguration, which is considered rigid, the basis states sand Tare obtained by diagonalizing the on-site Coulomb interaction between the two additional egelectrons, as embedded in the local part of the egHubbard model: 1 2∆ (nx−nz) + U(nx↑nx↓+nz↑nz↓) + (U−5 2JH)nxnz −2JHSxSz+JH(x† ↑x† ↓z↓z↑+z† ↑z† ↓x↓x↑).(1) Here we have introduced a shorthand notation by denoting the electron operators dx2−y2and d3z2−r2as xand z, respectively. Due to the pair hopping term ∝JH, the singlet ground state |s⟩is actually a linear combination |s⟩=cos θ z† ↑z† ↓−sin θ x† ↑x† ↓|t6 2g⟩(2) with θgiven by tan 2θ=JH/∆. In the regime of our interest, θis small and most of the weight is carried by the pair of electrons in 3z2−r2orbital as indicated in Fig. 1(c). The three members of the triplet Tread simply as |T+1⟩=x† ↑z† ↑|t6 2g⟩, |T0⟩=1 √2x† ↑z† ↓+x† ↓z† ↑|t6 2g⟩, |T−1⟩=x† ↓z† ↓|t6 2g⟩.(3) The energy splitting between the above states is an essential parameter of our model. Based solely on the ionic Hamiltonian, we get for the singlet-triplet splitting ET=E(|T⟩)−E(|s⟩) = q∆2+J2 H−3JH.(4) In reality, this splitting is influenced by various factors beyond the ionic Hamiltonian, for example covalency effects, and may be therefore regarded as a free parameter. However, for simplicity we will maintain the connection to the egHubbard model and use Eq. (4) hereafter. Note, that if the pair hopping term is omitted, we get the familiar expression ∆ −3JHfor ET. The exchange interactions among d8ions are derived in the usual way by perturbatively eliminating the intersite hopping to second order. In the case of egorbitals on the square lattice, the hopping at nearest-neighbor bond ij takes the following form fixed by Slater-Koster rules: −thx†x+1 3z†z∓1 √3(x†z+z†x)iij .(5) We parametrize the hopping by the larger amplitude t for the planar x2−y2orbitals. The doped holes to be introduced later will reside in 3z2−r2orbitals, having thus a significantly reduced bandwidth determined by t/3. The sign (∓) of the interorbital hopping is linked to the bond direction: −is applied at xbonds, + at y bonds. The exchange interactions are most intuitively depicted as bond processes involving s(scalar) and T(vector) hardcore bosons representing the singlet and triplet d8states. The hardcore bosons Tare hereafter referred to as triplons. Such a hardcore-boson representation is used in Fig. 2(a),(b), where we show a few examples of the rather rich possibilities to be discussed later. Fortunately, due to the inherent spin-isotropy imposed by the underlying Hubbard model, the exchange interactions have to obey spin conservation rules, which limits the number of independent model parameters. To express the resulting d8model Hamiltonian in a compact form, we utilize several auxiliary vector operators based on the 3 following Cartesian combinations of the triplet states: |Tx⟩=i √2(|T+1⟩−|T−1⟩), |Ty⟩=1 √2(|T+1⟩+|T−1⟩), |Tz⟩=−i|T0⟩.(6) The first vector operator e S= (e Sx,e Sy,e Sz) is local and is associated with the on-site s↔Ttransition: e Sα=−i(s†Tα−T† αs), α =x, y, z . (7) The second local operator S= (Sx, Sy, Sz) directly corresponds to the spin-1 carried by Tand is defined by Sα=−i ϵαβγ T† βTγ,i.e. S=−i(T†×T).(8) We also introduce analogous bond operators having the internal structure identical to e Sand S: e Sα ji =−i(s† jTαi −T† αjsi), Sα ji =−i ϵαβγT† βj Tγi .(9) The above operators enable us to express the d8model Hamiltonian in a manifestly spin-isotropic compact form Hd8=ETX i nT i +κX ⟨ij⟩e Sie Sj +X ⟨ij⟩hsgnij K(e SijSji +Sij e Sji) + JSiSji.(10) Here the first term counts the number of triplons via nT=Pα=x,y,z T† αTαand penalizes them by the energy ET, the other terms represent the exchange interactions being sorted as second, third, and fourth-order contributions in Toperators, i.e. according to their expected relevance at low triplon density. Due to egorbital symmetry, there is sensitivity to the bond direction inherited from the hopping (5). It is captured by the bond-dependent sign factor sgnij, which is equal to +1 or −1 for xand ybonds, respectively. Finally, when rewriting the operators in the Kterm using the elementary hardcore-boson operators sand T, these have to be taken in normal order with the creation operators moved to the left. We now briefly discuss the rather complex exchange processes contained in (10) with the help of the cartoon representations in Fig. 2(a),(b), capturing the exchange in the form of a generation, propagation, and mutual interactions of hardcore triplet particles Tin a background made of s. The κ-term e Sie Sjincludes two distinct types of processes depicted in Fig. 2(a). The first is a creation and annihilation of singlet pairs of Ton bonds, corresponding to a combination TxTx+TyTy+TzTz= T+1T−1−T0T0+T−1T+1. The other one is a simple T hopping preserving the Tlabel (either ±1, 0 or Cartesian x, y, z). The relative amplitude of these processes is fixed to −1 by the interaction term. Being of second order in Toperators, this is the most important exchange interaction at low triplon concentration of about nT< ∼0.2 per K J exchange [meV]  [eV] 15 20 25 30 1.6 1.8 2.0 2.2 K J exchange [meV] JH [eV] 15 20 25 30 0.4 0.5 0.6 0.7 0.8 (a) (b) (c) (d) s s κ κ s K J T+1 T s s TT T T TT T TT −1 0 +10 +1+1 +1 −1 0+1 FIG. 2. (a) Examples of exchange processes among the d8ions contributing to the κinteraction channel of Hd8 of Eq. (10). There are two distinct classes of processes— creation/annihilation of singlet pairs of triplons T(top) and Thopping in the sbackground (bottom). (b) Sample processes contained in the higher-order channels Kand Jof Hd8. (c) Exchange parameters entering Eq. (10) as obtained by second-order perturbation theory [see Eq. (11)] assuming U= 5 eV, JH= 0.6 eV, and t= 0.3 eV. (d) The same for fixed ∆ = 2 eV and varying JH. site. The third-order Kcontribution may be understood as an s↔Ttransition entangled with the spin-1 of another Tparticle, while the fourth-order Jterm is just the regular Heisenberg exchange between the spins-1 carried by T. Two examples of such higher-order processes are presented in Fig. 2(b). Through the perturbative calculation, one arrives at two more contributions not explicitly included in (10): (i) a small renormalization of ETlevel, (ii) repulsion between Tparticles. Both are marginal, opposing each other, and can be safely absorbed into ET(the latter one on a mean-field level), so that they do not alter the form of the d8model (10). The overall phase behavior of the d8model (10) is driven by the competition between the triplon cost ET and the strength of the exchange interactions that profit from the presence of the triplons and eventually lead to their condensation. The microscopic derivation provides the following values of the three exchange parameters: κ=t2 31 + cos 2θ E−1 +2 sin 2θ E0 +1−cos 2θ E+1 , K=2√2t2 3√3 1 E−1 +1 E0cos θ+1 E+1 +1 E0sin θ, J=t2 31 E−1 +10 3E0 +1 E+1 ,(11) 4 where the denominators contain virtual energies En= U+n∆ + JH+ 2ET(n= 0,±1) in the intermediate states. Plotted in Fig. 2(c),(d) are the exchange constants evaluated for representative values of the microscopic parameters. For the Ni ions, we assume U= 5 eV and JH= 0.6 eV throughout the paper, the crystal-field splitting ∆ is varied around 2 eV, leading to quasidegenerate sand Tlevels. The hopping tis set to a representative value of 0.3 eV in Fig. 2. The values of κ,K, and Jobserved in Fig. 2(c),(d) are relatively robust to variations in ∆ and JH. This is caused by the dominating values of U+n∆ in the denominators in (11), which are only little affected by changes ∆ and JHwithin the ranges shown. In the presented parameter regime we get the hierarchy J > K > κ, which, however, needs to be considered together with the different role of these interactions. Even though Jand Kare larger than κ, they are of higher order in Toperators and thus relevant only at larger concentration nTof triplons. In addition, to activate Jand Kinteractions, a prior presence of triplons is necessary, while κinteraction channel generates them itself. In some of the cases discussed later, it will be thus sufficient to study a simplified model including ETand κterms only. B. Phase diagram and excitations A basic exploration of the phase behavior of the singlet-triplet model and its excitation spectrum can be performed using the methods developed in the context of spin ladders or bilayer magnets, where the singlet-triplet basis is hosted by Heisenberg rungs connecting the rails or layers [17–20]. The possibility of triplon condensation at a sufficient exchange strength is addressed using the variational Ansatz |Ψ⟩=Y Rhp1−ρ s†+√ρ(u−iv)T†iR|vac⟩.(12) The energy obtained by averaging the Hamiltonian (10) in the trial state (12) Eavg =ETX i ρ+X ⟨ij⟩n4κρ (1 −ρ)vivj −sgnij Kρpρ(1 −ρ) [vi(u×v)j+ (u×v)ivj] + 4Jρ2(u×v)i(u×v)jo(13) is to be minimized with respect to the condensate density 0 ≤ρ≤1 and site-dependent real vectors u,vthat obey the condition u2+v2= 1. This approach does not include quantum fluctuations but still gives a semiquantitatively correct picture, as we will demonstrate later by cross checking with the exact diagonalization of the underlying Hubbard model. Presented in Fig. 3(a) is the variational phase diagram obtained by varying the microscopic parameters ∆ and tand evaluating the effective model parameters ET,κ, K,Jvia Eqs. (4) and (11). The variational approach identifies three phases described in the following. In the paramagnetic (PM) phase, the condensate is absent (ρ= 0) and triplons Tcorrespond to gapped elementary excitations. Their dispersion can be obtained by starting with Hd8of Eq. (10), making the replacement s, s†→√1−nTto account for the hardcore constraint, and expanding the result in Tα,T† αoperators. Afterwards, Tare treated as unconstrained bosons. At the same level of approximation as used when constructing the phase diagram itself, the expansion is performed to second order and the approximate Hamiltonian is readily solved by Fourier and Bogoliubov transformations. This standard procedure [19] gives a three-fold degenerate dispersion ωq=qET(ET+ 8κγq) (14) with γq=1 2(cos qx+ cos qy), visualized in Fig. 3(b) for a few sample points in the PM phase. As the exchange strength increases, following the increase of t, the excitations gradually soften at the momentum (π, π). When the critical point ET= 8κis reached, the excitation gap closes, signaling a condensation of triplons and the emergence of long-range order. The observation of the above excitations, e.g. by neutron scattering, would be an important experimental test of our model and would provide a valuable input for its quantification in a particular material. We also note that a significant broadening of the triplons due to their mutual interactions as well as interactions with doped carriers can be expected. In the new phase, labeled by FQ in Fig. 3(a), the optimization of the variational Ansatz (12) leads to a nonzero condensate density ρ. The condensate structure is given by staggered vR=veiQR with Q= (π, π) and zero u. The corresponding energy per site derived from (13) takes the value EFQ =−2κ(1 −ET/8κ)2and does not include Kand Jas these contribute at higher order via quantum fluctuations. Based on the above condensate structure, the phase is characterized by long-range correlations of the auxiliary quantity e S∝vwith the characteristic momentum Qand nonzero average ⟨e S⟩q=Qacting as an order parameter. In terms of the measurable spin-1 carried by T, one gets the average ⟨SαSβ⟩q=0 =ρ(δαβ −vαvβ) corresponding to a ferroquadrupolar (FQ) order, where vplays the role of the director. Note, however, that the primary order parameter is the above ⟨e S⟩q=Qand the FQ spin-1 correlations are reduced by having Tmixed with still prevailing s. Finally, if the triplon cost ETis low enough, the system switches into a fully saturated antiferromagnetic (AF) phase with ρ= 1, i.e. each site is occupied by spin-1 Tparticle. This transition may be thus interpreted as the spin-state crossover between the low-spin and highspin state of d8, realized in our approximation by a level crossing. Using Eq. (13), the AF phase gives the average energy per site EAF =ET−2J, which has to be com- 5 t [eV] Δ [eV] 0.0 0.1 0.2 0.3 0.4 0.5 1.6 1.8 2.0 2.2 0.0 0.2 0.4 0.6 0.8 1.0 nT one hole Δ [eV] 1.6 1.8 2.0 2.2 0.0 0.2 0.4 0.6 0.8 1.0 nT two holes t [eV] 0.0 0.1 0.2 0.3 0.4 0.5 0.0 0.2 0.4 0.6 0.8 1.0 nT no holes t [eV] 0.0 0.1 0.2 0.3 0.4 0.5 0 2 4 6 8 〈S S〉 (π,π) no holes 0 1 2 3 4 5 〈S S〉 (π,π) no holes ∼∼ ρ t [eV] Δ [eV] ET [eV] 0.0 0.1 0.2 0.3 0.4 0.5 1.6 1.8 2.0 2.2 0.0 0.2 0.4 0.0 0.2 0.4 0.6 0.8 1.0 PM FQ AF (a) (b) (c) (d) (e) (f) t=0.2 0.3 0.4eV ω [eV] ΓX M Γ 0.0 0.1 0.2 0.3 0.4 FIG. 3. (a) Phase diagram of the d8model in Eq. (10) obtained by a variational estimate based on the trial state (12). Shown is the condensate density ρwhich equals zero in the PM phase, saturates around 0.3 in the FQ phase and jumps to ≈1 upon entering the AF phase. The model parameters are given by Eqs. (4) and (11) evaluated for U= 5 eV, JH= 0.6 eV and variable ∆ and t. (b) Triplon Tdispersion in PM phase calculated for the three parameter points marked in (a) (∆ = 2 eV and t= 0.2, 0.3 or 0.4 eV) and plotted along the path connecting the high-symmetry points Γ = (0,0), X= (π, 0), and M= (π, π) in the Brillouin zone. At the PM/FQ boundary, the dispersion becomes gapless at the Mpoint. (c) Triplon number per site nTobtained by exact diagonalization of the egHubbard model for the 8-site cluster shown in (d). (e) Characteristic correlations in the cluster ground state revealing the tendencies towards longrange ordering—either AF (left) or FQ (right). (f) The same as in (c) but calculated for a reduced number of electrons on the cluster to effectively create 1 or 2 holes. pared with EFQ to find the transition line observed in Fig. 3(a). Following the AF/FQ phase boundary from the large tto small tlimit, we find the associated step in the Tdensity increasingly pronounced, reaching maximum at t= 0, where we hit the point corresponding to ET= 0. This point is a trivially exact point of the phase diagram—here ETchanges sign and in the absence of exchange interactions due to t= 0, the system immediately switches between the state consisting exclusively of either sor T. To assess the adequacy of the effective model and its simple variational treatment, we performed exact diagonalization of the two-orbital Hubbard model with eg-type hopping on a small square-lattice cluster. The first set of data to compare to the variational phase diagram is presented in Fig. 3(c). It shows the average triplet occupation nTper site calculated for the 8-site cluster of Fig. 3(d) by employing the same microscopic parameters ∆, t,U, and JHas those used to construct Fig. 3(a). Note that, in contrast to ρof Fig. 3(a), nTof Fig. 3(c) includes also the fluctuating part of Tdensity. While the map of nTitself gives already a good idea of the phase diagram, to reliably identify the phases, we check the tendency towards long-range order by calculating the characteristic correlations ⟨SS⟩q=Qand ⟨e Se S⟩q=Qassociated with the AF and FQ phases, respectively. Shown in Fig. 3(e) are the ground-state values of the correlations evaluated according to ⟨SS⟩q=N−1 site Pj=j′eiq(Rj−Rj′)⟨SjSj′⟩, i.e. excluding the on-site contribution. The overall agreement between Fig. 3(a) and Figs. 3(c),(e) is more than satisfactory, given that the variational phase diagram in fact involves two levels of approximation: (i) the effective model (10) was derived via perturbation theory to second order in electron hopping, (ii) the effective model was solved using a crude variational Ansatz. Moreover, the exact diagonalization for such a small cluster may be expected to overestimate the effect of quantum fluctuations, shifting the phase boundaries significantly. A comparison of the results for 4-, 6-, and 8-site cluster indeed shows a systematic trend of the PM/FQ boundary moving slightly up with increasing system size. Therefore, the Hubbard model calculation well confirms the topology of the phase diagram and a rough location of the phase boundaries observed in Fig. 3(a), suggesting additionally a smoother crossover between the high-spin AF phase and low-spin FQ phase. Finally, Fig. 3(f) demonstrates, that the overall structure of the phase diagram does not change upon hole doping that we simulate by removing one or two electrons from our Hubbard system, corresponding formally to 12.5% and 25% doping. Though the data presented in Fig. 3(f) are necessarily plagued by large finite-size effects, we may anticipate roughly the same location of the PM/FQ phase boundary and a growth of the FQ phase at the expense of the AF phase. This trend is in line with intuition, since the FQ phase characterized by the dynamical mixing of sand Tcan better accommodate the d7mobile carriers and their interplay with the d8background to be discussed in the next section. The FQ phase might be therefore energetically favorable to the AF one upon doping. On the other hand, from the same perspective there seems to be only a little difference between the FQ and PM phases. 6 III. DOPED CASE: d7HOLES WITHIN d8 BACKGROUND In this section, the effective model will be completed by introducing mobile carriers represented by holes corresponding to the d7configuration fshown in Fig. 1(c). We only include the d7states with the single egelectron occupying the 3z2−r2orbital. The participation of d7 configurations of x2−y2character is partially suppressed by their higher energy due to ∆, partially it is caused by on-site correlations preventing their motion even though the hopping amplitude associated with the x2−y2orbital is three times larger than that of 3z2−r2. We will discuss these issues at relevant points later. As a result, the x2−y2electrons “live” in the system mostly bound in the triplet d8configurations. In the following, we first discuss the derived interactions between the d7and d8objects in our model, then focus on the individual propagation of the doped holes, and finally examine their pairing tendencies and implications for possible superconductivity. A. Interaction between holes and background The interaction of the d7holes with the d8background consisting of sand Tparticles is derived in a straightforward way by projecting the electronic hopping on a d7–d8 bond onto the low-energy states constituting the basis of our model. It takes a form of fhopping accompanied by “counterflow” of d8objects sand T, which may involve changes in their state or even an s↔Ttransition. Several sample processes of this kind are illustrated in Fig. 4. To express the interaction in a compact form, we again utilize the bond d8operators e Sij and Sij defined by Eq. (9) and introduce new bond operators that capture the hopping of f. The first one is a regular spinindependent hopping nij =f† i↑fj↑+f† i↓fj↓,(15) the second one involves the fspin and is defined via σα ij =X s,s′=↑,↓ f† isσα ss′fjs′(16) with σα(α=x, y, z) denoting the Pauli matrices. In addition to transferring fbetween the sites, the latter operator performs either a flip of the fspin (σx,σy) or makes the hopping sensitive to the spin value (σz). The final d7–d8Hamiltonian reads as Hd7–d8=X ⟨ij⟩sgnij Aσij e Sji +BσijSji +nij(C0s† jsi+C1T† jTi)+ H.c.(17) As in the case Hd8of Eq. (10), the form of the Hamiltonian Hd7–d8is again explicitly spin-isotropic and in fact s f +1 T C0 A s A x2−y 2 (a) B fT f T f +1 0 (b) ss f f f TT f C1 +1 +1 (c) (d) f−1 T t/3 ff +1 T s FIG. 4. Cartoon representation of the fparticle motion through the d8background, i.e. examples of d7–d8bond processes contained in Eq. (17). (a) Hole hopping involving spin degrees of freedom: Within the process A, the hole hopping is entangled with a singlet-triplet transition; during process B, the hole exchanges its position with a triplon Tin a spinsensitive manner. (b) Spin-independent hole hopping: the hole exchanges its position with either a singlet s(process C0) or triplon T(process C1) without altering any spin state. (c) Process Adepicted using the representation of the spinorbital states as in Fig. 1(c). Its origin is the interorbital hopping of the electron from d8to d7site. The electron added to the x2−y2orbital gets bound into a triplon Tby virtue of Hund’s coupling. (d) Hopping process converting 3z2−r2 hole f↑into a virtual (x2−y2)↓fermion, at the energy cost of 3JHrequired to break S= 1 triplet Tstate at the neighboring site. dictated by this symmetry. Analogous forms of the interaction between holes and singlet-triplet background can be encountered e.g. in the context of hole-doped Heisenberg ladders or bilayers [21,22] or in the context of spincrossover cobaltates, where the triplons carry in addition an orbital degree of freedom [23,24]. Following the microscopic derivation, the individual parameters entering (17) are given by A=1 √6tcos θ , B =1 2t , C0=1 3tcos2θ , C1=1 2t . (18) A few examples of different types of processes contained in Hd7–d8are sketched in Fig. 4(a),(b). The most interesting one corresponds to the Aterm representing a hole hopping accompanied by s↔Ttransition. Ex- 7 panded in full, the relevant term 1 √2σij e Sji takes the form T† −1jf† i↑fj↓+T† 0j 1 √2(f† i↑fj↑−f† i↓fj↓)−T† +1jf† i↓fj↑+ H.c. (19) For brevity, we have omitted siin each contribution and used the 1 √2prefactor. The Aprocesses are essential to understand the effects arising due to the hole doping in the PM phase of the background. The dynamical generation and annihilation of triplet excitations Tby the doped hole has a profound impact on its propagation and leads to clear polaronic features discussed in Sec. III B. It may become also a source of hole pairing that is evaluated in Sec. III C. When focusing on the PM phase, the Aprocesses need to be considered together with the processes linked to the parameter C0, that just exchange the position of fand son the bonds and provide therefore a free motion of fthrough the predominantly sbackground. The processes of corresponding to the Band C1terms become important at larger Tdensity and are thus relevant in case of the ordered d8background that is beyond our scope here. Another possibility of visualizing the hopping processes in Hd7–d8—based on the schematic representation of the states of Fig. 1(c)—is used in Fig. 4(c). Starting with the fhole and its most frequent bond-neighbor s in the PM phase, the interorbital hopping can be easily associated with the Aprocess. We note that the hopping processes in the underlying egorbital Hubbard model allow a mixing of the fhole of 3z2−r2symmetry with the d7(x2−y2) states not included in our effective s-T-fmodel. Figure 4(d) shows such an example, where 3z2−r2intraorbital hopping converts f-fermion into a virtual d7(x2−y2) state. This process breaks S= 1 triplet bound state and costs a large energy of 3JH. As the latter far exceeds a mixing matrix element of t/3, admixture of x2−y2orbital into the f-fermion wave function is small and can be neglected. This is confirmed below by the ED analysis of the full egorbital Hubbard model. Thus in our effective s-T-fmodel, the x2−y2electrons are present only in a hidden form, i.e. bound in the triplon Tstates of d8configuration. B. Hole motion – renormalization within SCBA In the previous section we found that the motion of doped holes proceeds via rather complex bond processes, where the hopping of d7holes is intertwined in various ways with triplet excitations in the d8background. The aim of this section is to analyze, how the coupling to the triplet excitations affects the motion of doped holes. We will focus specifically on the PM phase, where it is sufficient to consider just ETand κterms from Hd8(10) describing the magnetic background and Aand C0processes from the Hd7–d8interaction Hamiltonian (17). The other terms that are higher-order in Toperators can be omitted in the PM phase due to the low triplon density. The bare propagation of the doped holes (fparticles) in the d8background, which is mostly composed of s singlets in the PM phase, is captured by the C0processes of (17) leading to the bare dispersion εk= 2C0(cos kx+ cos ky). The hopping amplitude C0≈1 3t[see Eq. (18)] is derived from the weak 3z2−r2intraorbital hopping, resulting in a small bare bandwidth. The holes couple to the triplet excitations represented as before by the triplons Tα(α=x, y, z) of hardcore bosonic character, carrying the elementary excitations of the d8background. The coupling is provided by the Aterm of (17), which transforms into the momentum representation as Hd7–d8→iA X kqαss′ ηk−qT† αqf† k−q,sσα ss′fk,s′+ H.c. (20) with ηq= 2(cos kx−cos ky). The approximate diagonalization of Hd8as described in Sec. II B leads to three degenerate eigenmodes αqwith the dispersion (14). The corresponding operators are constructed from Tx,y,z by Bogoliubov transformation Tαq=uqαq+vqα† −q, where the Bogoliubov factors take the form uq=1 √2sET+ 4κγq ωq + 1 , vq=2κγq ωquq .(21) Inserted into Eq. (20), this gives the linear coupling of f to the excitations αq: Hd7–d8→iA X kqαss′ Γkqα† qf† k−q,sσα ss′fk,s′+ H.c. ,(22) with the formfactor Γkq =ηk−quq−ηkvq. The propagation of fthrough the d8background including the coupling (22) to its excitations is treated within standard selfconsistent Born approximation (SCBA) scheme commonly used to address magnetic polarons (see e.g. [25–28]). For the case of a single hole, we get the selfenergy Σ(k, E) = 3A2X q Γ2 kq G(k−q, E −ωq),(23) which enters the single-hole propagator G(k, E)=[E+ i0+−εk−Σ(k, E)]−1to be determined selfconsistently. Presented in Fig. 5(a) is the corresponding spectral function A(k, E) = −π−1Im G(k, E) obtained for a representative parameter point that falls into the PM phase. The spectral function bears typical polaronic features. The bottom of the hole band around k=Mis only mildly affected, but the higher parts are strongly renormalized with the quasiparticle weight reduced by factor ≈2−3. Compared to the bare dispersion, the total bandwidth is about two times smaller. The energy position of the onset of the damping and dispersion flattening correlates with ET, in this case evaluated to ≈0.29 eV. The spectral function resulting from the effective model is well confirmed by exact diagonalization of the 8 E [eV] 0.0 0.5 1.0 1.5 ΓX M Γ E [eV] 0.0 0.5 1.0 1.5 ΓX M Γ (a) (b) (c) DOS [a.u.] E [eV] 3z2-r2 x2-y2 0 1 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 FIG. 5. (a) Single-hole spectral function obtained within SCBA applied to the effective s-T-fmodel. The parameter point corresponds to ∆ = 2 eV and t= 0.3 eV in Fig. 3(a). The spectral function is presented both in a form of a color map as well as profiles at selected high-symmetry points in the Brillouin zone. The dotted line shows the bare dispersion εk. (b) Spectral function of a 3z2−r2hole obtained by exact diagonalization of the egHubbard model at the 8-site cluster for the same parameter point. Both spectral functions were artificially broadened for an easier comparison and shifted to measure the energy from the bottom of the band. (c) Singleparticle density of states (DOS) for the 8-site cluster resolved according to the contributing egorbital. The DOS corresponds to an average of the spectral function over the eight accessible momenta in the Brillouin zone. underlying Hubbard model. In Fig. 5(b) we present the spectral function corresponding to an extraction of 3z2−r2electron AED(k, E) = −1 πIm ⟨GS|z† kσ 1 E+i0+−HHubb zkσ|GS⟩, (24) which may be directly compared to the fspectral function. The calculation was performed for the 8-site cluster of Fig. 3(d) used before. |GS⟩denotes the cluster ground state. Given the smallness of the cluster and the approximations made in the derivation and treatment of the effective model, the agreement is very good. This supports the picture that we have drawn based on the effective model, i.e. hole propagation along with intense generation and annihilation of the triplet excitations. In Fig. 5(c) we return to the question of the relevance of the x2−y2holes for the effective model, as discussed in the end of Sec. III A. Using the exact diagonalization of the Hubbard model, we can directly quantify the degree of their involvement. To this end, we study the orbital-resolved single-particle density of states, evaluated as an average of the corresponding spectral function over the cluster-compatible momenta k: DOS(E) = N−1 kPkA(k, E). Here A(k, E) is either the spectral function for the 3z2−r2hole presented for selected kpoints in Fig. 5(b), or analogous one for the x2−y2hole calculated as −π−1Im ⟨GS|x† kσ(E+i0+− HHubb)−1xkσ|GS⟩. As seen in Fig. 5(c), the DOS related to the x2−y2hole is indeed negligible compared to the 3z2−r2one and is concentrated in two energy regions. The first one below approximately 0.5 eV covers essentially the range of the low-energy band of propagating 3z2−r2hole. This feature indicates a slight mixing of the x2−y2and 3z2−r2holes stemming from the hole-type conversion processes discussed in the context of Fig. 4(d). The second region is located at energies higher by about ∆, which corresponds to the expected excitation energy of the x2−y2hole as compared to the 3z2−r2one. C. Hole pairing and BCS estimates Focusing on the PM phase again, we now consider the dynamic generation and annihilation of triplet excitations as a source of pairing among the holes. We will follow the usual procedures applied in the theory of conventional superconductivity based on electron-phonon coupling and derive the hole-hole interaction by perturbatively eliminating the coupling Hamiltonian (17) that connects the d7and d8sectors, obtaining thereby the effective interaction acting in the d7sector. Specifically, from Hd7–d8of Eq. (17) we again take the A-channel part, here in the form Pij sgnij Aσij e Sji. The sum now runs through all nearest-neighbor i,jand includes thus also the H.c. counterpart of the bond sum in (17). The d7hole operators are embedded in the bond operators σij with the components given by (16), while the d8sector is represented by e Sji defined by (9). By eliminating the above d7–d8coupling via second-order perturbation theory, we obtain the effective interaction among f, which can be expressed in the form Hf–f≈ −A2X ijα X klβ sgnij sgnkl Fαβ ji,lk σα ij σβ kl (25) with Fbeing the ⟨e Se S⟩propagator evaluated entirely within the d8sector. It is defined as the d8ground-state average Fαβ ji,lk =DGS e Sα ji 1 Hd8−EGS e Sβ lk GSE.(26) 9 Recalling the similarity to the electron-phonon problem, the propagator Fentering the hole-hole interaction is analogous to the static limit of the phonon propagator mediating effective electron-electron interaction. Thanks to the spin-space symmetry of the d8background, the quantity (26) is isotropic, Fαβ ji,lk =Fji,lkδαβ, and we need to consider just the diagonal elements Fji,lk in the following. Note that by construction, the pairs of indices i, jand k,lcorrespond to nearest neighbors (NN) and we thus deal with a bond-bond correlator. The evaluation of Fji,lk for the PM phase is performed using an approximate solution of the d8model limited to ETand κterms as in Sec. III B. One can proceed, for example, by converting (26) to the momentum representation and utilizing the previous description of the PM phase based on the Bogoliubov transformation of T. For simplicity and to enable a transparent interpretation of the result, we will limit ourselves to an expansion to first order in κ/ET, where Fji,lk can be also easily evaluated using real-space perturbation theory and the contributions visualized in a pictorial way as done in Fig. 6(a)- (c). In this expansion, we start with the ETpart of Hd8of (10) as the unperturbed Hamiltonian and add the κ-interaction playing the role of the perturbation. Up to first order in κ/ET, we get Fji,lk =1 ETδil −κ ET (δNN il +1 2δNN ik +1 2δNN jl ).(27) The leading term of zeroth order contains standard Kronecker δil which is equal to unity if i=l. The physical process behind this contribution to the effective f–f interaction (25) is a creation of triplon by hopping of f during the first Aprocess and its immediate absorption in the second Aprocess [see Fig. 6(a) for a cartoon representation]. The possibility of such a pair hopping and the related kinetic energy gain drives the pairing of the holes into short-range pairs. Higher-order processes, that involve some action of κ, give rise to a more extended attraction of holes. In the first order included in (27), there are two physically distinct contributions that can be associated with the processes depicted in Fig. 6(b),(c). The first contribution corresponds to the term δNN il that checks whether the sites iand lform a nearest-neighbor bond. The process behind is similar to that of Fig. 6(a), but there is an extra hopping of the triplon between sites iand lin the intermediate state, which is provided by the κinteraction in the d8sector. The second contribution, linked to the combination 1 2(δNN ik +δNN jl ) in (27), has a completely different physical origin. It exploits the fact, that the κinteraction preforms singlet pairs of triplons T, which may—through two correlated hoppings within the A-channel of Hd7–d8—resonate with singlet pairs of holes. The triplon pairs appear with amplitude ∝κ/ET, hence such processes appear as first-order in (27). The next step is to insert the approximate Fji,lk into (25) and collect all the possible process pathways. Here the real-space formulation proves quite helpful, as it enables to handle the various constraints explicitly. The T-pair T-hopping V1 / V0 ratio doping 0 1 2 3 4 0.0 0.1 0.2 0.3 0.4 0.5 BCS  ET [eV] 0.0 0.2 0.4 0.6 0.8 1.0 0.1 0.2 0.3 0.4 0.5 (a) (d) (e) (c) (b) κ A / ET ∼ κ A A ff sf f s A AfT+1 f A triplon singlet Cooper pair resonance FIG. 6. (a) Cartoon representation of the dominant pairing process taking place on two adjacent bonds. One of the fholes undergoes Ahopping leaving a triplon Texcitation behind, which is later absorbed by the other hole in the second Ahopping. The amplitude of this process is proportional to A2/ET. (b),(c) Examples of additional pairing processes arising as corrections to (a) in the first order in κ/ET. Process (b) is similar to (a) but with an extra hopping κof the intermediate triplon T. The process (c) corresponds to an absorption of a singlet pair of triplons T(created via the κ exchange channel and having an amplitude ∝κ/ET) by a pair of ffermions. The absorption occurs via two Ahoppings with the intermediate state containing one unpaired T, hence its amplitude is ∝A2/ET. Such a hopping chain may proceed in both directions, leading to a promotion of fsinglets by their resonance with the singlet triplon pairs fluctuating in the d8 background by virtue of the κexchange. (d) Relative contribution of first-order terms to s-wave pairing potential (31) averaged over a circular Fermi surface at various doping levels. Following Eq. (31), the numbers have to be multiplied by κ/ETwhen comparing the absolute pairing strength. (e) BCS λfor fixed t= 0.3 eV and varying ∆ presented as function of ETat doping level 0.25. Total value is shown as black solid line, contributions from V0,T-hopping, and T-pair firstorder terms are shown as dashed/red/blue lines, respectively. ETinterval corresponding to long-range order is indicated by shading. result is converted into momentum space and adapted to