## Persistence of small polarons into the superconducting phase of Ba 1 -x K BiO x 3

Muntaser Naamneh , 1 Eugenio Paris, 2 Daniel McNally, 2 Yi Tseng, 2 Wojciech R. Pudeglyph[suppress]ko, l 2 Dariusz J. Gawryluk, 3 J. Shamblin, 4 Eric O'Quinn, 4 Benjamin Cohen-Stead , 5, 6 Ming Shi, 2, 7 Milan Radovic, 2 M. Lang, 4 Thorsten Schmitt, 2 Steven Johnston , 5, 6 and Nicholas C. Plumb 2

1 Department of Physics, Ben-Gurion University of the Negev, Beer-Sheva, 84105, Israel

2 Center for Photon Science, Paul Scherrer Institut, CH-5232 Villigen PSI, Switzerland

3 Laboratory for Multiscale Materials Experiments,

Paul Scherrer Institut, CH-5232 Villigen PSI, Switzerland

4 Department of Nuclear Engineering, The University of Tennessee, Knoxville, TN 37996, USA

5

Department of Physics and Astronomy, The University of Tennessee, Knoxville, TN 37996, USA

Institute of Advanced Materials and Manufacturing, The University of Tennessee, Knoxville, TN 37996, USA

7 Center for Correlated Matter and School of Physics, Zhejiang University, 310058, Hangzhou, China

(Dated: August 2, 2024)

Bipolaronic superconductivity is an exotic pairing mechanism proposed for materials like Ba 1 -x K BiO x 3 (BKBO); however, conclusive experimental evidence for a (bi)polaron metallic state in this material remains elusive. Here, we combine resonant inelastic x-ray and neutron total scattering techniques with advanced modelling to study the local lattice distortions, electronic structure, and electron-phonon ( -ph) e coupling in BKBO as a function of doping. Data for the parent compound ( x = 0 ) indicates that the electronic gap opens in predominantly oxygen-derived states strongly coupled to a long-range ordered breathing distortion of the oxygen sublattice. Upon doping, short-range breathing distortions and sizable e -ph coupling persist into the superconducting regime ( x = 0 4 . ). Comparisons with exact diagonalization and determinant quantum Monte Carlo calculations further support this conclusion. Our results provide compelling evidence that BKBO's metallic phase hosts a liquid of small (bi)polarons derived from local breathing distortions of the lattice, with implications for understanding the low-temperature superconducting instability.

Unconventional superconductivity typically appears near a parent phase with an ordered insulating or poorly metallic ground state. Notable examples include antiferromagnetic Mott insulating order in the cuprates [1], spin-density-wave order in the Fe-based superconductors [2], charge order in kagome superconductors [3], and charge-density-wave (CDW) order in the bismuthates [4] and antimonates [5]. Understanding the parent compounds and their residual correlations upon doping are key steps toward identifying the pairing mechanisms in these materials. Addressing this question is challenging, however, as the parent compounds' correlations by themselves are complex, and various degrees of freedom can become intertwined as the materials are doped, leading to the emergence of competing or coexisting states [6] and novel quasi-particles.

[10-12]. The compounds do not exhibit any magnetic order [13], and their band structures closely match density functional theory (DFT) calculations, implying that short-range Hubbard-like electronic correlations play a secondary role, if any [14, 15]. This view is further supported by recent hybrid Monte Carlo (HMC) calculations for the parent compound BaBiO 3 (BBO), which reproduce quantitative details of its CDW state using an uncorrelated DFT-derived tight-binding model and coupling to the Bi-O bond-stretching modes [16]. Thus, in the bismuthates, much of the phenomenology associated with unconventional superconductors appears to arise from a relatively limited set of interactions, and Ba 1 -x K BiO x 3 (BKBO) presents a unique opportunity to investigate aspects related to electron-phonon ( e -ph) coupling, polaron formation, and superconductivity.

Perovskite bismuth oxides [4, 7] are a promising platform to disentangle the roles of various interactions in oxide superconductors. The superconductivity in these compounds can be regarded as 'unconventional' in the sense that it emerges out of an insulating parent phase, with a doping-dependent transition temperature, T c , that can reach over 30 K. Moreover, their normal state resistivity follows a 'strange metal'-like linear temperature dependence that becomes more Fermi liquid-like at high doping [8] - a common feature of several correlated electron systems, including the cuprates [9]. At the same time, superconductivity in the bismuth oxides appears to be phonon-mediated with an s -wave order parameter

The phase diagram of BKBO is shown in Fig. 1 a . The parent compound is a semiconductor with an indirect band gap of 0.25 eV [17]. It has a perovskite structure with a long-range 'breathing' structural distortion in which the BiO 6 octahedra are alternately expanded and collapsed following a Q = ( π, π, π ) /a wave vector ( a is the Bi-Bi bond distance), as sketched in Fig. 1 b . At x = 0, the system is half-filled, with charges coupled to the breathing distortion to form a checkerboard of trapped holes or (bi)polarons [16, 18]. The polarons are diluted upon hole doping, and superconductivity emerges at a doping level of x ≈ 0 3. . The maximum or optimal T c occurs near or slightly below x = 0 4. .

6

FIG. 1. Overview of Ba 1 -x K BiO . x 3 a Sketch of the doping-temperature ( x -T ) phase diagram. The superconducting phase is labeled SC. b Two-dimensional sketch of the long-range ordered oxygen breathing distortions in the insulating parent compound ( x = 0). c Scenario of short-range, disordered (bi)polaron breathing distortions persisting in the metallic phase.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00401v1-figure-b0fcf65e-44b1-4e79-9af2-5628f93c259d.jpg)

Neither the details of the insulator-superconductor transition nor the polarons' evolution with doping, along with their impact on the superconducting state, are clear. The current understanding suggests two possible scenarios to explain how superconductivity emerges in the bismuthates. The first scenario minimizes the significant influence of the breathing distortion of the parent compound state on the superconducting state. Instead, it embraces a conventional mechanism where pairing is mediated by e -ph coupling that is enhanced by long-range Coulomb interactions [15, 19]. An alternative perspective posits that the small polarons derived from local breathing distortions of the BiO 6 octahedra persist in the metallic state, as sketched in Fig. 1 c , creating a (bi)polaronic liquid that ultimately forms a superconducting condensate at low temperatures [20-23]. While some studies have found evidence supporting the existence of such a polaron liquid [24-26], others have cast doubt on this scenario [27, 28].

Here, we investigate the presence and evolution of small polarons in BKBO across its insulatorsuperconductor transition using two experimental techniques: neutron scattering with pair distribution function (PDF) analysis to probe the crystal structure and O K -edge resonant inelastic x-ray scattering (RIXS) to probe the e -ph coupling across the transition. The experiments are complemented by determinant quantum Monte Carlo (DQMC) and exact diagonalization (ED) calculations. In contrast to many other spectroscopic or diffraction techniques, the PDF and RIXS methods used here are sensitive to the local environments of the BiO 6 octahedra and the oxygen sublattice. The neutron PDF measures the distribution of bond lengths in the system; it is thus sensitive to the presence of local lattice distortions, even if they are highly disordered and hence 'invisible' to standard diffraction techniques. Meanwhile, RIXS, by virtue of being performed on the oxygen K -edge resonance, can be tuned to select specifically those sites where the oxygen states are most strongly coupled to lattice breathing mode distortions and provide information about the coupling strength [29-32]. Together, the experiments and simulations show a consistent picture in which breathing distortions persist across the insulator-superconducting transition and remain strongly coupled to the electrons. Our findings contrast with a recent claim that these distortions vanish at low doping before the emergence of superconductivity [28].

## Results

## Local structure characterization

The neutron scattering measurements were performed on polycrystalline BKBO samples, as described in the Methods section. We begin by analysing BKBO's average crystal structure when viewed over long length scales, which is provided by a traditional diffraction analysis based on Rietveld refinement. In this case, BKBO is expected to have different (average) crystal structures and symmetries depending on doping level x and temperature T [4]. The parent compound adopts monoclinic structures in space group P 2 1 /n or space group I 2 /m at low ( T &lt; 140 K) and high temperatures (140 K &lt; T &lt; 430 K), respectively, while orthorhombic structures in space group Ibmm appear for intermediate doping levels ( x = 0 2 0 3), and tetragonal structures in either space . , . group I 4 mm or I 4 /mcm form in the superconducting samples ( x = 0 4). . To confirm this, we performed Rietveld refinements on all neutron diffraction patterns, as shown in Fig. 2 a . As expected, the diffraction pattern of the parent compound is modeled best using the monoclinic structure ( P 2 1 /n ). The data in the intermediate doping region ( x = 0 2 and 0 3) are . . fit well with the orthorhombic structure ( Ibmm ), while the data for the superconducting sample ( x = 0 4) is . fit well with the higher-symmetry tetragonal structure ( I 4 /mcm ).

To investigate the evolution of the local crystal structure with doping, we converted the neutron scattering data into PDFs, following the Fourier transform

FIG. 2. Rietveld and small-box refinements of neutron diffraction data . a Rietveld refinement of neutron diffraction patterns. Data was taken at T = 100 K for Ba 1 -x K BiO x 3 samples with different doping levels. Black dots correspond to the measurements. Simulated diffraction patterns for each structure are shown in red. b c , Small-box refinements of neutron total scattering data. Results are collected from Ba 1 -x K BiO x 3 at 100 K and fit against b the structural model inferred from Rietveld refinement and c the monoclinic P 2 1 /n model, which contains two distinct Bi sites. All refinements were performed over the range 1 5 . ≤ r ≤ 10 ˚. A

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00401v1-figure-400c6a79-31f2-4124-9f0d-c1aef9b0ee6c.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00401v1-figure-7278374e-bc0f-41ba-a0da-e9179ce30b5a.jpg)

in Eq. (2). The results at T = 100 K are shown in Figs. 2 b and 2 c . (We obtain similar results at 300 K, see Supplementary Figure 1. Supplementary Figure 2 shows the same data but zoomed in on the nearestneighbor Bi-O bond distance.) While the orthorhombic and tetragonal structures inferred from the diffraction analysis reasonably describe the local structure for the intermediate-doped ( x = 0 2 0 3) and superconduct-. , . ing samples ( x = 0 4), . respectively (see Fig. 2 b ), the PDF can be better modeled using the monoclinic phase ( P 2 1 /n ) across the whole doping range, as quantified by the fit goodness R w (see Fig. 2 c ).

To better understand the nature of the short-range distortions, Fig. 3 investigates the refined local structure obtained from the different structural models in more detail. A key feature of the monoclinic local structure is that it allows for two inequivalent octahedral sites occupied by Bi (Wyckoff positions 2 a and 2 ) and thus for b any potential breathing distortion. In contrast, the longrange tetragonal and orthorhombic structure only allows for one octahedral site. Using these models, we analyzed the local Bi-O lengths (Fig. 3 a ) and octahedral volumes (Fig. 3 b ), which are indicative of the breathing distortion. We also examined the octahedral bond angle variance (Fig. 3 c ) and quadratic elongation (Fig. 3 d ), which are indicative of trigonal and tetragonal distortions, respectively.

Small-box refinements of the PDFs based on the monoclinic structure identify two sets of Bi-O bond lengths and BiO 6 volumes for each composition evaluated. We find a smaller octahedron (red diamonds) exhibiting very little change in volume across the series, while a larger

FIG. 3. Evolution of the local BiO 6 octahedra with doping. Each panel shows the evolution of different local structural parameters as inferred from Rietveld refinement (black) of the diffraction patterns and small-box refinement (red) of PDFs collected at 100 K. Results are shown for a the local Bi-O bond lengths, b BiO 6 octahedral volumes, c Bi-O bond angle variance, and d quadratic elongation of the BiO 6 octahedra.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00401v1-figure-eef1ab07-36d8-4502-971a-3e967b47f8ea.jpg)

octahedron (red squares) decreases in volume monotonically with x . The two octahedra also display distinct characteristics in terms of their bond angle variances and quadratic elongations. For all measured dopings, both of

these quantities are consistently near zero for the smaller octahedron and clearly finite for the larger octahedron. The smaller BiO 6 unit therefore maintains near-ideal octahedral shape while the larger one possesses both trigonal and tetragonal distortions.

The local structure obtained from the PDF analysis contrasts with the structure averaged over long distances, as determined by Rietveld refinement. In the Rietveld analysis (black squares), the x = 0 0 sample exhibits . two distinct octahedra that not only have different volumes but also their own unique, finite angle variances and quadratic elongations. At x ≥ 0 2, these merge into . a single octahedron whose volume, angle variance, and quadratic elongation all decrease as x further increases. The trigonal and tetragonal distortions signified by the angle variance and quadratic elongation eventually vanish at x = 0 4. .

The two diffraction analyses combine to tell a nuanced but self-consistent story: averaged over long length scales, the two inequivalent, trigonally, and tetragonally distorted Bi octahedra ( x = 0) merge into one for x &gt; 0 2 . and trend toward ideal geometry at x = 0 4. . However, across short length scales, two inequivalent Bi octahedra persist in all compositions: a smaller, nearly idealized octahedron and a larger one that consistently exhibits trigonal and tetragonal distortions. Therefore, the distinction between the two local Bi octahedra is not simply in the volume variance between the two but also the difference of distortions.

## DQMC simulations of the O displacements

We performed determinant quantum Monte Carlo (DQMC) simulations of a multi-orbital sp -model to study the role of e -ph coupling in establishing the persistence of local breathing distortions with doping. The model is identical to the one used in Ref. [16], which successfully describes the CDW insulating phase of BBO and its collapse with doping. Here, we consider cubic Bi 64 O 256 clusters with a unit cell that includes the Bi 6 s and O 2 p σ orbitals oriented along the three Bi-Bi bond directions. To model the O phonons and their coupling to the holes, we place a quantum harmonic oscillator on each O atom, which is polarized along the Bi-Bi bond direction and modulates the Bi-O hopping integral. We fix the corresponding phonon energy to ℏ Ω = 60 meV and set all remaining parameters to values obtained from DFT calculations [21]. For further details, see Ref. [16] and the Methods section.

Figures 4 a e -plot the distribution of O displacements generated during our simulations. Results are shown for T = 300 K and a range of hole concentrations from the parent compound ⟨ n ⟩ = 1 ( x = 0) to slightly over-doped ⟨ n ⟩ = 1 4 ( . x = 0 4) in units of holes/Bi. . Figs. 4 f j present the same data but projected on the eigenvector of a local breathing distortion with A 1g symmetry

$$X _ { A _ { i g } } = \frac { 1 } { \sqrt { 6 } } \left ( X _ { i, x } + X _ { i, y } + X _ { i, z } - X _ { i, - x } - X _ { i, - y } - X _ { i, - z } \right ). \quad \text{The c} \\ ( 1 ) \quad \text{rically $\varepsilon$}$$

FIG. 4. DQMC simulations of the Bi-O bond-lengths as a function of doping. a e -Probability distributions of the lattice displacements measured in our DQMC simulations of a cubic Bi 64 O 256 cluster at T = 300 K. Data for different excess hole concentrations x = ⟨ n ⟩ -1 (in holes/Bi) is shown, as indicated in each panel. f j -The same data but projected on the A 1g breathing mode. The gray shaded area in each panel is the measured distribution. The dashed black and solid red lines are one- and two-Gaussian fits to the data, respectively. The black dashed curves show the individual contributions of the peaks within the two-Gaussian fits. The bottom portions of panels d e , , i , and j plot the residuals between the fits and the data to better highlight the differences between the oneand two-Gaussian models (red and black, respectively).

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00401v1-figure-04eace2a-4dea-42c0-b00e-58932b9029ad.jpg)

Here, X i , ± δ denotes the displacement of the O atoms surrounding the Bi atom in unit cell i . In each panel, the solid red and black lines are one- and two-Gaussian fits to the data, respectively, with the dashed lines indicating the individual components of the latter. Panels d , e , i , and j also plot the difference between the histograms and the models in the lower half each panel.

The distributions for ⟨ n ⟩ = 1 and 1 1 each consist of . a pair of equally weighted Gaussians centered symmetrically about zero displacement, indicating that the O

FIG. 5. RIXS data from BBO thin film. a RIXS intensity map of BBO across the O K -edge as highlighted by blue in panel e . b A close-up view of the dashed rectangular region in the RIXS map of a , highlighting the low-energy excitations. The dashed line marks the resonance energy. c RIXS spectra near the elastic peak, plotted for various incident energies. The series displays an intensity enhancement of the low-energy excitations at the resonance energy of 527.9 eV (red curve), which corresponds to the peak of the O K -edge absorption curve illustrated in panel e . d Fluorescence measured by RIXS at the resonance energy (solid black line, left and bottom axes) plotted along with an angle-integrated photoemission spectrum (PES, red shaded area, right and top axes). The enhanced portions of the RIXS spectrum relative to the PES spectrum (yellow shading) indicate regions dominated by oxygen states, including the gap edge. e O K -edge XAS spectrum measured in partial fluorescence yield mode at T = 20 K. The blue shading highlights the photon energy range from which the RIXS map is acquired in panel a .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00401v1-figure-ae7ba36b-6ce7-4dae-ac47-a8c52cf5b969.jpg)

atoms fluctuate around two inequivalent equilibrium positions. The distribution of the A 1g projected displacements is also bifurcated, as expected for the long-range breathing distortion predicted by this model at this temperature [16]. Increasing the hole concentration suppresses the long-range order, causing the two distributions to begin merging. Nevertheless, for 1 2 . ≤ ⟨ n ⟩ ≤ 1 3, . the distributions are better described as a pair of Gaussian distributions, indicating some fraction of the BiO 6 octahedra remains compressed in the simulations. The persistence of this bifurcation is particularly evident in the distribution of A 1g projected displacements, which indicates that local breathing distortions persist at these doping levels. Once the system is heavily doped ( ⟨ n ⟩ = 1 4), it becomes more difficult to resolve the bifurcation . in ⟨ X ⟩ (Fig. 4 e ), but a slight splitting can still be resolved in the projected coupling (Fig. 4 ). j (This is particularly evident when one examines the differences between the fits and the data, as shown in the lower parts of Figs. 4 i and 4 .) j The fact that local structural distortions are more clearly detectable in the projected measurements for ⟨ n ⟩ ≥ 1 3 . suggests that they are fluctuating in na- ture.

The trends observed in the simulated data are entirely consistent with the PDF results, which strongly suggests that e -ph coupling is the likely mechanism for the persistence of the local structural distortions with doping. They are also consistent with earlier DQMC [22] and semi-classical Monte Carlo [23] simulations for related 2D models, but here for a realistic model for the 3D material. We stress that our model is derived from uncorrelated DFT calculations and does not include any potential enhancement of the e -ph coupling by the longrange Coulomb interaction [15, 19]. If we included these effects, we would likely obtain even clearer bifurcations at higher hole concentrations in the histograms. Given the relatively weak local symmetry breaking observed in Fig. 2 c , it may be possible to place an upper bound on the strength of the e -ph coupling, which would provide valuable insights into the nature of pairing at lower temperatures.

FIG. 6. Evolution of the electronic gap with potassium doping. a Cartoon sketch illustrating measurements of the unoccupied and occupied oxygen PDOS via XAS and RIXS, respectively. b Evolution of the electronic structure with doping as the indirect semiconducting gap between the occupied and the unoccupied state decreases with increasing doping. c Dependence of the energy gap on potassium doping estimated from XAS and RIXS data, compared with gap values extracted from optical absorption and magnetic susceptibility data [33, 34].

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00401v1-figure-b2dd63b2-96c5-48f5-aae0-029c2f4bfdd1.jpg)

## Signatures of bond disproportionation

RIXS experiments were performed on thin films of BKBO grown by pulsed laser deposition (PLD). By using thin film samples, we were able to overcome the challenge of synthesizing sufficiently large single crystals of BKBO in the under- to optimally-doped regime ( x = 0 to 0 4) and . avoid the spectral degradation that might otherwise occur due to the surface roughness of cleaved crystals. The large surfaces of the films also facilitate beam rastering to mitigate potential sample damage.

Figure 5 a maps the parent ( x = 0) compound's RIXS spectra as a function of the incident photon energy. The excitation energy is tuned over the lowest-lying O K -edge peak measured by x-ray absorption spectroscopy (XAS), as shown in Fig. 5 e . The RIXS spectra consist of two main components exhibiting distinct behaviors as a function of the excitation energy. The first and most obvious component is the prominent, broad peak centered roughly around -4 eV energy loss, with a weaker flat feature extending to an energy loss of about -1 eV. Increasing the incident photon energy directly shifts this portion of the RIXS spectrum to higher energy loss (as seen in Fig. 5 a ), indicating that these are valence band fluorescence features. The second key component of the spectrum, highlighted in Figs. 5 b c -, is an asymmetric broadening near the elastic line, which is a signature of low-energy excitations, such as excitations of the lattice due to e -ph coupling [29-32]. This feature is most evident at an excitation energy of 527.9 eV (red curve in Fig. 5 c ), where the inelastic portion of the spectrum appears largest compared to the elastic peak.

The fluorescence features provide a useful window into BKBO's electronic structure that affirms the relevance of negative charge transfer energies in the system by showing that the gap opens in predominantly oxygenderived states [35]. Fig. 5 d compares the RIXS spectrum of BBO with data from angle-integrated photoemission spectroscopy (PES). The spectra are aligned and scaled to roughly match features in the high energy loss/deep binding energy region. The PES data was acquired from an in situ BBO film, as described in Ref. [14]. The PES signal derives from the density of states (DOS), albeit modulated by photoexcitation matrix elements [36]. The RIXS fluorescence spectrum, on the other hand, is due to the decay of O 2 p states to fill the O 1 s core hole and thus relates to the oxygen partial density of states (PDOS). The overall similarities between these two spectra highlight the profound degree of oxygen hybridization throughout the band structure [14, 21, 35]. The differences in the spectra of Fig. 5 d are just as revealing. In particular, compared to the photoemission spectrum, the RIXS fluorescence exhibits an enhanced shelf-like protrusion in the region of about -1 to -3 eV energy loss. The band dispersion of these highest occupied states is seen in an angle-resolved view of the photoemission spectrum [14], but when the data is angle- (i.e., momentum-) integrated to obtain the PES spectrum, the intensity of this feature is almost totally drowned out by the height of the main peak near -3.5 eV binding energy. The pronounced appearance of the shelf-like feature in RIXS compared to PES implies a resonant enhancement, indicating that oxygen states dominate at the gap edge. A similar resonant enhancement can be seen on the right side of the main peak. Both enhancements are consistent with the calculated oxygen PDOS of BBO [35]. Hence, the RIXS fluorescence paints a picture of strong Bi-O hybridization accompanied by a gap opening in predominantly O-derived states, consistent with the model of a bond-

FIG. 7. Probing the e -ph coupling with oxygen K -edge RIXS a . Fitted phonon harmonics for the BBO ( x = 0) sample. The dashed lines show the contributions from each phonon harmonic while the solid line is the sum. b Calculated RIXS intensity obtained from our cluster calculations as a function of the e -ph coupling strength g . Increasing g results in a longer tail of lattice excitations extending to higher energy loss values. c Comparison of O K -edge RIXS spectra of the parent compound BBO and the predictions of our model for g = 200 meV. d Comparison of the measured O K -edge RIXS spectra in BKBO for different values of x . All spectra have a similar shape with a significant tail extending to higher energy loss, indicative of a strong e -ph coupling. The black dashed line in this panel corresponds to the same model data plotted in panel c . The elastic line has been removed from the experimental data points to focus on the phonon features.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00401v1-figure-4e6d6772-b8e8-49c9-ba72-45bf66c73feb.jpg)

disproportionated phase generated by a reverse transfer of holes onto the ligands [21, 35].

These insights are reinforced by a comparison of the XAS spectrum with the fluorescence emission in RIXS. As illustrated schematically in Fig. 6 a , the total fluorescence yield XAS intensity is due to the transition from the O 1 s core level to the O 2 p states, thus relating to the oxygen PDOS of the unoccupied bands. On the other hand, as noted earlier, the fluorescence intensity of the RIXS spectrum stems from electron decay from the 2 p bands to empty core level states, consequently tracking the oxygen PDOS of the occupied bands. Thus, the XAS and RIXS fluorescence spectra complement each other and provide a means to follow the evolution of the indirect bandgap as a function of doping.

Figure 6 b displays the indirect gap ∆ determined from the XAS and RIXS spectra for x = 0, 0 1, . and 0 4. . Upon hole doping, the dominant absorption peak shifts to lower photon energy, while the leading edge of the fluorescence in the RIXS spectra shifts to lower (i.e., closer to zero) energy loss. The estimated gap between the XAS and RIXS fluorescence spectra decreases with doping, as summarized in Fig. 6 c , and accounts for the change in the indirect gap of the oxygen bands. The changes in the band's filling and the gap between the XAS and RIXS fluorescence leading edge are reasonably close to the indirect gap estimates obtained from optical absorption and magnetic susceptibility measurements [17, 33, 34, 37]. The fact that the gap closure is observed via the oxygen PDOS further affirms the dominant role of O 2 p states at the gap edge, which closes due to the hole doping introduced by the potassium substitution.

## Probing the e -ph coupling with RIXS

Figure 7 analyzes the near-elastic component of the RIXS spectra to study the e -ph coupling as a function of doping. The data were acquired with incident photon energy tuned to the maximum of the lowest-lying O K -edge absorption, as was shown in Fig. 5. We can decompose the spectra for each doping level into an elastic component and a single phonon excitation at ℏ Ω ≈ 65 meV and its overtones. We illustrate this by fitting the parent compound spectrum with a series of weighted harmonics in Fig. 7 a . (Supplementary Figure 3 shows the decomposition for the other doping levels.)

To gain more insight into the nature of these harmonic phonon excitations, we use exact diagonalization to calculate the low-energy RIXS spectrum of a BiO 6 octahedron with open boundary conditions within the KramersHeisenberg formalism. To reduce the Hilbert space associated with orbital and phonon subspaces, our model includes the Bi 6 s orbital, the A 1g combination of ligand oxygen orbitals, and the corresponding local breathing mode, which directly modulates the Bi-O hybridization (see Methods). All model parameters are identical to those used in our DQMC calculations, except the phonon energy ℏ Ω and e -ph coupling g , which we now treat as fitting parameters. We further fix the inverse core hole lifetime parameter Γ = 150 meV, which is typical for the O K -edge [30, 38].

As shown in Fig. 7 b , the magnitude of the e -ph coupling in our model determines the relative intensities of the phonon harmonics. We can thus estimate the coupling strength's value and doping dependence from the overall shape of the spectrum's tail, with some important caveats. Because we are approximating the system with a small local cluster, our approach is very similar in spirit to the single-site framework introduced by Ament et al . [29] but extended to off-diagonal e -ph interactions. Small cluster approaches like these neglect charge and orbital fluctuations beyond the cluster, which can introduce additional momentum dependence to the spectra [32]. Our approach also does not account for coupling to multiple phonon modes [39] or any potential coupling between the lattice and the core-hole [40], which might affect the shape of the spectra and estimates for the e -ph coupling.

We obtain good agreement with the data for a microscopic coupling of g = 200 meV, as shown in

Fig. 7 c . In particular, the overall line shape of the phonon excitations is well reproduced by the model which only couples to the local breathing model, with slight deviations appearing in energy losses above the third harmonic. (We attribute this discrepancy to charge fluctuations beyond the length scale of the cluster.) The microscopic coupling extracted from the fits is a factor of two larger than the one inferred from first principles calculations, which we attribute to limitations in the model. Nevertheless, the results suggest that the oxygen holes in BBO are coupled strongly to the breathing motion of the oxygen atoms. Crucially, we observe nearly identical profiles for the lattice excitations in the doped samples. Fig. 7 d illustrates this by overlaying the RIXS spectra for x = 0, 0 1, . and 0 4 samples together . with the model data. Here, we have subtracted the contribution from the elastic line and normalized each spectrum to the intensity at the first phonon line to enable direct comparisons. Since the relative intensities of the phonon harmonics (and hence the overall shape of each spectrum) are determined by the strength of the e -ph coupling, this result suggests that substantial coupling to the breathing distortions persists with hole doping.

## Discussion

We have used resonant inelastic x-ray and neutron total scattering techniques to study the evolution of BKBO's local crystal structure, electronic structure, and e -ph coupling with doping. Analysis of the neutron PDF finds evidence that local BiO 6 breathing distortions - residuals of the parent insulating phase - persist deep into the metallic phase, even at the highest measured K concentration ( x = 0 4), . which is at or slightly beyond the optimal doping level for superconductivity. Details of the O K -edge RIXS measurements affirm the influence of an effective negative charge transfer energy in BKBO, which drives the doped holes to reside primarily on the O 2 p orbitals. Low-energy excitation features in the RIXS data indicate that significant coupling between the oxygen holes and the Bi-O bond-stretching phonons exists throughout the investigated doping range. Small cluster exact diagonalization (ED) and state-ofthe-art DQMC simulations on extended lattices demonstrate that our measurements can be explained using a multi-orbital model in which the oxygen vibrations directly modulate the Bi-O hybridization. Our results thus provide compelling evidence for the presence of small BiO 6 (bi)polarons in BKBO's metallic phase, which persist through optimal doping.

The findings from our neutron scattering experiments stand in contrast to a recent x-ray PDF study of BKBO [28], which reported that short-range breathing distortions already vanish at finite doping levels below the insulator-metal transition, and instead found evidence of local displacements of the Bi atoms along the Bi-O axis that break inversion symmetry. Such findings, if verified, would severely constrain the relevance of

(bi)polaron interactions for superconductivity in the bismuthates. At present, we cannot explain definitively why our findings diverge so dramatically from Ref. [28], but the experiments differ in some notable respects. Compared to neutrons, the sensitivity of hard x-rays to the crucial O atoms is vastly diminished, and the O ( Z = 8) scattering contribution is dwarfed by the signals from Ba ( Z = 56) and Bi ( Z = 83). Another difference is that hard x-ray scattering should have a sizeable inelastic component, while the thermal neutrons in our measurements scatter totally elastically. The use of hard x-rays also carries a risk of radiation-induced sample damage. Finally, while our neutron experiments were performed on powder samples obtained by solid state synthesis, the x-ray PDF measurements were performed on single crystals grown by an electrochemical method, which in some studies have exhibited significant doping inhomogeneity [41].

The bismuthates are classified as negative charge transfer systems with significant hybridization between the Bi 6 s orbitals and ligand oxygen orbitals [20, 21, 35]. However, other families of perovskites share similar characteristics, whereby a lattice-coupled transfer of hole pairs onto the ligands should occur [42-46]. Our results are likely relevant to this much broader class of materials. For example, the notion that bipolarons may condense into insulating or superconducting states provides a natural framework for understanding these systems' unusual metal-insulator transitions and superconductivity. Here, the nature of the broken symmetry state at low temperatures would depend on the relative strengths of the cation-ligand hopping, charge transfer energy, and density of holes. For instance, in rare earth nickelates, which undergo metalinsulator transitions, recent neutron total scattering experiments found similar evidence of local breathing distortions [47, 48], and dielectric spectroscopy [48] and transport [49] measurements find evidence for a polaronic liquid metallic phase. In another likely related system, the recently discovered superconductivity in Ba 1 -x K SbO x 3 arises out of an ordered insulator phase in a similar manner as in BKBO [5]. However, the insulating phase in the antimonate persists to larger hole concentrations, probably due to differences in the charge transfer energy and Sb-O hybridization. Our results thus call for systematic studies of these materials focusing on local structural properties and e -ph coupling.

## Methods

## Neutron scattering experiments

The neutron scattering measurements were performed on polycrystalline BKBO samples synthesized via solidstate reaction methods. X-ray diffraction (XRD) measurements at room temperature confirmed the samples' phase. The atomic ratios in the various BKBO samples were measured by x-ray fluorescence spectroscopy (XRF). Both XRD and XRF indicate successful and homogeneous Ba/K substitution in all the studied samples.

Neutron total scattering experiments were performed on the Nano-scale Ordered Materials Diffractometer (NOMAD) at the Spallation Neutron Source (SNS) at Oak Ridge National Laboratory (ORNL) [50]. All samples were placed in quartz capillaries (2 mm diameter, 0.01 mm wall thickness) and exposed to the neutron beam for 48 minutes. Calibration of NOMAD's timeof-flight detectors was performed by scattering from diamond powder. The structure function S Q ( ) was obtained by subtracting the background (an empty quartz capillary) from the sample measurement and normalizing to scattering from a vanadium rod to account for the neutron spectrum and detector effects. The experimental pair distribution G r ( ) was extracted via the Fourier transform

$$G ( r ) = \frac { 2 } { \pi } \int Q [ S ( Q ) - 1 ] \sin ( Q r ) \ d Q, \quad ( 2 ) \quad \text{RIX}$$

where the scattering vector Q ranged from 0.2 to 31.4 ˚ A -1 . The GSAS-II program [51] was used to perform Rietveld refinement on bank 3 of the neutron diffraction patterns. An instrument parameter file was created using neutron diffraction data collected from Si powder. The PDFGui program [52] was used for smallbox refinement of the pair distribution functions. The goodness-of-fit is calculated via the relation

$$R _ { w } ^ { 2 } = \frac { \sum _ { i = 1 } ^ { N } w ( r _ { i } ) [ G _ { \text{obs} } ( r _ { i } ) - G _ { \text{calc} } ( r _ { i } ) ] ^ { 2 } } { \sum _ { i = 1 } ^ { N } w ( r _ { i } ) [ G _ { \text{obs} } ( r _ { i } ) ] ^ { 2 } }, \quad ( 3 )$$

where w represents the weight of the data points, G obs is the experimental data point, and G calc is the simulated data point. All small-box refinements were performed over the range r = 1 5 to 10 ˚ . . A

## RIXS experiments

RIXS experiments were performed on BKBO thin films grown by PLD, as these samples provide sufficient size, quality, and surface flatness over the investigated doping range. BKBO films with a thickness of 12 nm were deposited on SrTiO (001) 3 substrates from ablation targets with doping levels of x = 0, 0 1, . and 0 4. . The growth conditions are described in previous work [14, 53]. The growth was epitaxial along the c -axis, as demonstrated by reflection high-energy electron diffraction (RHEED) measurements. The (pseudo)cubic lattice parameter was determined by room temperature XRD to be a = 4 31 ˚ and . A a = 4 28 ˚ for the doping levels . A of x = 0 and 0 4, respectively. . These values correspond to slight compressive strains of 1% ( x = 0) and 0.1% ( x = 0 4), indicating that the films are near bulklike in . their structure [54]. Transport measurements show a superconducting transition temperature T c = 20 K for the x = 0 4 sample. This value is typical of thin BKBO films . but lower than the maximal T c ≈ 30 K in bulk samples. The lower-doped samples are insulating.

The RIXS experiments were performed at the ADRESS beamline of the Swiss Light Source [55].

The samples were transferred from the film deposition chamber to the RIXS system via a vacuum suitcase with pressure better than 10 -7 mbar to minimize surface contamination and preserve the surface quality. The sample temperature during the experiments was 17 K unless otherwise noted. The energy resolution was determined to be 70 meV by measuring the full width at half maximum (FWHM) of the elastic line on carbon tape. All RIXS measurements were performed at the oxygen K -edge under fixed geometry corresponding to the scattering wavevector q = (0 24 0 0 17) . , , . r.l.u. The positioning of the sample was shifted after each measurement to minimize the risk of damage from irradiation.

## RIXS calculations

The lattice excitations in the RIXS spectra were modeled within the Kramers-Heisenberg formalism, where the initial | i ⟩ , intermediate | m ⟩ , and final states | f ⟩ of the scattering process are obtained from ED calculations. We approximate the system in the insulating phase of BBO with a single BiO 6 octahedron. The cluster Hamiltonian, written in hole-language, is H = H e + H ph + H e -ph , where

$$H _ { e } & = \epsilon _ { s } \sum _ { \sigma } s _ { \sigma } ^ { \dagger } s _ { \sigma } + \epsilon _ { p } \sum _ { \sigma, \nu } p _ { \nu, \sigma } ^ { \dagger } p _ { \nu, \sigma } \\ & + t _ { s p } \sum _ { \sigma, \nu } \left ( Q _ { \nu } s ^ { \dagger } p _ { \nu, \sigma } + h. c. \right )$$

describes the electronic sector,

$$H _ { \text{ph} } = \hbar { \Omega } \sum _ { \nu } \left ( b _ { \nu } ^ { \dagger } b _ { \nu } + \frac { 1 } { 2 } \right )$$

describes the phononic sector and

$$H _ { e - \text{ph} } = g \sum _ { \nu, \sigma } \left [ s ^ { \dagger } p _ { \nu, \sigma } ( b _ { \nu } ^ { \dagger } + b _ { \nu } ) + \text{h.c.} \right ], \quad ( 6 )$$

describes the e -ph coupling due to the modulation of the Bi-O hopping integrals. Here, s † σ creates a spinσ hole on the Bi 6s orbital; p † ν,σ creates a hole in each the O 2p ν orbital ( ν = ± x, ± y, ± z ) oriented toward the Bi atom; b † ν creates a phonon on the ligand oxygen site ν ; ϵ s and ϵ p are the on-site energies of the 6 s and 2 p orbitals, respectively; ℏ Ω is the phonon energy; t sp is the Bi-O hopping integral with phase factors Q ± y = Q ± y = Q ± z = ∓ 1 (see Fig. 1a of Ref. [16]). For simplicity, we have neglected the O-O hopping as it does not play a major role in the resulting model.

Next, we introduce a molecular orbital basis for the O

orbitals [20]

$$\text{orbitals [20]} & & \text{is the trial, i} \\ L _ { s } & = \frac { 1 } { \sqrt { 6 } } \left ( p _ { x } + p _ { y } + p _ { z } - p _ { - x } - p _ { - y } - p _ { - z } \right ), & \text{is the for} \\ L _ { x } & = \frac { 1 } { \sqrt { 2 } } \left ( p _ { x } + p _ { - x } \right ), & \text{is the in the} \\ L _ { y } & = \frac { 1 } { \sqrt { 2 } } \left ( p _ { z } + p _ { - z } \right ), & \text{the in} \\ L _ { z } & = \frac { 1 } { \sqrt { 2 } } \left ( p _ { z } + p _ { - z } \right ), & \text{The i} \\ L _ { x ^ { 2 } - y ^ { 2 } } & = \frac { 1 } { 2 } \left ( p _ { x } - p _ { y } - p _ { - x } + p _ { - y } \right ), \text{ and} & \text{obtain} \\ L _ { 3 z ^ { 2 } - r ^ { 2 } } & = \frac { 1 } { 2 \sqrt { 3 } } \left ( p _ { x } + p _ { y } - 2 p _ { z } - p _ { - x } - p _ { - y } + 2 p _ { - z } \right ), & \text{not only 1 at ou functi} \\ \text{where we have suppressed the spin index for brevity. We} \\ \text{can define analogous transformations for the phonon op-} \\ \text{arrow} \quad \text{$\mathcal{New}$} \quad \text{$\mathcal{new}$} \quad \text{$\mathcal{new}$} \quad \text{$\mathcal{new}$} \quad \text{$\mathcal{new}$}$$

where we have suppressed the spin index for brevity. We can define analogous transformations for the phonon operators. For example, the annihilation operator for the breathing mode is given by

$$B _ { s } = \frac { 1 } { \sqrt { 6 } } \left ( b _ { x } + b _ { y } + b _ { z } - b _ { - x } - b _ { - y } - b _ { - z } \right ).$$

This operator corresponds to the local breathing distortion with A 1g symmetry, as defined in Eq. (1). When written in this molecular basis, the cluster Hamiltonian simplifies to

$$\text{sphines to} \\ H & = \epsilon _ { s } \sum _ { \sigma } s _ { \sigma } ^ { \dagger } s _ { \sigma } + \epsilon _ { p } \sum _ { \sigma, \nu } L _ { \nu, \sigma } ^ { \dagger } L _ { \nu, \sigma } + \Omega \sum _ { \nu } \left ( B _ { \nu } ^ { \dagger } B _ { \nu } + 1 / 2 \right ) \, \text{smail} \, \\ & - \sqrt { 6 } t _ { s p } \sum _ { \sigma } \left ( s ^ { \dagger } L _ { s, \sigma } + h. c. \right ) \\ & + g \sum _ { \nu, \sigma } \left [ s _ { \sigma } ^ { \dagger } L _ { \nu, \sigma } \left ( B _ { \nu } ^ { \dagger } + B _ { \nu } \right ) + h. c. \right ] \,, \, \text{the} \, \\ & \quad \text{period} \, \\ \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \hat { \ \hat { \hat { \hat { \hat { \hat { \hat { \hat { \ \hat { \ \hat { \ \hat { \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \ \hat { \ \hat { \ \hat { \ \hat { \ \ \hat { \ \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \hat { \ \ \hat { \ \hat { \ \hat { \ \ \hat { \ \hat { \ \ \hat { \ \ \hat { \ \ \hat { \ \hat { \ \ \hat { \ \ \hat { \ \ \hat { \ \ \hat { \ \ \hat { \ \ \hat { \ \ \hat { \ \ \ \hat { \ \ \hat { \ \ \ \ \ \ \ \ \ \ } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } ) } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } { } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } \ } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } } { \text{sphinx} \, \, \,$$

(8)

where the index ν = s, x, y, z, x 2 -y 2 , and 3 z 2 -r 2 . We now further reduce the system to an effective two-level problem by neglecting all molecular orbitals that don't hybridize directly with the Bi orbitals

$$\text{hybridize directly with the $B$} \text{ orbitals} \\ H & = \epsilon _ { s } \sum _ { \sigma } s _ { \sigma } ^ { \dagger } s _ { \sigma } + \epsilon _ { p } \sum _ { \sigma } L _ { s, \sigma } ^ { \dagger } L _ { s, \sigma } + h \Omega \left ( B _ { s } ^ { \dagger } B _ { s } + \frac { 1 } { 2 } \right ) \text{                                                                                                                                                                                                        
  & + \sqrt { 6 } t _ { s p } \sum _ { \sigma } \left ( s ^ { \dagger } L _ { s, \sigma } + h. c. \right ) \text{                                                                                                                                                                                                 
  & \text{                                                                                                                                                                                                    
  & \text{$\text{dynamic}$} \text{dynamic}$} \text{the d} \text{in Rel} \\ & + g \sum _ { \sigma } \left [ s _ { \sigma } ^ { \dagger } L _ { s, \sigma } \left ( B _ { s } ^ { \dagger } + B _ { s } \right ) + h. c. \right ] \text{                                                                                                                                                                                                    

  & \text{Data} \\ \text{Eq. (9)} \text{ can now be easily diagonalized while retaining} \text{                                                                                                                                                                                                     The c}$$

Eq. (9) can now be easily diagonalized while retaining a large number of phonon quanta N max for the oxygen breathing mode.

We compute the RIXS intensity at the O K -edge using the Kramers-Heisenberg formula in this orbital basis

$$I ( \omega ) = \sum _ { f } | M _ { f, g } | ^ { 2 } \, \delta ( E _ { f } - E _ { i } - \omega ) \quad \ ( 1 0 ) \quad \text{Co} \\ \text{A } \, \text{J}$$

where ω is the energy transfer,

$$M _ { f g } = \sum _ { \sigma } \frac { \left \langle f \right | L _ { \sigma } ^ { \dagger } \left | n \right \rangle \left \langle n \right | L _ { \sigma } \left | g \right \rangle } { E _ { g } - E _ { n } + \omega _ { \text{in} } + \text{i} \Gamma } \quad \left ( 1 1 \right ) \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{  \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{$$

is the scattering amplitude, | g ⟩ , | n ⟩ , and | f ⟩ are the initial, intermediate, and final states of the scattering process, and Γ is the inverse core-hole lifetime.

For numerical calculations, we set (in units of eV) t sp = 2 31, . ϵ s = 6 23, . ϵ p = 4 14 based on DFT calcula-. tions [16, 21]. When evaluating Eq. (11), we tuned the incident photon energy to the location of the maximum in the computed XAS and adopted Γ = 150 meV [30]. The initial and final states are obtained by diagonalizing Eq. (9) in the two-hole, m z = 0 sector while retaining N max = 100 phonon quanta. The intermediate states are obtained by diagonalizing Eq. (9) in the one-hole sector and neglecting the influence of the O 1 s core hole, which only produces a trivial shift in the eigenvalue spectrum at our level of modeling. The energy-conserving δ -function in Eq. (10) was approximated with a Gaussian line shape with a standard deviation of σ = 30 meV, as set by the leading edge of the elastic line. Since the overall weight of the elastic line is determined by extrinsic factors like surface roughness, we re-scaled its contribution in our calculation by a factor of 0.17 to reproduce the experimental weight and normalized the spectra to their intensity at the first phonon line to compare to the data. The RIXS data is best reproduced for ℏ Ω = 65 meV and g = 0 2 eV, with the latter being . a factor of two larger than the value inferred from DFT calculations [16, 56]. This difference is not surprising; a small cluster suppresses band effects, which can impact the intensity of lattice excitations [32].

## DQMC Simulations

The DQMC simulations of the BKBO models were performed on L = 4 lattice at T = 300 K using the SmoQyDQMC.jl package [57, 58]. Simulations were run using 12 parallel walkers, each performing 5,000 hybrid Monte Carlo (HMC) updates to thermalize the system, followed by an additional 10,000 updates, after each of which measurements were made. The HMC update consisted of N t = 3 time-steps with corresponding step-size ∆ = t π/ (2 N t Ω). The chemical potential was dynamically updated during the simulation to achieve the desired hole density using the algorithm presented in Ref. [59].

## Data Availability

The data supporting this study has been deposited in a repository that will be made public once the final version of the paper is accepted for publication. Until then, data will be made available upon request.

## Code Availability

A Julia notebook for the RIXS calculations can be found at INSERT LINK. The SmoQyDQMC.jl package can be obtained at https://github.com/SmoQySuite/ SmoQyDQMC.jl .

## Acknowledgements

M. N. was supported by the ISRAEL SCIENCE FOUNDATION (grant No. 2509/20). E. O'Q., M. L., and S. J. were supported by the National Science Foundation Materials Research Science and Engineering Center program through the UT Knoxville Center for Advanced Materials and Manufacturing (DMR-2309083). M. N., W. R. P., and N. C. P. were supported by the Swiss National Science Foundation through Project Nos. 200021 159678 and 200021 185037. E. P., D. M., Y. T., and T. S. were supported by the Swiss National Science Foundation through the NCCR MARVEL (grant number 141828) and the Sinergia project Mott Physics Beyond the Heisenberg Model (MPBH) (grant number 160765). The synchrotron experiments were performed at the ADRESS beamline of the Swiss Light Source at the Paul Scherrer Institut (PSI). The neutron scattering experiments used resources at the Spallation Neutron Source, a DOE Office of Science User Facility operated by the Oak Ridge National Laboratory.

## Author Contributions

J. S., E. O'Q., and M. L. performed neutron scattering experiments. M. N., E. P., D. M., Y. T., W. P., M. R., N. C. P., and T. S. performed RIXS experiments. D. J. G. synthesized and characterized powder samples for the neutron experiments, as well as ablation target material for the thin films studied by RIXS. M. N. and M. R. grew and characterized the thin film samples. B. C-.S. performed DQMC simulations. S. J. performed RIXS calculations. M. N., S. J., and N. C. P. wrote the manuscript with input from all authors.

## Materials &amp; Correspondence

Requests for materials and correspondence should be addressed to S. J. (sjohn145@utk.edu) or N. C. P. (nicholas.plumb@psi.ch).

## Competing Interests

The authors declare no competing interests.

## References

- [1] Keimer, B., Kivelson, S. A., Norman, M. R., Uchida, S. &amp; Zaanen, J. From quantum matter to high-temperature superconductivity in copper oxides. Nature 518 , 179-186 (2015). URL https://doi.org/10.1038/nature14165 .
- [2] Stewart, G. R. Superconductivity in iron compounds. Rev. Mod. Phys. 83 , 1589-1652 (2011). URL https: //link.aps.org/doi/10.1103/RevModPhys.83.1589 .
- [3] Plokhikh, I. et al. Discovery of charge order above room-temperature in the prototypical kagome superconductor La(Ru 1 -x Fe x ) 3 Si 2 . Communications Physics 7 , 182 (2024). URL https://doi.org/10.1038/ s42005-024-01673-y .
- [4] Sleight, A. W. Bismuthates: BaBiO 3 and related superconducting phases. Physica C 514 , 152165 (2015). URL https://www.sciencedirect.com/ science/article/pii/S0921453415000398 .
- [5] Kim, M. et al. Superconductivity in (Ba,K)SbO . 3 Nat. Mater. 21 , 627-633 (2022). URL https://doi.org/10. 1038/s41563-022-01203-7 .
- [6] Fradkin, E., Kivelson, S. A. &amp; Tranquada, J. M. Colloquium: Theory of intertwined orders in high temperature superconductors. Rev. Mod. Phys. 87 , 457482 (2015). URL https://link.aps.org/doi/10.1103/ RevModPhys.87.457 .
- [7] Cava, R. J. et al. Superconductivity near 30 K without copper: the Ba 0 6 . K 0 4 . BiO 3 perovskite. Nature 332 , 814816 (1988). URL https://doi.org/10.1038/332814a0 .
- [8] Nagata, Y., Mishiro, A., Uchida, T., Ohtsuka, M. &amp; Samata, H. Normal-state transport properties of Ba 1 -x K BiO x 3 crystals. J. Phys. Chem. Solids 60 , 19331942 (1999). URL https://www.sciencedirect.com/ science/article/pii/S0022369799002176 .
- [9] Ayres, J. et al. Incoherent transport across the strange-metal regime of overdoped cuprates. Nature 595 , 661-666 (2021). URL https://doi.org/10.1038/ s41586-021-03622-z .
- [10] Kondoh, S., Sera, M., Ando, Y. &amp; Sato, M. Normal state properties and oxygen isotope effect of (Ba,K)BiO . 3 Physica C 157 , 469-477 (1989). URL https://www.sciencedirect.com/science/article/ pii/0921453489902724 .
- [11] Zhao, G.-m. et al. Large oxygen-isotope effect in Sr 0 4 . K 0 6 . BiO : 3 Evidence for phonon-mediated superconductivity. Phys. Rev. B 62 , R11977-R11980 (2000). URL https://link.aps.org/doi/10.1103/ PhysRevB.62.R11977 .
- [12] Snezhko, A. &amp; Prozorov, R. Temperature dependence of magnetic penetration depth in Bi 1 -x K BiO x 3 superconductor. Physica C 403 , 200-202 (2004). URL https://www.sciencedirect.com/science/article/ pii/S0921453403017040 .
- [13] Shang, T. et al. Strongto weak-coupling superconductivity in highT c bismuthates: Revisiting the phase diagram via µ SR. Phys. Rev. B 101 , 014508 (2020). URL https://link.aps.org/doi/10.1103/ PhysRevB.101.014508 .
- [14] Plumb, N. C. et al. Momentum-resolved electronic structure of the highT c superconductor parent compound BaBiO . 3 Phys. Rev. Lett. 117 , 037002 (2016). URL https://link.aps.org/doi/10.1103/ PhysRevLett.117.037002 .
- [15] Wen, C. H. P. et al. Unveiling the superconducting mechanism of Ba 0 51 . K 0 49 . BiO . 3 Phys. Rev. Lett. 121 , 117002 (2018). URL https://link.aps.org/doi/10. 1103/PhysRevLett.121.117002 .
- [16] Cohen-Stead, B., Barros, K., Scalettar, R. &amp; Johnston, S. A hybrid Monte Carlo study of bond-stretching electronphonon interactions and charge order in BaBiO . 3 npj Comput. Mater. 9 , 40 (2023). URL https://doi.org/ 10.1038/s41524-023-00998-6 .
- [17] Federici, J. F., Greene, B. I., Hartford, E. H. &amp; Hellman, E. S. Optical characterization of excited states in BaBiO . 3 Phys. Rev. B 42 , 923-926 (1990). URL https: //link.aps.org/doi/10.1103/PhysRevB.42.923 .
- [18] Franchini, C., Kresse, G. &amp; Podloucky, R. Polaronic hole trapping in doped BaBiO . 3 Phys. Rev. Lett. 102 , 256402 (2009). URL https://link.aps.org/doi/10. 1103/PhysRevLett.102.256402 .
- [19] Yin, Z. P., Kutepov, A. &amp; Kotliar, G. Correlationenhanced electron-phonon coupling: Applications of GW

and screened hybrid functional to bismuthates, chloronitrides, and other highT c superconductors. Phys. Rev. X 3 , 021011 (2013). URL https://link.aps.org/doi/10. 1103/PhysRevX.3.021011 .

- [20] Khazraie, A., Foyevtsova, K., Elfimov, I. &amp; Sawatzky, G. A. Bond versus charge disproportionation in the bismuth perovskites. Phys. Rev. B 98 , 205104 (2018). URL https://link.aps.org/doi/10.1103/ PhysRevB.98.205104 .
- [21] Khazraie, A., Foyevtsova, K., Elfimov, I. &amp; Sawatzky, G. A. Oxygen holes and hybridization in the bismuthates. Phys. Rev. B 97 , 075103 (2018). URL https://link. aps.org/doi/10.1103/PhysRevB.97.075103 .
- [22] Li, S. &amp; Johnston, S. Quantum Monte Carlo study of lattice polarons in the two-dimensional threeorbital Su-Schrieffer-Heeger model. npj Quantum Mater. 5 , 40 (2020). URL https://doi.org/10.1038/ s41535-020-0243-3 .
- [23] Jiang, M., Sawatzky, G. A., Berciu, M. &amp; Johnston, S. Polaron and bipolaron tendencies in a semiclassical model for hole-doped bismuthates. Phys. Rev. B 103 , 115129 (2021). URL https://link.aps.org/doi/10. 1103/PhysRevB.103.115129 .
- [24] Menushenkov, A. P., Troyan, I. A. &amp; Eremets, M. I. Resonant raman scattering in superconducting Ba 1 -x K BiO . x 3 JEPT Lett. 77 , 521-525 (2003).
- [25] Braden, M., Reichardt, W., Shiryaev, S. &amp; Barilo, S. Giant phonon anomalies in the bond-stretching modes in doped BaBiO : 3 comparison to cuprates manganites and nickelates. Physica C 378-381 , 8996 (2002). URL https://www.sciencedirect.com/ science/article/pii/S0921453402013886 .
- [26] Kim, B., Kim, Y., Kim, H.-T., Kang, K.-Y. &amp; Lee, J. EXAFS observation of two distinct Bi-O distances below T c for a Ba 0 6 . K 0 4 . BiO 3 single crystal. Physica C 392-396 , 286-290 (2003). URL https://www.sciencedirect. com/science/article/pii/S0921453403009353 .
- [27] Tajima, S., Yoshida, M., Koshizuka, N., Sato, H. &amp; Uchida, S. Raman-scattering study of the metalinsulator transition in Ba 1 -x K BiO . x 3 Phys. Rev. B 46 , 1232-1235 (1992). URL https://link.aps.org/doi/ 10.1103/PhysRevB.46.1232 .
- [28] Griffitt, S. et al. Local inversion-symmetry breaking in a bismuthate high-T c superconductor. Nat. Commun. 14 , 845 (2023). URL https://doi.org/10.1038/ s41467-023-36348-9 .
- [29] Ament, L. J. P., van Veenendaal, M. &amp; van den Brink, J. Determining the electron-phonon coupling strength from resonant inelastic x-ray scattering at transition metal L -edges. Europhys. Lett. 95 , 27008 (2011). URL https: //doi.org/10.1209/0295-5075/95/27008 .
- [30] Lee, W. S. et al. Role of lattice coupling in establishing electronic and magnetic properties in quasi-onedimensional cuprates. Phys. Rev. Lett. 110 , 265502 (2013). URL https://link.aps.org/doi/10.1103/ PhysRevLett.110.265502 .
- [31] Devereaux, T. P. et al. Directly characterizing the relative strength and momentum dependence of electronphonon coupling using resonant inelastic x-ray scattering. Phys. Rev. X 6 , 041019 (2016). URL https: //link.aps.org/doi/10.1103/PhysRevX.6.041019 .
- [32] Bieniasz, K., Johnston, S. &amp; Berciu, M. Beyond the single-site approximation modeling of electron-phonon coupling effects on resonant inelastic x-ray scattering
- spectra. SciPost Phys. 11 , 062 (2021). URL https: //scipost.org/10.21468/SciPostPhys.11.3.062 .
- [33] Kozlov, M. E., Ji, X., Minami, H. &amp; Uwe, H. Change of IR spectra of Ba 1 -x K BiO x 3 crystals at the lowtemperature phase transition. Phys. Rev. B 56 , 1221112216 (1997). URL https://link.aps.org/doi/10. 1103/PhysRevB.56.12211 .
- [34] Ahmad, J. &amp; Uwe, H. Small-polaron excitations in Ba 1 -x K BiO x 3 studied by optical reflectivity measurements. Phys. Rev. B 72 , 125103 (2005). URL https: //link.aps.org/doi/10.1103/PhysRevB.72.125103 .
- [35] Foyevtsova, K., Khazraie, A., Elfimov, I. &amp; Sawatzky, G. A. Hybridization effects and bond disproportionation in the bismuth perovskites. Phys. Rev. B 91 , 121114 (2015). URL https://link.aps.org/doi/10. 1103/PhysRevB.91.121114 .
- [36] Damascelli, A. Probing the electronic structure of complex systems by ARPES. Phys. Scr. 2004 , 61 (2004). URL https://dx.doi.org/10.1238/Physica.Topical. 109a00061 .
- [37] Uwe, H. &amp; Tachibana, K. Indirect excitons in threedimensional charge-density-wave system of BaBiO . 3 In Yamafuji, K. &amp; Morishita, T. (eds.) Advances in Superconductivity VII , 165-168 (Springer Japan, Tokyo, 1995). URL https://link.springer.com/chapter/10. 1007/978-4-431-68535-7\_36 .
- [38] Kotani, A. &amp; Toyozawa, Y. Theoretical Aspects of InnerLevel Spectroscopy , 169-229 (Springer Berlin Heidelberg, Berlin, Heidelberg, 1979). URL https://doi.org/10. 1007/978-3-642-81297-2\_4 .
- [39] Geondzhian, A. &amp; Gilmore, K. Generalization of the Franck-Condon model for phonon excitations by resonant inelastic x-ray scattering. Phys. Rev. B 101 , 214307 (2020). URL https://link.aps.org/doi/10. 1103/PhysRevB.101.214307 .
- [40] Geondzhian, A. &amp; Gilmore, K. Demonstration of resonant inelastic x-ray scattering as a probe of exciton-phonon coupling. Phys. Rev. B 98 , 214305 (2018). URL https://link.aps.org/doi/10.1103/ PhysRevB.98.214305 .
- [41] Minami, H., Uzumaki, T. &amp; Uwe, H. Characterization of Ba 1 -x K BiO x 3 -y synthesized by an electrochemical method with the influences of heat treatment. In Hayakawa, H. &amp; Enomoto, Y. (eds.) Advances in Superconductivity VIII , 429-432 (Springer Japan, Tokyo, 1996).
- [42] Park, H., Millis, A. J. &amp; Marianetti, C. A. Site-selective Mott transition in rare-earth-element nickelates. Phys. Rev. Lett. 109 , 156402 (2012). URL https://link.aps. org/doi/10.1103/PhysRevLett.109.156402 .
- [43] Johnston, S., Mukherjee, A., Elfimov, I., Berciu, M. &amp; 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). URL https:// link.aps.org/doi/10.1103/PhysRevLett.112.106404 .
- [44] Bisogni, V. et al. Ground-state oxygen holes and the metal-insulator transition in the negative charge-transfer rare-earth nickelates. Nat. Commun. 7 , 13017 (2016). URL https://doi.org/10.1038/ncomms13017 .
- [45] Dalpian, G. M., Liu, Q., Varignon, J., Bibes, M. &amp; Zunger, A. Bond disproportionation, charge selfregulation, and ligand holes in s -p and in d -electron ABX 3 perovskites by density functional theory. Phys.

Rev. B 98 , 075135 (2018). URL https://link.aps.org/ doi/10.1103/PhysRevB.98.075135 .

- [46] Benam, M. R., Foyevtsova, K., Khazraie, A., Elfimov, I. &amp; 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). URL https: //link.aps.org/doi/10.1103/PhysRevB.104.195141 .
- [47] Li, B. et al. Insulating pockets in metallic LaNiO . 3 Advanced Electronic Materials 2 , 1500261 (2016). URL https://onlinelibrary.wiley.com/doi/abs/10. 1002/aelm.201500261 .
- [48] Shamblin, J. et al. Experimental evidence for bipolaron condensation as a mechanism for the metalinsulator transition in rare-earth nickelates. Nat. Commun. 9 , 86 (2018). URL https://doi.org/10.1038/ s41467-017-02561-6 .
- [49] Tyunina, M., Savinov, M., Pacherova, O. &amp; Dejneka, A. Small-polaron transport in perovskite nickelates. Scientific Reports 13 , 12493 (2023). URL https://doi.org/ 10.1038/s41598-023-39821-z .
- [50] Neuefeind, J., Feygenson, M., Carruth, J., Hoffmann, R. &amp; Chipley, K. K. The Nanoscale Ordered MAterials Diffractometer NOMAD at the Spallation Neutron Source SNS. Nucl. Instrum. Meth. B 287 , 6875 (2012). URL https://www.sciencedirect.com/ science/article/pii/S0168583X12003291 .
- [51] Toby, B. H. &amp; Von Dreele, R. B. GSAS-II: the genesis of a modern open-source all purpose crystallography software package. J. Appl. Crystallogr. 46 , 544-549 (2013). URL https://onlinelibrary.wiley.com/doi/ abs/10.1107/S0021889813003531 .
- [52] Farrow, C. L. et al. PDFfit2 and PDFgui: computer programs for studying nanostructure in crystals. J. Phys. Condens. Matter 19 , 335219 (2007). URL https://dx.
- doi.org/10.1088/0953-8984/19/33/335219 .
- [53] Naamneh, M. et al. Cooling a polaronic liquid: Phase mixture and pseudogap-like spectra in superconducting Ba 1 -x K BiO . x 3 arXiv:1808.06135 (2018). URL https: //arxiv.org/abs/1808.06135 .
- [54] Pei, S. et al. Structural phase diagram of the Ba 1 -x K BiO x 3 system. Phys. Rev. B 41 , 4126-4141 (1990). URL https://link.aps.org/doi/10.1103/ PhysRevB.41.4126 .
- [55] Strocov, V. N. et al. High-resolution soft x-ray beamline ADRESS at the Swiss Light Source for resonant inelastic x-ray scattering and angle-resolved photoelectron spectroscopies. J. Synchrotron Radiat. 17 , 631-643 (2010). URL https://doi.org/10.1107/S0909049510019862 .
- [56] Meregalli, V. &amp; Savrasov, S. Y. Electron-phonon coupling and properties of doped BaBiO . 3 Phys. Rev. B 57 , 14453-14469 (1998). URL https://link.aps.org/doi/ 10.1103/PhysRevB.57.14453 .
- [57] Cohen-Stead, B. et al. SmoQyDQMC.jl: A flexible implementation of determinant quantum Monte Carlo for Hubbard and electron-phonon interactions. SciPost Phys. Codebases 29 (2024). URL https://scipost.org/10. 21468/SciPostPhysCodeb.29 .
- [58] Cohen-Stead, B. et al. Codebase release r0.3 for SmoQyDQMC.jl. SciPost Phys. Codebases 29r0.3 (2024). URL https://scipost.org/10.21468/ SciPostPhysCodeb.29-r0.3 .
- [59] 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). URL https://link.aps.org/doi/10. 1103/PhysRevE.105.045311 .
Muntaser Naamneh Eugenio Paris Daniel McNally Yi Tseng Wojciech R. Pudelko Dariusz J. Gawryluk J. Shamblin Eric OQuinn Benjamin Cohen-Stead Ming Shi Milan Radovic M. Lang Thorsten Schmitt Steven Johnston Nicholas C. Plumb
Persistence of small polarons into the superconducting phase of Ba$_{1-x}$K$_x$BiO$_3$
Bipolaronic superconductivity is an exotic pairing mechanism proposed for
materials like Ba$_{1-x}$K$_x$BiO$_3$ (BKBO); however, conclusive experimental
evidence for a (bi)polaron metallic state in this material remains elusive.
Here, we combine resonant inelastic x-ray and neutron total scattering
techniques with advanced modelling to study the local lattice distortions,
electronic structure, and electron-phonon coupling ($e$-ph) in BKBO as a
function of doping. Data for the parent compound ($x = 0$) indicates that the
electronic gap opens in predominantly oxygen-derived states strongly coupled to
a long-range ordered breathing distortion of the oxygen sublattice. Upon
doping, short-range breathing distortions and sizable ($e$-ph) coupling persist
into the superconducting regime ($x = 0.4$). Comparisons with exact
diagonalization and determinant quantum Monte Carlo calculations further
support this conclusion. Our results provide compelling evidence that BKBO's
metallic phase hosts a liquid of small (bi)polarons derived from local
breathing distortions of the lattice, with implications for understanding the
low-temperature superconducting instability