Quantum Monte Carlo study of lattice polarons in the two-dimensional multi-orbital Su-Schrieffer-Heeger model

We study a three-orbital Su-Schrieffer-Heeger model defined on a two-dimensional Lieb lattice and in the negative charge transfer regime using determinant quantum Monte Carlo. At half-filling (1 hole/unit cell), we observe a bipolaron insulating phase, where the ligand oxygen atoms collapse and expand about alternating cation atoms to produce a bond-disproportionated state. This phase is robust against moderate hole doping but is eventually suppressed at large hole concentrations, leading to a metallic polaron-liquid-like state with fluctuating patches of local distortions. Our results suggest that the polarons are highly disordered in the metallic state and freeze into a periodic array across the metal-to-insulator transition. We also find an $s$-wave superconducting state at finite doping that primarily appears on the oxygen sublattices. Our approach provides an efficient, non-perturbative way to treat bond phonons in higher dimensions and our results have implications for many materials where coupling to bond phonons is the dominant interaction.

There is an urgent need to address off-diagonal e-ph interactions in higher dimensions, because such couplings are not only relevant to many materialse.g. the organic charge-transfer solids [3,23,29], the rare-earth nickelates [30][31][32], and high-T c superconductors like the cuprates [33,34] and bismuthates [35] -but several recent studies in the few-particle limit have shown that the new physics can occur in such models. For example, strong off-diagonal interactions can produce highly mobile polarons with light effective masses [24], generate robust phonon-mediated pairing [25], and even stabilize and control the location of a type-II Dirac point [36]. Offdiagonal models have also gained attention in relation to 1D topological insulators [37] in the BDI class [38]. It is, therefore, imperative to study off-diagonal e-ph interactions in higher dimensions and at arbitrary fillings, as our intuition gained by studying diagonal models may not serve us well.
Another motivation for studying the SSH-type models in higher dimensions is to better understand its role in establishing the insulating and superconducting states of the high-T c superconducting bismuthates Ba 1−x K x BiO 3 (BKBO). BKBO is in the so-called "negative charge transfer" regime [39][40][41][42], where holes self-dope from the cation to the ligand oxygen atoms. The subsequent hybridization between the cation and the oxygen atoms then leads to a sizable e-ph interaction [32,35], which may be further enhanced by correlations [43], and is believed to drive a high-temperature metalto-insulator (MIT) transition. Here, the insulating state has a bond disproportionated structure with expanded and collapsed BiO 6 octahedra alternating through the material and pairs of holes condensed into the molecular orbitals formed from the ligand oxygen orbitals with A 1g symmetry [32,35,41,42,44,45]. The relevant model describing this case is a multiorbital SSH model; however, knowledge about such models is limited due to a lack of suitable approaches for solving it.
With this motivation, we developed a determinant quantum Monte Carlo (DQMC) method for simulating SSH-type interactions, which is applied to study a 2D multi-orbital model for the first time. At half-filling (one hole/unit cell), we find that the system is a bipolaronic insulator with a bond-disproportionated structure, similar to what is observed in BKBO [46] or the rare-earth nickelates [30]. Hole doping suppresses the insulating phase, giving way to state where the lattice distortions have short-range correlations suggestive of phase mixing and/or fluctuations. At high doping levels, we find evidence for a metallic phase where holes are strongly correlated with local structural distortions, forming a polaronliquid phase. Finally, at low temperatures, we find swave superconducting tendencies that form primarily on the oxygen sublattice and evidence for a superconducting dome. Our results are in qualitative agreement with the phase diagram of the bismuthate superconductors and provide theoretical support for a polaronic view of BKBO and other negative charge transfer oxides. The red and blue dots indicate the s and px,y orbitals, respectively, while the black arrow indicate the displacement pattern of each oxygen atom in the bond disproportionated structure. Panels (b) and (c) plot the lattice displacement correlation functions X r,xX0,x and X r,yX0,y as a function of distance r = nxa + nyb, respectively. Here, a and b are the primitive vectors along x-and y-directions, respectively. Panel (d) plots the real-space displacement correlation function X r,yX0,x indicating the two-sublattice structure of the bond disproportionated state. The distance between two nearest Bi atom in the undistorted square structure is a. whose orbital basis consists of a Bi 6s orbital and two O 2p orbitals, as shown in Fig. 1(a). We freeze the heavier Bi atoms into place and restrict lighter O atoms to move along the bond directions. The Hamiltonian is Here, . . . denotes a sum over nearest neighbor atoms, δ, δ = ±x, ±y index the oxygen atoms surrounding each Bi, and the operators s † r,σ s r,σ and p † r,δ,σ p r,δ,σ create (annihilate) spin σ holes on the Bi 6s and O 2p δ orbitals, respectively. The unit cells are indexed by r = n x a+n y b, The temperature dependence of the spectral weight at the Fermi level βG(r = 0, τ = β/2) and the direct current (dc) conductivity σ dc . (b) The temperature dependence of the charge-density-wave susceptibility χC (π, π). In both panels, the average filling is n = 1 corresponding to the "half-filled" case with one hole per unit cell.
where a = (a, 0), b = (0, a) are the primitive lattice vectors along x-and y-directions, respectively, and a is the Bi-Bi bond length (and our unit of length). To simplify the notation, we have introduced shorthand notation p r,−x,σ = p r−a,x,σ and p r,−y,σ = p r−b,y,σ . The operatorŝ n s r,σ = s † r,σ s r,σ andn pα r,σ = p † r,α,σ p r,α,σ are the number operators for s and p α (α = x, y) orbitals, respectively; s and p are the site energies; µ is the chemical potential; t sp and t pp are the Bi-O and O-O hopping integrals in the undistorted crystal; and α is the e-ph coupling constant. The phase factors are P x(y) = −P −x(−y) = 1, and P ±x,±y = P ±y,±x = −P ±x,∓y = −P ∓y,±x = 1. The motion of the O atoms described by the atomic displacement (momentum) operatorsX r,α (P r,α ). Here, M is the oxygen mass and K is the coefficient of elasticity between each Bi and O atom, and each O is linked by springs to the neighboring Bi atoms. Thus, bare phonon frequency is Ω = 2K/M . Finally, the atomic displacements modulate the hopping integral as t sp (P δ − αû r,δ ), where we have introduced the shorthandû r,x =X r,x , u r,−x =X r−a,x ,û r,y =X r,y , andû r,−y =X r−b,y . We study the model on a square lattice with N = 4 × 4 Bi atoms (48 orbitals in total) using DQMC. We stress that the model considered here is free of the Fermion sign problem. The details are provided in the supplementary materials [47], along with expressions for the standard quantities measured in this work, and supple- mentary exact diagonalization calculations. The details of all non-standard quantities are provided in the main text. Throughout, we adopt t sp = 2.08, t pp = 0.056, s = 6.42, and p = 2.42 (in units of eV), which are obtained from DFT calculations of BaBiO 3 [35]. We adopt a phonon energy Ω = √ 2t sp and e-ph coupling strength α = 4a −1 , which gives average displacement's squared N r X 2 r,y = 0.0356a 2 at half-filling, indicating that the oxygen atoms do not cross the bismuth atoms during the sampling. (Here, we are limited to large Ω by long autocorrelation times.) Results -Figures 1(b)-1(d) plots the lattice displacement correlation functions X r,xX0,x , X r,yX0,y , and X r,yX0,x , as a function of position at inverse temperature β = 10/t sp , which provides evidence for a bond disproportionated insulating state at n = 1. Both X r,xX0,x and X r,yX0,y alternate in sign following a checkerboard pattern while X r,yX0,x alternates in sign along x-and y-directions but is constant along the diagonal. This behavior reflects the breathing distortion sketched in Fig. 1(a), and is consistent with the bond disproportionation observed in the insulating phase of the BKBO [48][49][50]. Figure 2(a) plots the dc conductivity σ dc and orbitalresolved spectral weight βG γγ (r = 0, τ = β/2), where γ is the orbital index, for n = 1 [47,51]. The conductivity (black dots) initially increases as the temperature is lowered until reaching a maximum at β ≈ 5/t sp then it is suppressed. All three orbital spectral weights follow a similar trend, indicating a concomitant removal of spectral weight at the Fermi level. The insulating phase is characterized by a q = (π, π) charge order, as evidenced by the charge susceptibility χ C γγ (q) plotted in Fig. 2(b) as a function of temperature. Below 1/βt ps = 0.2, the charge correlations rapidly increase on the s orbital, while there is little change on the p orbitals. This observation implies that the charge density on the O sublattice is uniform, even in the bond disproportionated structure, while a charge modulation forms on the Bi sites in the insulating state. An examination of the real-space charge density, as shown in the inset of Fig. 2(b), confirms this. We stress, however, that the charge transfer between the Bi sites is on the order of 0.1 holes/Bi. From this analysis, it is clear that the model has a bond-disproportionated structure and a small charge modulation on the Bi atoms in the insulating state. This result supports the bond disproportionation scenario proposed for the bismuthates [41]. We now examine how this state evolves upon hole doping. Here, our focus is on the possible formation of lattice polarons, where holes are bound to local breathing distortions of the oxygen sublattice. These objects can be studied by considering the polaron number operatorp(r) =x r,Ls (n r,s +n r,Ls ), wheren r,Ls = σ L † r,s,σ L r,s,σ is the number operator for the A 1g combination of the ligand oxygen orbitals L r,s,σ = 1 2 (p r,x,σ + p r,y,σ − p r,−x,σ − p r,−y,σ ) [47] and x r,Ls = (X r,x +X r,y −X r,−x −X r,−y ). This operator measures the combined presense of holes in the A 1g molecular orbital surrounding a Bi site and a local contraction of those same orbitals, and can be used to trace the evolution of polarons with doping.
With increasing hole concentrations, we observe a MIT at β = 10/t sp . Figure 3(a) plots σ dc as a function of filling, where it increases upon hole doping until saturating at n ≈ 1.4, indicating metallic behavior. At the same time, the number of polarons 1 N r p(r) decreases as additional holes are introduced but remains nonzero even at the largest dopings [ Fig.3(b)], indicating that the free carriers have polaronic character. We also study polaron correlations in real space using the staggered polaron correlation function P (r) = (−1) nx+ny p(r)p(0) , which is plotted in Figs.3(c)-(f) for selected hole concentrations. At half filling, P (r) is positive for all r, indicating that the polarons are frozen into a long-range two-sublattice order, consistent with the patterns inferred from Figs. 1 and 2. With increasing hole concentrations, P (r) decreases at the larger distances, signalling an overall relaxation of the bond disproportionated on long length scales but the persistence of short-range correlations. Such behavior could reflect nanoscale phase separation [52]; however, studies on large clusters are likely needed to clarify this issue. Finally, in the high doping region, where the system is metallic (e.g. n > 1.44), the correlations become very short-ranged and extend up to at most one or two lattice constants.
We also examined the doping evolution of the bipolaron number, defined as 1 N r ĝ(r) , wherê g(r) =x r,Ls (n r,s,↑ +n r,Ls,↑ )(n r,s,↓ +n r,Ls,↓ ), and the staggered bipolaron correlation function BP (r) = (−1) rx+ry ĝ(r)ĝ(0) , as a function of doping. When computing the latter quantity, we considered the signal on the Bi site by keeping only the terms proportional tô n r,s,↑nr,s,↓ . This simplification is necessary due to the enormous number of terms generated by the Wick contraction of the product ofĝ(r) operators. The fact that we see excess charge density on the Bi sites at the center of a breathing distortion provides some justification for this simplification. Figure 3(b) plots the doping evolution of the bipolaron number operator. As with the polaron number, it is largest near half-filling and decreases slowly with doping. At large hole concentrations, however, it is still finite, suggesting that a significant amount of bipolarons are present in the system. The staggered bipolaron correlation function is plotted in Figs. 3(g)-(j). At n = 1, the bipolaron correlations are clear and long-ranged on the scale of the cluster. This result supports the interpretation that the insulating phase is a static bipolaron lattice. As the hole concentration increases, we find that the bipolaron correlations are suppressed at all length scales, while a finite number of bipolarons are present, as indicated in Fig. 3(b). These results can again be easily understood if the metallic phase is a polaron liquid.
Given the presence of lattice polarons in the metallic phase, we computed the s-wave orbital-resolved pair field susceptibility χ sc γ [47]. Figure 4(a) plots χ sc γ as a function of temperature at n = 1.59, and compares it against the dominant charge correlations χ C ss (π, π/2) at this doping. All three susceptibilities increase with decreasing temperature, but χ sc px = χ sc py ≡ χ sc p dominates below 1/βt sp ≈ 0.04. This observation implies pairing appears predominantly in the oxygen atoms. Extrapolating 1/χ sc p to zero (inset), yields an estimate for the superconducting β c ≈ 63.29/t sp ). This value is artificially high, due to the large value of Ω used in our calculations. Nevertheless, our results provide evidence that the bipolaronic rich metallic phase has a superconducting ground state. Fig. 4(b) plots χ sc p as a function of doping at 1/βt sp = 0.03, where we find that pairing susceptibility is suppressed in proximity to the insulating phase, suggesting the presence of a superconducting dome induced by competition with the insulating phase.
Summary -We have introduced a quantum Monte Carlo approach for studying bond phonons with SSHtype e-ph couplings in higher dimensions. While our approach has broad applications to many materials, we have used it to study a 2D three-orbital SSH model in the negative charge transfer regime for the first time. We obtained several results consistent with the observed properties of bismuthate high-T c superconductors. At half filling, we find a bond disproportionated state that can be viewed as a lattice of localized bipolarons. Upon hole-doping, this state gives way to a polaron-liquid-like state with short-range correlations, consistent with proposals for nano-scale phase separation or strongly fluctuating lattice polarons in doped BKBO [52][53][54][55][56][57]. We also find s-wave superconducting tendencies, which primarily form on the oxygen sublattice. It would be interesting to contrast our results with those obtained from an effective single band model to fully gauge the importance of the oxygen orbitals. Acknowledgments where Z = Tr e −βH is the partition function. We will use Z as an example to illustrate the method since generalizations to the operatorÔ are straightforward [S1]. which is controllable as ∆τ → 0.
Next, the phonon operators are treated by inserting a complete set of position and momentum eigenstates at each time slice. One then integrates out the phonon momenta analytically such that the partition function depends only on a trace over the continuous lattice displacements X r,x and X r,y and terms that are bilinear in the Fermion operators. The trace over the Fermion degrees of freedom can then be evaluated analytically and expressed as a product of matrix determinants [S3]. The final result is where dX x and dX y are shorthand for multidimensional integrals over the displacements 2 X r,x,l and X r,y,l defined at each oxygen site and time slice l. The matrix M σ is defined as The final step is to evaluate the displacement integrals using Metropolis sampling. In this work, we performed both single-site and block updates [S2], as described below.
Most observables can be expressed in terms of the single particle Green's function G σ (τ ).
For an electron propagating through field configurations {X r,x,l } and {X r,y,l }, the Green's function at time τ = l∆τ is given by where A σ (l) = B σ (l) · · · B σ (1)B σ (L) · · · B σ (l + 1),T τ is the time ordering operator, and i, j are combined orbital and site indicies. The determinant of M σ is related to the Green's function M σ = det[G σ (l)] −1 and is independent of l.

S1.2. Efficient single-site updates
Equation (S4) shows that the Green's function G σ (l + 1) can be obtained from G σ (l) using the identity After updating a field X r,α,l , the corresponding B σ (l) matrix must also be updated as where V contains the terms in H e−ph arising from the change in the phonon field. From the Hamiltonian, we can infer that V is a symmetric matrix with only four non-zero elements and it can be written in the form · · · 0 αt 0 sp ∆X r,l 0 · · · · · · αt 0 sp ∆X r,l 0 αt 0 sp ∆X r,l · · · · · · 0 αt 0 sp ∆X r,l 0 · · · . . . . . . . . . . . . . . .
To efficiently calculate the new B σ (l) matrix, we make the approximation which introduces an error on the order of the Trotter error and is valid when ∆τ is small.
The matrix e −∆τ V is then evaluated via e −∆τ V = P e −∆τ D P T , where P is the orthogonal transformation that diagonalizes V , and D is a diagonal matrix with only two non-zero Using these approximations, the Green's function can be efficiently updated after accepting a change in the phonon field using where Q = P T [I − G σ (l)]. Due to the sparsity of matrix ∆, ∆Q has only two non-zero rows where u and w are N × 2 and 2 × N matrices, respectively. Using the Woodbury matrix identity and Matrix determinant lemma, the updated Green's function is given by and the acceptance ratio is given by where I 2 is a 2 × 2 identity matrix. Evaluating these expressions involves O(N 2 ) operations, as opposed to computing the updated Green's function from scratch using Eq. (S4), which has a computational cost of O(N 3 ). Once updates have been performed for the fields on a given time slice l, G σ (l) is advanced to G σ (l + 1) using Eq. (S5) and the process repeated.
This update scheme is efficient but it relies on the approximation B σ (l) ≈ e −∆τ B σ (l). Section S1 S1.4 benchmarks this approach by comparing results obtained using our fast update method against results obtained by explicitly calculating G σ (l). We find that the approximate fast update reproduces the exact solution with comparable error bars for ∆τ 5 values that are typical of most DQMC calculations.

S1.3. Block updates
To reduce the autocorrelation time, we also periodically perform block updates of the phonon fields, where the lattice displacements for a given site are simultaneously updated on all imaginary time slices l. In other words, for a given (r, α), the fields are updated as X r,α,l → X r,α,l + ∆X r,α for all τ l ∈ [0, β]. This type of update efficiently moves phonon configurations out of false minima at low temperatures. There is, however, no fast method for updating Green's function following a block update; it must be calculated using Eq. (S4). As such, the inclusion of block updates slows down the simulations considerably. To strike a balance between the computational time and efficient sampling, we perform between two to four block updates at randomly selected sites for every full spacetime sweep of single-site updates.

S1.4. Reliability of the fast updates
To assess the reliability of the approximation underlying Eq. (S8), we performed two DQMC calculations. In the first calculation, the updated Green's function following each single-site update and the acceptance ratio is computed exactly via Eqs. (S4) and (S6), respectively.
In the second calculation, the fast update procedure is used following the single-site updates.
Both calculations were carried out on a 2 × 2 cluster with a hole density n = 1. We further set β = 6/t sp , ∆τ = 1/10t ps , while the other parameters are the same as in the main text. Figure S1 shows that the Green's functions for both calculations are the same, indicating ∆τ is small enough to reproduce the exact solution with comparable error bars.

S2. Measurements
The dc conductivity is approximated using σ dc = β 2 π Λ xx (q = 0, τ = β/2) [S4], where Λ xx (q, τ ) = r ĵ x (r, τ )ĵ x (0, 0) e iq·r is the current-current correlation function and  = (a, 0) and b = (0, a) are the primitive vectors along the xand y-directions, respectively, where a is the Bi-Bi distance. The first, second, third columns show results for G s,s , G px,px ,vand G py,py , respectively. The red circles and blue triangles represent results obtained using the exact and fast update procedures, respectively. The error bars for the two approaches are comparable.
is the current operator, with phase factors P δ (given in the main text) and Q ±x,±y = −Q ±y,±x = −Q ±x,∓y = Q ∓y,±x = 1.
A measure of the superconducting and charge ordering tendencies can be obtained from the orbitally-resolved charge χ C γ γ (q) and superconducting pair-field χ sc γ susceptibilities, where γ is an orbital index. The charge susceptibility is defined as where q is the momentum, τ is the imaginar ytime,n q,γ = i,σ e iq·r in r i ,γ,σ , and r i is the lattice vector. Similarily, the pair-field susceptibility in the s-wave channel is given by where ∆ s = r s r,↑ s r,↓ and ∆ p δ = r p r,δ,↑ p r,δ,↓ .

S3. A molecular orbital viewpoint
To help understand DQMC results, we also carried out a simplified molecular orbital analysis of a Bi 2 O 4 cluster, which provides a more transparent view of the physics. We refer the reader to Ref. [S5] for a similar discussion of the 3D case from an ab initio perspective.
Note, however, that Ref. [S5] uses electron language whereas we use hole language.
The first step of our analysis is to expand the simple square unit cell to allow for two distinct Bi 6s orbitals and four O 2p orbitals, as indicated by the black dashed frame in Fig.   S2(a). This expanded cell defines the cluster after we apply periodic boundary conditions.
The two Bi 6s orbitals are denoted as s 1 and s 2 .
(The L s and L d operators correspond to the A 1g and E g orbitals in Ref. [S5].) Similarily, we can introduce new phonon operatorŝ x r,Ls = 1 2 (û r,x +û r,y −û r,−x −û r,−y ) x r,L d = 1 2 (û r,x +û r,y −û r,−x −û r,−y ) Here, the sums on γ are taken over γ = s, d, x, y andn Lγ σ = L † γ,σ L γ,σ . We can glean several insights into the problem from this cluster model. In the atomic limit (t sp = t pp = 0) and in the negative charge transfer regime ( p < s in hole language), the four molecular orbitals are degenerate, as shown in Fig. S2(f). This degeneracy is lifted by the orbital overlaps: a nonzero t pp raises (lowers) the onsite energy of the L s (L d ) molecular orbital, while a nonzero t sp hybridizes the Bi s and molecular L s orbitals to form new bonding (sL s ), nonbonding (sL S ) 0 , and antibonding (sL s ) * states. Here, the bonding state's energy is lowered by 2t sp relative to the atomic values such that the two holes fill this state at half-filling, as shown in Fig. S2(g). This ground state charge distribution is analogous to the one inferred for 3D bismuthates in ab initio calculations [S5] and ARPES measurements [S6].
The impact of the e-ph coupling is also evident from this form of the Hamiltonian; holes hop between the L γ molecular orbital and the Bi sites while exciting phonon eigenmodes with the same symmetry. At half-filling, the holes in the (sL s ) bonding state will, therefore, excite the breathing phonon mode of the surrounding oxygen atoms. In an extensive system, this coupling can lead to a static breathing distortion of the lattice after a spontaneous symmetry breaking selects one of the Bi sublattices as the center of the compressed plaquettes. Upon doping, the additional holes will occupy the L d and L x,y orbitals, where they will couple to the orthogonal phonon modes. Since the superposition of the individual modes determines the total displacement of the oxygen atoms, the breathing distortion will relax as the other modes are excited, even though the (sL s ) holes remain coupled to the x Ls phonons.
To confirm this physical picture, we diagonalized the Hamiltonian H M on a Bi 2 O 4 cluster and evaluated several observables in the grand canonical ensemble, with β = 14.56/t sp , Ω = t sp , and µ was adjusted to set the particle number. When diagonalizing this model, we included up to N ph = 5 quanta for each phonon mode, which was sufficient to obtain converged results for our choice of parameters. Figure S3 summarizes the results of our exact diagonlization (ED) calculations. Figure   S3(a) and S3(b) plot the evolution of the hole density n Lγ on each molecular orbital, and the displacement fluctuations of each eigenmode δ(x δ ) = x 2 L δ − x L δ 2 , respectively, as a function of the total filling. As expected, the holes primarily occupy the L s orbital at  When additional holes are introduced they enter the L d , L x and L y molecular orbitals, as expected based on the level diagram shown in Fig. 1(h). In this case, the L d orbital has a larger hole occupation due to the finite value of t pp . At the same time, δ(x L d ), δ(x Lx ), and δ(x L d ) also increase linearly and the number phonon quanta for these modes grows. Both the displacement fluctuations and the number of excited phonon quanta are comparable for the four phonon modes once the hole doping reaches n = 2.4 − 2.6. Finally, the introduction of additional holes slightly suppresses the magnitude of δ(x Ls ) and the total number of x Ls modes.
Our ED calculations suggest that hole doping induces relaxation of the breathing distortion of the lattice, which is dominant at half-filling. However, it does so by exciting the orthogonal phonon modes rather than by suppressing the number of x Ls quanta in the system. In this context, it is interesting then to determine if the L s holes and x Ls modes can be viewed of as a composite object (i.e., a polaron). We checked this idea in our ED calculations by computing the expectation value of the polaron P and bipolaron BP number operators, defined as P = (n s1 +n Ls )x Ls − (n s2 +n Ls )x Ls (S18) and BP = (n s1 ↑ +n Ls ↑ )(n s1 ↓ +n Ls ↓ )x Ls − (n s2 ↑ +n Ls ↑ )(n s2 ↓ +n Ls ↓ )x Ls .  S3(d) plots the doping evolution of the P and BP . We find that the ground state has a significant amount of polaron and bipolaron character, which persists to higher doping levels. Our ED results strongly suggest that the system hosts polaronic carriers, where holes occupying the L s molecular orbitals are bound to local x Ls modes.
The molecular orbitals discussed here will of course form bands in the extended system.
Nevertheless, much of our analysis still applies in this case. To illustrate this, Fig. S4 plots the non-interacting band structure of our model in the insulating phase, where the bond disproportionated structure has been introduced by modifying the hopping integrals as t ij sp = t sp [1 + (−1) i+j × 0.3]. Fig. S4(a) and Fig. S4(b) provide fat band plots of the L s and L d molecular orbital weight, respectively, while Fig. S4(c) plots the total and orbitallyresolved density of states (DOS). As can be seen in Fig. S4(a), the occupied band below the Fermi level (E = 0) at half-filling is the bonding (sL s ) band, which couples to the breathing motion of the lattice. The first band above the Fermi level is mostly of L d and L x,y orbital character, such that doped holes will predominantly couple to the corresponding phonon