Abstract
The electron–lattice interaction gives rise to a rich set of phenomena in quantum materials. Microscopically, this interaction often arises from the modulation of orbital overlaps; however, many theoretical studies neglect such couplings. Here, we present an exact diagonalization and determinant quantum Monte Carlo study of a three-orbital Su–Schrieffer–Heeger (SSH) model, on a two-dimensional Lieb lattice and in the negative charge transfer regime. At half-filling (one hole/unit cell), we observe a bipolaron insulating phase with a bond-disproportionate lattice. This phase is robust against moderate hole doping but is suppressed at large hole concentrations, leading to a metallic polaron-liquid-like state with fluctuating patches of local distortions. We also find an s-wave superconducting state at large hole doping that primarily appears on the oxygen sublattice. Our work provides a non-perturbative view of SSH-type couplings in two dimensions with implications for materials where such couplings are dominant.
Similar content being viewed by others
Introduction
Lattice dimerization, referring to an alternation of two lattice constants, occurs in many quantum materials, including the organic charge-transfer solids1,2,3, and perovskites systems like the rare-earth nickelates RNiO34,5,6, and the high-temperature (high-Tc) superconducting bismuthates Bi1−xKxBiO3 (BKBO)7,8. Lattice dimerization in one-dimension was initially described by the Su–Schrieffer–Heeger (SSH) model1, in which the atomic motion modulates the overlap integrals between neighboring atoms. Recently, SSH-like models have attracted further attention due to the fact that they can produce highly mobile polarons with light effective masses9, generate robust phonon-mediated pairing10, and even stabilize and control the location of a type-II Dirac point11. SSH interactions are also believed to be dominant mechanism for electron–phonon (e–ph) interactions in materials like the high-Tc cuprates12,13.
The SSH model has primarily been studied in one dimension (1D)9,10,14,15,16,17,18,19,20,21,22. But it is also necessary to study this model in higher dimensions in the context of materials like RNiO3 and BKBO. In this work, we focus on BKBO and study how the SSH-type e–ph interaction produces both insulating and superconducting states as a function of doping. BKBO is in the so-called “negative charge transfer” regime23,24,25,26, where holes self-dope from the cation to the ligand oxygen atoms. The subsequent hybridization between the cation and the oxygen atoms then leads to a sizable e–ph interaction6,7, which may be further enhanced by correlations27, and is believed to drive a high-temperature metal-to-insulator (MIT) transition. Here, the insulating state has a dimerized (or “bond disproportionated”) structure with expanded and collapsed BiO6 octahedra alternating through the material, and pairs of holes condensed into the molecular orbitals formed from the ligand oxygen orbitals with A1g symmetry6,7,25,26,28,29 (see also “Methods” section). The relevant model describing this situation is a multiorbital SSH model; however, little is currently known about the physics of such models due to a lack of suitable approaches for studying them in dimensions higher than one.
To address this problem, we carried out a combined determinant quantum Monte Carlo (DQMC) and exact diagaonalization (ED) study of a two-dimensional (2D) multiorbital model with SSH-type interactions. In this case, our DQMC simulations are made possible by the development of a fast update procedure for local spacetime moves of the phonon fields (see “Methods” section). Having BKBO in mind, we then deployed this approach to study a three-orbital SSH model defined on a Lieb lattice whose orbital basis consists of a Bi 6s orbital and O 2px and 2py orbitals. In this case, we freeze the heavier Bi atoms into place and restrict lighter O atoms to move along the bond directions, and study the model using non-perturbative DQMC and ED. Further details about the model, and the changes necessary for the DQMC algorithm are provided in the “Methods” section and in the supplementary materials (see Supplementary Note 1 for more detailed results).
Results
An overview of the phase diagram
An overview of the phase diagram inferred here is presented in Fig. 1. Near half-filling (one hole/unit cell), we find that the system is a bipolaronic charge-density-wave (CDW) insulator with a bond-disproportionated structure, similar to what is observed in BKBO30. Hole doping suppresses the insulating phase, giving way to a state where the lattice distortions have short-range correlations. At high doping levels, we find evidence for a metallic phase where holes are strongly correlated with local structural distortions, forming a polaron-liquid-like phase. Finally, at low temperatures, we find s-wave superconducting tendencies that form primarily on the oxygen sublattice. Our results are in qualitative agreement with the phase diagram of the bismuthate superconductors and provide theoretical support for a polaronic view of BKBO and other negative charge transfer oxides. Our simulations are performed in 2D, and quantitative differences may appear when simulating larger systems for the same model in three-dimensions (3D). But we believe that the underlying physics of the system should be the same in both cases. In what follows, we discuss the numerical results leading to this phase diagram and what it means for our understanding of quantum materials dominated by SSH-type interactions.
A molecular orbital viewpoint
Before proceeding to our DQMC results, we present a simplified molecular orbital analysis of a Bi2O4 cluster, which provides a more transparent view of the physics. We refer the reader to ref. 25 for a similar discussion of the 3D case from an ab initio perspective. (Note that ref. 25 uses electron language whereas we use hole language.)
The first step of our analysis is to expand the simple square unit cell to allow for two distinct Bi 6s orbitals and four O 2p orbitals, as indicated by the black dashed frame in Fig. 2. This expanded cell defines the cluster after we apply periodic boundary conditions. The two Bi 6s orbitals are denoted as s1 and s2. Next, we transform the four ligand oxygen orbitals into a molecular orbital basis (see “Methods” section). We also perform an analogous transformation for the phonon operators. Fig. 3a–d indicate the phases of the ligand 2pδ (δ = x, y) orbitals for each molecular orbital using ± signs, and the black arrows indicate the displacement patterns of the transformed phonon eigenmodes. The bond disproportionated structure that forms in the model corresponds to a coherent state of the optical \({x}_{{\bf{r}},{L}_{s}}\) phonon modes in this representation, while the \({x}_{{\bf{r}},{L}_{x}}\) and \({x}_{{\bf{r}},{L}_{y}}\) modes form the basis for the acoustic phonon modes.
We can glean several insights into the problem from this cluster model. In the atomic limit (tsp = tpp = 0) and in the negative charge transfer regime (ϵp < ϵs in hole language), the four molecular orbitals are degenerate, as shown in Fig. 3e. This degeneracy is lifted by the orbital overlaps: a nonzero tpp raises (lowers) the onsite energy of the Ls (Ld) molecular orbital, while a nonzero tsp hybridizes the Bi s and molecular Ls orbitals to form bonding (sLs), nonbonding \({(s{L}_{s})}^{0}\), and antibonding \({(s{L}_{s})}^{* }\) states. Here, the energy of the bonding state is lowered by 2tsp relative to the atomic values such that the two holes fill this state at half-filling, as shown in Fig. 3f. This ground state charge distribution is analogous to the one inferred for 3D bismuthates in ab initio calculations25 and ARPES measurements26.
The impact of the e–ph coupling is also evident from the form of the transformed Hamiltonian; holes hop between the Lγ molecular orbital and the Bi sites while exciting phonon eigenmodes with the same symmetry. At half-filling, the holes in the (sLs) bonding state will, therefore, excite the breathing phonon mode of the surrounding oxygen atoms. In an extensive system, this coupling can lead to a static breathing distortion of the lattice after a spontaneous symmetry breaking selects one of the Bi sublattices as the center of the compressed plaquettes. Upon doping, the additional holes will occupy the Ld and Lx,y orbitals, where they couple to the orthogonal phonon modes. Since the superposition of the individual modes determines the total displacement of the oxygen atoms, the breathing distortion will relax as the other modes are excited, even though the (sLs) holes remain coupled to the \({x}_{{L}_{s}}\) phonons.
To confirm this physical picture, we diagonalized the transformed Hamiltonian HM on a Bi2O4 cluster and evaluated several observables in the grand canonical ensemble with β = 14.56/tsp, Ω = tsp, and μ was adjusted to set the particle number. When diagonalizing this model, we included up to Nph = 5 quanta for each phonon mode, which was sufficient to obtain converged results for our choice of parameters.
Figure 4 summarizes the ED results. Figure 4a and b plot the evolution of the hole density \(\langle {\hat{n}}_{{L}_{\gamma }}\rangle\) on each molecular orbital, and the displacement fluctuations of each eigenmode \(\delta ({x}_{\delta })=\langle {\hat{x}}_{{L}_{\delta }}^{2}\rangle -{\langle {\hat{x}}_{{L}_{\delta }}\rangle }^{2}\), respectively, as a function of the hole concentration. As expected, holes primarily occupy the Ls orbital at half-filling. (The missing hole weight is split between the two Bi sites and is not shown.) At the same time, the displacement of the \({x}_{{L}_{s}}\) mode fluctuates significantly, while the remaining eigenmodes have fluctuations consistent with zero point motion. This behavior is also reflected in the expectation value of the phonon numbers (Fig. 4c), where the \({x}_{{L}_{s}}\) mode is excited while the remaining phonon modes are in their ground state. Note that we do not observe a distortion \(\langle {\hat{x}}_{{L}_{s}}\rangle\, \ne\, 0\) due to the absence of any symmetry breaking in the cluster; however, our DQMC simulations performed on larger lattices do find such a state.
When additional holes are introduced they occupy the Ld, Lx, and Ly molecular orbitals, as expected based on the level diagram shown in Fig. 3g. In this case, the Ld orbital has a larger hole occupation due to the finite value of tpp. At the same time, \(\delta ({x}_{{L}_{d}})\), \(\delta ({x}_{{L}_{x}})\), and \(\delta ({x}_{{L}_{d}})\) also increase linearly and the expectation value of the number of phonon quanta for these modes grows. Both the displacement fluctuations and the number of excited phonon quanta are comparable for the four phonon modes once the hole doping reaches \(\langle \hat{n}\rangle\) = 2.4 – 2.6. Finally, the introduction of additional holes slightly suppresses the magnitude of \(\delta ({x}_{{L}_{s}})\) and the total number of \({x}_{{L}_{s}}\) modes.
Our ED calculations suggest that hole doping induces a relaxation of the breathing distortion of the lattice that is dominant at half-filling. However, it does so by exciting the orthogonal phonon modes rather than by suppressing the number of \({x}_{{L}_{s}}\) quanta in the system. In this context, it is interesting then to determine if the Ls holes and \({x}_{{L}_{s}}\) modes can be viewed of as a composite object (i.e., a polaron). We checked this idea in our ED calculations by computing the expectation value of the polaron \(\langle {\hat{N}}_{{\rm{p}}}\rangle\) and bipolaron \(\langle {\hat{N}}_{{\rm{bp}}}\rangle\) number operators (see the “Methods” section). Figure 4d plots the doping evolution of the \(\langle {\hat{N}}_{{\rm{p}}}\rangle\) and \(\langle {\hat{N}}_{{\rm{bp}}}\rangle\). We find that the ground state has a significant amount of polaron and bipolaron character, which persists to higher doping levels. Our ED results strongly suggest that the system hosts polaronic carriers, where holes occupying the Ls molecular orbitals are bound to local \({x}_{{L}_{s}}\) modes.
The molecular orbitals discussed here will of course form bands in the extended system. Nevertheless, much of our analysis still applies in this case. To illustrate this, Fig. 5 plots the non-interacting band structure of our model in the insulating phase, where the dimerized structure has been introduced by modifying the hopping integrals as \({t}_{sp}^{ij}={t}_{sp}[1+{(-1)}^{i+j}\times 0.3]\). Figure 5a and b provide fat band plots of the Ls and Ld molecular orbital weight, respectively, while Fig. 5c plots the total and orbitally resolved density of states (DOS). As can be seen in Fig. 5a, the occupied band below the Fermi level (E = 0) at half-filling is the bonding (sLs) band, which couples to the breathing motion of the lattice. The first band above the Fermi level is mostly of Ld and Lx,y orbital character, such that doped holes will predominantly couple to the corresponding phonon modes.
DQMC simulations of an extended lattice
The molecular orbital picture presented in the previous section provides an intuitive way of understanding the physics of the model. With this in mind, we now turn to DQMC simulations on an extended cluster with N = 4 × 4 Bi atoms (48 orbitals in total).
We first examine the insulating phase that forms at \(\langle \hat{n}\rangle =1\). Figure 6b–d show the lattice displacement correlation functions \(\langle {\hat{X}}_{{\bf{r}},x}{\hat{X}}_{0,x}\rangle ,\langle {\hat{X}}_{{\bf{r}},y}{\hat{X}}_{0,y}\rangle\), and \(\langle {\hat{X}}_{{\bf{r}},y}{\hat{X}}_{0,x}\rangle\), as a function of position at temperature \({(\beta {t}_{sp})}^{-1}=0.1\), which shows evidence of a static bond disproportionated structure. For example, both \(\langle {\hat{X}}_{{\bf{r}},x}{\hat{X}}_{0,x}\rangle\) and \(\langle {\hat{X}}_{{\bf{r}},y}{\hat{X}}_{0,y}\rangle\) alternate in sign following a checkerboard pattern while \(\langle {\hat{X}}_{{\bf{r}},y}{\hat{X}}_{0,x}\rangle\) alternates in sign along x- and y-directions but is constant along the diagonal. These results are consistent with the breathing distortion pattern sketched in Fig. 6a, as well as the observed lattice distortion that appears in the insulating phase of the bismuthates31,32,33.
Remaining at \(\langle \hat{n}\rangle =1\), we now examine the temperature evolution of this phase. Fig. 7a plots the conductivity weight σ(β/2) and orbital-resolved spectral weight as a function of temperature \({(\beta {t}_{sp})}^{-1}\) (see section “Methods”). At high-temperature σ(β/2) (black dots) initially increases as the temperature is lowered until reaching a maximum at \({(\beta {t}_{sp})}^{-1}\approx 0.2\) before it rapidly falls off signaling the formation of an insulating state. The orbital-resolved spectral weight βGγ,γ(r = 0, τ = β/2), where γ is the orbital index, also reflects this behavior. Above (βtsp)−1 = 0.2, βGs,s(r = 0, β/2) increases as temperature decreases while \(\beta {G}_{{p}_{x/y},{p}_{x/y}}({\bf{r}}=0,\beta /2)\) remains relatively flat. Below \({(\beta {t}_{sp})}^{-1}=0.2\), however, the spectral weights of all three orbitals decrease rapidly as the insulating state forms, signaling the removal of spectral weight from the Fermi level.
The formation of charge order in the insulating phase can also be observed in the charge susceptibility \({\chi }_{\gamma ,\gamma }^{{\rm{C}}}({\bf{q}})\). Figure 7b plots the temperature evolution of \({\chi }_{\gamma ,\gamma }^{{\rm{C}}}({\bf{q}})\) at q = (π, π), corresponding to the real space ordering inferred from Fig. 6. Below \({(\beta {t}_{sp})}^{-1}=0.2\), the charge correlations rapidly increase on the s orbital, while there is little change in the signal on the p orbitals. The transition temperature for the CDW can be estimated using the temperature at which the correlation length equals the lattice size (see Supplementary Note 3 for more detailed results). At half filling, the transition temperature is about \({(\beta {t}_{sp})}^{-1}=0.15\). We stress that such an estimate is likely an overestimate, as finite size effects can be quite strong on small clusters34.
The behavior observed in Fig. 7 implies that the average density on the O sublattice remains uniform above and below the MIT transition, while a density modulation forms on the Bi sites. We confirm this picture in the inset of Fig. 7b, which plots the equal time Bi–Bi charge density along the high symmetry directions of the cluster, where we find a clear (π, π)-density modulation. The fact that the charge order signal appears in the Bi orbital component can be understood once one recognizes that all of the oxygen orbitals in the system are equivalent, even in the bond disproportionated structure. In this case, the breathing distortions increase the hybridization between the Ls orbital and one of the Bi s orbitals at the expense of the other, creating a density modulation on the two Bi sites. From a charge disproportionation point of view, however, the density modulation shown in the inset of Fig. 7b corresponds to a charge transfer of only ~0.1 holes between the two Bi sites.
Next, we study effects of hole doping on the MIT at fixed temperature \({(\beta {t}_{sp})}^{-1}=0.1\). Figure 8a plots σ(β/2) as a function of filling, where it increases upon hole doping until \(\langle \hat{n}\rangle \approx 1.5\), and after which it levels off within error bars, indicating metallic behavior. Moreover, we also find evidence for the formation of mobile polarons in this region, where holes are bound to local breathing distortions of the oxygen sublattice. This behavior can be seen by examining the polaron operator \(\hat{p}({\bf{r}})={\hat{x}}_{{\bf{r}},{L}_{s}}({\hat{n}}_{{\bf{r}},s}+{\hat{n}}_{{\bf{r}},{L}_{s}})\), which is analogous to the operator considered during our ED analysis. Figure 8b plots the doping evolution of the number of polarons given by \(\frac{1}{N}{\sum }_{{\bf{r}}}\langle \hat{p}({\bf{r}})\rangle\), where it decreases as additional holes are introduced but remains nonzero even at the largest dopings. We conclude that a finite number of polarons are present in the system at all dopings, which are formed from holes in the Ls orbitals and local breathing phonons.
Considering our DQMC and ED results, we suggest that the suppression of long-range polaron and bipolaron correlations with hole doping should be induced by introducing other phonon modes rather than directly suppressing the breathing phonon mode. We also note that the orbital character and electronic structure of BaBiO3 computed with ab initio methods share many similarities to the plot shown in Fig. 5. This fact suggests that our analysis can be extended to the bulk 3D material.
To explore the doping evolution of polarons in real space, we measured the staggered polaron correlation function \(\langle {\chi }_{{\rm{P}}}({\bf{r}})\rangle ={(-1)}^{{r}_{x}+{r}_{y}}\langle \hat{p}({\bf{r}})\hat{p}(0)\rangle\), which is plotted in Figs. 9a–d for selected hole concentrations. At half filling, 〈χP(r)〉 is positive for all r, indicating that the polarons are frozen into a long-range (on the length scale of the cluster) two-sublattice order, consistent with the patterns inferred from Figs. 6 and 7. With increasing hole concentrations, 〈χP(r)〉 decreases at the larger distances, indicating an overall relaxation of the bond disproportionated structure but the persistence of short-range correlations. Such behavior may be indicative of nanoscale phase separation35,36; however, studies on large clusters will be needed to clarify this issue. Finally, in the high doping region, where the system is metallic (e.g. \(\langle \hat{n}\rangle >1.44\)), the correlations become very short-ranged and extend up to one or two lattice constants at most. The observation of a finite number of polarons at high doping, but with short-range structure in real space is consistent with the proposal that a significant number of holes remain in the Ls orbitals, where they locally excite the breathing modes of the lattice, creating a polaronic metallic phase with fluctuation local breathing distortions.
We also examined the doping evolution of bipolarons in the system by measuring the bipolaron number, defined as \(\frac{1}{N}{\sum }_{{\bf{r}}}\langle \hat{g}({\bf{r}})\rangle\), where \(\hat{g}({\bf{r}})={\hat{x}}_{{\bf{r}},{L}_{s}}({\hat{n}}_{{\bf{r}},s,\uparrow }+{\hat{n}}_{{\bf{r}},{L}_{s},\uparrow })({\hat{n}}_{{\bf{r}},s,\downarrow }+{\hat{n}}_{{\bf{r}},{L}_{s},\downarrow })\), and the staggered bipolaron correlation function \(\langle {\chi }_{{\rm{BP}}}({\bf{r}})\rangle ={(-1)}^{{r}_{x}+{r}_{y}}\langle \hat{g}({\bf{r}})\hat{g}(0)\rangle\), as a function of doping. When computing the latter quantity, we considered the signal on the Bi site by keeping only the terms proportional to \({\hat{n}}_{{\bf{r}},s,\uparrow }{\hat{n}}_{{\bf{r}},s,\downarrow }\). (This simplification is necessary due to the many number of terms generated by the Wick contraction of the product of \(\hat{g}({\bf{r}})\) operators. The fact that we see excess charge density on the Bi sites at the center of a breathing distortion provides some justification for this simplification).
Figure 8b plots the doping evolution of the bipolaron number operator. As with the polaron number, it is largest near half-filling and decreases slowly with doping. At large hole concentrations, however, it is still finite, suggesting that a significant number of bipolarons are present in the system. The staggered bipolaron correlation function is plotted in Fig. 9e–h. At \(\langle \hat{n}\rangle =1\), the bipolaron correlations are clear and long-ranged on the scale of the cluster. This result supports the interpretation that the insulating phase is a static bipolaron lattice. As the hole concentration increases, we find that the bipolaron correlations are entirely suppressed at all length scales, while a finite number of bipolarons are present in the system, as indicated in Fig. 8b. These results can again be easily understood if the metallic phase is a polaron liquid.
This scenario raises questions regarding possible superconductivity. We, therefore, computed the pair field susceptibility \({\chi }_{\gamma }^{{\rm{sc}}}\). Since \({\chi }_{{p}_{x}}^{{\rm{sc}}}\) and \({\chi }_{{p}_{y}}^{{\rm{sc}}}\) are the same, we use \({\chi }_{p}^{{\rm{sc}}}\) to denote pairing on the O atoms. Here, we find that the s-wave pairing is dominant, which is perhaps not surprising given that pairing is mediated solely by the e–ph interaction. Figure 10 plots \({\chi }_{\gamma }^{{\rm{sc}}}\) as a function of temperature at \(\langle \hat{n}\rangle =1.59\), and compares it against the dominant charge correlations χC(π, π/2) at this doping. All three susceptibilities increase with decreasing temperature, but \({\chi }_{p}^{{\rm{sc}}}\) dominates below \({(\beta {t}_{sp})}^{-1}\approx 0.04\). The superconducting Tc can be crudely estimated by fitting the inverse pairing susceptibility data with a functional form \({({\chi }^{{\rm{sc}}})}^{-1}=a\sqrt{T-{T}_{{\rm{c}}}}\). This functional form captures the rapid variation in \({({\chi }^{{\rm{sc}}})}^{-1}\) expected as T → Tc34, while other functional forms such as linear extrapolation would not. Extrapolating \({({\chi }_{p}^{{\rm{sc}}})}^{-1}\) to zero (as shown in the inset), yields Tc ≈ 382 K [\({(\beta {t}_{sp})}^{-1}\approx 0.0158\)]. As with the CDW, we stress that this value is artificially high, due to the large value of Ω used in our calculations and likely finite size effects. Nevertheless, our results provide evidence that the bipolaronic-rich metallic phase has a superconducting ground state. The pairing susceptibility is suppressed slightly at lower hole dopings, suggesting the presence of a superconducting dome.
Discussion
We have introduced a quantum Monte Carlo approach for studying bond phonons with SSH-type e–ph couplings in higher dimensions. While our approach has broad applications to many materials, we have used it to study a 2D three-orbital SSH model in the negative charge transfer regime. We obtained a phase diagram consistent with the bismuthate high-Tc superconductors. At half filling, we find a bond disproportionated state. Upon hole-doping, this state gives way to a polaron-liquid-like state with short-range correlations, consistent with proposals for nano-scale phase separation or strongly fluctuating lattice polarons in doped BKBO35,36,37,38,39,40. We also find s-wave superconducting tendencies, which primarily form on the oxygen sublattice. It would be interesting to contrast our results with those obtained from an effective single band model to fully gauge the importance of the oxygen orbitals.
Methods
Model
The unit cell for our multiorbital SSH model contains one Bi 6s orbital and the oxygen 2px,y orbitals oriented along each of the Bi–Bi bond directions, as shown in Fig. 1a. In what follows, the unit cells are indexed by r = nxa + nyb, where a = (a, 0), b = (0, a) are the primitive lattice vectors along x- and y-directions, respectively, a is the Bi–Bi bond length (and our unit of length). In addition, \(\delta ,{\delta }^{\prime}=\pm\! x\), ±y will be used to identify the oxygen atoms surrounding each Bi.
In the SSH model, the atomic displacements modulate the hopping integrals between the Bi and O atoms as \({t}_{sp}({Q}_{\delta }-\alpha {\hat{u}}_{{\bf{r}},\delta })\), where we have introduced the shorthand \({\hat{u}}_{{\bf{r}},x}={\hat{X}}_{{\bf{r}},x}\), \({\hat{u}}_{{\bf{r}},-x}={\hat{X}}_{{\bf{r}}-{\bf{a}},x}\), \({\hat{u}}_{{\bf{r}},y}={\hat{X}}_{{\bf{r}},y}\), and \({\hat{u}}_{{\bf{r}},-y}={\hat{X}}_{{\bf{r}}-{\bf{b}},y}\). The Hamiltonian is written as H = Hel + Hlat + Hint, where
Here, 〈…〉 denotes a sum over nearest-neighbor atoms, and the operators \({s}_{{\bf{r}},\sigma }^{\dagger }({s}_{{\bf{r}},\sigma })\) and \({p}_{{\bf{r}},\delta ,\sigma }^{\dagger }({p}_{{\bf{r}},\delta ,\sigma })\) create (annihilate) spin σ holes on the Bi 6s and O 2pδ orbitals, respectively. To simplify the notation, we have introduced shorthand notation pr,−x,σ = pr−a,x,σ and pr,−y,σ = pr−b,y,σ. The operators \({\hat{n}}_{{\bf{r}},\sigma }^{s}={s}_{{\bf{r}},\sigma }^{\dagger }{s}_{{\bf{r}},\sigma }\) and \({\hat{n}}_{{\bf{r}},\sigma }^{{p}_{\alpha }}={p}_{{\bf{r}},\alpha ,\sigma }^{\dagger }{p}_{{\bf{r}},\alpha ,\sigma }\) are the number operators for s and pα (α = x, y) orbitals, respectively. Finally, ϵs and ϵp are the onsite energies; μ is the chemical potential; tsp and tpp are the Bi–O and O–O hopping integrals in the cubic crystal; and α is the e–ph coupling constant. The phase factors are Qx = Qy = −Q−x−Q−y = 1, and Q±x,±y = Q±y,±x = −Q±x,∓y = −Q∓y,±x = 1. The motion of the O atoms is described by the atomic displacement (momentum) operators \({\hat{X}}_{{\bf{r}},\alpha }\) (\({\hat{P}}_{{\bf{r}},\alpha }\)). Here, M is the oxygen mass and K is the coefficient of elasticity between each Bi and O atom, and each O is linked by springs to its neighboring Bi atoms. Thus, bare phonon frequency is \(\Omega =\sqrt{2K/M}\).
Parameters
We adopted tsp = 2.08, tpp = 0.056, ϵs = 6.42, and ϵp = 2.42 (in units of eV), which were obtained from DFT calculations of BaBiO37. We also adopted \(\Omega =\sqrt{2}{t}_{sp}\) and α = 4a−1. For these parameters, the DQMC calculations give \(\frac{1}{N}{\sum }_{{\bf{r}}}\langle {\hat{X}}_{{\bf{r}},x}\rangle =\frac{1}{N}{\sum }_{{\bf{r}}}\langle {\hat{X}}_{{\bf{r}},y}\rangle =0\) and \(\frac{1}{N}{\sum }_{{\bf{r}}}\langle {\hat{X}}_{{\bf{r}},x}^{2}\rangle =\frac{1}{N}{\sum }_{{\bf{r}}}\langle {\hat{X}}_{{\bf{r}},y}^{2}\rangle =0.57{(a/4)}^{2}=0.0356{a}^{2}\) at half-filling, indicating that the oxygen atoms do not cross the bismuth atoms during the DQMC sampling. Our parameters are therefore in a physically reasonable region. (Here, we are limited to large Ω by long autocorrelation times41). Finally, our DQMC calculations are carried out on a square lattice with N = 4 × 4 Bi atoms (48 orbitals in total).
Transformation to the molecular orbital basis
The molecular orbital basis can be obtained after expanding the unit cell to allow for two distinct Bi 6s orbitals and four O 2p orbitals, as shown in Fig. 2. The molecular orbital basis is then obtained by defining
(The Ls and Ld operators correspond to the A1g and Eg orbitals in ref. 25.) It is also useful to define oxygen displacement operators in the new basis as
with analogous definitions for the momentum operators.
After introducing the new basis and applying periodic boundary conditions, the Hamiltonian HM for the Bi2O4 cluster is \({H}^{M}={H}_{{\rm{el}}}^{M}+{H}_{{\rm{lat}}}^{M}+{H}_{{\rm{int}}}^{M}\), where
Here, the sums on γ are taken over γ = s, d, x, y and \({\hat{n}}_{\sigma }^{{L}_{\gamma }}={L}_{\gamma ,\sigma }^{\dagger }{L}_{\gamma ,\sigma }\).
The polaron \({\hat{N}}_{{\rm{p}}}\) and bipolaron \({\hat{N}}_{{\rm{bp}}}\) number operators for the transformed cluster Hamiltonian are defined as
and
These operators measure the combined presence of holes in the (sLs) bonding orbital together with a compression of the ligand oxygen. (The same quantities are also measured in our DQMC calculations [see Fig. 8b].) The minus sign in front of the second terms accounts for the fact that a compression of the O atoms around the second Bi site corresponds to a negative displacement of the \({x}_{{L}_{s}}\) mode as we have defined it.
DQMC simulations
DQMC computes the expectation value of an observable \(\hat{O}\) in the grand canonical ensemble \(\langle \hat{O}\rangle ={Z}^{-1}{\rm{Tr}}\left[\hat{O}{\mathrm {{e}}}^{-\beta H}\right]\), where \(Z={\rm{Tr}}\left[{\mathrm {{e}}}^{-\beta H}\right]\) is the partition function42. Details of the DQMC algorithm can be found in refs. 42,43,44. Here, we outline the parts of the algorithm that need to be modified for the SSH model.
DQMC divides the imaginary time interval [0, β] into L discrete steps of length Δτ = β/L such that the partition function can be rewritten using the Trotter formula \(Z={\rm{Tr}}\left({\mathrm {{e}}}^{-\Delta \tau LH}\right)\approx {\rm{Tr}}{\left({\mathrm {{e}}}^{-\Delta \tau {H}_{{\rm{int}}}}{\mathrm {{e}}}^{-\Delta \tau K}\right)}^{L}\), where K = Hel + Hlat contains the non-interacting terms of the Hamiltonian and Hint contains the e–ph interaction. The Trotter approximation neglects terms of order \({\mathcal{O}}{(\Delta \tau )}^{2}\), which is controllable as Δτ → 0. Next, the trace over the phonon momenta and bilinear Fermion operators can be performed analytically to yield a result expressed as a product of matrix determinants44 and an integral over the remaining lattice displacements
Here, ∫dXx and ∫dXy are shorthand for multidimensional integrals over the displacements Xr,x,l and Xr,y,l defined at each oxygen site and time slice l.
The matrix Mσ is defined as Mσ = I + Bσ(L)Bσ(L − 1) ⋯ Bσ(1), where I is an N × N identity matrix and \({B}_{\sigma }(l)={\mathrm {{e}}}^{-\Delta \tau {H}_{{\rm{int}}}}{\mathrm {{e}}}^{-\Delta \tau {H}_{{\rm{el}}}}\). (Each Bσ(l) is an N × N matrix and independent of spin.) Note that Hint has off-diagonal terms in orbital space due to the nature of the SSH interaction, unlike the case of the Holstein model where it is a diagonal matrix. The lattice contribution to the action is defined as
The integrals over the displacements are evaluated using the Metropolis–Hastings algorithm.
Most observable quantities can be expressed in terms of the single particle equal time Green’s function Gσ(τ). For an electron propagating through field configurations {Xr,x,l} and {Xr,y,l}, the Green’s function at time τ = lΔτ is given by
where Aσ(l) = Bσ(l) ⋯ Bσ(1)Bσ(L) ⋯ Bσ(l + 1), \({\hat{T}}_{\tau }\) is the time ordering operator, and i, j are combined orbital and site indices. The determinant of Mσ is related to the Green’s function \({M}_{\sigma }=\det \ {[{G}_{\sigma }(l)]}^{-1}\) and is independent of l.
Efficient updates of the phonon fields
Equation (19) shows that the Green’s function Gσ(l + 1) can be obtained from Gσ(l) using the identity \({G}_{\sigma }(l+1)={B}_{\sigma }(l+1){G}_{\sigma }(l){B}_{\sigma }^{-1}(l+1)\). This observation forms the basis for an efficient single-site Sherman–Morris updating scheme42. The DQMC algorithm starts by computing the Green’s function on time slice l = 0 using Eq. (19). A series of individual updates are then proposed by sweeping through the sites (r, α = x, y) proposing updates \({X}_{{\bf{r}},\alpha ,l}\to {X}_{{\bf{r}},\alpha ,l}^{\prime}={X}_{{\bf{r}},\alpha ,l}+\Delta {X}_{{\bf{r}},\alpha ,l}\) while holding the other phonon fields {\({X}_{{{\bf{r}}}^{\prime}\ne {\bf{r}},{\alpha }^{\prime}\ne \alpha ,l}\)} fixed. These updates are accepted with probability \(p=\min (1,R)\), where
and \({M}_{\sigma }^{\prime}\) and \({M}_{\sigma }^{\prime}\) correspond to matrices computed with the new and old phonon field configurations, respectively. Note that the product \(\det [{M}_{\uparrow }]\det [{M}_{\downarrow }]\) is positive definite for the model considered here, and there is no Fermion sign problem.
After updating a field, the corresponding Bσ(l) matrix must also be updated as \({B}_{\sigma }(l)\to {B}_{\sigma }^{\prime}(l)={\mathrm {{e}}}^{-\Delta \tau {H}_{{\rm{int}}}^{\prime}}{\mathrm {{e}}}^{-\Delta \tau {H}_{{\rm{el}}}}={\mathrm {{e}}}^{-\Delta \tau ({H}_{{\rm{int}}}+V)}{\mathrm {{e}}}^{-\Delta \tau {H}_{{\rm{el}}}}\), where V contains the terms in Hint arising from the change in the phonon field. V is a symmetric matrix with only four non-zero elements of the form
To efficiently calculate the new \({B}_{\sigma }^{\prime}(l)\) matrix, we make the approximation \({B}_{\sigma }^{\prime}(l)\approx {\mathrm {{e}}}^{-\Delta \tau V}{B}_{\sigma }(l)\), which introduces an error on the order of the Trotter error and is valid when Δτ is small. The matrix e−ΔτV is then computed as e−ΔτV = Pe−ΔτDPT, where P is the orthogonal transformation that diagonalizes V, and D is a diagonal matrix with only two non-zero elements \({[D]}_{11}=-\sqrt{2}\alpha {t}_{sp}^{0}\Delta {X}_{{\bf{r}},l}\) and \({[D]}_{NN}=\sqrt{2}\alpha {t}_{sp}^{0}\Delta {X}_{{\bf{r}},l}\). The \({B}_{\sigma }^{\prime}(l)\) matrix can then be written as \({B}_{\sigma }^{\prime}(l)=P{\mathrm {{e}}}^{-\Delta \tau D}{P}^{\mathrm {{T}}}{B}_{\sigma }(l)=P(I+\Delta ){P}^{\mathrm {{T}}}{B}_{\sigma }(l),\) where [Δ]ij = 0 except \({[\Delta ]}_{11}={\mathrm {{e}}}^{\sqrt{2}\alpha {t}_{sp}^{0}\Delta {X}_{{\bf{r}},l}}-1\) and \({[\Delta ]}_{NN}={\mathrm {{e}}}^{-\sqrt{2}\alpha {t}_{sp}^{0}\Delta {X}_{{\bf{r}},l}}-1\).
Using these approximations, the Green’s function can be efficiently updated after accepting a change in the phonon field using
where Q = PT[I −Gσ(l)]. Due to the sparsity of matrix Δ, ΔQ has only two non-zero rows
where u and w are N × 2 and 2 × N matrices, respectively. Using the Woodbury matrix identity and Matrix determinant lemma, the updated Green’s function is given by
and the acceptance ratio is given by
where I2 is a 2 × 2 identity matrix. Evaluating these expressions involves O(N2) operations, as opposed to computing the updated Green’s function from scratch using Eq. (19), which has a computational cost of O(N3). This update scheme is efficient but it relies on the approximation \({B}_{\sigma }^{\prime}(l)\approx {\mathrm {{e}}}^{-\Delta \tau V}{B}_{\sigma }(l)\). Benchmarks for this approximation are provided in the supplementary materials (see Supplementary Note 1) (see Supplementary Material for more detailed results). In general, we find that the approximate fast update reproduces the exact solution with comparable error bars for Δτ values that are typical of most DQMC calculations.
Block updates
We periodically performed block updates of the phonon fields to reduce the autocorrelation time. In this case, we follow the procedure given in ref. 43.
Measurements
The conductivity weight is defined as \(\sigma (\beta /2)=\frac{{\beta }^{2}}{\pi }{\Lambda }_{xx}({\bf{q}}=0,\tau =\beta /2)\)45, where \({\Lambda }_{xx}({\bf{q}},\tau )={\sum }_{{\bf{r}}}\langle {\hat{j}}_{x}({\bf{r}},\tau ){\hat{j}}^\dagger_{x}(0,0)\rangle {e}^{{\rm{i}}{\bf{q}}\cdot {\bf{r}}}\) is the current–current correlation function and
is the current operator, with phase factors Qδ (as before) and Q±x,±y = −Q±y,±x = −Q±x,∓y = Q∓y,±x = 1. At low temperature, the conductivity weight is a proxy of the conductivity. More relevant results are shown in the Supplementary Note 2 (see Supplementary Material for more detailed results).
A measure of the superconducting and charge ordering tendencies can be obtained from the orbitally resolved charge \({\chi }_{{\gamma }^{\prime}\gamma }^{{\rm{C}}}({\bf{q}})\) and superconducting pair-field \({\chi }_{\gamma }^{{\rm{sc}}}\) susceptibilities, where γ is an orbital index. The charge susceptibility is defined as
where q is the momentum, τ is the imaginary time, \({\hat{n}}_{{\bf{q}},\gamma }={\sum }_{i,\sigma }{\mathrm {{e}}}^{{\rm{i}}{\bf{q}}\cdot {{\bf{r}}}_{i}}{\hat{n}}_{{{\bf{r}}}_{i},\gamma ,\sigma }\), and ri is the lattice vector. Similarly, the pair-field susceptibility in the s-wave channel is given by
where Δs = ∑rsr,↑sr,↓ and \({\Delta }_{{p}_{\delta }}={\sum }_{{\bf{r}}}{p}_{{\bf{r}},\delta ,\uparrow }{p}_{{\bf{r}},\delta ,\downarrow }\).
Data availability
The data that support the findings of this study will be made available upon reasonable requests to the corresponding author.
Code availability
The computer codes used in this study will be made available upon reasonable requests to the corresponding author.
References
Su, W. P., Schrieffer, J. R. & Heeger, A. J. Solitons in polyacetylene. Phys. Rev. Lett. 42, 1698–1701 (1979).
Li, S., Dong, X., Yi, D. & Xie, S. Theoretical investigation on magnetic field effect in organic devices with asymmetrical molecules. Org. Electron. 14, 2216–2222 (2013).
Clay, R. T. & Mazumdar, S. From charge- and spin-ordering to superconductivity in the organic charge-transfer solids. Phys. Rep. 788, 1–89 (2019).
Medarde, M. L. Structural, magnetic and electronic properties of RNiO3 perovskites (R = rare earth). J. Phys.: Condens. Matter 9, 1679–1707 (1996).
Shamblin, J. et al. Experimental evidence for bipolaron condensation as a mechanism for the metal–insulator transition in rare-earth nickelates. Nat. Commun. 9, 86 (2018).
Johnston, S., Mukherjee, A., Elfimov, I., Berciu, M. & Sawatzky, G. A. Charge disproportionation without charge transfer in the rare-earth-element nickelates as a possible mechanism for the metal–insulator transition. Phys. Rev. Lett. 112, 106404 (2014).
Khazraie, A., Foyevtsova, K., Elfimov, I. & Sawatzky, G. A. Oxygen holes and hybridization in the bismuthates. Phys. Rev. B 97, 075103 (2018).
Wen, C. H. P. et al. Unveiling the superconducting mechanism of Ba0.51K0.49BiO3. Phys. Rev. Lett. 121, 117002 (2018).
Marchand, D. J. J. et al. Sharp transition for single polarons in the one-dimensional Su–Schrieffer–Heeger model. Phys. Rev. Lett. 105, 266605 (2010).
Sous, J., Chakraborty, M., Krems, R. V. & Berciu, M. Light bipolarons stabilized by Peierls electron–phonon coupling. Phys. Rev. Lett. 121, 247001 (2018).
Möller, M. M., Sawatzky, G. A., Franz, M. & Berciu, M. Type-II Dirac semimetal stabilized by electron–phonon coupling. Nat. Commun. 8, 2267 (2017).
Weber, W. Electron–phonon interaction in the new superconductors La2−x(Ba,Sr)xCuO4. Phys. Rev. Lett. 58, 1371–1374 (1987).
Lanzara, A. et al. Evidence for ubiquitous strong electron–phonon coupling in high-temperature superconductors. Nature 412, 510-514 (2001).
Fradkin, E. & Hirsch, J. E. Phase diagram of one-dimensional electron–phonon systems. I. The Su–Schrieffer–Heeger model. Phys. Rev. B 27, 1680–1697 (1983).
Caron, L. G. & Moukouri, S. Density matrix renormalization group applied to the ground state of the XY spin-Peierls system. Phys. Rev. Lett. 76, 4050–4053 (1996).
Bakrim, H. & Bourbonnais, C. Quantum vs classical aspects of one dimensional electron–phonon systems revisited by the renormalization group method. Phys. Rev. B 76, 195115 (2007).
Weber, M., Assaad, F. F. & Hohenadler, M. Excitation spectra and correlation functions of quantum Su–Schrieffer–Heeger models. Phys. Rev. B 91, 245147 (2015).
Bakrim, H. & Bourbonnais, C. Nature of ground states in one-dimensional electron–phonon Hubbard models at half filling. Phys. Rev. B 91, 085114 (2015).
Weber, M., Toldin, F. P. & Hohenadler, M. Competing orders and unconventional criticality in the Su–Schrieffer–Heeger model. Phys. Rev. Res. 2, 023013 (2020).
Sengupta, P., Sandvik, A. W. & Campbell, D. K. Peierls transition in the presence of finite-frequency phonons in the one-dimensional extended Peierls–Hubbard model at half-filling. Phys. Rev. B 67, 245103 (2003).
Hohenadler, M. Interplay of site and bond electron–phonon coupling in one dimension. Phys. Rev. Lett. 117, 206404 (2016).
Tang, S. & Hirsch, J. E. Peierls instability in the two-dimensional half-filled Hubbard model. Phys. Rev. B 37, 9546–9558 (1988).
Mizokawa, T. et al. Origin of the band gap in the negative charge-transfer-energy compound NaCuO2. Phys. Rev. Lett. 67, 1638–1641 (1991).
Zaanen, J., Sawatzky, G. A. & Allen, J. W. Band gaps and electronic structure of transition-metal compounds. Phys. Rev. Lett. 55, 418–421 (1985).
Foyevtsova, K., Khazraie, A., Elfimov, I. & Sawatzky, G. A. Hybridization effects and bond disproportionation in the bismuth perovskites. Phys. Rev. B 91, 121114 (2015).
Plumb, N. C. et al. Momentum-resolved electronic structure of the high-Tc superconductor parent compound BaBiO3. Phys. Rev. Lett. 117, 037002 (2016).
Yin, Z. P., Kutepov, A. & Kotliar, G. Correlation-enhanced electron–phonon coupling: applications of GW and screened hybrid functional to bismuthates, chloronitrides, and other high-Tc superconductors. Phys. Rev. X 3, 021011 (2013).
Park, H., Millis, A. J. & Marianetti, C. A. Site-selective Mott transition in rare-earth-element nickelates. Phys. Rev. Lett. 109, 156402 (2012).
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).
Sleight, A. W. Bismuthates: BaBiO3 and related superconducting phases. Physica C 214, 152–165 (2015).
Cox, D. E. & Sleight, A. W. Mixed-valent Ba2Bi3+Bi5+O6: structure and properties vs temperature. Acta Crystallogr. Sect. B 35, 1–10 (1979).
Rice, T. M. & Sneddon, L. Real-space and \(\overrightarrow{{\rm{k}}}\)-space electron pairing in BaPb1−xBiO3. Phys. Rev. Lett. 47, 689–692 (1981).
Kim, G. et al. Suppression of three-dimensional charge density wave ordering via thickness control. Phys. Rev. Lett. 115, 226402 (2015).
Dee, P. M., Nakatsukasa, K., Wang, Y. & Johnston, S. Temperature-filling phase diagram of the two-dimensional Holstein model in the thermodynamic limit by self-consistent Migdal approximation. Phys. Rev. B 99, 024514 (2019).
Giraldo-Gallo, P. et al. Stripe-like nanoscale structural phase separation in superconducting BaPb1−xBixO3. Nat. Commun. 6, 8231 (2015).
Naamneh, M. et al. Cooling a polaronic liquid: phase mixture and pseudogap-like spectra in superconducting Ba1−xKxBiO3. https://arxiv.org/abs/1808.06135.
Bischofs, I. B., Allen, P. B., Kostur, V. N. & Bhargava, R. Topological doping of a three-dimensional Peierls system: predicted structure of doped BaBiO3. Phys. Rev. B 66, 174108 (2002).
Climent-Pascual, E., Ni, N., Jia, S., Huang, Q. & Cava, R. J. Polymorphism in BaPb1−xBixO3 at the superconducting composition. Phys. Rev. B 83, 174512 (2011).
Tajima, S. et al. Optical study of the metal–semiconductor transition in BaPb1−xBixO3. Phys. Rev. B 32, 6302–6311 (1985).
Nagata, Y., Mishiro, A., Uchida, T., Ohtsuka, M. & Samata, H. Normal-state transport properties of Ba1−xKxBiO3 crystals. J. Phys. Chem. Solids 60, 1933–1942 (1999).
Hohenadler, M. & Lang, T. Computational Many-Particle Physics, 357–366 (Springer, Berlin, 2008).
White, S. R. et al. Numerical study of the two-dimensional Hubbard model. Phys. Rev. B 40, 506–516 (1989).
Johnston, S. et al. Determinant quantum Monte Carlo study of the two-dimensional single-band Hubbard–Holstein model. Phys. Rev. B 87, 235133 (2013).
Blankenbecler, R., Scalapino, D. J. & Sugar, R. L. Monte Carlo calculations of coupled boson–fermion systems. I. Phys. Rev. D 24, 2278–2286 (1981).
Trivedi, N., Scalettar, R. T. & Randeria, M. Superconductor–insulator transition in a disordered electronic system. Phys. Rev. B 54, R3756–R3759 (1996).
Acknowledgements
We thank M. Berciu, N.C. Plumb, G.A. Sawatzky, and R.T. Scalettar for useful discussions. This work was supported by the Scientific Discovery through Advanced Computing (SciDAC) program funded by the U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research and Basic Energy Sciences, Division of Materials Sciences and Engineering. An award of computer time was provided by the INCITE program. This research also used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725.
Author information
Authors and Affiliations
Contributions
S.L. performed the calculations. S.L. and S.J. developed the code, discussed the results, and wrote manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Li, S., Johnston, S. Quantum Monte Carlo study of lattice polarons in the two-dimensional three-orbital Su–Schrieffer–Heeger model. npj Quantum Mater. 5, 40 (2020). https://doi.org/10.1038/s41535-020-0243-3
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41535-020-0243-3
This article is cited by
-
Local inversion-symmetry breaking in a bismuthate high-Tc superconductor
Nature Communications (2023)
-
A hybrid Monte Carlo study of bond-stretching electron–phonon interactions and charge order in BaBiO3
npj Computational Materials (2023)
-
Stripe correlations in the two-dimensional Hubbard-Holstein model
Communications Physics (2022)
-
Phase diagram of the two-dimensional Hubbard-Holstein model
Communications Physics (2020)