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

The electron–lattice interaction gives rise to a rich set of phenomena in quantum materials. Microscopically, this interaction often arises from the modulation of orbital overlaps; however, many theoretical studies neglect such couplings. Here, we present an exact diagonalization and determinant quantum Monte Carlo study of a three-orbital Su–Schrieffer–Heeger (SSH) model, on a two-dimensional Lieb lattice and in the negative charge transfer regime. At half-filling (one hole/unit cell), we observe a bipolaron insulating phase with a bond-disproportionate lattice. This phase is robust against moderate hole doping but is suppressed at large hole concentrations, leading to a metallic polaron-liquid-like state with fluctuating patches of local distortions. We also find an s-wave superconducting state at large hole doping that primarily appears on the oxygen sublattice. Our work provides a non-perturbative view of SSH-type couplings in two dimensions with implications for materials where such couplings are dominant.


INTRODUCTION
Lattice dimerization, referring to an alternation of two lattice constants, occurs in many quantum materials, including the organic charge-transfer solids [1][2][3] , and perovskites systems like the rare-earth nickelates RNiO 3 [4][5][6] , and the high-temperature (high-T c ) superconducting bismuthates Bi 1−x K x BiO 3 (BKBO) 7,8 . Lattice dimerization in one-dimension was initially described by the Su-Schrieffer-Heeger (SSH) model 1 , in which the atomic motion modulates the overlap integrals between neighboring atoms. Recently, SSH-like models have attracted further attention due to the fact that they can produce highly mobile polarons with light effective masses 9 , generate robust phonon-mediated pairing 10 , and even stabilize and control the location of a type-II Dirac point 11 . SSH interactions are also believed to be dominant mechanism for electron-phonon (e-ph) interactions in materials like the high-T c cuprates 12,13 .
The SSH model has primarily been studied in one dimension (1D) 9,10,[14][15][16][17][18][19][20][21][22] . But it is also necessary to study this model in higher dimensions in the context of materials like RNiO 3 and BKBO. In this work, we focus on BKBO and study how the SSH-type e-ph interaction produces both insulating and superconducting states as a function of doping. BKBO is in the so-called "negative charge transfer" regime [23][24][25][26] , 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 6,7 , which may be further enhanced by correlations 27 , and is believed to drive a high-temperature metal-to-insulator (MIT) transition. Here, the insulating state has a dimerized (or "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 6,7,25,26,28,29 (see also "Methods" section). The relevant model describing this situation is a multiorbital SSH model; however, little is currently known about the physics of such models due to a lack of suitable approaches for studying them in dimensions higher than one.
To address this problem, we carried out a combined determinant quantum Monte Carlo (DQMC) and exact diagaonalization (ED) study of a two-dimensional (2D) multiorbital model with SSH-type interactions. In this case, our DQMC simulations are made possible by the development of a fast update procedure for local spacetime moves of the phonon fields (see "Methods" section). Having BKBO in mind, we then deployed this approach to study a three-orbital SSH model defined on a Lieb lattice whose orbital basis consists of a Bi 6s orbital and O 2p x and 2p y orbitals. In this case, we freeze the heavier Bi atoms into place and restrict lighter O atoms to move along the bond directions, and study the model using non-perturbative DQMC and ED. Further details about the model, and the changes necessary for the DQMC algorithm are provided in the "Methods" section and in the supplementary materials (see Supplementary Note 1 for more detailed results).

RESULTS
An overview of the phase diagram An overview of the phase diagram inferred here is presented in Fig. 1. Near half-filling (one hole/unit cell), we find that the system is a bipolaronic charge-density-wave (CDW) insulator with a bonddisproportionated structure, similar to what is observed in BKBO 30 . Hole doping suppresses the insulating phase, giving way to a state where the lattice distortions have short-range correlations. At high doping levels, we find evidence for a metallic phase where holes are strongly correlated with local structural distortions, forming a polaron-liquid-like phase. Finally, at low temperatures, we find s-wave superconducting tendencies that form primarily on the oxygen sublattice. 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. Our simulations are performed in 2D, and quantitative differences may appear when simulating larger systems for the same model in three-dimensions (3D). But we believe that the underlying physics of the system 1 should be the same in both cases. In what follows, we discuss the numerical results leading to this phase diagram and what it means for our understanding of quantum materials dominated by SSHtype interactions.
A molecular orbital viewpoint Before proceeding to our DQMC results, we present 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. 25 for a similar discussion of the 3D case from an ab initio perspective. (Note that ref. 25 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. 2. 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 . Next, we transform the four ligand oxygen orbitals into a molecular orbital basis (see "Methods" section). We also perform an analogous transformation for the phonon operators. Fig. 3a-d indicate the phases of the ligand 2p δ (δ = x, y) orbitals for each molecular orbital using ± signs, and the black arrows indicate the displacement patterns of the transformed phonon eigenmodes. The bond disproportionated structure that forms in the model corresponds to a coherent state of the optical x r;Ls phonon modes in this representation, while the x r;Lx and x r;Ly modes form the basis for the acoustic phonon modes.
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. 3e. 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 bonding (sL s ), nonbonding ðsL s Þ 0 , and antibonding ðsL s Þ Ã states. Here, the energy of the bonding state is lowered by 2t sp relative to the atomic values such that the two holes fill this state at halffilling, as shown in Fig. 3f. This ground state charge distribution is analogous to the one inferred for 3D bismuthates in ab initio calculations 25 and ARPES measurements 26 .
The impact of the e-ph coupling is also evident from the form of the transformed 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 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 transformed 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 4 summarizes the ED results. Figure 4a and b plot the evolution of the hole density hn Lγ i on each molecular orbital, and the displacement fluctuations of each eigenmode δðx δ Þ ¼ hx 2 L δ i À hx L δ i 2 , respectively, as a function of the hole concentration. As expected, holes primarily occupy the L s orbital at half-filling. (The missing hole weight is split between the two Bi sites and is not shown.) At the same time, the displacement of the x Ls mode fluctuates significantly, while the remaining eigenmodes have fluctuations consistent with zero point motion. This behavior is also reflected in the expectation value of the phonon numbers ( Fig. 4c), where the x Ls mode is excited while the remaining phonon modes are in their ground state. Note that we do not observe a distortion hx Ls i ≠ 0 due to the absence of any symmetry breaking in the cluster; however, our DQMC simulations performed on larger lattices do find such a state.
When additional holes are introduced they occupy the L d , L x , and L y molecular orbitals, as expected based on the level diagram shown in Fig. 3g. 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 expectation value of the number of phonon quanta for these modes grows. Both the  Fig. 1 The temperature-hole concentration (T − n) phase diagram of the 2D three-orbital SSH model inferred from our results. We observe charge-density-wave (CDW), superconducting (SC), and polaronic metal phases. The solid square and triangular symbols indicate the transition temperatures estimated for the CDW and SC phases (see main text). The open circles indicate all points where explicit DQMC calculations were performed. Here, we have scaled temperature axis by a factor of 4.5 to be consistent with the phase diagram of BKBO 26 . Our ED calculations suggest that hole doping induces a relaxation of the breathing distortion of the lattice that 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 hN p i and bipolaron hN bp i number operators (see the "Methods" section). Figure 4d plots the doping evolution of the hN p i and hN bp i. 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   zero point motion  where the dimerized structure has been introduced by modifying the hopping integrals as t ij sp ¼ t sp ½1 þ ðÀ1Þ iþj 0:3. Figure 5a and b provide fat band plots of the L s and L d molecular orbital weight, respectively, while Fig. 5c plots the total and orbitally resolved density of states (DOS). As can be seen in Fig. 5a, 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 modes.

DQMC simulations of an extended lattice
The molecular orbital picture presented in the previous section provides an intuitive way of understanding the physics of the model. With this in mind, we now turn to DQMC simulations on an extended cluster with N = 4 × 4 Bi atoms (48 orbitals in total).
We first examine the insulating phase that forms at hni ¼ 1. Figure 6b-d show the lattice displacement correlation functions hX r;xX0;x i; hX r;yX0;y i, and hX r;yX 0;x i, as a function of position at temperature ðβt sp Þ À1 ¼ 0:1, which shows evidence of a static bond disproportionated structure. For example, both hX r;xX 0;x i and hX r;yX 0;y i alternate in sign following a checkerboard pattern while hX r;yX 0;x i alternates in sign along x-and y-directions but is constant along the diagonal. These results are consistent with the breathing distortion pattern sketched in Fig. 6a, as well as the observed lattice distortion that appears in the insulating phase of the bismuthates [31][32][33] .
The formation of charge order in the insulating phase can also be observed in the charge susceptibility χ C γ;γ ðqÞ. Figure 7b plots the temperature evolution of χ C γ;γ ðqÞ at q = (π, π), corresponding to the real space ordering inferred from Fig. 6. Below ðβt sp Þ À1 ¼ 0:2, the charge correlations rapidly increase on the s orbital, while there is little change in the signal on the p orbitals. The transition temperature for the CDW can be estimated using the temperature at which the correlation length equals the lattice size (see Supplementary Note 3 for more detailed results). At half filling, the transition temperature is about ðβt sp Þ À1 ¼ 0:15. We stress that such an estimate is likely an overestimate, as finite size effects can be quite strong on small clusters 34 .
The behavior observed in Fig. 7 implies that the average density on the O sublattice remains uniform above and below the MIT transition, while a density modulation forms on the Bi sites. We confirm this picture in the inset of Fig. 7b, which plots the equal time Bi-Bi charge density along the high symmetry directions of 1 2 3 Fig. 6 The dimerized structure at half filling. a A sketch of bond disproportionated lattice structure. The red and blue dots indicate the s and p x,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 hX r;xX 0;x i and hX r;yX 0;y i as a function of distance r = n x a + n y b, respectively. Here, a = (a, 0) and b = (0, a) are the primitive vectors along x-and y-directions, respectively. Panel d plots the real-space displacement correlation function hX r;yX 0;x i indicating the two-sublattice structure of the bond disproportionated state. The distance between two nearest Bi atom in the undistorted square structure is a. The inverse temperature is ðβt sp Þ À1 ¼ 0:1. The error bars represent 1σ statistical errors, which were obtained from the sample variances measured during Monte Carlo sampling.
the cluster, where we find a clear (π, π)-density modulation. The fact that the charge order signal appears in the Bi orbital component can be understood once one recognizes that all of the oxygen orbitals in the system are equivalent, even in the bond disproportionated structure. In this case, the breathing distortions increase the hybridization between the L s orbital and one of the Bi s orbitals at the expense of the other, creating a density modulation on the two Bi sites. From a charge disproportionation point of view, however, the density modulation shown in the inset of Fig. 7b corresponds to a charge transfer of only~0.1 holes between the two Bi sites. Next, we study effects of hole doping on the MIT at fixed temperature ðβt sp Þ À1 ¼ 0:1. Figure 8a plots σ(β/2) as a function of filling, where it increases upon hole doping until hni % 1:5, and after which it levels off within error bars, indicating metallic behavior. Moreover, we also find evidence for the formation of mobile polarons in this region, where holes are bound to local breathing distortions of the oxygen sublattice. This behavior can be seen by examining the polaron operatorpðrÞ ¼x r;Ls ðn r;s þn r;Ls Þ, which is analogous to the operator considered during our ED analysis. Figure 8b plots the doping evolution of the number of polarons given by 1 N P r hpðrÞi, where it decreases as additional holes are introduced but remains nonzero even at the largest dopings. We conclude that a finite number of polarons are present in the system at all dopings, which are formed from holes in the L s orbitals and local breathing phonons.
Considering our DQMC and ED results, we suggest that the suppression of long-range polaron and bipolaron correlations with hole doping should be induced by introducing other phonon modes rather than directly suppressing the breathing phonon mode. We also note that the orbital character and electronic structure of BaBiO 3 computed with ab initio methods share many similarities to the plot shown in Fig. 5. This fact suggests that our analysis can be extended to the bulk 3D material.
To explore the doping evolution of polarons in real space, we measured the staggered polaron correlation function hχ P ðrÞi ¼ ðÀ1Þ rx þry hpðrÞpð0Þi, which is plotted in Figs. 9a-d for selected hole concentrations. At half filling, 〈χ P (r)〉 is positive for all r, indicating that the polarons are frozen into a long-range (on the length scale of the cluster) two-sublattice order, consistent with the patterns inferred from Figs. 6 and 7. With increasing hole concentrations, 〈χ P (r)〉 decreases at the larger distances, indicating an overall relaxation of the bond disproportionated structure but the persistence of short-range correlations. Such behavior may be indicative of nanoscale phase separation 35,36 ; however, studies on large clusters will be needed to clarify this issue. Finally, in the high doping region, where the system is metallic (e.g. hni>1:44), the correlations become very short-ranged and extend up to one or two lattice constants at most. The observation of a finite number of polarons at high doping, but with short-range structure in real space is consistent with the proposal that a significant number of holes remain in the L s orbitals, where they locally excite the breathing modes of the lattice, creating a polaronic metallic phase with fluctuation local breathing distortions.
We also examined the doping evolution of bipolarons in the system by measuring the bipolaron number, defined as 1 N P r hĝðrÞi, whereĝðrÞ ¼x r;Ls ðn r;s;" þn r;Ls;" Þðn r;s;# þn r;Ls;# Þ, and the staggered bipolaron correlation function hχ BP ðrÞi ¼ ðÀ1Þ rx þry hĝðrÞĝð0Þi, as a function of doping. When computing the latter quantity, we considered the signal on the Bi site by  Fig. 7 The evolution of the CDW phase as a function of temperature. a The temperature dependence of the spectral weight at the Fermi level βG(r = 0, τ = β/2) and the conductivity weight σ(β/2). b The temperature dependence of the charge-density-wave susceptibility χ C (π, π). In both panels, the average filling is hni ¼ 1 corresponding to the "half-filled" case with one hole per unit cell. The error bars represent 1σ statistical errors, which were obtained from the sample variances measured during Monte Carlo sampling.  The error bars represent 1σ statistical errors, which were obtained from the sample variances measured during Monte Carlo sampling, and error bars smaller than the marker size have been suppressed for clarity.
keeping only the terms proportional ton r;s;"nr;s;# . (This simplification is necessary due to the many 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 8b 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 number of bipolarons are present in the system. The staggered bipolaron correlation function is plotted in Fig. 9e-h. At hni ¼ 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 entirely suppressed at all length scales, while a finite number of bipolarons are present in the system, as indicated in Fig. 8b. These results can again be easily understood if the metallic phase is a polaron liquid. This scenario raises questions regarding possible superconductivity. We, therefore, computed the pair field susceptibility χ sc γ . Since χ sc p x and χ sc p y are the same, we use χ sc p to denote pairing on the O atoms. Here, we find that the s-wave pairing is dominant, which is perhaps not surprising given that pairing is mediated solely by the e-ph interaction. Figure 10 plots χ sc γ as a function of temperature at hni ¼ 1:59, and compares it against the dominant charge correlations χ C (π, π/2) at this doping. All three susceptibilities increase with decreasing temperature, but χ sc p dominates below ðβt sp Þ À1 % 0:04. The superconducting T c can be crudely estimated by fitting the inverse pairing susceptibility data with a functional form ðχ sc Þ À1 ¼ a ffiffiffiffiffiffiffiffiffiffiffiffiffi T À T c p . This functional form captures the rapid variation in ðχ sc Þ À1 expected as T → T c 34 , while other functional forms such as linear extrapolation would not. Extrapolating ðχ sc p Þ À1 to zero (as shown in the inset), yields T c ≈ 382 K [ðβt sp Þ À1 % 0:0158]. As with the CDW, we stress that this value is artificially high, due to the large value of Ω used in our calculations and likely finite size effects. Nevertheless, our results provide evidence that the bipolaronic-rich metallic phase has a superconducting ground state. The pairing susceptibility is suppressed slightly at lower hole dopings, suggesting the presence of a superconducting dome.

DISCUSSION
We have introduced a quantum Monte Carlo approach for studying bond phonons with SSH-type 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. We obtained a phase diagram consistent with the bismuthate high-T c superconductors. At half filling, we find a bond disproportionated state. Upon    Fig. 10 The evolution of the superconducting correlations as a function of temperature. Results are shown for the charge χ C (π/2, π/ 2) and pair-field χ sc susceptibilities as a function of temperature ðβt sp Þ À1 at hni ¼ 1:59. The inset plots ðχ sc p Þ À1 as a function of temperature ðβt sp Þ À1 . The black dashed line is the fitting result. The error bars represent 1σ statistical errors, which were obtained from the sample variances measured during Monte Carlo sampling.
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 [35][36][37][38][39][40] . 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.

METHODS Model
The unit cell for our multiorbital SSH model contains one Bi 6s orbital and the oxygen 2p x,y orbitals oriented along each of the Bi-Bi bond directions, as shown in Fig. 1a. In what follows, the unit cells are indexed by r = n x a + n y b, where a = (a, 0), b = (0, a) are the primitive lattice vectors along x-and y-directions, respectively, a is the Bi-Bi bond length (and our unit of length). In addition, δ; δ 0 ¼ ±x, ±y will be used to identify the oxygen atoms surrounding each Bi.
In the SSH model, the atomic displacements modulate the hopping integrals between the Bi and O atoms as t sp ðQ δ À αû r;δ Þ, where we have introduced the shorthandû r;x ¼X r;x ,û r;Àx ¼X rÀa;x ,û r;y ¼X r;y , and u r;Ày ¼X rÀb;y . (1) H int ¼ αt sp X hr;δi;σû r;δ s y r;σ p r;δ;σ þ h:c: : Here, 〈…〉 denotes a sum over nearest-neighbor atoms, and the operators s y r;σ ðs r;σ Þ and p y r;δ;σ ðp r;δ;σ Þ create (annihilate) spin σ holes on the Bi 6s and O 2p δ orbitals, respectively. 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 y r;σ s r;σ andn p α r;σ ¼ p y r;α;σ p r;α;σ are the number operators for s and p α (α = x, y) orbitals, respectively. Finally, ϵ s and ϵ p are the onsite energies; μ is the chemical potential; t sp and t pp are the Bi-O and O-O hopping integrals in the cubic crystal; and α is the e-ph coupling constant. The phase factors are Q x = Q y = −Q −x −Q −y = 1, and Q ±x,±y = Q ±y,±x = −Q ±x,∓y = −Q ∓y,±x = 1. The motion of the O atoms is 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 its neighboring Bi atoms. Thus, bare phonon frequency is

Parameters
We adopted t sp = 2.08, t pp = 0.056, ϵ s = 6.42, and ϵ p = 2.42 (in units of eV), which were obtained from DFT calculations of BaBiO 3 7 . We also adopted Ω ¼ ffiffi ffi 2 p t sp and α = 4a −1 . For these parameters, the DQMC calculations r;y i ¼ 0:57ða=4Þ 2 ¼ 0:0356a 2 at half-filling, indicating that the oxygen atoms do not cross the bismuth atoms during the DQMC sampling. Our parameters are therefore in a physically reasonable region. (Here, we are limited to large Ω by long autocorrelation times 41 ). Finally, our DQMC calculations are carried out on a square lattice with N = 4 × 4 Bi atoms (48 orbitals in total).

Transformation to the molecular orbital basis
The molecular orbital basis can be obtained after expanding the unit cell to allow for two distinct Bi 6s orbitals and four O 2p orbitals, as shown in Fig. 2.
DQMC divides the imaginary time interval [0, β] into L discrete steps of length Δτ = β/L such that the partition function can be rewritten using the Trotter formula Z ¼ Tr e ÀΔτLH ð Þ%Tr e ÀΔτHint e ÀΔτK ð Þ L , where K = H el + H lat contains the non-interacting terms of the Hamiltonian and H int contains the e-ph interaction. The Trotter approximation neglects terms of order OðΔτÞ 2 , which is controllable as Δτ → 0. Next, the trace over the phonon momenta and bilinear Fermion operators can be performed analytically to yield a result expressed as a product of matrix determinants 44 and an integral over the remaining lattice displacements Here, ∫dX x and ∫dX y are shorthand for multidimensional integrals over the displacements X r,x,l and X r,y,l defined at each oxygen site and time slice l. The matrix M σ is defined as M σ = I + B σ (L)B σ (L − 1) ⋯ B σ (1), where I is an N × N identity matrix and B σ ðlÞ ¼ e ÀΔτHint e ÀΔτHel . (Each B σ (l) is an N × N matrix and independent of spin.) Note that H int has off-diagonal terms in orbital space due to the nature of the SSH interaction, unlike the case of the Holstein model where it is a diagonal matrix. The lattice contribution to the action is defined as S lat ¼ KX 2 r;x;l þ KX 2 r;y;l þ M 2 X r;x;lþ1 À X r;x;l Δτ 2 þ M 2 X r;y;lþ1 À X r;y;l Δτ 2 : The integrals over the displacements are evaluated using the Metropolis-Hastings algorithm.
Most observable quantities can be expressed in terms of the single particle equal time 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 indices. The determinant of M σ is related to the Green's function M σ ¼ det ½G σ ðlÞ À1 and is independent of l.

Efficient updates of the phonon fields
Equation (19) shows that the Green's function G σ (l + 1) can be obtained from G σ (l) using the identity G σ ðl þ 1Þ ¼ B σ ðl þ 1ÞG σ ðlÞB À1 σ ðl þ 1Þ. This observation forms the basis for an efficient single-site Sherman-Morris updating scheme 42 . The DQMC algorithm starts by computing the Green's function on time slice l = 0 using Eq. (19). A series of individual updates are then proposed by sweeping through the sites (r, α = x, y) proposing updates X r;α;l ! X 0 r;α;l ¼ X r;α;l þ ΔX r;α;l while holding the other phonon fields {X r 0 ≠r;α 0 ≠α;l } fixed. These updates are accepted with probability p ¼ minð1; RÞ, where R ¼ e ÀΔτðSlat½fX 0 r;α;l gÀSlat½fXr;α;l gÞ det½M 0 " det½M 0 # det½M " det½M # ; (20) and M 0 σ and M 0 σ correspond to matrices computed with the new and old phonon field configurations, respectively. Note that the product det½M " det½M # is positive definite for the model considered here, and there is no Fermion sign problem.
After updating a field, the corresponding B σ (l) matrix must also be updated as B σ ðlÞ ! B 0 σ ðlÞ ¼ e ÀΔτH 0 int e ÀΔτHel ¼ e ÀΔτðHintþVÞ e ÀΔτHel , where V contains the terms in H int arising from the change in the phonon field. V is a symmetric matrix with only four non-zero elements of the form ; 0Þie iqÁr is the current-current correlation function and j x ðr; τÞ ¼ Àit sp X δ;σ Q δ À αû r;δ À Á s y r;σ p r;δ;σ À h:c: þ it pp X δ;δ 0 ;σ Q δ;δ 0 p y r;δ;σ p r;δ 0 ;σ ; is the current operator, with phase factors Q δ (as before) and Q ±x,±y = −Q ±y,±x = −Q ±x,∓y = Q ∓y,±x = 1. At low temperature, the conductivity weight is a proxy of the conductivity. More relevant results are shown in the Supplementary Note 2 (see Supplementary Material for more detailed results). A measure of the superconducting and charge ordering tendencies can be obtained from the orbitally resolved charge χ C γ 0 γ ð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 imaginary time,n q;γ ¼ P i;σ e iqÁrin ri ;γ;σ , and r i is the lattice vector. Similarly, the pair-field susceptibility in the swave channel is given by where Δ s = ∑ r s r,↑ s r,↓ and Δ p δ ¼ P r p r;δ;" p r;δ;# .

DATA AVAILABILITY
The data that support the findings of this study will be made available upon reasonable requests to the corresponding author.