# Intersite Coulomb Interactions in Charge Ordered Systems

Bo Gyu Jang  Present Address: Theoretical Division, Los Alamos National
Laboratory, Los Alamos, New Mexico 87545, USA  Korea Institute for Advanced
Study, Seoul 02455, Republic of Korea  Minjae Kim<EMAIL_ADDRESS>Korea Institute for Advanced Study, Seoul 02455, Republic of Korea  Sang-Hoon
Lee  Korea Institute for Advanced Study, Seoul 02455, Republic of Korea  Wooil
Yang  Department of Physics, Pohang University of Science and Technology,
Pohang 37673, Republic of Korea  Korea Institute for Advanced Study, Seoul
02455, Republic of Korea  Seung-Hoon Jhi  Department of Physics, Pohang
University of Science and Technology, Pohang 37673, Republic of Korea  Young-
Woo Son<EMAIL_ADDRESS>Korea Institute for Advanced Study, Seoul 02455,
Republic of Korea

###### Abstract

Using ab initio approaches for extended Hubbard interactions coupled to
phonons, we reveal that the intersite Coulomb interaction plays important
roles in determining various distinctive phases of the paradigmatic charge
ordered materials of Ba1-xK${}_{x}A$O3 ($A=$ Bi and Sb). We demonstrated that
all their salient doping dependent experiment features such as breathing
instabilities, anomalous phonon dispersions, and transition between charge-
density wave and superconducting states can be accounted very well if self-
consistently obtained nearest neighbor Hubbard interaction are included, thus
establishing a minimal criterion for reliable descriptions of spontaneous
charge orders in solids.

Since Verwey found the metal-to-insulator transition (MIT) in magnetite
(Fe3O4) owing to long-range order of alternating Fe2+ and Fe3+ ions [1], the
charge ordered state has been one of the central issues in condensed matter
physics. It often occurs near MITs, superconducting (SC), or charge density
wave (CDW) states [2, 3, 4, 5, 6, 7, 8, 9]. These local charge ordering (CO)
can lead to colossal magnetoresistance, ferroelectricity or multiferroicity
[10, 11]. In addition, the CO in high-temperature cuprate superconductors has
also been studied intensively to understand its roles as a leading competitor
of SC state [4, 5, 6, 7].

Charge ordered materials host atoms with disparate charging states that are
closely placed, invoking the strong Coulomb interactions and distorting
lattices to relieve their energetic cost [12]. Hence, theories for the systems
should treat the interaction as well as its coupling to lattices on an equal
footing. Typical approaches based on density functional theory with the local
density approximation (DFT-LDA) [13] or generalized gradient approximation
(GGA) [14] fail to describe their properties [15, 16]. The addition of on-site
Hubbard interaction ($U$) within DFT (DFT$+U$) [17, 18] captures a correct CO
state when the intersite interaction is screened [15]. Except few cases,
however, it is insufficient in obtaining the ground states [19, 20, 21].
Therefore, advanced methods beyond the local corrections are needed to figure
the role of the interactions in COs.

Perovskite BaBiO3 (BBO) and BaSbO3 (BSO) are prototypical CO materials that
are not well understood by DFT-GGA and DFT+$U$ [22, 23, 24, 25, 26, 27, 28,
29, 30, 31, 32, 33, 34]. They show CDW states characterized by breathing
distortion of oxygen octahedra with the charge disproportionations of Bi ions
[22, 23, 24]. The CDW is suppressed by substituting Ba with K and SC phase
occurs at a relatively higher transition temperature ($T_{c}$) [8, 9]. A few
recent studies based on the Heyd-Scuseria-Ernzerhof (HSE) hybrid functional
and $GW$ approximation (GWA) can capture the insulating ground state of BBO
and obtain enhanced electron-phonon ($e$-$ph$) interactions, providing a clue
to understanding the observed $T_{c}$ [32, 31, 33, 34, 35]. However, realistic
full phonon spectra to understand the transition between SC and CDW states for
experimentally accessible doping levels are hardly available due to demanding
computational resources and computed frequency is usually overestimated for
strongly coupled phonons [26, 33, 34].

On the other hand, recent developments on self-consistent evaluations of $U$
[36, 37] and intersite Hubbard interactions ($V$) [38, 39] within DFT
(DFT+$U$+$V$) [40] successfully describe various properties of solids [38, 39,
40, 41, 42, 43, 44, 45, 46, 47, 48, 49]. Owing to their low computational cost
comparable to DFT-LDA and improved accuracy to GWA [38, 46], the new method
provides an opportunity to study the correlated solids in large scale
structures and full phase space of interests. Motivated by these developments
as well as the works on BBO and BSO, we have carried out ab initio study to
explore the role of the interactions for interplay between their electronic
and structural properties.

In this Letter, we theoretically demonstrate that the evolution of charge
ordered states in bismuthates and antimonates with potassium doping is
essentially controlled by the doping-dependent non-local Coulomb interaction.
Our new parameter-free ab initio method for extended Hubbard interactions [38,
46] can compute electronic energy bands as well as full phonon dispersions of
Ba1-xK${}_{x}A$O3 ($A=$ Bi and Sb) for the whole phase space with structural
phase transitions, agreeing well with all the key measurements only when
intersite Coulomb interactions are included. Our establishments imply that the
explicit treatment of the non-local interactions is critical for the
description of the intertwined charge and lattice degrees of freedom in charge
ordered materials.

Our DFT+$U$+$V$ method uses the total energy functional of
$E_{\text{tot}}=E_{\text{DFT}}+E_{\text{Hub}}$ that can be decomposed into
(semi)local density functional of $E_{\textrm{DFT}}$ and Hubbard functional
with double counting corrections [40],
$E_{\textrm{Hub}}=\frac{1}{2}\sum_{I}U_{I}\sum_{m,m^{\prime},\sigma}(\delta_{mm^{\prime}}-n^{II\sigma}_{mm^{\prime}})n^{II\sigma}_{m^{\prime}m}$
$-\frac{1}{2}\sum_{\\{I,J\\}}V_{IJ}\sum_{m,m^{\prime},\sigma}n^{IJ\sigma}_{mm^{\prime}}n^{JI\sigma}_{m^{\prime}m},$
where the generalized occupation matrix is
$n_{mm^{\prime}}^{IJ\sigma}=\sum_{{\bf k}\nu}f_{{\bf k}\nu}^{\sigma}$
$\bra{\psi_{{\bf k}\nu}^{\sigma}}\ket{\phi_{m^{\prime}}^{J}}$
$\bra{\phi_{m}^{I}}\ket{\psi_{{\bf k}\nu}^{\sigma}}$, $f_{{\bf
k}\nu}^{\sigma}$ the Fermi-Dirac function of Kohn-Sham orbital of $\psi_{{\bf
k}\nu}^{\sigma}$ of the $\nu$th band with spin $\sigma$ at momentum ${\bf k}$,
and $\phi_{m}^{I}$ the localized orbitals with angular quantum number $m$.
Here, $I$ and $J$ are abbreviated indexes for atomic positions and principal
and azimuthal quantum numbers together and $\\{I,J\\}$ denotes a pair of atoms
within the nearest neighboring distances. We obtain self-consistent $U_{I}$
and $V_{IJ}$ using new pseudohybrid functionals for Hubbard interactions [36,
37, 38, 39]. We also consider rotationally invariant interactions [40, 38] so
that $U$ and $V$ for valence $s$ and $p$ orbitals of Bi (Sb) and 2$p$ orbital
of O are computed as shown in Table 1. Detailed parameters are in
Supplementary Materials (SM) [50].

Table 1: Calculated $U$ and $V$ (in eV) for Ba$A$O3 ($A=$ Bi and Sb). $U_{s(p)}^{A}$ ($U_{p}^{O})$ is on-site Hubbard parameters of valence $s(p)$ orbital of $A(O)$ atom. $V_{sp(pp)}$ is the intersite Hubbard parameters between $s(p)$ orbital of $A$ and $p$ orbital of O.  $A$ | $U_{s}^{A}$ | $U_{p}^{A}$ | $U_{p}^{O}$ | $V_{sp}$ | $V_{pp}$
---|---|---|---|---|---
Bi | 1.03 | 0.11 | 8.18 | 1.78 | 1.59
Sb | 0.92 | 0.13 | 8.18 | 1.86 | 1.61
Figure 1: (a) Atomic model for distorted octahedra in CDW state of BB(S)O.
Breathing distortion of $\delta_{B}$ and tilting angle of $\theta_{T}$ are
defined as an averaged length difference between Bi(Sb)-O bonds and as an
angle between vertical axes belong to adjacent octahedra, respectively. (b)
Double well potential as a function of $\delta_{B}$. Blue, red, and grey lines
indicate DFT-GGA, DFT+$U$+$V$ and experiments [51, 24], respectively. The DFT-
GGA results for $\delta_{B}>0$ are enlarged by ten times to show the minimum
for breathing distortion clearly.

At low-temperature, BBO is the monoclinic structure having oxygen octahedra
with breathing distortion of displacement by $\delta_{B}$ and a tilting angle
of $\theta_{T}$ between them as shown in Fig. 1(a), while BSO has $\delta_{B}$
only, resulting in fcc structure (Fm$\bar{3}$m) [24]. We first investigate the
artificial double-well potential induced by $\delta_{B}$ without tilting
[Fig.1(b)]. For BBO, our DFT-GGA calculation underestimates $\delta_{B}$ with
a very shallow potential well [26, 28]. The breathing distortion, however, is
significantly enhanced in DFT+$U$+$V$calculation. Unlike BBO, BSO has a
relatively deep potential with DFT-GGA. Nonetheless, the energy gain from the
breathing distortion in BSO also becomes larger with the extended Hubbard
interactions.

Table 2: Calculated structural and electronic properties of BBO along with computational and experimental data from previous studies. $v$: volume, $\beta$: monoclinic angle, $\delta_{B}$: breathing distortion, $\theta_{T}$: tilting distortion, and $E_{g}$: band gap  | Experiments 111References 22, 51, 52 | HSE222Reference 31 | DFT-GGA | DFT+$U$+$V$
---|---|---|---|---
$v$ (Å3) | 81.80 $\sim$ 82.54 | 82.10 | 85.03 | 82.94
$\beta$ (deg) | 90.16 $\sim$ 90.27 | 90.24 | 90.39 | 90.34
$\delta_{B}$ (Å) | 0.08 $\sim$ 0.09 | 0.09 | 0.08 | 0.10
$\theta_{T}$ (∘) | 10.12 $\sim$ 10.72 | 11.9 | 11.75 | 10.46
$E_{g}$ (eV) | 0.8 $\sim$ 1.1 | 0.84 | 0.0 | 0.99

As shown in the model potential wells in Fig. 1 (b), DFT+$U$+$V$ well captures
a part of CDW states of the undoped cases. Our results on a fully relaxed
monoclinic BBO are summarized and compared with the previous studies (Table
2). If tilting of $\theta_{T}$ is allowed, DFT-GGA accidentally reproduces
$\delta_{B}$ owing to its overestimation of volume by $\sim$ 4% and
$\theta_{T}$ by $\sim$ 14% but still cannot describe the insulating gap of CDW
state like a previous work [31]. On the other hand, DFT+$U$+$V$ calculation
well describes all the critical experimental parameters.

To investigate the effects of $U$ and $V$ in metallic BBO above the CDW
transition temperature of 800 K, we compute the energy bands of BBO in the
perfect cubic perovskite phase with experimental volume, using DFT-GGA,
DFT+$U$+$V$, and GWA as shown in Figs. 2 (a) and (b). The width of the energy
band crossing the Fermi level ($E_{F}$) is enhanced with both DFT+$U$+$V$and
GWA compared with one with DFT-GGA. We also note that the DFT-GGA bands for
fully occupied states associated with oxygen $p$ orbitals shifted down by
including $U^{O}_{p}$ in the DFT+$U$+$V$ and by the self-energy corrections in
the GWA method.

Figure 2: Energy bands of BBO of cubic perovskite structure obtained from (a)
DFT (blue) and DFT+$U$+$V$ (red) and (b) GWA method. Here we use $GW_{0}$
approximation for the GWA. (c) Density of states for monoclinic CDW structure.
Blue, red, and green line indicate DFT, DFT+$U$+$V$ and previous GWA results
[31], respectively. Empty circles are from photoemission spectroscopy
measurements [53, 54].

Effects of intersite interaction are also conspicuous for electronic
structures of the CDW state. Figure 2 (c) shows the calculated density of
states (DOS) of BBO with the fully relaxed monoclinic CDW structure obtained
from each method (Table 2) together with previous GWA results and
photoemission spectroscopy (PES) data [53, 54]. The two experiments show the
quite different positions of the highest PES peak below $E_{F}$ ($-3$ eV and
$-5$ eV, respectively), that may originate from the different substrate
conditions [54]. As shown in Fig. 2(c), the DOS peak position for oxygen $2p$
orbitals from DFT-GGA is significantly off the experimental one. Furthermore,
the energy gap of CDW state is absent despite the correct $\delta_{B}$. Unlike
DFT-GGA results, the DOS using DFT+$U$+$V$ well agrees with the experiments
and the previous GWA calculation [31]. The O 2$p$ peak position is located at
around $-4$ eV, which is in between two experimental results. The band gap is
about 1 eV, consistent with the experimental results. Although $U^{O}_{p}$ is
critical for the down shifted bands for O 2$p$ orbitals, we note that the
DFT+$U$ without $V$ still substantially underestimates the CDW band gap (Table
S1 in the SM).

The results so far demonstrate that our new method describes correctly the
electronic and structural properties of undoped bismuthate and antimonates
without serious computational costs. Thus, it enables the study of the non-
local interaction effects on phonons with various K doping levels thoroughly,
which could not be done with HSE or GWA method easily. From now on, we present
the comprehensive phonon dispersions using the frozen phonon techniques [55]
with extended Hubbard interactions to examine the structural instability of
potassium doped Ba1-xKxBiO3 and Ba1-xKxSbO3 (BKBO and BKSO). Related
computational details are in SM [50].

Figure 3: (a) Phonon dispersion and DOS of Ba0.6K0.4BiO3. Open diamonds and
blue (red) lines indicate experiments [56, 57] and the DFT-GGA (DFT+$U$+$V$)
results, respectively. Frequencies of the four selected modes are marked by
open and filled triangles and circles at $M$ and $R$ points, respectively.
Phonon frequency variations of the four modes as increasing K doping level in
(b) BKBO and (c) BKSO. Top and bottom panels are for the modes obtained by
DFT+$U$+$V$ and DFT-GGA methods, respectively. Experimental [51, 24] and our
theoretical phase boundaries between CDW and SC states are shown on the upper
abscissa and denoted by background color changes in (b) and (c), respectively.

Figure 3 (a) shows the phonon dispersion and DOS of Ba0.6K0.4BiO3 together
with experimental data [56, 57]. Here we focus four representative phonon
modes related with CDW states in BKBO and BKSO systems; stretching and in-
phase tilting modes at $M$ point and breathing and anti-phase titling modes at
$R$ point. For stretching mode, Bi-O bond length changes along only two axes
of oxgen octahedron while it changes along all three axes in the breathing
mode. In Fig. 3(a), it is immediately noticeable that the unstable
in(anti)-phase tilting mode at M(R) point obtained by DFT-GGA hardens enough
to be stable when $V$ is included, agreeing with experiments. Not only low
frequency modes, but the high frequency optical branches obtained with $V$
also agree with experiments. Specifically, our result well matches the
anomalous dispersion of LO mode along $\Gamma$X [56, 58, 59] related with the
instability toward the charge ordering [58, 59].

We present doping-dependent evolution of frequencies for the selected phonon
modes of BKBO and BSBO in Figs. 3(b) and (c). Our calculations with the
extended Hubbard interactions fruitfully reflect the measured trends of
structural distortions for the both systems. In the case of BKBO shown in Fig.
3(b), the breathing mode computed using DFT-GGA becomes to be stable when
$x\geq 0.15$ while it does only when $x\geq 0.4$ in the experiments [60, 61,
59] and our results with $V$. BKBO with $x\geq 0.4$ show SC states without
structural distortions [9, 62, 51, 56] while the tilting instabilities still
remain in DFT-GGA as shown in Figs. 3 (a) and (b) [34]. This is in sharp
contrast to complete absence of instability in DFT+$U$+$V$results when $x\geq
0.4$. Thus, the non-local Coulomb interaction is decisive in capturing the
transition between CDW and SC states of BKBO.

As already shown in Fig. 1 (b), the effect of $V$ is not as crucial in BKSO as
it is in BKBO. Nevertheless, our method makes improvements in describing
experimental phase diagram as shown in Fig. 3 (c). For undoped case, DFT-GGA
shows the tilting instability in addition to the breathing instability while
only the latter is observed in the experiment. Our calculation with $V$
results in a perfect fcc structure [24]. In addition, CDW phase survives up to
higher doping of $x\simeq 0.4$ in DFT+$U$+$V$ calculation, which is closer to
the experimental phase boundary of $x\simeq 0.5$, due to the enhanced
stretching instability. Comprehensive comparisons between results from DFT+$U$
and DFT+$U$+$V$ and detailed discussions are in Sec. D, Figs. S2, S5 and S7 of
the SM [50].

Figure 4:  Reduced $e$-$ph$ matrix element ($D^{L}_{R}$) as a function of
doping $x$. Diamonds (circles) denote BKB(S)O and blue, orange, and red color
indicate DFT-GGA, DFT+$U$, and DFT+$U$+$V$results, respectively. Open and
right-half diamonds are HSE results for the element from Ref. 35 and Ref. 33,
respectively.

Finally, we estimate the effects of non-local interactions on $e$-$ph$
coupling constant ($\lambda$) for the breathing mode that is an important
factor for the SC state. Instead of calculating $\lambda$ explicitly, we
compute ‘reduced $e$-$ph$ matrix element’ [33, 30, 63], $D^{\mu{\bf
k}}_{\alpha{\bf q}}=\partial\epsilon_{\mu{\bf k}}/{\partial{\bf
u}^{\alpha}_{\bf q}}$ to compare with previous studies based on HSE directly
[33, 35]. Here, ${\bf u}^{\alpha}_{\bf q}$ is a displacement vector of phonon
mode $\alpha$ with the wave-vector q and $\epsilon_{\mu\bf k}$ energy of
$\mu$th band at ${\bf k}$. Specifically, following previous works [33, 35], we
obtain the reduced element of $D^{L}_{R}$ for the energy band crossing the
$E_{F}$ at ${\bf k}$ of L point by the breathing mode at ${\bf q}$ of R point
(see Fig. S3 in SM [50]).

For both systems, $D^{L}_{R}$ increases with doping within DFT-GGA and DFT+$U$
as shown in Fig. 4. In sharp contrast to this, the matrix element with $V$
decreases with increasing $x$ for BKBO while it remains more or less the same
for BKSO. So, the matrix element of the latter becomes larger than that of the
former for $x>0.3$. Considering that the coupling can be roughly estimated as
$\lambda\sim(D^{\mu{\bf k}}_{\alpha{\bf q}}/\omega^{\alpha}_{\bf q})^{2}$ for
phonon frequency $\omega_{\bf q}^{\alpha}$ [33, 30, 63], Fig. 4 implies a
higher $T_{c}$ of BKSO than BKBO at larger doping levels with $V$, agreeing
with a recent experiment [24]. We also note that our estimations are
consistent with two previous HSE results for a few selected dopings [35, 33].

The opposite trends of $D^{L}_{R}$’s for different computational methods can
be understood by considering intersite interactions under lattice distortions.
The paired octahedra with the breathing modes in Fig. 1 have an elongated long
bond (LB) and a shrinking short bond (SB) between Bi(Sb) and oxygens [24].
Based on a perturbation theory, $D^{L}_{R}$ can be written as $D^{L}_{R}\delta
u_{B}\simeq\varepsilon_{\textrm{SB}}-\varepsilon_{\textrm{LB}}-V(n_{\textrm{SB}}-n_{\textrm{LB}})$
where $\varepsilon_{\textrm{SB(LB)}}$ and $n_{\textrm{SB(LB)}}$ are energy
level and density matrix for SB(LB), respectively, and $\delta u_{B}$ is the
perturbation amplitude of the breathing mode (See SM [50] for detailed
derivations). Since $n_{\textrm{LB}}>n_{\textrm{SB}}$ for lower dopings, the
intersite interaction should enhance the matrix elements over those from DFT-
GGA and DFT+$U$. As doping level is increased, the difference between
$n_{\textrm{LB}}$ and $n_{\textrm{SB}}$ diminishes so that the effect of
nonlocal interactions vanishes as shown in Fig. 4, thus highlighting again a
critical role of $V$ for doping dependent electron-phonon interactions.

In summary, we present a comprehensive study on the doping dependent
electronic and structural properties of prototypical charge ordered materials,
BKBO and BKSO using a newly developed ab initio computational method. We
demonstrated that the non-local Coulomb interactions between the nearest
neighbours are essential physical parameters in determining the doping
dependent evolution of CDW and SC states, highlighting nontrivial
relationships between the nonlocal interaction, charge order and lattice
distortion in correlated materials.

###### Acknowledgements.

B.G.J. M. K. and W. I. were supported by KIAS individual Grants (No. QP081301,
CG083501, QP090101). S.-H.J. was supported by NRF of Korea (Grant No.
2022R1A2C1006530) funded by the Korea government (MSIT). Y.-W.S. was supported
by NRF of Korea (Grant No. 2017R1A5A1014862, SRC program: vdWMRC center) and
KIAS individual Grant (No. CG031509). Computations were supported by the CAC
of KIAS.

## References

  * Verwey [1939] E. Verwey, Electronic conduction of magnetite (Fe3O4) and its transition point at low temperatures, Nature 144, 327 (1939).
  * Yamada  _et al._ [1996] Y. Yamada, O. Hino, S. Nohdo, R. Kanao, T. Inami, and S. Katano, Polaron Ordering in Low-Doping ${\mathrm{La}}_{1-\mathit{x}}{\mathrm{Sr}}_{\mathit{x}}{\mathrm{MnO}}_{3}$, Phys. Rev. Lett. 77, 904 (1996).
  * Van Aken  _et al._ [2003] B. B. Van Aken, O. D. Jurchescu, A. Meetsma, Y. Tomioka, Y. Tokura, and T. T. M. Palstra, Orbital-order-induced metal-insulator transition in La1-xCaxMnO3, Phys. Rev. Lett. 90, 066403 (2003).
  * Tranquada  _et al._ [1995] J. Tranquada, B. Sternlieb, J. Axe, Y. Nakamura, and S. Uchida, Evidence for stripe correlations of spins and holes in copper oxide superconductors, Nature 375, 561 (1995).
  * da Silva Neto  _et al._ [2014] E. H. da Silva Neto, P. Aynajian, A. Frano, R. Comin, E. Schierle, E. Weschke, A. Gyenis, J. Wen, J. Schneeloch, Z. Xu, S. Ono, G. Gu, M. L. Tacon, and A. Yazdani, Ubiquitous interplay between charge ordering and high-temperature superconductivity in cuprates, Science 343, 393 (2014).
  * da Silva Neto  _et al._ [2015] E. H. da Silva Neto, R. Comin, F. He, R. Sutarto, Y. Jiang, R. L. Greene, G. A. Sawatzky, and A. Damascelli, Charge ordering in the electron-doped superconductor Nd2-xCexCuO4, Science 347, 282 (2015).
  * Frano  _et al._ [2020] A. Frano, S. Blanco-Canosa, B. Keimer, and R. J. Birgeneau, Charge ordering in superconducting copper oxides, J. Phys.: Condens. Matt. 32, 374005 (2020).
  * Mattheiss  _et al._ [1988] L. Mattheiss, E. Gyorgy, and D. Johnson Jr, Superconductivity above 20 K in the Ba-K-Bi-O system, Phys. Rev. B 37, 3745 (1988).
  * Cava  _et al._ [1988] R. J. Cava, B. Batlogg, J. J. Krajewski, R. Farrow, L. W. Rupp, A. E. White, K. Short, W. F. Peck, and T. Kometani, Superconductivity near 30 K without copper: the Ba0.6K0.4BiO3 perovskite, Nature 332, 814 (1988).
  * Ramirez [1997] A. Ramirez, Colossal magnetoresistance, J. Phys.: Condens. Matt. 9, 8171 (1997).
  * Van Den Brink and Khomskii [2008] J. Van Den Brink and D. I. Khomskii, Multiferroicity due to charge ordering, J. Phys.: Condens. Matt. 20, 434217 (2008).
  * Attfield [2006] J. P. Attfield, Charge ordering in transition metal oxides, Sol. State Sci. 8, 861 (2006).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. 140, A1133 (1965).
  * Perdew  _et al._ [1992] J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh, and C. Fiolhais, Atoms, molecules, solids, and surfaces: Applications of the generalized gradient approximation for exchange and correlation, Phys. Rev. B 46, 6671 (1992).
  * Anisimov  _et al._ [1996] V. Anisimov, I. Elfimov, N. Hamada, and K. Terakura, Charge-ordered insulating state of Fe3O4 from first-principles electronic structure calculations, Phys. Rev. B 54, 4387 (1996).
  * Tsumuraya  _et al._ [2020] T. Tsumuraya, H. Seo, and T. Miyazaki, First-principles study of the charge ordered phase in $\kappa-{\mathrm{D}}_{3}{(\text{Cat-EDT-TTF/ST})}_{2}$: Stability of $\pi$-electron deuterium coupled ordering in hydrogen-bonded molecular conductors, Phys. Rev. B 101, 045114 (2020).
  * Hubbard [1963] J. Hubbard, Electron correlations in narrow energy bands, Proc. Roy. Soc. A 276, 238 (1963).
  * Schüler  _et al._ [2013] M. Schüler, M. Rösner, T. Wehling, A. Lichtenstein, and M. Katsnelson, Optimal Hubbard models for materials with nonlocal coulomb interactions: graphene, silicene, and benzene, Phys. Rev. Lett. 111, 036601 (2013).
  * Seo [2000] H. Seo, Charge ordering in organic ET compounds, J. Phys. Soc. Jpn. 69, 805 (2000).
  * Yamamoto  _et al._ [2007] S. Yamamoto, T. Fujiwara, and Y. Hatsugai, Electronic structure of charge and spin stripe order in ${\mathrm{La}}_{2-x}{\mathrm{Sr}}_{x}\mathrm{Ni}{\mathrm{O}}_{4}$ $(x=\frac{1}{3},\frac{1}{2})$, Phys. Rev. B 76, 165114 (2007).
  * Terletska  _et al._ [2017] H. Terletska, T. Chen, and E. Gull, Charge ordering and correlation effects in the extended Hubbard model, Phys. Rev. B 95, 115149 (2017).
  * Cox and Sleight [1976] D. Cox and A. Sleight, Crystal structure of Ba2Bi3+Bi5+O6, Solid State Commun. 19, 969 (1976).
  * Cox and Sleight [1979] D. E. Cox and A. W. Sleight, Mixed-Valent Ba2Bi3+Bi5+O6: Structure and Properties, Acta. Cryst. B35, 1 (1979).
  * Kim  _et al._ [2022] M. Kim, G. M. McNally, H.-H. Kim, M. Oudah, A. S. Gibbs, P. Manuel, R. J. Green, R. Sutarto, T. Takayama, A. Yaresko, U. Wedig, M. Isobe, R. K. Kremer, D. A. Bonn, B. Keimer, and H. Takagi, Superconductivity in (Ba,K)SbO3, Nature Materials 10.1038/s41563-022-01203-7 (2022).
  * Liechtenstein  _et al._ [1991] A. I. Liechtenstein, I. I. Mazin, C. O. Rodriguez, O. Jepsen, O. K. Andersen, and M. Methfessel, Structural phase diagram and electron-phonon interaction in Ba1-xKxBiO3 , Phys. Rev. B 44, 5388 (1991).
  * Meregalli and Savrasov [1998] V. Meregalli and S. Y. Savrasov, Electron-phonon coupling and properties of doped BaBiO3, Phys. Rev. B 57, 14453 (1998).
  * Meregalli and Savrasov [1999] V. Meregalli and S. Y. Savrasov, Electron-phonon coupling and properties of doped BaBiO3, J. Supercond. Nov. Magn. 12, 185 (1999).
  * Thonhauser and Rabe [2006] T. Thonhauser and K. M. Rabe, Fcc breathing instability in BaBiO3 from first principles, Phys. Rev. B 73, 212106 (2006).
  * Korotin  _et al._ [2012] D. Korotin, V. Kukolev, A. V. Kozhevnikov, D. Novoselov, and V. I. Anisimov, Electronic correlations and crystal structure distortions in BaBiO3, J. Phys.: Condens. Matter 24, 415603 (2012).
  * Korotin  _et al._ [2014] D. M. Korotin, D. Novoselov, and V. I. Anisimov, Correlation effects and phonon modes softening with doping in Ba1-xKxBiO3, J. Phys.: Condense. Matt. 26, 195602 (2014).
  * Franchini  _et al._ [2010] C. Franchini, A. Sanna, M. Marsman, and G. Kresse, Structural, vibrational, and quasiparticle properties of the Peierls semiconductor BaBiO3: A hybrid functional and self-consistent $GW+$vertex-corrections study, Phys. Rev. B 81, 1 (2010).
  * Franchini  _et al._ [2009] C. Franchini, G. Kresse, and R. Podloucky, Polaronic hole trapping in doped BaBiO3, Phys. Rev. Lett. 102, 256402 (2009).
  * Yin  _et al._ [2013] Z. P. Yin, A. Kutepov, and G. Kotliar, Correlation-enhanced electron-phonon coupling: Applications of $GW$ and screened hybrid functional to bismuthates, chloronitrides, and other high-Tc superconductors, Phys. Rev. X 3, 021011 (2013).
  * Li  _et al._ [2019] Z. Li, G. Antonius, M. Wu, F. H. Da Jornada, and S. G. Louie, Electron-Phonon Coupling from Ab Initio Linear-Response Theory within the GW Method: Correlation-Enhanced Interactions and Superconductivity in Ba1-xKxBiO3, Phys. Rev. Lett. 122, 186402 (2019).
  * Yuan  _et al._ [2022] Z. Yuan, P. Zheng, Y. Peng, R. Liu, X. Ma, G. Wang, T. Yu, and Z. Yin, Correlation-enhanced electron-phonon coupling and superconductivity in (Ba, K)SbO3 superconductors, Phys. Rev. B 105, 014517 (2022).
  * Agapito  _et al._ [2015] L. A. Agapito, S. Curtarolo, and M. Buongiorno Nardelli, Reformulation of $\mathrm{DFT}+U$ as a Pseudohybrid Hubbard Density Functional for Accelerated Materials Discovery, Phys. Rev. X 5, 011006 (2015).
  * Tancogne-Dejean  _et al._ [2018] N. Tancogne-Dejean, M. A. Sentef, and A. Rubio, Ultrafast Modification of Hubbard $U$ in a Strongly Correlated Material: Ab initio High-Harmonic Generation in NiO, Phys. Rev. Lett. 121, 097402 (2018).
  * Lee and Son [2020] S.-H. Lee and Y.-W. Son, First-principles approach with a pseudohybrid density functional for extended hubbard interactions, Phys. Rev. Research 2, 043410 (2020).
  * Tancogne-Dejean and Rubio [2020] N. Tancogne-Dejean and A. Rubio, Parameter-free hybridlike functional based on an extended Hubbard model: $\mathrm{DFT}+U+V$, Phys. Rev. B 102, 155117 (2020).
  * Campo Jr and Cococcioni [2010] V. L. Campo Jr and M. Cococcioni, Extended DFT$+U+V$ method with on-site and inter-site electronic interactions, J. Phys.: Condens. Matt. 22, 055602 (2010).
  * Huang  _et al._ [2020] J. Huang, S.-H. Lee, Y.-W. Son, A. Supka, and S. Liu, First-principles study of two-dimensional ferroelectrics using self-consistent Hubbard parameters, Phys. Rev. B 102, 165157 (2020).
  * Kulik and Marzari [2011] H. J. Kulik and N. Marzari, Transition-metal dioxides: A case for the intersite term in hubbard-model functionals, J. Chem. Phys. 134, 094103 (2011).
  * Cococcioni and Marzari [2019] M. Cococcioni and N. Marzari, Energetics and cathode voltages of $\mathrm{Li}M{\mathrm{PO}}_{4}$ olivines ($M=\mathrm{Fe}$, Mn) from extended Hubbard functionals, Phys. Rev. Materials 3, 033801 (2019).
  * Ricca  _et al._ [2020] C. Ricca, I. Timrov, M. Cococcioni, N. Marzari, and U. Aschauer, Self-consistent $\mathrm{DFT}+U+V$ study of oxygen vacancies in ${\mathrm{SrTiO}}_{3}$, Phys. Rev. Research 2, 023313 (2020).
  * Timrov  _et al._ [2020] I. Timrov, F. Aquilante, L. Binci, M. Cococcioni, and N. Marzari, Pulay forces in density-functional theory with extended Hubbard functionals: From nonorthogonalized to orthogonalized manifolds, Phys. Rev. B 102, 235159 (2020).
  * Yang  _et al._ [2021] W. Yang, S.-H. Jhi, S.-H. Lee, and Y.-W. Son, Ab initio study of lattice dynamics of group IV semiconductors using pseudohybrid functionals for extended Hubbard interactions, Phys. Rev. B 104, 104313 (2021).
  * Timrov  _et al._ [2021] I. Timrov, N. Marzari, and M. Cococcioni, Self-consistent hubbard parameters from density-functional perturbation theory in the ultrasoft and projector-augmented wave formulations, Phys. Rev. B 103, 045141 (2021).
  * Timrov  _et al._ [2022] I. Timrov, F. Aquilante, M. Cococcioni, and N. Marzari, Accurate electronic properties and intercalation voltages of olivine-type Li-ion cathode materials from extended hubbard functionals, PRX Energy 1, 033003 (2022).
  * Yang  _et al._ [2022] W. Yang, B. G. Jang, Y.-W. Son, and S.-H. Jhi, Lattice dynamical properties of antiferromagnetic oxides calculated using self-consistent extended hubbard functional method, J. of Phys.: Condens. Matter 34, 295601 (2022).
  * [50] Supplementary Materials for detailed computational methods, comprehensive phonon dispersions and self-consistent extended Hubbard parameters, which includes Ref. [64, 65, 66, 67, 68].
  * Pei  _et al._ [1990] S. Pei, J. D. Jorgensen, B. Dabrowski, D. G. Hinks, D. R. Richards, A. W. Mitchell, J. M. Newsam, S. K. Sinha, D. Vaknin, and A. J. Jacobson, Structural phase diagram of the Ba1-xKxBiO3, Phys. Rev. B 41, 4126 (1990).
  * Kennedy  _et al._ [2006] B. J. Kennedy, C. J. Howard, K. S. Knight, Z. Zhang, and Q. Zhou, Structures and phase transitions in the ordered double perovskites Ba2Bi${}^{\textrm{III}}$Bi${}^{\textrm{V}}$O6 and Ba2Bi${{}^{\textrm{III}}}$Sb${{}^{\textrm{V}}}$O6, Acta Cryst. B 62, 537 (2006).
  * Plumb  _et al._ [2016] N. C. Plumb, D. J. Gawryluk, Y. Wang, Z. Ristić, J. Park, B. Q. Lv, Z. Wang, C. E. Matt, N. Xu, T. Shang, K. Conder, J. Mesot, S. Johnston, M. Shi, and M. Radović, Momentum-Resolved Electronic Structure of the High-Tc Superconductor Parent Compound BaBiO3, Phys. Rev. Lett. 117, 1 (2016), 1603.01745 .
  * Oh  _et al._ [2018] J. S. Oh, M. Kim, G. Kim, H. G. Lee, H. K. Yoo, S. Sinn, Y. J. Chang, M. Han, C. Jozwiak, A. Bostwick, E. Rotenberg, H. D. Kim, and T. W. Noh, Evidence for absence of metallic surface states in BiO2-terminated BaBiO3 thin films, Curr. Appl. Phys. 18, 658 (2018).
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, First principles phonon calculations in materials science, Scr. Mater. 108, 1 (2015).
  * Braden  _et al._ [1995] M. Braden, W. Reichardt, W. Schmidbauer, A. S. Ivanov, and A. Y. Rumiantsev, Lattice dynamics of (Ba/K)BiO3, J. Supercond. 8, 595 (1995).
  * Loong  _et al._ [1992] C. K. Loong, P. Vashishta, R. K. Kalia, W. Jin, M. H. Degani, D. G. Hinks, D. L. Price, J. D. Jorgensen, B. Dabrowski, A. W. Mitchell, D. R. Richards, and Y. Zheng, Phonon density of states and oxygen-isotope effect in Ba1-xKxBiO3, Phys. Rev. B 45, 8052 (1992).
  * Braden  _et al._ [1996] M. Braden, W. Reichardt, A. S. Ivanov, and A. Y. Rumiantsev, Anomalous dispersion of LO phonon branches in Ba0.6K0.4BiO3, Europhys. Lett. 34, 531 (1996).
  * Braden  _et al._ [2002] M. Braden, W. Reichardt, S. Shiryaev, and S. N. Barilo, Giant phonon anomalies in the bond-stretching modes in doped BaBiO3: Comparison to cuprates manganites and nickelates, Phys. C Supercond. its Appl. 378-381, 89 (2002).
  * Tajima  _et al._ [1992] S. Tajima, M. Yoshida, N. Koshizuka, H. Sato, and S. Uchida, Raman-scattering study of the metal-insulator transition in Ba1-xKxBiO3, Phys. Rev. B 46, 1232 (1992).
  * Braden  _et al._ [2000] M. Braden, W. Reichardt, E. Elkaim, J. P. Lauriat, S. Shiryaev, and S. N. Barilo, Structural distortion in superconducting Ba1-xKxBiO3, Phys. Rev. B 62, 6708 (2000).
  * Fleming  _et al._ [1988] R. M. Fleming, P. Marsh, R. J. Cava, and J. J. Krajewski, Temperature dependence of the lattice parameters in the 30-K superconductor Ba0.6K0.4BiO3, Phys. Rev. B 38, 7026 (1988).
  * Kang and Kotliar [2019] C.-J. Kang and G. Kotliar, Material design of indium-based compounds: Possible candidates for charge, valence, and bond disproportionation and superconductivity, Phys. Rev. Materials 3, 015001 (2019).
  * Mosey and Carter [2007] N. J. Mosey and E. A. Carter, Ab initio evaluation of coulomb and exchange parameters for $\mathrm{DFT}+\mathrm{U}$ calculations, Phys. Rev. B 76, 155123 (2007).
  * Mosey  _et al._ [2008] N. J. Mosey, P. Liao, and E. A. Carter, Rotationally invariant ab initio evaluation of coulomb and exchange parameters for dft+u calculations, J. Chem. Phys. 129, 014103 (2008).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, _et al._ , QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials, J. Phys.: Condens. Matt. 21, 395502 (2009).
  * van Setten  _et al._ [2018] M. J. van Setten, M. Giantomassi, E. Bousquet, M. J. Verstraete, D. R. Hamann, X. Gonze, and G.-M. Rignanese, The PseudoDojo: Training and grading a 85 element optimized norm-conserving pseudopotential table, Comp. Phys. Comm. 226, 39 (2018).
  * Kim  _et al._ [2015] G. Kim, M. Neumann, M. Kim, M. D. Le, T. D. Kang, and T. W. Noh, Suppression of Three-Dimensional Charge Density Wave Ordering via Thickness Control, Phys. Rev. Lett. 115, 1 (2015).

## Supplementary Information

###  .1 A. Pseudohybrid density functional for extended Hubbard interactions

Figure S1: (a) Energy bands of BBO of cubic perovskite structure obtained from
DFT+$U$+$V$ with different $V$ cutoff length. For 3.2 Åcutoff case, only the
nearest Bi-O pairs are considered while up to 3rd nearest Bi-O pairs are
considered for 7.0 Åcutoff case. (b) Comparison between energy band obtained
from DFT+$U$ ($U^{Bi}_{s},U^{Bi}_{p}$, and $U^{O}_{p}$) and DFT+$U$+$V$. Here,
$U^{A}_{o}$ for $o$-orbital of atom $A$ shown in Table S2

The total energy formula with the extended Hubbard interactions including both
on-site and intersite Coulomb interaction ($U$ and $V$) is given by [40]

$\displaystyle E_{\textrm{DFT+U+V}}=E_{\textrm{DFT}}+E_{\textrm{UV}},$ (S1)

where $E_{\textrm{DFT}}$ is the (semi)local density functional and
$E_{\textrm{UV}}$ the extended Hubbard funcitonal. Within a rotationally
invariant form with the fully localized limit (FLL) double counting
correction, $E_{\textrm{UV}}$ is written as [40]

$\displaystyle E_{\textrm{UV}}$ $\displaystyle=$
$\displaystyle\sum_{I}\sum_{mm^{\prime}\sigma}\frac{U^{I}}{2}(\delta_{mm^{\prime}}-n^{II\sigma}_{mm^{\prime}})n^{II\sigma}_{m^{\prime}m}$
(S2)
$\displaystyle-\sum_{\\{I,J\\}}\sum_{mm^{\prime}\sigma}\frac{V^{IJ}}{2}n^{IJ\sigma}_{mm^{\prime}}n^{JI\sigma}_{m^{\prime}m},$

where the generalized occupation matrix is
$n_{mm^{\prime}}^{IJ\sigma}=\sum_{{\bf k}\nu}f_{{\bf
k}\nu}^{\sigma}\bra{\psi_{{\bf
k}\nu}^{\sigma}}\ket{\phi_{m^{\prime}}^{J}}\bra{\phi_{m}^{I}}\ket{\psi_{{\bf
k}\nu}^{\sigma}}$, $f_{{\bf k}\nu}^{\sigma}$ the Fermi-Dirac function of Kohn-
Sham orbital of $\psi_{{\bf k}\nu}^{\sigma}$ of the $\nu$th band with spin
$\sigma$ at momentum ${\bf k}$. The Löwdin orthonormalized atomic wave
function, $\phi_{m}^{I}$ is used as a projector for the localized atomic
orbital where $m$ is angular quatum number. Here, $I$ and $J$ are abbreviated
indexes for atomic positions and principal and azimuthal quantum numbers
together. $\\{I,J\\}$ in Eq. S2 denotes a pair of atoms within the nearest
neighboring distances.

Table S1: Calculated structural and electronic properties of BBO along with computational and experimental data from previous studies. $v$: volume, $\beta$: monoclinic angle, $\delta_{B}$: breathing distortion, $\theta_{T}$: tilting distortion, and $E_{g}$: band gap  | DFT-GGA | DFT+$U$ (w/o O $2p$) | DFT+$U$ | DFT+$U$+$V$
---|---|---|---|---
$v$ (Å3) | 85.03 | 85.03 | 82.19 | 82.94
$\beta$ (deg) | 90.39 | 90.40 | 90.27 | 90.34
$\delta_{B}$ (Å) | 0.08 | 0.08 | 0.10 | 0.10
$\theta_{T}$ (∘) | 11.75 | 11.76 | 9.71 | 10.46
$E_{g}$ (eV) | 0.00 | 0.00 | 0.71 | 0.99
$U^{Bi}_{p}$ | - | 0.14 | 0.09 | 0.11
$U^{O}_{p}$ | - | - | 8.38 | 8.23

###  .2 B. Calculation details

Following the ansatz by Mosey $et$ $al$. [64, 65], we extend the ACBN0
pseudohybrid functional for $U$ [36] into $V$ [38] by considering the
renormalized occupation number $N^{IJ\sigma}_{\psi_{\boldsymbol{k}v}}$ and
density matrix $P^{IJ\sigma}_{mm^{\prime}}$ for a pair of different atoms $I$
and $J$:

$\displaystyle N^{IJ\sigma}_{\psi_{{\bf k}v}}$ $\displaystyle=$
$\displaystyle\sum_{{I}}\sum_{m}\bra{\psi_{{\bf
k}v}^{\sigma}}\ket{\phi_{m}^{I}}\bra{\phi_{m}^{I}}\ket{\psi_{{\bf
k}v}^{\sigma}}$ (S3) $\displaystyle+\sum_{{J}}\sum_{m^{\prime}}\bra{\psi_{{\bf
k}v}^{\sigma}}\ket{\phi_{m^{\prime}}^{J}}\bra{\phi_{m^{\prime}}^{J}}\ket{\psi_{{\bf
k}v}^{\sigma}},$ $\displaystyle P^{IJ\sigma}_{mm^{\prime}}$ $\displaystyle=$
$\displaystyle\sum_{{\bf k}v}w_{{\bf k}}f_{{\bf k}v}N^{IJ\sigma}_{\psi_{{\bf
k}v}}\bra{\psi_{{\bf
k}v}^{\sigma}}\ket{\phi_{m}^{I}}\bra{\phi_{m^{\prime}}^{J}}\ket{\psi_{{\bf
k}v}^{\sigma}},$ (S4)

where $w_{\boldsymbol{k}}$ is $\bf{k}$-point weight and $f_{{\bf k}v}$ is the
Fermi-Dirac function of the Bloch state $\ket{\psi_{{\bf k}v}^{\sigma}}$. By
using Eq. S3 and S4, $U^{I}$, $J^{I}$, and $V^{IJ}$ can be defined as [36, 38,
39]:

$\scriptstyle
U^{I}=\frac{\sum_{m,\cdots,m^{\prime\prime\prime}}\sum_{\sigma\sigma^{\prime}}P^{II\sigma}_{mm^{\prime}}P^{II\sigma^{\prime}}_{m^{\prime\prime}m^{\prime\prime\prime}}(mm^{\prime}|m^{\prime\prime}m^{\prime\prime\prime})}{\sum_{m\neq
m^{\prime}}\sum_{\sigma}n^{II\sigma}_{mm}n^{II\sigma}_{m^{\prime}m^{\prime}}+\sum_{\\{m\\}}\sum_{\sigma}n^{II\sigma}_{mm}n^{II-\sigma}_{m^{\prime}m^{\prime}}},$
(S5) $\scriptstyle
J^{I}=\frac{\sum_{m,\cdots,m^{\prime\prime\prime}}\sum_{\sigma}P^{II\sigma}_{mm^{\prime}}P^{II\sigma}_{m^{\prime\prime}m^{\prime\prime\prime}}(mm^{\prime}|m^{\prime\prime}m^{\prime\prime\prime})}{\sum_{m\neq
m^{\prime}}\sum_{\sigma}n^{II\sigma}_{mm}n^{II\sigma}_{m^{\prime}m^{\prime}}},$
(S6) $\scriptstyle
V^{IJ}=\frac{1}{2}\frac{\sum_{m,\cdots,m^{\prime\prime\prime}}\sum_{\sigma\sigma^{\prime}}[P^{II\sigma}_{mm}P^{JJ\sigma^{\prime}}_{m^{\prime}m^{\prime}}-\delta_{\sigma\sigma^{\prime}}P^{IJ\sigma}_{mm^{\prime}}P^{JI\sigma^{\prime}}_{m^{\prime}m}](mm^{\prime}|m^{\prime\prime}m^{\prime\prime\prime})}{\sum_{\\{m\\}}\sum_{\sigma\sigma^{\prime}}[n^{II\sigma}_{mm}n^{JJ\sigma^{\prime}}_{m^{\prime}m^{\prime}}-\delta_{\sigma\sigma^{\prime}}n^{IJ\sigma}_{mm^{\prime}}n^{JI\sigma^{\prime}}_{m^{\prime}m}]},$
(S7)

where $(mm^{\prime}|m^{\prime\prime}m^{\prime\prime\prime})\equiv\int
d\textbf{r}_{1}d\textbf{r}_{2}\phi_{m}^{I*}(\textbf{r}_{1})\phi_{m^{\prime}}^{I}(\textbf{r}_{1})|\textbf{r}_{1}-\textbf{r}_{2}|^{-1}\phi_{m^{\prime\prime}}^{J*}(\textbf{r}_{2})\phi_{m^{\prime\prime\prime}}^{J}(\textbf{r}_{2})$.

Figure S2: Double well potential as a function of $\delta_{B}$. Blue, orange
closed diamond, orange open diamond, and red indicate DFT-GGA, DFT+$U$
($U^{Bi}_{s},U^{Bi}_{p}$), DFT+$U$ ($U^{Bi}_{s},U^{Bi}_{p}$, and $U^{O}_{p}$),
and DFT+$U$+$V$ respectively.

We implemented the pseudohybrid density functionals for extended Hubbard
interactions shown in Eqs. S5, S6, and S7 in QUANTUM ESPRESSO package [66, 38,
46]. We used the norm-conserving pseudopotentials provided by the PseudoDojo
project [67]. The energy cutoff for charge density was set to 480 Ry. A
15$\times$15$\times$15 $k$-point mesh was used for self-consistent calculation
of perfect perovskite structure (5 atoms in the unitcell). The $k$-mesh
density was kept the same for all other calculations. The Hubbard $U$ and $V$
parameters are obtained during the self-consistent calculation from Eqs. S5,
S6 and S7 as shown in Table S2.

Figure S3: Experimental pseudocubic lattice parameter $a_{P}$ as a function of
K doping of $x$ for BKBO and BKSO [51, 24]. Cubic supercell structures up to
40 atoms are used to model realistic K substitution structures. DFT-GGA
calculations overestimate the lattice constant while DFT+$U$+$V$ results well
agree with experimental observations.

The potassium (K) substitution effects are simulated by reducing the number of
electrons in the unitcell and inserting a compensating background charge. For
lattice optimization with doping, we used large cubic supercell structures
containing 40 atoms to model realistic doping. As shown in Fig. S3, our
optimized pseudocubic lattice parameter of $a_{P}$ using DFT+$U$+$V$ formalism
agrees with experimental data [51, 24] very well while DFT-GGA fails. With
these optimized structures, the phonon dispersions were calculated using the
finite displacement method by Phonopy [55]. We used fixed $U$ and $V$ values,
which are obtained from the equilibrium structure, to calculate the supercell
structures with atomic displacement. Full phonon band structures of BKBO and
BKSO are shown in Fig. S4.

For the direct comparison with the previous HSE result (Fig.4 of main text),
the reduced matrix element for estimating $e$-$ph$ interaction is obtained by
using the same scheme used in Ref. 33. For example, the breathing distortion
induces band splitting as shown Fig. S6. The energy splitting at L point
($\Delta E$) is used to determine the matrix element with the given oxygen
displacement.

Figure S4: Phonon band structures of BKBO and BKSO. Blue and red lines
indicate DFT and DFT+$U$+$V$ results, respectively.  Figure S5: Phonon
frequency variations of the four modes as increasing K doping level in (a)
BKBO and (b) BKSO. Top, middle, and bottom panels are for the modes obtained
by DFT-GGA, DFT+$U$ ($U^{Bi}_{s},U^{Bi}_{p}$), and DFT+$U$
($U^{Bi}_{s},U^{Bi}_{p}$, and $U^{O}_{p}$) methods, respectively. Experimental
[51, 24] and our theoretical phase boundaries between CDW and SC states are
shown on the upper abscissa and denoted by background color changes in (a) and
(b), respectively.

###  .3 C. Doping dependence of the reduced el-ph matrix elements

Table S2: Calculated $U$ and $V$ (in eV) for Ba$A_{1-x}$KxO3 ($A=$ Bi and Sb). $U_{s(p)}^{A}$ ($U_{p}^{O})$ is on-site Hubbard parameters of valence $s(p)$ orbital of $A(O)$ atom. $V_{sp(pp)}$ is the inter-site Hubbard parameters between $s(p)$ orbital of $A$ and $p$ orbital of O. $V_{sp}^{Ba}$ is the inter-site Hubbard parameters of between $s$ orbital of $Ba$ and $p$ orbital of O.  $x$ | $U_{s}^{Bi}$ | $U_{p}^{Bi}$ | $U_{p}^{O}$ | $V_{sp}$ | $V_{pp}$ | $V_{sp}^{Ba}$ |  | $U_{s}^{Sb}$ | $U_{p}^{Sb}$ | $U_{p}^{O}$ | $V_{sp}$ | $V_{pp}$ | $V_{sp}^{Ba}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---
$0.00$ | 1.03 | 0.11 | 8.18 | 1.78 | 1.59 | 0.77 |  | 0.92 | 0.13 | 8.18 | 1.86 | 1.61 | 0.75
$0.05$ | 1.00 | 0.11 | 8.17 | 1.76 | 1.58 | 0.76 |  | 0.95 | 0.13 | 8.15 | 1.82 | 1.59 | 0.74
$0.10$ | 0.98 | 0.11 | 8.15 | 1.74 | 1.58 | 0.75 |  | 0.91 | 0.13 | 8.12 | 1.78 | 1.58 | 0.74
$0.15$ | 0.96 | 0.11 | 8.14 | 1.72 | 1.57 | 0.75 |  | 0.88 | 0.13 | 8.08 | 1.74 | 1.56 | 0.73
$0.20$ | 0.94 | 0.11 | 8.13 | 1.70 | 1.56 | 0.74 |  | 0.85 | 0.13 | 8.05 | 1.71 | 1.55 | 0.72
$0.25$ | 0.92 | 0.11 | 8.11 | 1.68 | 1.56 | 0.74 |  | 0.83 | 0.12 | 8.01 | 1.68 | 1.54 | 0.71
$0.30$ | 0.90 | 0.11 | 8.09 | 1.66 | 1.55 | 0.73 |  | 0.80 | 0.12 | 7.97 | 1.65 | 1.52 | 0.71
$0.35$ | 0.89 | 0.11 | 8.08 | 1.64 | 1.54 | 0.72 |  | 0.78 | 0.12 | 7.93 | 1.63 | 1.51 | 0.69
$0.40$ | 0.87 | 0.11 | 8.06 | 1.63 | 1.54 | 0.72 |  | 0.76 | 0.12 | 7.89 | 1.60 | 1.49 | 0.68
$0.45$ | 0.85 | 0.10 | 8.04 | 1.61 | 1.52 | 0.72 |  | 0.74 | 0.12 | 7.84 | 1.58 | 1.48 | 0.68
$0.50$ | 0.83 | 0.10 | 8.02 | 1.59 | 1.52 | 0.70 |  | 0.72 | 0.12 | 7.79 | 1.55 | 1.46 | 0.67
$0.55$ | 0.81 | 0.10 | 8.00 | 1.57 | 1.51 | 0.70 |  | 0.70 | 0.12 | 7.74 | 1.53 | 1.44 | 0.66
$0.60$ | 0.80 | 0.10 | 7.98 | 1.56 | 1.50 | 0.69 |  | 0.68 | 0.12 | 7.69 | 1.51 | 1.43 | 0.65

The essential trend of the K doping ($x$) dependent variation of the reduced
el-ph matrix ($D^{L}_{R}$) of Ba1-xKx$A$O3 ($A$=Bi, Sb) shown in Fig. 4 of the
main text can be understood by considering simplified perturbative expansions
of interactions within DFT+$U$+$V$ scheme. The extended Hubbard interactions
are treated as perturbations to the Kohn-Sham (KS) equation of (semi)local
functional such that

${\mathcal{H}}_{\textrm{total}}={\mathcal{H}}_{\textrm{DFT}}+{\mathcal{H}}_{\textrm{UV}}$
(S8)

where ${\mathcal{H}}_{\textrm{total}}\psi^{\sigma}_{{\bf
k}\nu}=\epsilon^{\sigma}_{{\bf k}\nu}\psi^{\sigma}_{{\bf k}\nu}$ and
${\mathcal{H}}_{\textrm{DFT}}\varphi^{\sigma}_{{\bf
k}\nu}\equiv\left[-\nabla^{2}+V_{\textrm{DFT}}\right]\varphi^{\sigma}_{{\bf
k}\nu}=\varepsilon^{\sigma}_{{\bf k}\nu}\varphi^{\sigma}_{{\bf k}\nu}$. Here,
$V_{\textrm{DFT}}$ is the KS potential corresponding to (semi)local energy
functional ($E_{\textrm{DFT}}$) and $\varphi^{\sigma}_{{\bf k}\nu}$ the KS
wavefunction with eigenvalue of $\varepsilon^{\sigma}_{{\bf k}\nu}$. The
perturbation can be written as ${\mathcal{H}}_{\textrm{UV}}\psi^{\sigma}_{{\bf
k}\nu}={\delta E_{\textrm{UV}}}/\delta(\psi^{\sigma}_{{\bf k}\nu})^{*}$. As
discussed above, we only consider the alternation of energy bands crossing the
Fermi level at the $L$ point in the presence of breathing distortion at the
$R$ point as shown in Fig S6 so that the momentum, spin and band indexes will
be dropped hereafter.

The energy band considered here mainly originate from the anti-bonding state
between $s$ orbtial of Bi (Sb) at the center of octahedra and $p$ orbitals of
oxygen at their vertices [68]. Thus, the pair of distorted perovskites (see
Fig. 1a of the main text) has elongated and shortened bonds for the anti-
bonding states. We will call them the long and short bonding (LB and SB),
respectively. From these considerations, the split energies of the band by the
breathing mode in Fig. S6 can be approximated as energy eigenvalues of
$\epsilon_{\textrm{LB}}$ and $\epsilon_{\textrm{SB}}$ (or
$\varepsilon_{\textrm{LB}}$ and $\varepsilon_{\textrm{SB}}$ without
${\mathcal{H}}_{\textrm{UV}}$) for the static atomic configurations associated
with LB and SB, respectively. The resulting difference of
$\epsilon_{\textrm{SB}}-\epsilon_{\textrm{LB}}$ (or
$\varepsilon_{\textrm{SB}}-\varepsilon_{\textrm{LB}}$ without
${\mathcal{H}}_{\textrm{UV}}$) can be assigned as $D^{L}_{R}$ because the
energy level of the anti-bonding SB is higher than that of the LB.

Figure S6: Electronic band structure of Ba0.6K0.4BiO3 with and with out the
oxygen-breathing displacement obtained from DFT+$U$+$V$. The breathing
distortion induces band splitting dipicted by orange dotted lines. The band
splitting of $\Delta E$ indicated by the arrow are used to calculate the
reduced $e$-$ph$ matrix element of $D^{L}_{R}$.

With the amplitude $\delta u_{B}$ of the breathing mode, the
$E_{\textrm{DFT}}$ contribution is simply written as $D^{L}_{R}\cdot\delta
u_{B}=\varepsilon_{\textrm{SB}}-\varepsilon_{\textrm{LB}}$. In case of a
pristine octahedron, we can assign an overlap integral ($t_{sp}$) between $p$
orbital of oxygen and $s$ orbital of Bi(Sb) and the orbital energies for them
as $\omega_{p}$ and $\omega_{s}$ respectively. Then, with distorted structures
by the breathing mode, the energy splitting can be expressed up to the second
order of $t_{sp}$ as follows,

$D^{L}_{R}\delta
u_{B}=\varepsilon_{\textrm{SB}}-\varepsilon_{\textrm{LB}}\sim\frac{4t^{2}_{sp}\delta\tau_{sp}}{|\omega_{p}-\omega_{s}|},$
(S9)

where $t_{sp}$ changes to $t_{sp}(1\pm\delta\tau_{sp})$ for SB (LB) and
$|\omega_{p}-\omega_{s}|\gg t_{sp}$ is assumed. We note that the K doping
($x$) enhances the hole contribution to the $p$ orbital, thereby reducing the
$|\omega_{p}-\omega_{s}|$. So, Eq. S9 implies that $D^{L}_{R}$ based on DFT-
GGA increases with $x$, explaining the increasing matrix elements shown in
Fig. 4. We also note that a small contraction of lattice constant with doping
in Fig. S3 also could increases $t_{sp}$, thus corroborating the trend.

With onsite and intersite Hubbard interaction for $s$ orbital of Bi(Sb) and
$p$ orbital oxygen ($U_{p}$ and $V_{sp}$), the perturbative Hamiltonian for
the Hubbard interactions can be written as

${\mathcal{H}}_{\textrm{UV}}={\mathcal{H}}_{\textrm{U}}+{\mathcal{H}}_{\textrm{V}}\equiv
U_{p}(1-2n_{pp})-V_{sp}n_{sp},$ (S10)

where $n_{pp}$ is the density matrix for the $p$ orbital of oxygen and
$n_{sp}$ the generalized one for the neighboring $s$ orbital of Bi(Sb) and
oxygen $p$ orbitals. Here, the spin index is suppressed and the significantly
smaller onsite energy for $s$ orbital as shown in Table. S2 is neglected. For
$n_{pp}\neq 0$, the onsite Hubbard energy shifts the oxygen level of
$\omega_{p}$ toward $\omega_{s}$ so that, according to Eq. S9, the resulting
$D^{L}_{R}$ should enhance compared with those from DFT-GGA and increases as
the amount of doped potassium accumulates. This tendency is indeed confirmed
by our calculation only with $U$ as shown in Fig. 4.

From Eq. S10, we can estimate the contribution from the nonlocal interactions
for the matrix element as follows,

$D^{L}_{R}\delta
u_{B}=\epsilon_{\textrm{SB}}-\epsilon_{\textrm{LB}}\simeq\tilde{\varepsilon}_{\textrm{SB}}-\tilde{\varepsilon}_{\textrm{LB}}-V_{sp}(n_{sp}^{\textrm{SB}}-n_{sp}^{\textrm{LB}}),$
(S11)

where $\tilde{\varepsilon}_{\textrm{SB(LB)}}$ is the renormalized energy level
of Eq. S9 by the on-site energy of $U_{p}$ in Eq. S10 and
$n_{sp}^{\textrm{SB(LB)}}$ is the generalized density matrix for the SB (LB).
Here we note that we neglect the perturbation to the wave function in
considering Eq. S11 that seems to be smaller than the variations of density
matrix under perturbations. The density matrix for the lower energy state
should be larger than that for the higher, i.e.,
$n_{sp}^{\textrm{LB}}>n_{sp}^{\textrm{SB}}$ so that Eq. S11 immediately
implies the enhanced $D^{L}_{R}$ with $V$ as shown in Fig. 4 of the main text.
Furthermore, as doping increases, the charge disproportionation or difference
between the SB and LB should decreases and then disappear as $x$ approaches 1.
Therefore, $n_{sp}^{\textrm{SB}}-n_{sp}^{\textrm{LB}}$ in Eq. S11 decreases to
zero as doping increases so that the effect of intersite Hubbard interactions
to $D^{L}_{R}$ diminishes as $x$ increases as shown in Fig. 4.

###  .4 D. Comparison between DFT+$U$ and DFT+$U$+$V$

Figure S7: Reduced $e$-$ph$ matrix element ($D^{L}_{R}$) as a function of
doping $x$. Diamonds (circles) denote BKB(S)O and blue, orange open, orange
filled and red color indicate DFT-GGA, DFT+$U$ ($U^{Bi}_{s},U^{Bi}_{p}$),
DFT+$U$ ($U^{Bi}_{s},U^{Bi}_{p}$, and $U^{O}_{p}$) and DFT+$U$+$V$results,
respectively

Like in the previous study [29], DFT+$U$ on Bi $s$ (and $p$) atomic orbital
does not give any improvement compared to DFT-GGA method. The optimized
crystal structure (Table S1), double well potential (Fig. S2), phonon
frequency as a function of K doping level (Fig. S5), and reduced $e$-$ph$
matrix element (Fig. S7) obtained from DFT+$U$ ($U^{Bi}_{s},U^{Bi}_{p}$) [or
labeled as DFT+$U$ (w/o O $2p$)] are consistent with those from DFT-GGA.

When $U^{O}_{p}$ is included in DFT+$U$ calculation, it seemingly gives
improvement for some physical properties. Unlike a typical partially filled
case, the $U^{O}_{p}$ in BB(S)O changes the energetic position of O $p$
orbitals with respect to the Bi(Sb) $s$ orbitals as explained in the main
text. The O $p$ level is shifted down so that the energy difference between
bands for O $p$ orbitals and the Bi(Sb) $s$ orbital becomes smaller and the
covalency between them is enhanced (Fig. 2). As a result, the breathing
distortion and the band gap are enhanced compared to DFT-GGA as shown in Table
S1 and Fig. S2.

However, it is still not enough to describe charge ordered phase of BB(S)O.
The band gap of the monoclinic CDW phase obtained from DFT+$U$
($U^{Bi}_{s},U^{Bi}_{p}$, and $U^{O}_{p}$) is about 70% of that from
DFT+$U$+$V$and experimental observation. The computed phonon stabilities also
fail to explain experimental phase boundaries as shown in Fig. S5. The most
notable failure can be found in the reduced $e$-$ph$ matrix element
($D^{L}_{R}$) as shown in Fig. S7. As already discussed in the previous
section, the inclusion of $U^{O}_{p}$ affects the energy level and the
covalency so that $D^{L}_{R}$ is enhanced compared to DFT-GGA results.
However, $D^{L}_{R}$ keeps increasing as the K doping increases like DFT-GGA
result. Without intersite Coulomb interaction $V$, one cannot explain the
$e$-$ph$ coupling depending on K doping correctly as explained in the previous
section. The decreasing trend of $D^{L}_{R}$ depending K doping can be
captured only in DFT+$U$+$V$calculations and it well agrees with the previous
HSE results.