Emergence of superconductivity in doped multiorbital Hubbard chains

We introduce a variational state for one-dimensional two-orbital Hubbard models that intuitively explains the recent computational discovery of pairing in these systems when hole doped. Our ansatz is an optimized linear superposition of Affleck–Kennedy–Lieb–Tasaki valence-bond states, rendering the combination a valence-bond liquid dubbed orbital resonant valence bond. We show that the undoped (one-electron/orbital) quantum state of two sites coupled into a global spin singlet is exactly written employing only spin-1/2 singlets linking orbitals at nearest-neighbor sites. Generalizing to longer chains defines our variational state visualized geometrically expressing our chain as a two-leg ladder, with one orbital per leg. As in Anderson’s resonating valence-bond state, our undoped variational state contains preformed singlet pairs that via doping become mobile, leading to superconductivity. Doped real materials with one-dimensional substructures, two near-degenerate orbitals, and intermediate Hubbard U/W strengths—W the carrier’s bandwidth—could realize spin-singlet pairing if on-site anisotropies are small. If these anisotropies are robust, spin-triplet pairing emerges.

Within quantum materials, quasi-one-dimensional (1D) ladders and chains are attractive because powerful computational techniques, such as the density matrix renormalization group (DMRG) 39,40 and Lanczos 41 , allow for the study of model Hamiltonians with accuracy. This removes the veil of theoretical uncertainty in higher dimensions that complicates the comparison theory vs. experiment. In particular, this 1D avenue may allow for the challenging study of systems where both topology and correlations are simultaneously relevant.
In this context, there are few studies of the effects of hole doping and magnetic-based hole pairing on topological interacting systems. Early work in the t-J limit (no double occupancy) for doped S = 1 chains, indicated a narrow region of pairing, suppressed by competing ferromagnetism 42 . In related work, triplet superconductivity was also analyzed in 1D 43,44 . More recent efforts using multiorbital Hubbard models unveiled robust tendencies to spin-singlet pairing, an exciting result 45,46 . However, these valuable computational efforts did not provide a concrete mechanism as an explanation. In particular, we lack a simple intuitive picture connecting the topological properties of Haldane chains and the emergence of hole pairs in Hubbard models. Developing such a simple "cartoon" may allow generalizations to other systems and also facilitate the experimental search for realizations in particular materials.
Here we fill this conceptual gap. Our main conclusion is illustrated in Fig. 1. We consider the two-orbital Hubbard model on a chain (Fig. 1a), using the two orbitals a and b as legs of a mathematically equivalent two-leg ladder (Fig. 1b). We rely on a hereby proposed variational state: the orbital generalization of the resonating valence-bond concepts 47 . We employ preformed spin-1/2 singlets as in the original formulation, but now in the enlarged space spanned by the real chain in one direction and the orbital index in another (Fig. 1c). More simply, our state-the orbital resonant valence bond (ORVB)-is an optimized linear combination of Affleck-Kennedy-Lieb-Tasaki (AKLT) valence-bond solids 3,48,49 , rendering the proposed state a liquid. Doping this state with two holes in principle could break two singlets. But when holes are close to one another, they break only one singlet minimizing the energy and leading to an effective singlet hole pairing (Fig. 1d), in agreement with computational results 45,46 . Our undoped and doped states are variational, not exact, but they capture the essence of the problem, as shown below.
Our conclusions are not obvious: naively there are preformed triplets at each site because of the robust J H /U. Actually, we found that triplet pairing-a rarity 50,51 -becomes stable when easyplane anisotropies are not negligible. However, doping particular quasi-1D materials-with two active fairly equivalent orbitals and weak spin anisotropy-should lead instead to singlet pairing.
The present effort for two-orbital chains has qualitative implications for other multiorbital systems, such as iron-based superconductors. More specifically, our results, and the hole binding found in multiorbital ladders 33 , show that magnetic fluctuations induce pairing in repulsive Hubbard models. In this framework, these efforts are as important as the theoretical studies of Cu ladders in the 1990s 13,15 : if pairing occurs convincingly in 1D systems, the same Hamiltonian may induce analogous tendencies in higher dimensions where many-body techniques are not as accurate.

Model
We use a canonical two-orbital Hubbard model with kinetic energy and interaction terms written as H = H K + H I + H D . The tight-binding portion is where c y iγs (c iγσ ) creates (destroys) an electron at site i of a chain, orbital γ (a and b in our case, although our Hamiltonian notation is generic for arbitrary number of orbitals), and spin projection σ. The nearest-neighbor (NN) electron hopping is here a 2 × 2 orbital-space unit matrix, that is, t γγ′ = t δ γγ′ , with t the energy unit throughout the publication. The non-interacting bandwidth is W = 4.0t. The hopping symmetry between the two orbitals, and the absence of crystal-field splitting, prevents the appearance of the orbital-selective Mott physics recently studied in related multiorbital models [52][53][54] .
The electronic interaction is standard for multiorbital fermionic systems 55 : The first term is the intra-orbital Hubbard repulsion U. The second contains the inter-orbital repulsion at different orbitals, with the usual relation U′ = U -2J H due to rotational invariance. The third term involves the Hund's coupling J H , and the last term represents the on-site inter-orbital electron-pair hopping Later it will also be important to incorporate an easy-plane anisotropy component (D > 0): The spin-1/2 operators (S x , S y , S z ) are defined as S α iγ ¼ 1 2 À Á P σσ 0 c y iγσ σ α σ;σ 0 c iγσ 0 via Pauli matrices. For our results, we used the Lanczos method as well as DMRG, with up to m = 1800 states and truncation errors below 10 −6 as in previous investigations 45 .
Undoped two-orbital Hubbard model at intermediate U/W vs. Haldane state We focus on multiorbital models in iron-based superconductors where ladders and chains can be synthesized, but our results are valid for other transition metal compounds. Iron superconductors are "intermediate" between weak and strong coupling, and U/W ≈ 1 is considered realistic [4][5][6] . Because the iron family is not at U/W » 1, a pure spin model is not appropriate and interacting itinerant fermions must be used.
Consider first whether the model discussed here-with mobile electrons, intermediate U/W, and hopping unit matrix-is smoothly connected to the Haldane limit. At one particle per orbital and U/W » 1-with concomitant growth of J H fixed at the often used ratio J H /U = 1/4 5,56 -our model certainly develops S = 1 states at every site, antiferromagnetically Heisenberg coupled. To analyze if intermediate U/W ≈ 1 and strong coupling U/W » 1 (with S = 1 spins on-site) are qualitatively similar, we compute with DMRG the entanglement spectra (ES) 57 . For example, at U/W = 1.6 where hole binding is maximized (see below), Fig. 2a, b indicates that increasing J H /U the Hubbard ES clearly resembles the S = 1 chain ES 58 . However, our model is not merely a S = 1 chain: the inset of Fig.  2c indicates that the von Neumann entropy 59-61 S VN converges to ln(2) (S = 1 chain result) only at U/W ≈ 5 and beyond. At typical couplings of iron compounds, S VN is approximately double the U/ W » 1 limit. Thus, the two-orbital Hubbard model qualitatively resembles the Haldane chain, but at U/W ≈ 1 there are quantitative differences likely caused by non-negligible charge fluctuations.
Consider now the evolution increasing D/t. Recent work found a transition between the gapped Haldane region and a gapped state with trivial topology 46 . In Fig. 2c indeed S VN at fixed U/W = 1.6 and J H /U = 0.25 does not evolve smoothly from D/t = 0connected to the large U/W Haldane limit-to the anisotropic large D/t "XY" limit. The ground state in this limit has a spin triplet with zero z-projection at every site, and no edge states. At 0.1 < D/ t < 0.2, an abrupt change occurs and eventually S VN → ln(1) as D/t grows, compatible with a product state of zero z-projection triplets [see discussion below, Eq. (4)].
In summary, although with quantitative differences, the undoped Hubbard model qualitatively resembles the Haldane chain as long as D/t does not cross a threshold beyond which edge states disappear and a topologically trivial regime develops.
Pairing in the doped two-orbital Hubbard model Our main focus is why pairing occurs and why in the channel it occurs. However, before addressing these issues, let us review and extend recent studies about hole-pair formation and pair-pair correlations in the doped two-orbital Hubbard model. This analysis will provide hints for the intuitive explanation. In Fig. 3a the twohole binding energy vs. U/W is shown, parametric with J H /U. This binding energy is defined as ΔE = E(2) − E(0) − 2[E(1) − E(0)], with E(M) the ground-state energy with M holes (zero holes refers to the half-filled state with one electron per orbital). When ΔE becomes negative, it signals a two-hole bound state. Clearly, Fig. 3a indicates pair formation with maximum |ΔE| at 1<U/W<2, as in ref. 45 , and growing with increasing J H /U (note J H /U should Fig. 1 Summary main results. a Sketch of a chain with two active orbitals a and b. b Representation of a splitting the orbitals into legs forming a fictitious two-leg ladder, with legs only connected by the Hund coupling J H . c One component of the variational state proposed in the text. Arrows indicate spin-1/2 singlets linking nearest-neighbor sites. Although their spin is zero, they are oriented objects because singlets are antisymmetric under the exchange of spins. The full ORVB state is an optimized linear combination of all possible arrangements of these singlets, that is, a linear combination of AKLT valence-bond solids. d Doped state: holes "h" are effectively paired when a spin singlet is removed. be <1/3 to remain smaller than U′/U due to the constraint U′ = U − 2J H ). At U/W » 1, ferromagnetism for 1 and 2 holes-see "Discussion" below-prevents pairing, suggesting that directly doping the S = 1 chain is not the proper theoretical approach. In Fig. 3b we show new results, now increasing D/t at fixed J H /U = 0.25. Robust pairing is observed again. However, while the binding curves are almost identical at D/t = 0.0 and 0.2, at larger D/t they rapidly increase in magnitude. This reflects qualitative differences in pairing, compatible with the von Neumann analysis increasing D/t in Fig. 2c that indicated a topological change in the same D/t range.
The qualitative transition in Fig. 3b also occurs in Fig. 3c, d where pairing correlations are shown. At D/t = 0 and J H /U = 0.25, that is, doping a region smoothly connected to the Haldane chain, spin-singlet pairing dominates (triplet is exponentially suppressed). With increasing D/t at fixed hole density x, a transition from singlet to triplet dominance is observed. For example, in Fig. 3d we observe that spin-triplet pairing, heavily suppressed at D/t = 0, instead dominates as D/t increases (while singlet is exponentially suppressed).
The "global summary" in Fig. 4 is based on a large set of DMRG data. It contains a phase diagram varying D/t and doping x, with only a few representative points displayed. The red region near x = 0 and D/t = 0 is where the model resembles the Haldane state according to the entropy entanglement. Here, at light doping x singlet pairing dominates at intermediate U/W, but eventually other non-superconducting channels [spin density wave (SDW) and charge density wave (CDW)] take over as x grows. Increasing D/t, at small x a transition from singlet-to tripletdominated pairing occurs. In the singlet regime, holes are primarily located at NN sites and different orbitals, while in the triplet regime they are primarily at the same site in different orbitals.
The DMRG results unveiled a parameter space region (small D/t, low hole doping, intermediate U/W, and robust J H /U) where superconducting spin-singlet correlations dominate. These results are surprising. First, the connection with the S = 1 chain suggests that antiferromagnetic (AFM) fluctuations are short range and perhaps not sufficiently strong for pairing. Second, at every site and at intermediate-strong U/W, a nonzero magnetic moment develops due to J H /U. Naively, these same-site electrons can be considered as local preformed triplets. Then, after hole doping the resulting ground state could be envisioned as these triplets becoming mobile. For D/t > 0.2, this naive perspective is compatible with numerical results in Fig. 4. However, better understanding the on-site spin-triplet pairing will require further work because in principle a Hund coupling J H /U = 1/4 is not sufficient to overcome the inter-orbital repulsion. The anisotropy D, which influences on the energy, seems important to stabilize the triplet pairing as our computational results indicate. The product-state wave function Eq. (4) is a good variational approximation (undoped system), exact as D diverges: Then, why singlets dominate at small D/t? Although in a Haldane regime all triplet correlations must decay exponentially, such reasoning does not explain why singlet pairing is enhanced. Hints for the variational state presented below arise from the AKLT exact solution 3,48,49 , where a S = 1 spin model was considered employing two auxiliary idealized S = 1/2 degrees of freedom at every site. These auxiliary states form spin singlets with other S = 1/2 auxiliary states at the next site. Below, we show that the two-orbital Hubbard model shares properties similar to this intuitive idea.
Variational state for the undoped two-orbital Hubbard chain We now introduce a variational state for both the undoped and lightly doped two-orbital Hubbard chain, at intermediate and strong U/W. We argue that these ground states can be qualitatively described in terms of S = 1/2 spin singlets involving NN sites, connecting the same or different orbitals. Our main result is that the small D/t region of the two-orbital Hubbard model has hidden "preformed" singlets that become mobile with doping.
Knowing what type of Hubbard model develops pairing, and with what type of hoppings, allow us to predict what characteristics a material must display to realize this physics. We propose a variational state inspired by an exact equality. Consider first only two sites, say 1 and 2, and construct the quantum global spin-zero state using one electron per orbital. With only one orbital this has the canonical expression: where |1/2,1/2〉 means total spin-1/2 and z-projection ↑, and so on. For two spins 1, the global spin-zero state is still relatively simple Because in our case each site S = 1 arises from two real S = 1/2 electrons at each orbital and same site, we use the 1; 1 notation and an analogous expression at site 2. Then, simple algebra leads to What is remarkable is that this last expression can be exactly rewritten as a combination of S = 1/2 singlets, involving either different or the same orbitals: Intuition is gained when this result is represented visually (Fig. 5), where we have rewritten exactly Eq. (7) doubling the number of valence-bond states for an easier extrapolation to more sites (using the total spin at each site S i = S ia + S ib , the identity (S ia + S ib ) · (S ia + S ib ) = S ia · S ia + S ib · S ib + S ia · S ib + S ib · S ia also helps in this context).
Our exact result is counterintuitive: with perfect S = 1 states at each site, the total spin-zero state of the two-orbital two-site Hubbard model can be represented exactly as a linear combination of S = 1/2 singlets. This resembles the original AKLT perspective 3 , although here applied to a fermionic system. Fig. 5 is as in the views of Anderson 47 and Affleck et al. 3 , but with orbitals as legs of a two-leg ladder, with these legs only connected by J H (no inter-leg hopping).
Then, intuitively, as in the AKLT states, the undoped state has preformed S = 1/2 singlet pairs in all possible arrangements that upon doping should become mobile, leading to spin-singlet pairing dominance. Thus, we predict that doping real quasi-one-dimensional materials with two dominant nearly degenerate orbitals should lead to superconductivity in the spin-singlet channel if anisotropies are not large. We need two "similar" orbitals because we used a 2 × 2 unit hopping matrix.
How is this generalized to more sites? The two-site exact result in Fig. 5 assuming periodic boundary conditions (PBCs) establishes a rule: at each elementary 2 × 2 plaquette only one singlet can be used, either along a diagonal or along a leg. Each ladder site S = 1/ 2 can be used only once: after forming a singlet they disappear from the picture. The two-site example has other properties common to a longer chain. The global singlet state is even under the exchange of orbitals a and b and also even under a reflection with respect to the middle of the plaquette. Now extending to more sites becomes natural. For example, in Fig. 6 we show the 16 valence-bond states needed for four sites using PBCs, as well as the "representative" of each class (i.e., applying to a representative translations and orbital exchange, the full original class can be reconstructed). We remark again that, by construction, all states have perfect S = 1 spins at every site, as in the AKLT setup, even using spin-1/2 singlets as building blocks.
Lanczos overlaps How accurate is this state? Using Lanczos, we calculated the normalized exact ground-state GS of the two-orbital Hubbard model in short chains and computed the overlap with the ORVB linear combination of the individual AKLT-like states of Fig. 6. The coefficients for each class were optimized to maximize the global overlap, arriving to a final normalized-to-one state-dubbed ORVB. Care must be taken because the individual AKLT components do not form an orthogonal set. By this procedure, at U/W = 20 and x is the number of holes divided by 96 (48 sites, 2 orbitals). For completeness, we repeated several ΔE/t calculations in a, b removing the electron-pair hopping term from Eq. (2) to avoid the impression that this term, arising from Coulomb energy matrix levels 55 , may cause the binding. In all cases studied, the binding curves were barely affected by removing the electron-pair hopping. 2 sites, |〈ORVB|GS〉| is virtually 1, because double occupancy is much suppressed, as in Eq. (7). Using four sites, the binding energy is now optimized at U/W = 4. At this coupling and size, | 〈ORVB|GS〉| = 0.95, indicating that ORVB state is a good representation of the ground state.
We extended to more sites using PBCs. As already explained, below when "classes" are mentioned for the undoped case, they represent groups of valence-bond states related by applying translations and orbital exchange to a particular representative. For six sites, there are eight classes and at U/W = 2, where binding is maximized, the overlap is 0.79. As the lattice size grows, other configurations involving longer S = 1/2 singlets and especially doubly occupied orbitals will contribute to the ground state because the optimal U/W where binding is maximized is reduced towards the intermediate range. But finding a robust 0.79 overlap with six sites indicates that ORVB is a good variational state. The same occurs for eight sites: here the number of classes is 16, and when U/W = 1.5 is chosen because it optimizes the binding energy, the overlap |〈ORVB|GS〉| remains robust at 0.61. Increasing the system size, the optimal binding converges to intermediate U/W (see Fig. 3).
Should we worry about a reducing overlap with increasing size? As example consider the simple (π,π) spin staggered state stagg j i ¼ "#"# ::: j i of the S = 1/2 Heisenberg model twodimensional square lattice. For a 2 × 2 cluster, its overlap with the true ground state is 0.58, but for the 4 × 4 cluster, it decreases to 0.29. However, stagg j i is certainly a good variational state. A reducing overlap is natural because with increasing cluster size, the fraction of the total Hilbert space spanned by simplified variational states-such as proposed here for two orbitals or the spin staggered state for Heisenberg models-rapidly decreases. Thus, after confirming the overlap is robust for small clusters, what matters more is whether the proposed state captures the essence of the ground state, as shown in the next.
Doped variational state and superconductivity Let us generalize our variational state to the doped case, a topic barely addressed in the AKLT context. In our DMRG studies in Fig.  3 and in previous efforts 45 , we found that in the two-hole ground state, the largest-weight configuration occurs when holes are placed at NN sites and in different orbitals. Having the two holes in the same leg is not optimal because they collide: with one hole per leg they can move without obstacles, while taking advantage of the effective attraction in the variational state arrangement. For this reason, and to reduce complexity, our proposed doped state will have only one hole per orbital and will be obtained primarily from the undoped state by removing one diagonal singlet. This is exemplified for four-site PBCs in Fig. 7 (left) where classes are shown. Note that the diagonal character of the 2 × 2 hopping matrix establishes that the number of holes per orbital is conserved.
Additional remarks: (1) here we use the unit hopping matrix, but in most realistic situations a nonzero crystal field among the orbitals (rendering them non-equivalent even after a change of basis) as well as inter-orbital hoppings and non-equal diagonal hoppings could be present. What occurs in these conditions remains to be studied. (2) While at very large U/W the AKLT guidance should work well, at intermediate U/W the form of the orbital hopping matrix influences on the energy. Then, in the ORVB Ansatz, the states with all inter-orbital singlets will not have the same weight as the states with all intra-orbital singlets.
Two extra ingredients are needed. First, quantum mechanically each hole in each orbital is "oscillating" (zero-point motion) via the intra-orbital hopping because these are not frozen holes. Then, the configuration with two holes in the same rung must be included because it is generated by oscillations within the bound state. Second, to avoid unpaired S = 1/2 electrons left and right of that two-hole rung, a singlet across is required, as in the three classes in Fig. 7 (right). Indeed, a π-shift across-the-hole develops in the spin correlations of the two-orbital Hubbard ground state 45 . This also occurs in one-orbital t-J models 62,63 . For completeness, singlets across-the-hole were also added for diagonal hole configurations as in the bottom left class of Fig. 7.
This procedure resembles qualitatively the exact solution of the infinite U single-orbital Hubbard chain 64 : holes and spins are independent in this limit. Our mobile holes can be visualized as effectively inserted in between the original singlets of the undoped valence-bond state.
The ORVB state from Fig. 7, but now for six sites and two holes, 23 classes are required [for two holes, to generate all states not only translational symmetry and exchange of orbitals are needed but also reflection (parity) with respect to the middle]. The overlap with the two-hole Lanczos exact GS at U/W = 2 is 0.59. For eight sites and two holes, 84 classes are needed, and the overlap at  U/W = 1.5 is 0.48. These are good numbers, but more important is how qualitatively these states capture the essence of the problem. For instance, from the six-site exact GS the spin-spin correlations can be measured using special projection operators 45,65 for when the mobile holes are at their highest-probability ground-state position 33,45 . The Lanczos and our variational results are contrasted in Fig. 8a, b. The agreement is remarkable. While away from the holes, the pattern resembles the undoped case, with ferromagnetic (FM) rungs and AFM legs 45 , near the hole the AFM correlation "across-the-hole," typical of carriers in an AFM background, is reproduced. Moreover, a puzzling FM link diagonally placed opposite to the diagonal of holes is also observed. Naively this may suggest that triplets are needed in the undoped variational state along diagonals. However, this effective FM correlation is merely a consequence of the mobility of the holes that displace the original on-site FM triplet-contained in our variational state by construction because of the AKLT projection to spin-1 at every site-from rung to diagonal. Figure 8c also shows that at very large U/W our proposed picture breaks down. In this regime, the effective AFM superexchange weakens because it scales as 1/U. As a consequence, alternative tendencies such as "double exchange" as in the manganite context 55,66,67 are enhanced, leading to ferromagnetism to improve the kinetic energy of the now unbounded holes. This suggests that simply doping the S = 1 Haldane chain may not be sufficient, but U/W must be limited to intermediate values to Fig. 6 Four-site variational. Individual states used for the four-site, half-filled two-orbital Hubbard model with PBCs. The color-framed states represent four distinct "classes." A representative of each class is shown at the bottom with the same color convention. By applying translation and orbital exchange for each representative (reflection is not needed for the undoped case), we recover all the states in the upper frames. Note that the 16 AKLT states displayed are not orthogonal to one another. Our proposed variational state is an energy-optimized linear combination of these 16 states (i.e., the weight of the four classes is different). Blue arrows are spin singlets, see Fig. 5.  Here the probability of single occupancy of one orbital is 96%, indicating that local moments S ≈ 1 are well formed. However, the Heisenberg limit is only reached at U/W > 5 (Fig. 2c, inset). This suggests that other terms in the strong coupling expansion are of relevance in the regime of binding.
avoid the ferromagnetic competition. However, in the reported hole-binding range the total spin quantum numbers are qualitatively compatible with those of our variational picture. These optimal values of U/W are also compatible with our ladder work 33 and with alternative explanations focused on pairing amplitudes and effective exchange interaction optimizations 68,69 .
Varying D/t When easy-plane anisotropies are included, the spin-triplets product state at each site TPS, Eq. (4), becomes increasingly a better approximation as D/t grows (at D/t diverges, TPS is the exact ground state). This evolution is illustrated in Fig. 9. In Fig. 9a, the half-filling overlaps |〈ORVB|GS〉| and |〈TPS|GS〉| are shown. The ORVB (TPS) overlap decreases (increases) with increasing D/t, as expected. Note that the TPS state used here is crude: larger overlaps in the D/t range shown could be obtained adding fluctuations, but this is irrelevant for our main focus, that is, the origin of the spin-singlet pairing at small D/t.
In addition, we observed that the spin-spin correlations involving the z components within the two orbitals of the same rung correctly evolve with increasing D/t. At small D/t, they are FM because a spin close to S = 1 forms at each site. At large D/t, they become AFM because only the spin-zero z-projection component survives in the XY product state. Also note that the ORVB and TPS states are not orthogonal to one another, but their overlap is very small (0.02 for N = 8, 0.07 for N = 6, both at D = 0 PBCs).
The case of two holes (Fig. 9b) is more interesting. Here a level crossing occurs (Fig. 9c): at small D/t the ORVB two-hole-doped state has a nonzero overlap with the two-hole Lanczos state because both states have quantum number (−1) under orbital exchange (orbital antisymmetric). However, at D/t ≈ 0.125, where the von Neumann entropy (Fig. 2c) signaled a qualitative transition in the undoped case, a level crossing occurs in the Lanczos ground state. At large D/t the quantum number under orbital exchange becomes +1, leading to a nonzero overlap with the orbital symmetric two-hole TPS state.

DISCUSSION
We introduced a variational state for the undoped and holedoped two-orbital Hubbard chain, verified its accuracy, and explained the development of spin-singlet pairing upon doping. Our analysis relies on valence-bond states defined in an extended ladder-like geometry spanned by the real chain in the long axis and the number of orbitals in the short axis. Our variational state is an optimized linear combination of AKLT singlets. Using DMRG and Lanczos, we found excellent agreement with our variational predictions at intermediate U/W. The ES and von Neuman entropy indicate that the undoped intermediate U/W regime is connected to the Haldane limit at U/W » 1. However, in the realm of spin models at U/W » 1, a strong competition with ferromagnetism upon doping prevents pairing from occurring.
Our variational state relies on a mathematical expression involving two S = 1 spins at NN sites, with one electron per orbital. The global spin-zero state of these two sites is written exactly exclusively using spin-1/2 singlets linking electrons at those NN sites, involving the same or different orbitals. When extended to more sites, the proposed variational state is an optimized linear combination of spin-1/2 singlets in all possible NN arrangements. Note that same-rung singlets are excluded because of the large ferromagnetic Hund coupling. Using Lanczos, our variational state was shown to be a good approximation to the true ground state for small easy-plane anisotropy D/t, robust Hund coupling J H /U, intermediate U/W, and light hole doping.
The preformed NN spin singlets-Cooper pairs-become mobile upon hole doping, and will form a coherent superconducting state if employing the canonical BCS (Bardeen, Cooper, and Schrieffer)-like product-state construction, at least within the limitations of one dimensionality that only allow for power-law decays. A weak coupling among chains will render the state truly superconducting with long-range order, as in two-leg Cu-or Fe-based ladders. Because of the small size of the Cooper pairs, the coherent state is likely in the Bose-Einstein condensation class.
What occurs if more orbitals are used? The pioneering work of Haldane established that integer and half-integer spin chains are intrinsically different. Thus, if we use three orbitals and still a unit-matrix hopping, a generalized ORVB state with a spin gap due to the presence of spin singlets should not be a good variational state. However, using four orbitals we should return to the class of two. The situation becomes more complicated, and difficult to predict, if in addition to modifying the number of orbitals we also add crystal fields, inter-orbital hoppings, or assign different values in the diagonal for different orbitals. In this case, the subject is totally open. First indications 33 suggest that binding is possible for ladders with a non-trivial hopping matrix, but probably the undoped state is not topological as in the Haldane chain. These many open issues will be addressed in the future.
In summary, our study combines topological concepts with correlation effects. Pairing emerges with hole doping via the liberation of preformed spin-1/2 singlets already contained in the undoped limit. For experimentalists to realize our model, the challenge is to find quasi-1D materials with two nearly degenerate dominant active orbitals, and with similar overlaps along the chain so that the hopping matrix is nearly the unit matrix, as in our calculations. How robust the Haldane regime is with regards to small deviations from this hopping symmetric case, as well as the introduction of an orbital small crystal-field splitting, remains to be investigated (the orbital-selective Mott phase is close in parameter space [52][53][54]. Doping the existing physical realizations of undoped Haldane chains 7,10 can provide a starting point towards the Fig. 9 Overlaps. Overlaps |〈ORVB|GS〉| and |〈TPS|GS〉| vs. D/t for: a zero hole (half-filling) and b two holes. The discontinuity in b indicates a first-order transition due to a level crossing and associated change in the ground-state quantum numbers under orbital exchange; c ground-state chain energies, as well as their derivatives to emphasize sudden slope changes, for two holes E N − 2 . All results (a-c) obtained with Lanczos using N = 6 and 8 sites and PBCs.
predicted superconductivity, but, again, intermediate U/W is a more attractive parameter region than U/W » 1. For this reason, the iron superconductor family provides a natural starting point, although realizations could also be found in other multiorbitalactive compounds.

METHODS Operators
To address pairing, we defined a general pair creation operator as Δ γγ0 y ði;jÞ ± ¼ 1 ffiffi ffi 2 p c y iγ" c y jγ0# ± c y iγ# c y jγ0" ; where i,j are sites, γ,γ′ are orbitals (a or b), and ± sign represents a spin singlet or triplet. We only focused on two different pair operators: (1) NN pair that is odd under spin (singlet) and under orbital exchange (S y nn (i) below), (2) on-site inter-orbital pair that is even under spin (triplet) and under orbital exchange (T y on (i) below). These are defined as In the main text, we often refer to S y nn (i) as odd diagonal singlet and to T y on (i) as on-site triplet, where "diagonal" and "rung" refers to the ladder representation of a two-orbital chain Fig. 1. Using these pair creation operators, we study the decay of the pair-pair correlations Þi; (10) where N R represents the number of total neighbors at distance R with respect to site i, summed over all sites.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon request.