# A hybrid Monte Carlo study of bond-stretching electron-phonon interactions
and charge order in BaBiO3

Benjamin Cohen-Stead  Department of Physics and Astronomy, The University of
Tennessee, Knoxville, Tennessee 37996, USA  Institute for Advanced Materials
and Manufacturing, The University of Tennessee, Knoxville, Tennessee 37996,
USA  Kipton Barros  Theoretical Division and CNLS, Los Alamos National
Laboratory, Los Alamos, New Mexico 87545, USA  Richard Scalettar  Department
of Physics and Astronomy, University of California, Davis, California 95616,
USA  Steven Johnston<EMAIL_ADDRESS>Department of Physics and Astronomy,
The University of Tennessee, Knoxville, Tennessee 37996, USA  Institute for
Advanced Materials and Manufacturing, The University of Tennessee, Knoxville,
Tennessee 37996, USA

###### Abstract

The role of electron-phonon ($e$-ph) interactions and their relationship to
charge-density-wave (CDW) order in the bismuthate family of high-temperature
superconductors remains unresolved. We address this question using state-of-
the-art nonperturbative hybrid quantum Monte Carlo calculations for the parent
compound BaBiO3. Our model includes the Bi $6s$ and O $2p_{\sigma}$ orbitals
and coupling to the Bi-O bond-stretching branch of optical phonons via the
modulation of the Bi-O hopping integral. Using recently proposed algorithms
that are fast and scalable, we perform simulations on large three-dimensional
clusters of up to 4000 orbitals, with input model parameters taken directly
from ab initio electronic structure calculations and a physically motivated
phonon energy $\hbar\Omega_{0}=60$ meV. Our results demonstrate that the
coupling to the bond-stretching modes is sufficient to account for the CDW
transition in this system, despite the relatively small dimensionless coupling
to this branch predicted by density functional theory. We also find that the
transition deviates from the weak-coupling Peierls picture, despite the
noninteracting Fermi surface being well-nested at the CDW ordering vector. Our
results highlight the richness of $e$-ph interactions that are off-diagonal in
orbital space and suggest that such interactions play a crucial role in
establishing the bismuthate phase diagram. redThe 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 BaBiO3. 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 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 Ba1-xKxBiO3 and BaBi1-xPbxO3
exhibit high-temperature (high-$T_{\mathrm{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
Ba1-xKxBiO3 family with $T_{\mathrm{c}}\approx 32$ K for $x\approx 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 has also received new urgency redis
pertinent with the discovery of superconductivity in related (Ba,K)SbO3 [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 spin-fluctuation-based pairing mechanisms
proposed for high-$T_{\mathrm{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 an excellent reda compelling platform for studying
high-$T_{\mathrm{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, 13, 12], and theoretical [14, 15, 10, 16, 17, 18, 19]
techniques.

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. Bi4+ is a skipped oxidation state that does not
typically appear in nature. Early interpretations of the insulating phase,
therefore, proposed that the would-be Bi4+ ions undergo charge
disproportionation to (Bi4+δBi4-δ) [20, 21, 22]. In the low-temperature mixed
valant phase, Bi-O bond distances then relax in response to Bi charge,
resulting in a $\mathbf{q}=(\pi,\pi,\pi)/a$
red$\mathbf{q}=(\pi,\pi,\pi)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 [25,
26, 9]. This picture is also at odds with the significant hybridization
between the Bi $6s$ and O $2p$ orbitals and ab initio electronic structure
calculations [14, 10, 9], 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_{\mathrm{1g}}$ symmetry
[29, 10, 30]. Recent ARPES data for the parent compound support this view of
the electronic structure [9]. In this scenario, the $A_{\mathrm{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 [30, 18, 19].

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”
BCS-like mechanism [31, 12]. The standard BCS scenario was originally ruled
out because a naïve estimate for $T_{\mathrm{c}}$ based on the McMillian-
Allen-Dynes formula suggests that $\lambda\approx 1$ is needed to reproduce
the optimal $T_{\mathrm{c}}=32$ K in Ba1-xKxBiO3. Early theoretical estimates
based on density functional theory (DFT) [32] instead placed $\lambda\approx
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 $\lambda$ 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
$\mathbf{q}=(\pi,\pi,\pi)/a$ red$\mathbf{q}=(\pi,\pi,\pi)a^{-1}$ mode of the
optical bond-stretching phonon branch and estimated $\lambda\approx 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 $\lambda$ 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 BaBiO3, 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_{\sigma}$ 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 $\hbar\Omega_{0}=60$ meV [36,
32]. 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.

##  II Results

The model. We consider a four-orbital $sp$-model for BaBiO3 defined on an
$N=L^{3}$ cubic perovskite lattice (a total of $4L^{3}$ orbitals). The orbital
basis includes the Bi $6s$ and O $2p_{\sigma}$ orbitals polarized along the
Bi-Bi bond directions, as shown in Fig. 1(a). The Hamiltonian, written in hole
language, is
$\hat{H}=\hat{H}_{e}+\hat{H}_{\textrm{ph}}+\hat{H}_{e\textrm{-ph}}$, where

$\displaystyle\hat{H}_{e}=$
$\displaystyle(\epsilon_{s}-\mu)\sum_{\mathbf{i},\sigma}\hat{n}_{\mathbf{i},s,\sigma}+(\epsilon_{p}-\mu)\sum_{\mathbf{i},|\delta|,\sigma}\hat{n}_{\mathbf{i},|\delta|,\sigma}$
(1)
$\displaystyle+\sum_{\mathbf{i},\delta,\sigma}\mathcal{Q}^{\phantom{\dagger}}_{\delta}t^{\phantom{\dagger}}_{sp}[\hat{s}_{\mathbf{i},\sigma}^{\dagger}\hat{p}_{\mathbf{i},\delta,\sigma}^{\phantom{\dagger}}+\textrm{h.c.}]$
$\displaystyle+\sum_{\mathbf{i},\langle\delta,\delta^{\prime}\rangle,\sigma}\mathcal{Q}^{\phantom{\dagger}}_{\delta,\delta^{\prime}}t^{\phantom{\dagger}}_{pp}[\hat{p}_{\mathbf{i},\delta,\sigma}^{\dagger}\hat{p}_{\mathbf{i},\delta^{\prime},\sigma}^{\phantom{\dagger}}+\textrm{h.c.}]$

describes the noninteracting electron degrees of freedom,

$\hat{H}_{\textrm{ph}}=\sum_{\mathbf{i},|\delta|}\bigg{[}\frac{1}{2M}\hat{P}_{\mathbf{i},|\delta|}^{2}+\frac{1}{2}M\Omega_{0}^{2}\hat{X}_{\mathbf{i},|\delta|}^{2}\bigg{]}$
(2)

describes the noninteracting lattice degrees of freedom, and

$\hat{H}_{e\textrm{-ph}}=\sum_{\mathbf{i},\delta,\sigma}\alpha^{\phantom{\dagger}}_{sp}\hat{X}^{\phantom{\dagger}}_{\mathbf{i},\delta}[\hat{s}_{\mathbf{i},\sigma}^{\dagger}\hat{p}_{\mathbf{i},\delta,\sigma}^{\phantom{\dagger}}+\textrm{h.c.}]$
(3)

describes the $e$-ph interaction. Here, $\hat{s}^{\dagger}_{\bf{i},\sigma}$
and $\hat{p}^{\dagger}_{\bf{i},|\delta|,\sigma}$ create a spin-$\sigma$ hole
on a Bi $6s$ and O $2p_{\delta}$ orbital in the unit cell containing a Bi atom
located at $\mathbf{i}$ in the lattice, respectively. A sum over $\delta\
(|\delta|)$ runs over the set $\delta\in\\{\pm x,\pm y,\pm z\\}\
(|\delta|\in\\{x,y,z\\})$, where we have introduced the shorthand
$\hat{p}_{\mathbf{i},-|\delta|,\sigma}=\hat{p}_{\mathbf{i}-\tfrac{a}{2}|\bm{\delta}|,|\delta|,\sigma}$
and
$\hat{X}_{\mathbf{i},-|\delta|}=\hat{X}_{\mathbf{i}-\tfrac{a}{2}|\bm{\delta}|,|\delta|}$,
such that $a$ is the nearest neighbor Bi-Bi distance. A sum over
$\langle\delta,\delta^{\prime}\rangle$ denotes a sum over nearest-neighbor O
$2p$ orbital pairs that surround a central Bi $6s$ orbital. Finally,
$\hat{n}^{\phantom{\dagger}}_{\mathbf{i},s,\sigma}=\hat{s}_{\mathbf{i},\sigma}^{\dagger}\hat{s}^{\phantom{\dagger}}_{\mathbf{i},\sigma}$
and
$\hat{n}^{\phantom{\dagger}}_{\mathbf{i},|\delta|,\sigma}=\hat{p}_{\mathbf{i},|\delta|,\sigma}^{\dagger}\hat{p}^{\phantom{\dagger}}_{\mathbf{i},|\delta|,\sigma}$
are the spin-$\sigma$ hole number operators for a Bi $6s$ and O $2p$ orbital
in unit cell $\mathbf{i}$, respectively.

The on-site energies for the Bi $6s$ and O $2p$ orbitals are $\epsilon_{s}$
and $\epsilon_{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
$\mathcal{Q}_{\delta}$ and $\mathcal{Q}_{\delta,\delta^{\prime}}$, such that
$\mathcal{Q}_{\pm\delta}=\mp 1$,
$\mathcal{Q}_{\pm|\delta|,\pm|\delta^{\prime}|}=1$ and
$\mathcal{Q}_{\pm|\delta|,\mp|\delta^{\prime}|}=-1$, where
$\delta\neq\delta^{\prime}$. The sign convention enforced by these phase
factors is shown in Fig. 1(a). $\hat{X}_{\mathbf{i},|\delta|}$ and
$\hat{P}_{\mathbf{i},|\delta|}$ are the position and momentum operators for
each oxygen atom, where $M$ is the mass of an oxygen atom, $\Omega_{0}$ is the
bare phonon frequency, and $\alpha_{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\ \textrm{eV}$. The on-site
energies are set to $\epsilon_{s}=6.23$ and $\epsilon_{p}=4.14\ \textrm{eV}$
unless otherwise stated, consistent with the negative charge-transfer nature
of these materials [10, 9, 30]. We adjust the chemical potential to obtain an
average filling of $\langle n\rangle=1$ hole/Bi using a recently developed
$\mu$-tuning algorithm [38]. Fig. 1(b) displays the corresponding hole band
structure for this model with $\alpha_{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 off-diagonal in
orbital space, necessitates a multi-orbital description of the problem.

To estimate the $e$-ph coupling $\alpha_{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}\big{[}1+\tfrac{X}{a/2}\big{]}^{-2}$, where
$a=4.34\ \textrm{\AA}$ [39]. By expanding $t_{sp}(X)$ for small $X$ and
assuming a linear coupling, we obtain $\alpha_{sp}=\frac{4t_{sp}}{a}=2.13\
\mathrm{eV/\AA}$. This value corresponds to a bare dimensionless coupling
$\lambda=0.18$ in the standard Migdal-Eliashberg formulation, where

$\lambda=\frac{2}{\hbar\Omega_{0}N(0)N^{2}}\sum_{{\bf k},{\bf q}}|g({\bf
k},{\bf q})|^{2}\delta(\epsilon_{\bf k}-\mu)\delta(\epsilon_{{\bf k}+{\bf
q}}-\mu).$ (4)

Here, $N(0)$ is the density of states per spin at the Fermi energy and $g({\bf
k},{\bf q})$ is the $e$-ph coupling constant in momentum space. This value is
consistent with early DFT-based 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 $\hbar\Omega_{0}=60$ meV, which is typical for the bond-
stretching optical oxygen phonons in most oxides [43, 32]. 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_{\gamma\textrm{-}\nu}(\mathbf{q})=S_{\gamma\textrm{-}\nu}(\mathbf{q},\tau=0)$
given by

$\displaystyle
S_{\gamma\textrm{-}\nu}(\mathbf{q},\tau)=\frac{1}{N}\sum_{\mathbf{r},\mathbf{i}}e^{-\rm{i}\mathbf{q}\cdot\mathbf{r}}\langle\hat{n}_{\mathbf{i}+\mathbf{r},\gamma}(\tau)\hat{n}_{\mathbf{i},\nu}(0)\rangle,$
(5)

and the corresponding charge susceptibility

$\displaystyle\chi_{\gamma\textrm{-}\nu}(\mathbf{q})=\int_{0}^{\beta}S_{\gamma\textrm{-}\nu}(\mathbf{q},\tau)\
d\tau,$ (6)

where
$\hat{n}_{\mathbf{i},\nu}=(\hat{n}_{\mathbf{i},\nu,\uparrow}+\hat{n}_{\mathbf{i},\nu,\downarrow})$,
and $\gamma$ and $\nu$ specify the orbital species.

We also calculate the renormalized phonon energy

$\displaystyle\hbar\Omega(\mathbf{q},0)=\sqrt{\hbar^{2}\Omega_{0}^{2}+\Pi(\mathbf{q},0)},$
(7)

where $\Pi(\mathbf{q},\nu_{n})$ is a function related to the phonon self-
energy and $\nu_{n}=2\pi nT$ is a bosonic Matsubara frequency. We obtain
$\hbar\Omega(\mathbf{q},0)$ by measuring the phonon Green’s function

$\displaystyle
D(\mathbf{r},\tau)=-\frac{2M\Omega_{0}}{3N\hbar}\sum_{\mathbf{i},|\delta|}\langle\hat{X}_{\mathbf{r}+\mathbf{i},|\delta|}(\tau)\hat{X}_{\mathbf{i},|\delta|}(0)\rangle,$
(8)

where the Fourier transform in both $\mathbf{r}$ and $\tau$ is related to
$\Pi({\bf q},\nu_{n})$ by [33]

$D(\mathbf{q},\rm{i}\nu_{n})=\frac{2\hbar\Omega_{0}}{(\rm{i}\nu_{n})^{2}-\hbar^{2}\Omega_{0}^{2}-\Pi(\mathbf{q},\rm{i}\nu_{n})}.$
(9)

Charge-density-wave phase transition in BaBiO3. Both formal valence counting
and electronic structure calculations [32, 46, 10] indicate that BaBiO3, in
its high-temperature 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 $\langle n\rangle=1$ hole/Bi atom.

According to high-resolution neutron powder diffraction measurements, BaBiO3
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 $\mathbf{q}_{\textrm{cdw}}=(\pi,\pi,\pi)/a$
ordering vector and expansion and contraction of alternating BiO6 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
$\chi_{\textrm{6s-6s}}(\mathbf{q})$ versus $\mathbf{q}$, as shown in Fig.
2(a), we see a narrow peak form at the expected ordering wave-vector
$\mathbf{q}_{\textrm{cdw}}$. Fig. 2(b), therefore, focuses on the Bi charge
structure factor $S_{\textrm{6s-6s}}(\mathbf{q}_{\textrm{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_{\textrm{6s-6s}}(\mathbf{q}_{\textrm{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 $\mathbb{Z}_{2}$ symmetry between the
two otherwise equivalent Bi $6s$ sub-lattices, this phase transition falls in
the three-dimensional Ising universality class. Therefore, we use the
corresponding critical exponents $\gamma=1.237$ and $\nu=0.630$ to perform a
finite-size scaling (FSS) analysis to determine the transition temperature
[49]. Figure 2(c) displays the relevant
$S_{\textrm{6s-6s}}(\mathbf{q}_{\textrm{cdw}})L^{-\gamma/\nu}$ versus $\beta$
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_{\textrm{cdw}}\approx 10~{}\mathrm{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_{\textrm{6s-6s}}(\mathbf{q}_{\textrm{cdw}})L^{-\gamma/\nu}$ is plotted
against $L^{1/\nu}(T-T_{\textrm{cdw}})/T_{\textrm{cdw}}$, demonstrating the
quality of the scaling analysis. In this case, we have omitted the $L=4$ and
$6$ results from the collapse. Note that a similar 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].

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_{\textrm{cdw}}$, the average hole density on the Bi $6s$ sub-lattices are
equal, $\langle n_{6s}^{\textrm{A}}\rangle\approx\langle
n_{6s}^{\textrm{B}}\rangle$. For $T<T_{\textrm{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 $R$NiO3 [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 $\langle
n_{6s}^{\textrm{A}}\rangle-\langle n_{6s}^{\textrm{B}}\rangle\approx 0.47$ at
$\beta=17\ \textrm{eV}^{-1}$. It has been proposed that the formation of CDW
order in BaBiO3 is driven by charge disproportionation and the formation of a
mixed valence state $\textrm{Bi}^{4+\delta}\textrm{Bi}^{4-\delta}$,
corresponding to $\langle n_{6s}^{\textrm{A}}\rangle-\langle
n_{6s}^{\textrm{B}}\rangle=2\delta$. In the extreme case, the charge
disproportionation is complete, with $\delta=1$. Our model, however, predicts
a much smaller charge transfer between neighboring Bi orbitals, which is more
in line with experimental observations [26, 9].

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_{\mathrm{1g}}$ molecular orbital formed from the six surrounding O
$2p_{\sigma}$ 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_{\textrm{rms}}=\sqrt{\langle X^{2}\rangle}$, as well as the analytic value
$X^{0}_{\textrm{rms}}$ in the $\alpha_{sp}=0\ {\rm eV/\AA}$ limit, with the
main panel displaying the difference
$X^{\phantom{0}}_{\textrm{rms}}-X^{0}_{\textrm{rms}}$. For
$T>T_{\textrm{cdw}}$, this difference initially decreases as the temperature
is lowered, reflecting the reduction of thermal fluctuations of the lattice.
However, at $T\approx T_{\textrm{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_{\textrm{rms}}$ that is
nearly temperature independent and consistent with the experimental values
obtained from high-resolution neutron powder diffraction. For example, at
$\beta=13$ eV-1 (893 K), we obtained $X_{\mathrm{rms}}=0.100\pm 0.008$ Å,
which agrees with the experimental value $0.094~{}\mathrm{\AA}$ measured at
this temperature [47]. These results explain why the measured Bi-O bond
lengths in BaBiO3 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}\pm\alpha_{sp}\langle X\rangle$.

Phonon Renormalizations. To gain further insight into the nature of the CDW
transition, we calculated the renormalized phonon energy
$\Omega(\mathbf{q},0)$ as a function of scattering momentum $\mathbf{q}$. In
the $\alpha_{sp}=0\ \textrm{eV/\AA}$ limit, the vibrations of the oxygen atoms
in the Bi-Bi bond direction are described by a dispersionless optical phonon
branch $\Omega(\mathbf{q},0)=\Omega_{0}$. Figure 4 shows that this behavior is
approximately recovered at temperatures well above $T_{\textrm{cdw}}$, with
only a relatively small dip in $\Omega(\mathbf{q},0)$ occurring at the
ordering wave-vector $\mathbf{q}_{\textrm{cdw}}$. As the temperature is
lowered, however, $\Omega(\mathbf{q},0)$ is significantly renormalized, with a
pronounced softening of the phonon branch occurring at
$\mathbf{q}_{\textrm{cdw}}$. Moreover, this softening is quite sharp; at
$\beta=15.0\ \textrm{eV}^{-1}$ we see that
$\Omega(\mathbf{q}\neq\mathbf{q}_{\textrm{cdw}},0)\approx\Omega_{0}$, while
$\Omega(\mathbf{q}_{\textrm{cdw}},0)$ has softened to a value close to zero.
The inset of Fig. 4(a) shows that for $T<T_{\textrm{cdw}}$, the frequency of
the $\mathbf{q}_{\textrm{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 $\chi({\bf
q}_{\mathrm{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 BaBiO3 can be captured using the bond-
disproportionation 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 $\Delta\epsilon=(\epsilon_{p}-\epsilon_{s})$ has
on the character of the CDW order. This analysis will also reveal some subtle
aspects of this transition.

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
$\Delta\epsilon$. In particular, for $\Delta\epsilon\approx-1.7$ eV, the
average total density of holes on the oxygen atoms equals the hole density on
bismuth atoms, $\langle n_{6s}\rangle\approx\langle
n_{2p_{x}}+n_{2p_{y}}+n_{2p_{z}}\rangle$. Increasing (decreasing)
$\Delta\epsilon$ transfers holes to the bismuth (oxygen) orbitals.

Figure 5(b) examines how the distribution of holes affects the charge
correlations. Increasing $\Delta\epsilon$ not only increases $\langle
n_{6s}\rangle$ but $S_{\textrm{6s-6s}}(\mathbf{q}_{\textrm{cdw}})$ as well.
Referring back to Figs. 1(c-d), we see that the non-interacting charge
susceptibility $\chi_{0}(\mathbf{q})$ has a peak at
$\mathbf{q}_{\textrm{cdw}}$ that grows with $\beta$. However, the height of
this peak also decreases monotonically with increasing $\Delta\epsilon$. We
can therefore conclude that the growth in
$S_{\textrm{6s-6s}}(\mathbf{q}_{\textrm{cdw}})$ with $\Delta\epsilon$ is not
the result of changes in the bare electronic structure and the degree of Fermi
surface nesting. Rather, the increase in
$S_{\textrm{6s-6s}}(\mathbf{q}_{\textrm{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_{\textrm{6s-6s}}(\mathbf{q}_{\textrm{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_{\textrm{6s-6s}}(\mathbf{q}_{\textrm{cdw}})$
does not appear to correlate to the degree of phonon softening strictly. For
example, if you were to only look at
$S_{\textrm{6s-6s}}(\mathbf{q}_{\textrm{cdw}})$, you might conclude that the
CDW correlations grow by nearly 33% as $\Delta\epsilon$ is changed from
$-3~{}\mathrm{eV}$ to $1~{}\mathrm{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 $\mathbf{q}_{\textrm{cdw}}$ mode. In
Fig. 5(c) we see that $\Omega(\mathbf{q}_{\textrm{cdw}},0)$ is a nonmonotonic
function of $\Delta\epsilon$ that does not track the dependence of
$S_{\textrm{6s-6s}}(\mathbf{q}_{\textrm{cdw}})$ on the charge transfer energy.
Moreover, $\Omega(\mathbf{q}_{\textrm{cdw}},0)$ is minimized near
$\Delta\epsilon\approx-1.7$ eV, the same point at which $\langle
n_{6s}\rangle\approx\langle n_{2p_{x}}+n_{2p_{y}}+n_{2p_{z}}\rangle$. This
observation highlights the crucial role of the hybridization between the Bi
$6s$ and surrounding O $2p_{\sigma}$ orbitals in establishing the CDW order.

Doping dependence of the charge correlations. Our results demonstrate that
coupling to the bond-stretching motion of the oxygen atoms in BaBiO3 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 $\chi_{\mathrm{6s-6s}}({\bf
q}_{\mathrm{cdw}})$ as a function of doping, as shown in Fig. 6. Here, we
focus on an $L=4$ cluster and an inverse temperature of
$\beta=18~{}\mathrm{eV}^{-1}$ ($\approx 644~{}\mathrm{K}$). While the strength
of the CDW correlations decreases with doping, the correlations remain strong
to carrier concentrations as large as $\langle n\rangle=1.25$ holes/Bi. The
strength of these correlations is consistent with the Ba1-xKxBiO3 phase
diagram, where the CDW phase extends out to $\langle n\rangle\approx 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
$\lambda$ can account for many aspects of the bismuthate temperature-doping
phase diagram.

##  III Discussion

Our model has a relatively weak bare dimensionless coupling $\lambda=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
$\chi_{0}({\bf q})$ shown in Fig. 1(c) indicates that Fermi surface for the
noninteracting band structure is well nested for the ordering wave-vector
${\bf q}_{\mathrm{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_{\mathrm{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 $\lambda_{\mathrm{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
$\lambda^{\text{eff}}=0.22$ and $\lambda^{\text{eff}}_{{\bf
q}_{\text{cdw}}}\approx 9.92$, respectively, at $\beta=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 $\Delta\epsilon$. Our HMC results show
that $\Omega({\bf q}_{\mathrm{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 $\Delta\epsilon$, contrary to the exact
numerics.

Instead, our results point towards a bipolaronic scenario. In this picture,
the ligand $A_{\mathrm{1g}}$ molecular orbital strongly hybridizes with the Bi
$6s$ orbital at its center with an effective hopping
$t_{\mathrm{eff}}=\sqrt{6}t_{sp}$ [30]. This hybridization creates bonding and
antibonding states that shift $\pm t_{\mathrm{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_{\mathrm{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 $\Omega({\bf
q}_{\mathrm{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 ($\approx
10~{}\mathrm{eV}/\mathrm{\AA}$) between the $A_{\mathrm{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
$\lambda^{\text{eff}}_{{\bf q}_{\mathrm{cdw}}}/\lambda_{{\bf
q}_{\mathrm{cdw}}}\approx 31\approx\Omega_{0}/\Omega({{\bf
q}_{\mathrm{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\times 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 $\lambda=0.18$, which is below the (bare) value $\lambda_{c}=0.5$
where polaronic effects become significant in the 2D Holstein model [33].
These results thus underscore the notion that intuition derived from Holstein-
and 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\sim\alpha^{2}/(\hbar\Omega_{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}$ red$T_{\rm 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 $\hbar\Omega_{0}/t\approx 0.5-1$
[23, 55, 33, 18, 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 $R$NiO3 [51, 52] or the various families of $sp$-$ABX_{3}$
perovskites [5].

##  IV Methods

The non-interacting band structure. In the absence of $e$-ph coupling, the
Hamiltonian can be diagonalized exactly in momentum space. Introducing Fourier
transforms $s^{\dagger}_{{\bf i},\sigma}=\frac{1}{\sqrt{N}}\sum_{\bf
k}e^{-\mathrm{i}{\bf k}\cdot{\bf R}_{i}}s^{\dagger}_{{\bf k},\sigma}$ and
$p^{\dagger}_{{\bf i},\delta,\sigma}=\frac{1}{\sqrt{N}}\sum_{\bf
k}e^{-\mathrm{i}{\bf k}\cdot({\bf R}_{i}+\hat{\delta}a/2)}p^{\dagger}_{{\bf
k},\delta,\sigma}$ for the Bi $6s$ and O $2p_{\gamma}$ orbitals, we have
$H_{0}=\sum_{{\bf k},\sigma}\psi^{\dagger}_{{\bf
k},\sigma}H^{\phantom{\dagger}}_{\bf k}\psi^{\phantom{\dagger}}_{{\bf
k},\sigma}$, where $\psi^{\dagger}_{{\bf k},\sigma}=[s^{\dagger}_{{\bf
k},\sigma},p^{\dagger}_{{\bf k},x,\sigma},p^{\dagger}_{{\bf
k},y,\sigma},p^{\dagger}_{{\bf k},z,\sigma}]$ and

$H_{\bf k}=\left[\begin{array}[]{cccc}\epsilon_{s}-\mu&\mathrm{i}t_{x}({\bf
k})&\mathrm{i}t_{y}({\bf k})&\mathrm{i}t_{z}({\bf k})\\\ -\mathrm{i}t_{x}({\bf
k})&\epsilon_{p}-\mu&t^{\prime}_{x,y}({\bf k})&t^{\prime}_{x,z}({\bf k})\\\
-\mathrm{i}t_{y}({\bf k})&t^{\prime}_{x,y}({\bf
k})&\epsilon_{p}-\mu&t^{\prime}_{y,z}({\bf k})\\\ -\mathrm{i}t_{z}({\bf
k})&t^{\prime}_{x,z}({\bf k})&t^{\prime}_{y,z}({\bf
k})&\epsilon_{p}-\mu\end{array}\right],$ (10)

with $t_{\gamma}({\bf k})=-2t_{sp}\sin(k_{\gamma}a/2)$ and
$t^{\prime}_{\delta,\delta^{\prime}}({\bf
k})=4t_{pp}\sin(k_{\delta}a/2)\sin(k_{\delta^{\prime}}a/2)$.

The dimensionless electron-phonon coupling. The dimensionless electron-phonon
coupling constant for an arbitrary phonon branch $\Omega({\bf q})$ that is
relevant in theories of superconductivity is given by
$\lambda=\frac{1}{N}\sum_{\bf q}\lambda_{\bf q}$, where $\lambda_{\bf q}$ is a
mode-resolved coupling defined as

$\displaystyle\lambda_{\bf q}=\frac{2}{N(0)N}\sum_{\bf k}\frac{|g({\bf k},{\bf
q})|^{2}}{\hbar\Omega(\bf q)}\delta(\epsilon_{\bf k}-\mu)\delta(\epsilon_{{\bf
k}+{\bf q}}-\mu).$ (11)

Here, $N(0)=\frac{1}{N}\sum_{\mathbf{k}}\delta(\epsilon_{\mathbf{k}})$ is
density of states at the Fermi level per spin species, and $g({\bf k},{\bf
q})$ is the $e$-ph vertex for scattering within the band crossing the Fermi
level
$g(\mathbf{k},\mathbf{q})=\bm{\phi}_{0,\mathbf{k}+\mathbf{q}}^{\dagger}G({\bf
k},{\bf q})\bm{\phi}^{\phantom{\dagger}}_{0,\mathbf{k}}$, where

$G({\bf k},{\bf q})=\left[\begin{array}[]{cccc}0&g_{x}^{*}({\bf k}+{\bf
q})&g_{y}^{*}({\bf k}+{\bf q})&g_{z}^{*}({\bf k}+{\bf q})\\\
g_{x}(\mathbf{k})&0&0&0\\\ g_{y}(\mathbf{k})&0&0&0\\\
g_{z}(\mathbf{k})&0&0&0\end{array}\right].$ (12)

Here $g_{\delta}(\mathbf{k})=2g\cos\big{(}\tfrac{k_{\delta}a}{2}\big{)}$, and
$g=\sqrt{(\hbar\alpha_{sp}^{2})/(2M\Omega_{0})}=99.4$ meV. In the above,
$\bm{\phi}_{0,\mathbf{k}^{\prime}}$ is the vector appearing in the Bloch wave
function
$\bm{\psi}_{0,\mathbf{k}^{\prime}}(\mathbf{r})=\bm{\phi}_{0,\mathbf{k}^{\prime}}e^{{\rm
i}\mathbf{k}^{\prime}\cdot\mathbf{r}}$ for the band crossing the Fermi level,
and is solved for by numerically diagonalizing $H_{\mathbf{k}}$.

To calculate $\lambda$ in the noninteracting limit, we set $\Omega({\bf
q})=\Omega_{0}$ and approximate the delta functions appearing in Eq. (11) with
Lorentz distributions $\delta(x)\approx L_{\gamma}(x)$ and half-width at half-
maximum set to $\gamma=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 $\lambda=0.18$. The
corresponding value of the mode-resolved coupling $\lambda_{\bf
q_{\mathrm{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(\mathbf{k},\mathbf{q})\approx g$ and $N(0)\approx
1/W$, where $W=5.62$ eV is the bandwidth of the lowest energy band, the result
is much smaller, with $\lambda\approx 0.0586$.

To estimate the effective total $\lambda^{\text{eff}}$ and mode-resolved
$\lambda_{\bf q}^{\text{eff}}$ dimensionless coupling constants, we
approximated $\Omega_{\bf q}$ in Eq. (11) with the renormalized dispersion
$\Omega({\bf 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
$\lambda^{\text{eff}}=0.22$ and $\lambda^{\text{eff}}_{\mathbf{q}_{\rm
cdw}}\approx 9.92$, respectively, at $\beta=17$ eV-1. The enhancement in the
mode-resolved coupling at $\mathbf{q}_{\rm cdw}$ can be attributed to the
phonon softening, with $\lambda^{\text{eff}}_{\mathbf{q}_{\rm
cdw}}/\lambda_{\mathbf{q}_{\rm cdw}}\approx[\Omega(\mathbf{q}_{\rm
cdw},0)/\Omega_{0}]^{-1}.$

An RPA estimate for $T_{\mathrm{cdw}}$. The transition temperature to the CDW
is signaled by a divergence in the charge susceptibility $\chi({\bf q})$ in
the thermodynamic limit. For weak coupling, $\chi({\bf q})$ can be
approximated using the RPA and calculated from the infinite sum polarization
bubble diagrams to yield [63]

$\chi_{\mathrm{RPA}}({\bf q})=\frac{\chi_{0}({\bf q})}{1+D_{0}({\bf
q},0)\chi_{g,0}({\bf q})}.$ (13)

Here, $\chi_{0}({\bf q})$ is the bare (Lindhard) susceptibility defined in the
main text, $D_{0}({\bf q},0)=-2/\hbar\Omega_{0}$ is the noninteracting phonon
Green’s function at zero frequency, and $\chi_{g,0}({\bf q})$ is shorthand for

$\displaystyle\chi_{g,0}({\bf q})=-\frac{2}{N}\sum_{{\bf k}}|g({\bf k},{\bf
q})|^{2}\frac{f(\epsilon_{{\bf k}+{\bf q}})-f(\epsilon_{\bf
k})}{\epsilon_{{\bf k}+{\bf q}}-\epsilon_{\bf k}},$ (14)

with $f(x)=[1+\exp(\beta 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, $\chi_{\rm RPA}(\mathbf{q}_{\rm cdw})$ diverges at
$T_{\rm cdw}=2630$ K, over a factor of two larger than the value obtained from
HMC.

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 $\Delta\tau=0.05\ {\rm eV}^{-1}$. With this choice, global
discretization errors of order $\mathcal{O}(\Delta\tau^{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 $\epsilon_{\rm max}=10^{-5}$. For the Metropolis
accept/reject step, we reduce the threshold to $\epsilon_{\rm 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 $\Delta 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 $\Delta
t^{\prime}=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_{\rm reg}=0.5$) that modifies the
effective relaxation time scale according to the imaginary-time wavelength.
Systems were initially thermalized using $N_{\rm therm}=1000$ HMC trial
updates. After thermalization, an additional $N_{\rm sim}=4000$ updates were
performed. During this simulation stage, measurements were taken after every
HMC trial update. Each measurement employed $N_{\rm 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.

## Acknowledgements

The authors thank M. Berciu, S. Li, and G. A. Sawatzky for useful discussions.
This work was supported by the U.S. Department of Energy, Office of Science,
Office of Basic Energy Sciences, under Award Number DE-SC0022311.

## Author Contributions

S.J. conceived of the project and supervised the research. B.C-S., K.B. and
R.S. developed the code. B.C-S. collected the results and performed the
calculations. B.C-S., S.J., K.B. and R.S. discussed the results. B.C-S. and
S.J. wrote the manuscript.

## Competing Interests

The authors declare no competing financial or non-financial interests.

## References

  * [1] Mattheiss, L. F., Gyorgy, E. M. & Johnson, D. W. Superconductivity above 20 K in the Ba-K-Bi-O system. Phys. Rev. B 37, 3745–3746 (1988).
  * [2] Cava, R. J. et al. Superconductivity near 30 K without copper: the Ba0.6K0.4BiO3 perovskite. Nature 332, 814–816 (1988).
  * [3] Sleight, A. W. Bismuthates: BaBiO3 and related superconducting phases. Phys. C: Supercond. Appl. 514, 152–165 (2015).
  * [4] Minu, K. et al. Superconductivity in (Ba,K)SbO3. Nat. Mater. 21, 627–633 (2022).
  * [5] Benan, M. R., Foyevtsova, K., Khazraie, A., Elfimov, I. & Sawatzky, G. A. Bond versus charge disproportionation and nature of the holes in $s-p~{}AB{X}_{3}$ perovskites. Phys. Rev. B 104, 195141 (2021).
  * [6] Scalapino, D. J. A common thread: The pairing interaction for unconventional superconductors. Rev. Mod. Phys. 84, 1383–1417 (2012).
  * [7] Mattheiss, L. F. & Hamann, D. R. Electronic structure of $\mathrm{Ba}{\mathrm{Pb}}_{1-x}{\mathrm{Bi}}_{x}{\mathrm{O}}_{3}$. Phys. Rev. B 28, 4227–4241 (1983).
  * [8] Harrison, W. A. Valence-skipping compounds as positive-U electronic systems. Phys. Rev. B 74, 245128 (2006).
  * [9] Plumb, N.C. el al. Momentum-Resolved Electronic Structure of the High-${T}_{c}$ Superconductor Parent Compound ${\mathrm{BaBiO}}_{3}$. Phys. Rev. Lett. 117, 037002 (2016).
  * [10] Foyevtsova, K., Khazraie, A., Elfimov, I. & Sawatzky, G. Hybridization effects and bond disproportionation in the bismuth perovskites. Phys. Rev. B 91, 121114 (2015).
  * [11] Jiao, Y., Cheng, W., Deng, Q., Yang, H. & Wen, H. Collective vortex pinning and merging of the irreversibility line and second peak effect in optimally doped Ba1-xKxBiO3 single crystals. Phys. C: Supercond. Appl. 545, 43-49 (2018).
  * [12] Wen, C.H.P. et al. Unveiling the Superconducting Mechanism of ${\mathrm{Ba}}_{0.51}{\mathrm{K}}_{0.49}{\mathrm{BiO}}_{3}$. Phys. Rev. Lett. 121, 117002 (2018).
  * [13] Nicoletti, D. el al. Anomalous relaxation kinetics and charge-density-wave correlations in underdoped BaPb1-xBixO3. Proc. Natl. Acad. Sci. U.S.A. 114, 9020-9025 (2017).
  * [14] Franchini, C., Kresse, G. & Podloucky, R. Polaronic Hole Trapping in Doped $\mathrm{BaBiO}_{3}$. Phys. Rev. Lett. 102, 256402 (2009).
  * [15] Nourafkan, R., Marsiglio, F. & Kotliar, G. Model of the Electron-Phonon Interaction and Optical Conductivity of ${\mathrm{Ba}}_{1-x}{\mathrm{K}}_{x}{\mathrm{BiO}}_{3}$ Superconductors. Phys. Rev. Lett. 109, 017001 (2012).
  * [16] Naamneh, M. et al. Cooling a polaronic liquid: Phase mixture and pseudogap-like spectra in superconducting Ba1-xKxBiO3. Preprint at https://arxiv.org/abs/1808.06135 (2018).
  * [17] Li, Z., Antonius, G., Wu, M., da Jornada, F.H. & Louie, S.G. Electron-Phonon Coupling from Ab Initio Linear-Response Theory within the $GW$ Method: Correlation-Enhanced Interactions and Superconductivity in ${\mathrm{Ba}}_{1-x}{\mathrm{K}}_{x}{\mathrm{BiO}}_{3}$. Phys. Rev. Lett. 122, 186402 (2019).
  * [18] Li, S. & Johnston, S. Quantum Monte Carlo study of lattice polarons in the two-dimensional three-orbital Su–Schrieffer–Heeger model. npj Quantum Mater. 5, 40 (2020).
  * [19] Jiang, M., Sawatzky, G.A., Berciu, M. & Johnston, S. Polaron and bipolaron tendencies in a semiclassical model for hole-doped bismuthates. Phys. Rev. B 115129, 115129 (2021).
  * [20] Cox, D.E. & Sleight, A.W. Crystal structure of Ba2Bi3+Bi5+O6, Solid State Commun. 19, 969-973 (1976).
  * [21] Rice, T. M. & Sneddon, L. Real-Space and $\mathbf{k}$-Space Electron Pairing in $\mathrm{Ba}{\mathrm{Pb}}_{1-x}{\mathrm{Bi}}_{x}{\mathrm{O}}_{3}$. Phys. Rev. Lett. 47, 689–692 (1981).
  * [22] Varma, C. M. Missing valence states, diamagnetic insulators, and superconductors. Phys. Rev. Lett. 61, 2713–2716 (1988).
  * [23] Scalettar, R. T. et al. Phase diagram of the two-dimensional negative-$U$ Hubbard model. Phys. Rev. Lett. 62, 1407–1410 (1989).
  * [24] Fontenele, R. A., Costa, N. C., dos Santos, R. R. & Paiva, T. Two-dimensional attractive Hubbard model and the BCS-BEC crossover. Phys. Rev. B 105, 184502 (2022).
  * [25] Orchard, A. F. & Thornton, G. Evidence of ‘mixed-valency’ character in the X-ray photoelectron spectra of $\alpha$-diantimony tetraoxide and barium bismuth trioxide. J. Chem. Soc., Dalton Trans. 13, 1238-1240 (1977).
  * [26] Wertheim, G. K., Remeika, J. P. & Buchanan, D. N. E. Electronic structure of $\mathrm{Ba}{\mathrm{Pb}}_{1-x}{\mathrm{Bi}}_{x}{\mathrm{O}}_{3}$. Phys. Rev. B 26, 2120–2123 (1982)
  * [27] Zaanen, J., Sawatzky, G. A. & Allen, J. W. Band gaps and electronic structure of transition-metal compounds. Phys. Rev. Lett. 55, 418–421 (1985).
  * [28] Khomskii, D. I. Unusual valence, negative charge-transfer gaps and self-doping in transition-metal compounds. Lithuanian J. Phys. 37, 65 (1997).
  * [29] Ignatov, A.Yu. Relationship between the electronic and local structure in BaPbxBi1-xO3 and Ba1-xKxBiO3 perovskites. Nucl. Instrum. Methods. Phys. Res. B 448, 332-339 (2000)
  * [30] Khazraie, A., Foyevtsova, K., Elfimov, I., Sawatzky, G. A. Oxygen holes and hybridization in the bismuthates. Phys. Rev. B 97, 075103 (2018).
  * [31] Yin, Z. P., Kutepov, A., & Kotliar, G. Correlation-Enhanced Electron-Phonon Coupling: Applications of $GW$ and Screened Hybrid Functional to Bismuthates, Chloronitrides, and Other High-${T}_{\rm c}$ Superconductors. Phys. Rev. X 3, 021011 (2013).
  * [32] Meregalli, V. & Savrasov, S. Y. Electron-phonon coupling and properties of doped ${\mathrm{BaBiO}}_{3}$. Phys. Rev. B 57, 14453–14469 (1998).
  * [33] Esterlis, I. et al. Breakdown of the Migdal-Eliashberg theory: A determinant quantum Monte Carlo study. Phys. Rev. B 97, 140501 (2018).
  * [34] Cohen-Stead, B. et al. Fast and scalable quantum Monte Carlo simulations of electron-phonon models. Phys. Rev. E 105, 065302 (2022).
  * [35] ElPhDynamics, an open source code. Available online, https://github.com/cohensbw/ElPhDynamics.
  * [36] Braden, M., Reichardt, W., Ivanov, A. S. & Rumiantsev, A. Y. Anomalous dispersion of LO phonon branches in Ba0.6K0.4BiO3. Europhys. Lett. 34, 531–536 (1996).
  * [37] White, S. R. et al. Numerical study of the two-dimensional Hubbard model. Phys. Rev. B 40, 506–516 (1989).
  * [38] Miles, C. et al. Dynamical tuning of the chemical potential to achieve a target particle number in grand canonical Monte Carlo simulations. Phys. Rev. E 105, 045311 (2022).
  * [39] Harrison, W. A., Elementary electronic structure (World Scientific Publishing Company, Singapore, 2004).
  * [40] Duane, S., Kennedy, A. D., Pendleton, B. J. & Duncan, R. Hybrid Monte Carlo. Phys. Let. B 195, 216-222 (1987).
  * [41] Beyl, S., Goth, F. & Assaad, F. F. Revisiting the hybrid quantum Monte Carlo method for Hubbard and electron-phonon models. Phys. Rev. B 97, 085144 (2018).
  * [42] Batrouni, G. G. & Scalettar, R. T. Langevin simulations of a long-range electron-phonon model. Phys. Rev. B 99, 035114 (2019).
  * [43] Khosroabadi, H. et al. Softening of bond-stretching phonon mode in Ba1-xKxBiO3 at the metal-insulator transition. Phys. Rev. B 83, 224525 (2011).
  * [44] Prokof’ev, N. V. & Svistunov, B. V. Phonon modulated hopping polarons: $x$-representation technique. Phys. Rev. B 106, L041117 (2022).
  * [45] Li, S., & Johnston, S. Suppressed superexchange interactions in the cuprates by bond-stretching oxygen phonons. Preprint at https://arxiv.org/abs/2205.12678 (2022).
  * [46] Liechtenstein, A. I. et al. Structural phase diagram and electron-phonon interaction in ${\mathrm{Ba}}_{1\mathrm{-}\mathit{x}}$${\mathrm{K}}_{\mathit{x}}$${\mathrm{BiO}}_{3}$. Phys. Rev. B 44, 5388–5391 (1991).
  * [47] Kennedy, B. J., Howard, C. J., Knight, K. S., Zhang, Z. & Zhou, Q. Structures and phase transitions in the ordered double perovskites Ba2BiIIIBiVO6 and Ba2BiIIISbVO6. Acta Cryst. B 62, 537–546 (2006).
  * [48] Klinkova, L. A., Nikolaichik, V. I., Barkovskii, N. V. & Fedotov V. K. Thermal Stability of the Perovskite BaBiO3. J. Solid State Chem. 146, 439-447 (1999).
  * [49] Hasenbusch, M. Finite size scaling study of lattice models in the three-dimensional Ising universality class. Phys. Rev. B 82, 174433 2010.
  * [50] Cohen-Stead, B. et al. Langevin simulations of the half-filled cubic Holstein model. Phys. Rev. B 102, 161108 (2020).
  * [51] Park, H., Millis, A. J. & Marianetti, C. A. Site-Selective Mott Transition in Rare-Earth-Element Nickelates. Phys. Rev. Lett. 109, 156402 (2012).
  * [52] Johnston, S., Mukherjee, A., Elfimov, I., Berciu, M. & Sawatzky, G. A. Charge Disproportionation without Charge Transfer in the Rare-Earth-Element Nickelates as a Possible Mechanism for the Metal-Insulator Transition. Phys. Rev. Lett. 112, 106404 (2014).
  * [53] Khazraie, A., Foyevtsova, K., Elfimov, I., & Sawatzky, G. A. Bond versus charge disproportionation in the bismuth perovskites. Phys. Rev. B 98, 205104 (2018).
  * [54] Bickers, N. E., Scalapino, D. J. & Scalettar, R. T. CDW and SDW mediated pairing interactions. Int. J. Mod. Phys. B 1, 687–695 (1987).
  * [55] Weber, M., Assaad, F. F., Hohenadler, M. Phonon spectral function of the one-dimensional Holstein-Hubbard model. Phys. Rev. B 91, 235150 (2015).
  * [56] Nosarzewski, B. et al. Superconductivity, charge density waves, and bipolarons in the Holstein model. Phys. Rev. B 103, 235156 (2021).
  * [57] Cai, X., Li, Z. & Yao, H. Antiferromagnetism Induced by Bond Su-Schrieffer-Heeger Electron-Phonon Coupling: A Quantum Monte Carlo Study. Phys. Rev. Lett. 127, 247203 (2021).
  * [58] Costa, N. C., Seki, K. & Sorella, S. Magnetism and Charge Order in the Honeycomb Lattice. Phys. Rev. Lett. 126, 107205 (2021).
  * [59] Xing, B., Chiu, W., Poletti, D., Scalettar, R. T. & Batrouni, G. Quantum Monte Carlo Simulations of the 2D Su-Schrieffer-Heeger Model. Phys. Rev. Lett. 126, 017601 (2021).
  * [60] Graebner, J. E., Schneemeyer, L. F. & Thomas, J. K. Heat capacity of superconducting ${\mathrm{Ba}}_{0.6}$${\mathrm{K}}_{0.4}$${\mathrm{BiO}}_{3}$ near ${T}_{\rm c}$. Phys. Rev. B 39, 9682–9684 (1989).
  * [61] Marsiglio, F., Carbotte, J. P., Puchkov, A. & Timusk, T. Imaginary part of the optical conductivity of ${\mathrm{Ba}}_{1-x}{\mathrm{K}}_{x}\mathrm{Bi}{\mathrm{O}}_{3}$. Phys. Rev. B 53, 9433–9441 (1996).
  * [62] Wechsler, F. & Heintzmann, R. FourierTools.jl - Efficiently working with Fourier space. Available online at https://github.com/bionanoimaging/FourierTools.jl.
  * [63] Eiter, H.-M. et al. Alternative route to charge density wave formation in multiband systems. Proc. Natl. Acad. Sci. U.S.A. 110, 64-69 (2013).

## Figure Captions

Figure 1: Model definition and electronic properties in the noninteracting
limit. (a) Lattice structure of the three-dimensional, four-orbital model. A
Bi $6s$ orbital is surrounding by six O $2p_{\sigma}$ orbitals. Red arrows
indicate the relevant hopping integrals, with the phase factor included in
hole language. (b) The band structure in the noninteracting limit,
$\alpha_{sp}=0$. The horizontal red line indicates the Fermi energy at a
target filling of 1 hole per Bi atom. (c) For the target filling, the
noninteracting charge susceptibility $\chi_{0}(\mathbf{q})$ is computed using
the lowest energy band (black) in panel (b) for three different temperatures.
As the temperature is lowered, the peak at $\mathbf{q}=(\pi,\pi,\pi)/a$
red$\mathbf{q}=(\pi,\pi,\pi)a^{-1}$ grows. (d) At the fixed temperature
$\beta=10\ \textrm{eV}^{-1}$, $\chi_{0}(\mathbf{q})$ for four values of the
charge transfer energy $\epsilon_{2p}-\epsilon_{6s}$, the difference in the
on-site energy of the O $2p_{\sigma}$ and Bi $6s$ orbitals. The peak at
$\mathbf{q}=(\pi,\pi,\pi)/a$ red$\mathbf{q}=(\pi,\pi,\pi)a^{-1}$ is
approximately independent of the charge transfer energy. Figure 2: Finite size
scaling to determine transition temperature $T_{\textrm{cdw}}$. (a) The charge
susceptibility $\chi_{\textrm{6s-6s}}(\mathbf{q})$ versus $\mathbf{q}$ for
various $\beta$, plotted on a semilog scale for an $L=10$ cluster. As the
temperature is lowered, a peak forms at $\mathbf{q}=(\pi,\pi,\pi)/a$
red$\mathbf{q}=(\pi,\pi,\pi)a^{-1}$. (b) The charge structure factor
$S_{\mathrm{6s-6s}}(\mathbf{q}_{\textrm{cdw}})$, sensitive to CDW order on
Bi-6s orbitals, versus $\beta$. (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_{\textrm{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. Figure 3: Local quantities
versus $\beta$ for the $L=10$ system. (a) The total density of holes on the O
$2p_{\sigma}$ orbitals (black circles) remains approximately independent of
$\beta$ as the temperature is lowered below $T_{\mathrm{cdw}}$. The hole
densities on the Bi $6s$ orbital $A$ and $B$ sub-lattices, $\langle
n_{\textrm{6s}}^{\textrm{A}}\rangle$ and $\langle
n_{\textrm{6s}}^{\textrm{B}}\rangle$ respectively, bifurcate near
$T_{\mathrm{cdw}}$, indicating the onset of charge order. (b) The inset shows
the root mean square fluctuations in the phonon position as measured in the
presence (black circles) and the absence of interactions, $\alpha_{sp}=0$
(dashed line). The main figure displays the difference between these two
quantities, with an inflection point occurring near $T_{\textrm{cdw}}$. The
error bars in this figure are the standard deviation of the mean of our
measurements. Figure 4: The renormalized phonon dispersion. (a)
$\Omega(\mathbf{q},\mathrm{i}\nu_{n}=0)$ along the high-symmetry cuts of the
first Brillouin zone, computed on an $L=10$ lattice. As the temperature is
lowered $\Omega(\mathbf{q},\mathrm{i}\nu_{n}=0)$ softens at the ordering wave
vector $\mathbf{q}=(\pi,\pi,\pi)/a$ red$\mathbf{q}=(\pi,\pi,\pi)a^{-1}$, with
other $\mathbf{q}$ remaining essentially pinned at $\Omega_{0}$. The inset
shows the softening of the $\Omega(\mathbf{q}_{\textrm{cdw}},0)$ mode
specifically versus $\beta$. (b) $\Omega(\mathbf{q},\mathrm{i}\nu_{n}=0)$ for
different lattice sizes, measured at a fixed inverse temperature
$\beta=12.5~{}{\rm eV}^{-1}$ below the transition temperature. Here,
$\Omega(\mathbf{q},0)$ softens at $\mathbf{q}_{\rm cdw}$ for all lattice sizes
$L=6,8,\textrm{ and }10$. The inset shows $\Omega(\mathbf{q}_{\rm cdw},0)$
decreasing with increasing lattice size, with the dashed line acting as a
guide for the eye. The error bars in this figure are the standard deviation of
the mean of our measurements. Figure 5: Sweep through on-site energy
difference $\Delta\epsilon=(\epsilon_{p}-\epsilon_{s})$. Results were
collected on an $L=4$ lattice at inverse temperature $\beta=18\
\textrm{eV}^{-1}$. The vertical dashed lines indicates the ab initio value of
$\Delta\epsilon$. (a) The total hole density on the O $2p_{\sigma}$ orbitals
(cyan) and the Bi $6s$ orbitals versus $\Delta\epsilon$. Near
$\Delta\epsilon\approx-1.7$ eV, these two hole densities are commensurate. (b)
The structure factor $S_{\textrm{6s-6s}}(\mathbf{q}_{\textrm{cdw}})$
monotonically increases with $\Delta\epsilon$. (c) The renormalized phonon
energy of the $\Omega(\mathbf{q}_{\textrm{cdw}},0)$ mode has a local minimum
near $\Delta\epsilon\approx-1.7$ eV. The black line is a guide for the eye.
The error bars in this figure are the standard deviation of the mean of our
measurements. Figure 6: Change in susceptibilities upon doping. Results were
collected on an $L=4$ lattice at inverse temperature
$\beta=18~{}\textrm{eV}^{-1}$. The charge susceptibility
$\chi_{\textrm{6s,6s}}(\mathbf{q}_{\textrm{cdw}})$ monotonically decreases
upon doping. The error bars in this figure are the standard deviation of the
mean of our measurements.