A hybrid Monte Carlo study of bond-stretching electron-phonon interactions and charge order in BaBiO$_3$

The relationship between electron-phonon ($e$-ph) interactions and charge-density-wave (CDW) order in the bismuthate family of high-temperature superconductors remains unresolved. We address this question using nonperturbative hybrid Monte Carlo calculations for the parent compound BaBiO$_3$. Our model includes the Bi $6s$ and O $2p_\sigma$ orbitals and coupling to the Bi-O bond-stretching branch of optical phonons via modulations of the Bi-O hopping integral. We simulate three-dimensional clusters of up to 4000 orbitals, with input model parameters taken from {\it ab initio} electronic structure calculations and a phonon energy $\hbar\Omega_0 = 60$~meV. Our results demonstrate that the coupling to the bond-stretching modes is sufficient to reproduce the CDW transition in this system, despite a relatively small dimensionless coupling. We also find that the transition deviates from the weak-coupling Peierls' picture. This work demonstrates that off-diagonal $e$-ph interactions in orbital space are vital in establishing the bismuthate phase diagram.


I. INTRODUCTION
The barium bismuthate families of superconductors Ba 1−x K x BiO 3 and BaBi 1−x Pb x O 3 exhibit hightemperature (high-T c ) superconductivity in proximity to a charge-density-wave (CDW) insulating phase [1][2][3]. The superconducting transition temperature in these materials is largest in the Ba 1−x K x BiO 3 family with T c ≈ 32 K for x ≈ 0.35, which is a record among non-cuprate oxides. These materials remain of fundamental interest due to possible analogies with spin-density-wave-driven exotic superconductivity in unconventional superconductors; nevertheless, the underlying pairing mechanism and its relationship to the proximate CDW phase remains unsettled. Resolving this question is pertinent with the discovery of superconductivity in related (Ba,K)SbO 3 [4] and theoretical proposals for superconductivity in other s-p ABX 3 perovskites [5].
The bismuthates have no magnetic order in their phase diagram [3], which speaks against unconventional spinfluctuation-based pairing mechanisms proposed for high-T c superconductors [6]. Moreover, the extended nature of the Bi 6s orbitals and their significant overlap with the neighboring oxygen atoms suggests that strong electronic correlations play a minor role in these materials [7][8][9][10]. These aspects make the bismuthates a compelling platform for studying high-T c superconductivity in proximity to a charge-ordered state, using uncorrelated models free from issues like the fermion sign problem. Indeed, there has been a recent resurgence of interest in these compounds, using a host of modern thin film syntheses [9,11,12], spectroscopic [9,12,13], and theoretical [10,[14][15][16][17][18][19] techniques. * sjohn145@utk.edu In analogy to the doped Mott insulator picture for the cuprates, several leading hypotheses for the pairing mechanism draw insight from the nature of the neighboring CDW phase. Bi 4+ is a skipped oxidation state that does not typically appear in nature. Early interpretations of the insulating phase, therefore, proposed that the would-be Bi 4+ ions undergo charge disproportionation to (Bi 4+δ Bi 4−δ ) [20][21][22]. In the low-temperature mixed valant phase, Bi-O bond distances then relax in response to Bi charge, resulting in a q = (π, π, π)a −1 breathing distortion of the oxygen atoms and an insulating state. In this case, it has been proposed that the displacement of the highly polarizable O atoms towards the Bi atoms generates a −U center on the Bi site [8]. This scenario would then naturally explain the superconducting state upon doping [21][22][23][24]. Spectroscopic measurements, however, do not support this charge disproportionation view [9,25,26]. This picture is also at odds with the significant hybridization between the Bi 6s and O 2p orbitals and ab initio electronic structure calculations [9,10,14], which find significant oxygen character for the states near the Fermi level in the metallic phase. As an alternative, it has been argued that the bismuthates are best viewed as negative charge transfer materials [27,28], where the system "self-dopes" and transfers one of the Bi holes to a ligand oxygen molecular orbital with A 1g symmetry [10,29,30]. Recent ARPES data for the parent compound support this view of the electronic structure [9]. In this scenario, the A 1g ligand holes couple strongly to the breathing motion of the lattice [30], which induces the insulating phase. Upon doping, the e-ph coupling and potential polaronic effects generate the subsequent pairing [18,19,30].
Finally, it has been argued that the long-range Coulomb interaction in these materials can enhance the e-ph interaction, leading to a more "conventional" BCSlike mechanism [12,31]. The standard BCS scenario was originally ruled out because a naïve estimate for T c based on the McMillian-Allen-Dynes formula suggests that λ ≈ 1 is needed to reproduce the optimal T c = 32 K in Ba 1−x K x BiO 3 . Early theoretical estimates based on density functional theory (DFT) [32] instead placed λ ≈ 0.33 (integrated over all branches), with the strongest coupling occurring to an optical phonon branch dispersing between 66−74 meV. However, recent calculations using more sophisticated hybrid functionals found that the inclusion of the long-range Coulomb interaction could enhance the e-ph coupling, leading to λ between 1 and 1.3 [31]. Similarly, DFT calculations emphasizing the negative charge transfer character of the bismuthates also derive strong e-ph coupling to the q = (π, π, π)a −1 mode of the optical bond-stretching phonon branch and estimated λ ≈ 0.89 [30].
These results call for a re-evaluation of the e-ph interaction in these materials and its role in establishing the CDW phase at low doping. However, given the large λ values estimated above, which are well outside the range where Migdal-Eliashberg theory is expected to be valid [33], it is crucial to perform such an analysis using numerically exact methods capable of treating possible polaronic effects.
Here we study the CDW transition in BaBiO 3 , the parent compound of the bismuthate superconductors, using a scalable hybrid Monte Carlo (HMC) algorithm [34]. We consider a multi-orbital model that includes the Bi 6s and O 2p σ orbitals and coupling to the Bi-O bond-stretching branch of optical phonons via the Su-Schrieffer-Heeger (SSH)-like (or Peierls-like) modulation of the Bi-O hopping integral. A highly efficient numerical code [35] makes possible simulations of large cubic clusters with input model parameters taken directly from ab initio electronic structure calculations [30], and a physically motivated phonon energy Ω 0 = 60 meV [32,36]. Our results demonstrate that an SSH coupling mechanism to the bond-stretching optical oxygen modes accounts for all aspects of the CDW phase in this compound without invoking any additional enhancement of the e-ph coupling constant beyond the early DFT-derived values once the polaronic effects are accounted for. Just as a complete understanding of the antiferromagnetic phase of the cuprates, provided by numerical solutions of the half-filled Hubbard Hamiltonian [37], has been an essential foundation for the rich physics created by doping [6], our work provides for a similar computationally rigorous solution of a model appropriate to the CDW phase of the bismuthates.
The on-site energies for the Bi 6s and O 2p orbitals are s and p , respectively; t sp and t pp are the amplitude of the nearest neighbor Bi-O and O-O hopping integrals, respectively. The overall sign of the hopping integral is given by the phase factors Q δ and Q δ,δ , such that Q ±δ = ∓1, Q ±|δ|,±|δ | = 1 and Q ±|δ|,∓|δ | = −1, where δ = δ . The sign convention enforced by these phase factors is shown in Fig. 1(a).X i,|δ| andP i,|δ| are the position and momentum operators for each oxygen atom, where M is the mass of an oxygen atom, Ω 0 is the bare phonon frequency, and α sp is the e-ph coupling strength.
Throughout, we use tight-binding parameters obtained from ab initio electronic structure calculations [30]. Accordingly, we set the amplitude of the hopping integrals to t sp = 2.31 and t pp = 0.335 eV. The on-site energies are set to s = 6.23 and p = 4.14 eV unless otherwise stated, consistent with the negative charge-transfer nature of these materials [9,10,30]. We adjust the chemical potential to obtain an average filling of n = 1 hole/Bi using a recently developed µ-tuning algorithm [38]. Fig. 1 [eV] As the temperature is lowered, the peak at q = (π, π, π)a −1 grows. (d) At the fixed temperature β = 10 eV −1 , χ0(q) for four values of the charge transfer energy 2p − 6s, the difference in the on-site energy of the O 2pσ and Bi 6s orbitals. The peak at q = (π, π, π)a −1 is approximately independent of the charge transfer energy.
model with α sp = 0. The horizontal dashed red line indicates the location of the Fermi energy at the target filling, which only intersects the lowest energy band. Based on this electronic structure, one might be tempted to adopt a single-band description of the system; however, we will demonstrate that the SSH-like interaction, which is offdiagonal in orbital space, necessitates a multi-orbital description of the problem.
To estimate the e-ph coupling α sp , we assume that the sp hopping integrals depend on the oxygen displacement along the Bi-Bi bond direction X according to t sp (X) = t sp 1 + X a/2 −2 , where a = 4.34Å [39]. By expanding t sp (X) for small X and assuming a linear coupling, we obtain α sp = 4tsp a = 2.13 eV/Å. This value corresponds to a bare dimensionless coupling λ = 0.18 in the standard Migdal-Eliashberg formulation, where Here, N (0) is the density of states per spin at the Fermi energy and g(k, q) is the e-ph coupling constant in momentum space. This value is consistent with early DFTbased estimates [32] without correlation-driven enhancements [31] (see Methods for further details). We solve the model using a recently proposed set of algorithms [34]. Specifically, HMC [40,41] and Fourier acceleration are used to efficiently sample decorrelated phonon fields [42]. A physics-inspired preconditioner reduces the cost of each integration time-step. Efficient measurement of quantum observables is achieved using stochastic techniques that leverage the fast Fourier transform. See Methods for more information. This approach allows us to simulate system sizes up to L = 10 (4L 3 = 4000 orbitals) and consider optical phonons with energies comparable to those in the real system. To this end, we take Ω 0 = 60 meV, which is typical for the bond-stretching optical oxygen phonons in most oxides [32,43]. Note that in all of our simulations, we have confirmed that the displacements of the oxygen atoms remain small enough that the effective hopping never changes sign, ensuring the validity of a linear SSH coupling [44,45].
To detect charge order we measure the static charge structure factor S γ-ν (q) = S γ-ν (q, τ = 0) given by and the corresponding charge susceptibility wheren i,ν = (n i,ν,↑ +n i,ν,↓ ), and γ and ν specify the orbital species. We also calculate the renormalized phonon energy where Π(q, ν n ) is a function related to the phonon selfenergy and ν n = 2πnT is a bosonic Matsubara frequency. versus q for various β, plotted on a semilog scale for an L = 10 cluster. As the temperature is lowered, a peak forms at q = (π, π, π)a −1 . (b) The charge structure factor S6s−6s(q cdw ), sensitive to CDW order on Bi-6s orbitals, versus β. (c) Crossing plot where the vertical axis is scaled using the three-dimensional Ising critical exponents. The vertical dashed black line indicates the transition temperature T cdw = 1160 K. The inset displays the corresponding collapse of the charge structure factor for the L = 8 and 10 data. The error bars in this figure are the standard deviation of the mean of our measurements.
We obtain Ω(q, 0) by measuring the phonon Green's function where the Fourier transform in both r and τ is related to Π(q, ν n ) by [33] B. Charge-density-wave phase transition in BaBiO3 Both formal valence counting and electronic structure calculations [10,32,46] indicate that BaBiO 3 , in its hightemperature metallic phase, has on average one hole/Bi. In this case, the large t sp hopping integral results in a half-filled sp hybridized antibonding band (in electron language), which disperses through the Fermi level, as shown in Fig. 1(b) (hole language). Throughout this paper, we focus on this case and fix n = 1 hole/Bi atom.
According to high-resolution neutron powder diffraction measurements, BaBiO 3 is a CDW insulator over a wide temperature range (4 − 973 K) [47], with the upper bound approaching the material's melting point [48]. The insulating phase is characterized by a q cdw = (π, π, π)/a ordering vector and expansion and contraction of alternating BiO 6 octahedra. This ordering leads to a double perovskite structure, with two inequivalent Bi sites whose average Bi-O bond lengths are largely independent of temperature. We, therefore, begin by studying the CDW transition in our model. In all of our calculations, we find that the O-O correlations remain weak at all temperatures, so no charge modulations develop on the oxygen sites at any temperature. This result is consistent with a bond-disproportionation scenario, where all oxygen orbitals are equivalent at all temperatures [10,30]. Looking at the dressed charge susceptibility on the Bi sites χ 6s-6s (q) versus q, as shown in Fig. 2(a), we see a narrow peak form at the expected ordering wave-vector q cdw . Fig. 2(b), therefore, focuses on the Bi charge structure factor S 6s-6s (q cdw ), which reflects the formation of a charge modulation in which there are distinct local densities on two sublattices of the Bi sites. We find that S 6s-6s (q cdw ) increases with decreasing temperature, with the low-temperature structure factor also growing with system size L at lower temperatures. This behavior is indicative of the formation of long-range order.
As this CDW phase spontaneously breaks a Z 2 symmetry between the two otherwise equivalent Bi 6s sub-lattices, this phase transition falls in the threedimensional Ising universality class. Therefore, we use the corresponding critical exponents γ = 1.237 and ν = 0.630 to perform a finite-size scaling (FSS) analysis to determine the transition temperature [49]. Figure 2(c) displays the relevant S 6s-6s (q cdw )L −γ/ν versus β curves. Here, the locus of crossing points of the data for each cluster size provides an estimate of the critical temperature, provided the clusters are large enough to capture the relevant spatial correlations. Using this data for the crossing of the L = 8 and 10 data sets, we obtain T cdw ≈ 10 eV −1 = 1160 K. This value is above the material's melting point [48] and thus explains why this material remains a CDW insulator at all temperatures. The inset panel displays the corresponding collapse that occurs when S 6s-6s (q cdw )L −γ/ν is plotted against analysis for the CDW transition in the 3D Holstein model also found that L = 4 and 6 clusters lie outside the scaling regime and need to be excluded from a finite size scaling analysis [50].

C. Local quantities
A signature of the phase transition is also present in the local densities shown in Fig. 3(a). At temperatures above the transition, T > T cdw , the average hole density on the Bi 6s sub-lattices are equal, n A 6s ≈ n B 6s . For T < T cdw these two local densities bifurcate with the onset of CDW order. Notably, the hole density on the oxygen atoms remains approximately constant across the transition. In other words, there is no net transfer of charge to the ligand oxygen atoms as charge modulations form on the Bi 6s orbitals. This behavior is reminiscent of the behavior seen in models for the rare-earth nickelates RNiO 3 [51,52], which are also negative charge transfer systems.
At low temperatures the difference in the hole density on each of the two Bi-6s sub-lattices begins to saturate, with n A 6s − n B 6s ≈ 0.47 at β = 17 eV −1 . It has been proposed that the formation of CDW order in BaBiO 3 is driven by charge disproportionation and the formation of a mixed valence state Bi 4+δ Bi 4−δ , corresponding to n A 6s − n B 6s = 2δ. In the extreme case, the charge disproportionation is complete, with δ = 1. Our model, however, predicts a much smaller charge transfer between neighboring Bi orbitals, which is more in line with experimental observations [9,26].
These results indicate that the holes in the compressed octahedra are not localized on a single Bi atom but instead occupy a more spatially distributed state involving the ligand oxygen orbitals. Previous investigations suggest this state is well described by hybridization between a Bi 6s orbital and an A 1g molecular orbital formed from the six surrounding O 2p σ orbitals [10,18]. The formation of this spatially extended hybridized state coincides with the oxygen octahedron surrounding a central bismuth atom contracting, thus forming a bipolaron. Conversely, the oxygen octahedron surrounding the six adjacent bismuth atoms must expand. There is a clear signature of this behavior in Fig. 3(b). The inset shows both the measured root mean square of the fluctuations in the phonon position, X rms = X 2 , as well as the analytic value X 0 rms in the α sp = 0 eV/Å limit, with the main panel displaying the difference X rms − X 0 rms . For T > T cdw , this difference initially decreases as the temperature is lowered, reflecting the reduction of thermal fluctuations of the lattice. However, at T ≈ T cdw it increases rapidly, then gradually levels off as the temperature is lowered further. This behavior can be attributed to the formation of an alternating pattern of expanded and contracted oxygen octahedron associated with the onset of CDW order. At temperatures well into the CDW phase, we obtain an X rms that is nearly temperature independent and consistent with the experimental values obtained from high-resolution neutron powder diffraction. For example, at β = 13 eV −1 (893 K), we obtained X rms = 0.100 ± 0.008Å, which agrees with the experimental value 0.094Å measured at this temperature [47]. These results explain why the measured Bi-O bond lengths in BaBiO 3 are essentially independent of temperature in the CDW insulating phase [47]. They also confirm that our model yields physically reasonable displacements and that no unphysical sign changes occur in the effective hopping integrals t sp ± α sp X .

D. Phonon renormalizations
To gain further insight into the nature of the CDW transition, we calculated the renormalized phonon energy Ω(q, 0) as a function of scattering momentum q. In the α sp = 0 eV/Å limit, the vibrations of the oxy- gen atoms in the Bi-Bi bond direction are described by a dispersionless optical phonon branch Ω(q, 0) = Ω 0 . Figure 4 shows that this behavior is approximately recovered at temperatures well above T cdw , with only a relatively small dip in Ω(q, 0) occurring at the ordering wave-vector q cdw . As the temperature is lowered, however, Ω(q, 0) is significantly renormalized, with a pronounced softening of the phonon branch occurring at q cdw . Moreover, this softening is quite sharp; at β = 15.0 eV −1 we see that Ω(q = q cdw , 0) ≈ Ω 0 , while Ω(q cdw , 0) has softened to a value close to zero. The inset of Fig. 4(a) shows that for T < T cdw , the frequency of the q cdw phonon mode levels off at a fixed value close to zero. The fact that it does not reach zero exactly is an expected finite size effect [see Fig. 4(b)], with a true divergence in χ(q cdw ) and complete softening only expected in the thermodynamic limit.
Our results for the average Bi 6s occupations and Bi-O bond distances demonstrate that the CDW transition in BaBiO 3 can be captured using the bonddisproportionation picture associated with the negative charge transfer nature of its electronic structure. However, Ref. [53] has pointed out that the charge and bond-disproportionation are just two ends of a continuous spectrum of possibilities. With this in mind, we next consider what effect modifying the on-site energy difference ∆ = ( p − s ) has on the character of the CDW order. This analysis will also reveal some subtle aspects of this transition.
E. The role of the charge-transfer energy Figure 5(a) shows how the average density of holes on the oxygen and bismuth atoms changes with ∆ . In particular, for ∆ ≈ −1.7 eV, the average total density of holes on the oxygen atoms equals the hole density on bismuth atoms, n 6s ≈ n 2px + n 2py + n 2pz . Increasing (decreasing) ∆ transfers holes to the bismuth (oxygen) orbitals. Figure 5(b) examines how the distribution of holes affects the charge correlations. Increasing ∆ not only increases n 6s but S 6s-6s (q cdw ) as well. Referring back to  Figs. 1(c-d), we see that the non-interacting charge susceptibility χ 0 (q) has a peak at q cdw that grows with β. However, the height of this peak also decreases monotonically with increasing ∆ . We can therefore conclude that the growth in S 6s-6s (q cdw ) with ∆ is not the result of changes in the bare electronic structure and the degree of Fermi surface nesting. Rather, the increase in S 6s-6s (q cdw ) can strictly be understood as resulting from the increasing hole density on the Bi-6s orbitals.
At face value, the observed increase in S 6s-6s (q cdw ) is not surprising as a large Bi hole density can support a stronger charge modulation on these orbitals. What is surprising is that S 6s-6s (q cdw ) does not appear to correlate to the degree of phonon softening strictly. For example, if you were to only look at S 6s-6s (q cdw ), you might conclude that the CDW correlations grow by nearly 33% as ∆ is changed from −3 eV to 1 eV. However, the CDW phase is characterized by both the hole density modulations and the periodic disproportionation of the Bi-O bond distances. As previously discussed, this periodic displacement of the oxygen atoms away from the midpoint between two neighboring bismuth atoms softens the q cdw mode. In Fig. 5(c) we see that Ω(q cdw , 0) is a nonmonotonic function of ∆ that does not track the dependence of S 6s-6s (q cdw ) on the charge transfer energy. Moreover, Ω(q cdw , 0) is minimized near ∆ ≈ −1.7 eV, the same point at which n 6s ≈ n 2px + n 2py + n 2pz . This observation highlights the crucial role of the hybridization between the Bi 6s and surrounding O 2p σ orbitals in establishing the CDW order.

F. Doping dependence of the charge correlations
Our results demonstrate that coupling to the bondstretching motion of the oxygen atoms in BaBiO 3 is sufficient to explain the CDW phase transition in this system at half-filling. To assess the relevance of this interaction for the doped system, we also calculated the charge susceptibility χ 6s−6s (q cdw ) as a function of doping, as shown in Fig. 6. Here, we focus on an L = 4 cluster and an inverse temperature of β = 18 eV −1 (≈ 644 K). While the strength of the CDW correlations decreases with doping, the correlations remain strong to carrier concentrations as large as n = 1.25 holes/Bi. The strength of these correlations is consistent with the Ba 1−x K x BiO 3 phase diagram, where the CDW phase extends out to n ≈ 1.2 holes/Bi at this temperature [3]. These results further support the notion that the SSH interaction with a relatively weak bare dimensionless coupling λ can account for many aspects of the bismuthate temperature-doping phase diagram.

III. DISCUSSION
Our model has a relatively weak bare dimensionless coupling λ = 0.18 when the e-ph coupling is characterized in the standard way [see Eq. (4) and Methods]. Moreover, the prominent peak in the bare charge susceptibility χ 0 (q) shown in Fig. 1(c) indicates that Fermi surface for the noninteracting band structure is well nested for the ordering wave-vector q cdw . Given this strong nesting condition, weak e-ph coupling, and the pronounced softening of the phonon dispersion predicted in Fig. 4, one might be tempted to conclude that the CDW transition is driven by a weak-coupling Peierls-like scenario. This picture, however, is inconsistent with many of our results. For example, a naïve RPA calculation over predicts T cdw by a factor of two (see Methods). This discrepancy suggests that electron and phonon self-energies must be included to describe the CDW transition in this system. Such effects would also result in an effective coupling λ eff that is much larger than the bare coupling computed from the noninteracting phononic and electronic structure. In fact, if we use the renormalized phonon energies to calculate the effective total and mode-resolved dimensionless e-ph couplings (see Methods), we obtain λ eff = 0.22 and λ eff q cdw ≈ 9.92, respectively, at β = 17 eV −1 . But what is perhaps more important is the fact that the RPA cannot account for the non-monotonic dependence of the soft phonon mode energy as a function of the 6s-2p (Bi-O) site energy difference ∆ . Our HMC results show that Ω(q cdw ) is minimized when the hybridization between the Bi 6s and O 2p orbitals is the largest, which is a consequence of the off-diagonal character of the SSH coupling and the multi-orbital nature of the problem. In contrast, a single-band RPA description predicts a monotonic behavior as a function of ∆ , contrary to the exact numerics.
Instead, our results point towards a bipolaronic scenario. In this picture, the ligand A 1g molecular orbital strongly hybridizes with the Bi 6s orbital at its center with an effective hopping t eff = √ 6t sp [30]. This hybridization creates bonding and antibonding states that shift ±t eff above and below the Fermi energy in hole language. A subsequent breathing distortion of the lattice creates a two sublattice structure and increases t eff . Therefore, the pair of holes originally residing on opposite sublattices in the undistorted structure can lower their kinetic energies significantly by condensing into the hybridized state below the Fermi energy. In this picture, the condensed holes and corresponding breathing distortions should be viewed as a bipolaron, and the CDW state as a frozen bipolaron lattice [18]. This scenario naturally explains why the largest softening of the Ω(q cdw ) mode occurs at the point of maximal hybridization between the Bi 6s and O 2p orbitals since this state would also produce the largest shift in the antibonding states and, subsequently, the largest binding energy between the holes. DFT calculations have found that such a polaronic mechanism results in a significant coupling (≈ 10 eV/Å) between the A 1g holes and the breathing lattice distortion [30], which far exceeds the bare e-ph coupling strength estimated using more conventional DFT methods. This is consistent with our estimated enhancement of the mode-resolved coupling with λ eff q cdw /λ q cdw ≈ 31 ≈ Ω 0 /Ω(q cdw , 0). We note that this bipolaron scenario was also advocated for in Ref. 18 based on determinant quantum Monte Carlo simulations of small 4 × 4 clusters in 2D, and for an unphysically large phonon energy. Our results demonstrate that this picture survives for the 3D lattice and a more realistic material-specific model. The bipolaronic description of the CDW phase holds for this system, despite the relatively small value of the bare dimensionless e-ph coupling λ = 0.18, which is below the (bare) value λ c = 0.5 where polaronic effects become significant in the 2D Holstein model [33]. These results thus underscore the notion that intuition derived from Holsteinand Fröhlich-like models, which are diagonal in orbital space, may not serve us well when studying off-diagonal e-ph interactions.
Our work also demonstrates that the early estimate for the bare coupling [32] is sufficient to account for the CDW phase in the bismuthates and, therefore, calls for a re-evaluation of e-ph interactions in the doped metallic phase. Previous attempts to explore superconductivity mediated by CDW exchange within the random phase approximation (RPA) [21,54] have examined an effective extended Hubbard Hamiltonian in which an intersite electron-electron repulsion replaces the electron-phonon interaction V ∼ α 2 /( Ω 0 ) 2 obtained by integrating out the phonons. While the resulting s-wave pairing interaction strength and spectral weight could account for superconductivity at high temperatures, the calculations were always considered problematic given the tendency of the RPA to overestimate T c and also the absence of retardation with instantaneous e-ph interactions. Our treatment resolves these issues by providing a nonperturbative solution to the original e-ph model at a physically motivated phonon frequency allowing for retardation. What remains now is to push these methods into the low-temperature metallic phase of the doped system.
In summary, we have focused on the parent compound of the bismuthate families of superconductors and obtained nonperturbative results for a model derived from ab initio electronic structure calculations. This is in contrast to more traditional QMC studies of these systems, which have traditionally focused on model Hamiltonians with large phonon energies Ω 0 /t ≈ 0.5−1 [18,23,33,[55][56][57][58][59]. Our work, therefore, represents a crucial bridge between ab initio nonpertubative numerical simulations of quantum materials with strong e-ph interactions. In this spirit, our results mark an important initial step to understanding a much broader range of perovskite materials involving the bond-stretching motion of their ligand anions. These include other negative charge-transfer materials like the rare-earth nickelates RNiO 3 [51,52] or the various families of sp-ABX 3 perovskites [5].

A. The non-interacting band structure
In the absence of e-ph coupling, the Hamiltonian can be diagonalized exactly in momentum space.

B. The dimensionless electron-phonon coupling
The dimensionless electron-phonon coupling constant for an arbitrary phonon branch Ω(q) that is relevant in theories of superconductivity is given by λ = 1 N q λ q , where λ q is a mode-resolved coupling defined as Here, N (0) = 1 N k δ( k ) is density of states at the Fermi level per spin species, and g(k, q) is the e-ph ver-tex for scattering within the band crossing the Fermi level g(k, q) = φ † 0,k+q G(k, q)φ 0,k , where Here g δ (k) = 2g cos k δ a 2 , and g = ( α 2 sp )/(2M Ω 0 ) = 99.4 meV. In the above, φ 0,k is the vector appearing in the Bloch wave function ψ 0,k (r) = φ 0,k e ik ·r for the band crossing the Fermi level, and is solved for by numerically diagonalizing H k .
To calculate λ in the noninteracting limit, we set Ω(q) = Ω 0 and approximate the delta functions appearing in Eq. (11) with Lorentz distributions δ(x) ≈ L γ (x) and half-width at half-maximum set to γ = 10 meV.We then calculate the average over the FS using a 50 3 grid of momentum points. (We have checked that our results are well converged for these values.) This procedure results in λ = 0.18. The corresponding value of the moderesolved coupling λ q cdw = 0.335. Interestingly, these values are comparable to early DFT estimates for coupling to the bond-stretching modes [32], as well as the value for the total coupling extracted from experiments on the bismuthates [60,61]. If we instead neglect the momentum dependence in the vertex and apply the approximations g(k, q) ≈ g and N (0) ≈ 1/W , where W = 5.62 eV is the bandwidth of the lowest energy band, the result is much smaller, with λ ≈ 0.0586.
To estimate the effective total λ eff and mode-resolved λ eff q dimensionless coupling constants, we approximated Ω q in Eq. (11) with the renormalized dispersion Ω(q, 0) [Eq. 7]. In this case, we interpolated the results for the L = 10 lattice onto the 50 3 grid of momentum points using a sinc interpolation scheme [62]. Using this approach, we obtain λ eff = 0.22 and λ eff q cdw ≈ 9.92, respectively, at β = 17 eV −1 . The enhancement in the mode-resolved coupling at q cdw can be attributed to the phonon softening, with λ eff q cdw /λ q cdw ≈ [Ω(q cdw , 0)/Ω 0 ] −1 .
C. An RPA estimate for T cdw The transition temperature to the CDW is signaled by a divergence in the charge susceptibility χ(q) in the thermodynamic limit. For weak coupling, χ(q) can be approximated using the RPA and calculated from the infinite sum polarization bubble diagrams to yield [63] χ RPA (q) = χ 0 (q) 1 + D 0 (q, 0)χ g,0 (q) .
Here, χ 0 (q) is the bare (Lindhard) susceptibility defined in the main text, D 0 (q, 0) = −2/ Ω 0 is the noninteracting phonon Green's function at zero frequency, and χ g,0 (q) is shorthand for with f (x) = [1 + exp(βx)] −1 the usual Fermi factor. The CDW transition thus occurs when the denominator of Eq. (13) goes to zero. Using a set of 50 3 momentum points, χ RPA (q cdw ) diverges at T cdw = 2630 K, over a factor of two larger than the value obtained from HMC.

D. HMC simulation details
As previously discussed, we used a recently introduced collection of algorithms to perform the numerical simulations described in this work. Below we highlight some of the relevant parameter values used in the simulation and refer the reader to Ref. [34] for more information on their definitions.
The HMC approach employed in this work is based on a path-integral formulation of the quantum partition function. For this, we used an imaginary-time discretization of ∆τ = 0.05 eV −1 . With this choice, global discretization errors of order O(∆τ 2 ) are well controlled. Updates to the phonon field were proposed using HMC trajectories of N t = 10 time-steps, and accepted or rejected according to a Metropolis criterion that ensures detailed balance. We solve for the force in our effective dynamics using the conjugate gradient method with a relative residual error threshold of max = 10 −5 . For the Metropolis accept/reject step, we reduce the threshold to max = 10 −10 .
Forces in the HMC dynamics can be decomposed into fermionic and bosonic parts. The former are more expensive to calculate, but also of smaller magnitude, so we use them to take relatively large time-steps of ∆t = 2. Within each of these large time-steps, we employed a time-step splitting algorithm that involves n t = 100 sub time-steps. The much smaller time-step value ∆t = 0.02 accounts for the rapidly evolving bosonic forces. To reduce autocorrelation times further, we apply Fourier acceleration via a dynamical mass matrix (regularization parameter m reg = 0.5) that modifies the effective relaxation time scale according to the imaginary-time wavelength. Systems were initially thermalized using N therm = 1000 HMC trial updates. After thermalization, an additional N sim = 4000 updates were performed. During this simulation stage, measurements were taken after every HMC trial update. Each measurement employed N rv = 10 random vectors to estimate the Green's function and higher-order correlation functions.

DATA AVAILABILITY
The data that support the findings of this study are accessible as a public repository at https://doi.org/ 10.5281/zenodo.7516925.

CODE AVAILABILITY
The primary Julia based simulation code used to generate the results presented in this paper is accessible as a public repository, https://github.com/cohensbw/ ElPhDynamics.git. Instruction on how to use the code will be provided upon reasonable request.