Quantum Monte Carlo study of lattice polarons in the two-dimensional multi-orbital Su-Schrieffer-Heeger model

TL;DR

This study uses quantum Monte Carlo simulations to explore lattice polarons in a two-dimensional multi-orbital model, revealing a bipolaron insulator, a metallic polaron-liquid, and an s-wave superconducting state, with implications for materials dominated by bond phonon interactions.

Contribution

First non-perturbative quantum Monte Carlo analysis of lattice polarons in a 2D multi-orbital model, uncovering phase transitions and superconductivity related to bond phonons.

Findings

Bipolaron insulating phase at half-filling with bond-disproportionation.

Metallic polaron-liquid state with fluctuating local distortions at high doping.

Emergence of s-wave superconductivity on oxygen sublattices at finite doping.

Abstract

We study a three-orbital Su-Schrieffer-Heeger model defined on a two-dimensional Lieb lattice and in the negative charge transfer regime using determinant quantum Monte Carlo. At half-filling (1 hole/unit cell), we observe a bipolaron insulating phase, where the ligand oxygen atoms collapse and expand about alternating cation atoms to produce a bond-disproportionated state. This phase is robust against moderate hole doping but is eventually suppressed at large hole concentrations, leading to a metallic polaron-liquid-like state with fluctuating patches of local distortions. Our results suggest that the polarons are highly disordered in the metallic state and freeze into a periodic array across the metal-to-insulator transition. We also find an $s$-wave superconducting state at finite doping that primarily appears on the oxygen sublattices. Our approach provides an efficient, non-perturbative way to treat bond phonons in higher dimensions and our results have implications for many materials where coupling to bond phonons is the dominant interaction.

Equations

$$H_0$$

$$H_{\mathrm{lat}}$$

$$H_{e-ph}$$

Full text

Introduction — Model Hamiltonians for electron-phonon (e-ph) interactions are broadly divided into two categories based on whether the coupling is diagonal or off-diagonal in orbital space. Diagonal e-ph interactions (e.g., Holstein Holstein (1959) or Fröhlich Fröhlich (1954) models) couple the atomic displacements to the charge density while off-diagonal e-ph couplings (e.g., the Su-Schrieffer-Heeger [SSH] modelSu et al. (1979)) modulate the carrier’s kinetic energy via the overlap integrals. To date, diagonal e-ph interactions have received the most attention Sangiovanni et al. (2006); Macridin et al. (2006); Bauer (2010); Bauer and Hewson (2010); Nowadnick et al. (2012); Murakami et al. (2013); Li et al. (2017); Berciu (2006); Noack et al. (1991); Costa et al. (2018); Bonča et al. (1999); Ku et al. (2002); Hague et al. (2006); Romero et al. (1999); Li et al. (2015); Chakraverty et al. (1998); Prokof’ev and Svistunov (1998); Weber and Hohenadler (2018); Dee et al. (2019), while studies of off-diagonal models have mainly been restricted to one-dimension [1D] Su et al. (1979); Li et al. (2013); Marchand et al. (2010); Sous et al. (2018); Sengupta et al. (2003); Hohenadler (2016); Tang and Hirsch (1988).

There is an urgent need to address off-diagonal e-ph interactions in higher dimensions, because such couplings are not only relevant to many materials – e.g. the organic charge-transfer solids Su et al. (1979); Li et al. (2013); Clay and Mazumdar (2018), the rare-earth nickelates Medarde (1996); Shamblin et al. (2018); Johnston et al. (2014), and high-$\mathrm{T}_{\mathrm{c}}$ superconductors like the cuprates Lanzara et al. (2001); Weber (1987) and bismuthates Khazraie et al. (2018) – but several recent studies in the few-particle limit have shown that the new physics can occur in such models. For example, strong off-diagonal interactions can produce highly mobile polarons with light effective masses Marchand et al. (2010), generate robust phonon-mediated pairing Sous et al. (2018), and even stabilize and control the location of a type-II Dirac point Möller et al. (2017). Off-diagonal models have also gained attention in relation to 1D topological insulators Qi and Zhang (2011) in the BDI class Schnyder et al. (2008). It is, therefore, imperative to study off-diagonal e-ph interactions in higher dimensions and at arbitrary fillings, as our intuition gained by studying diagonal models may not serve us well.

Another motivation for studying the SSH-type models in higher dimensions is to better understand its role in establishing the insulating and superconducting states of the high-Tc superconducting bismuthates $\mathrm{Ba}_{1-x}\mathrm{K}_{x}\mathrm{BiO}_{3}$ (BKBO). BKBO is in the so-called “negative charge transfer” regime Mizokawa et al. (1991); Zaanen et al. (1985); Foyevtsova et al. (2015); Plumb et al. (2016), 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 interaction Johnston et al. (2014); Khazraie et al. (2018), which may be further enhanced by correlations Yin et al. (2013), and is believed to drive a high-temperature metal-to-insulator (MIT) transition. Here, the insulating state has a 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 $A_{1g}$ symmetry Foyevtsova et al. (2015); Plumb et al. (2016); Khazraie et al. (2018); Johnston et al. (2014); Park et al. (2012); Bisogni et al. (2016). The relevant model describing this case is a multiorbital SSH model; however, knowledge about such models is limited due to a lack of suitable approaches for solving it.

With this motivation, we developed a determinant quantum Monte Carlo (DQMC) method for simulating SSH-type interactions, which is applied to study a 2D multi-orbital model for the first time. At half-filling (one hole/unit cell), we find that the system is a bipolaronic insulator with a bond-disproportionated structure, similar to what is observed in BKBO Sleight (2015) or the rare-earth nickelates Medarde (1996). Hole doping suppresses the insulating phase, giving way to state where the lattice distortions have short-range correlations suggestive of phase mixing and/or fluctuations. At high doping levels, we find evidence for a metallic phase where holes are strongly correlated with local structural distortions, forming a polaron-liquid phase. Finally, at low temperatures, we find $s$-wave superconducting tendencies that form primarily on the oxygen sublattice and evidence for a superconducting dome. 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.

Model and Methods — Keeping BKBO in mind, we adopt a three-orbital SSH model defined on a Lieb lattice whose orbital basis consists of a Bi $6s$ orbital and two O $2p$ orbitals, as shown in Fig. 1(a). We freeze the heavier Bi atoms into place and restrict lighter O atoms to move along the bond directions. The Hamiltonian is $H=H_{0}+H_{\mathrm{lat}}+H_{{\it e}-ph}$, where

[TABLE]

Here, $\langle\dots\rangle$ denotes a sum over nearest neighbor atoms, $\delta,\delta^{\prime}=\pm x$, $\pm y$ index the oxygen atoms surrounding each Bi, and the operators $s_{{\bf r},\sigma}^{\dagger}$ $\left(s_{{\bf r},\sigma}^{\phantom{\dagger}}\right)$ and $p_{{\bf r},\delta,\sigma}^{\dagger}$ $\left(p^{\phantom{\dagger}}_{{\bf r},\delta,\sigma}\right)$ create (annihilate) spin $\sigma$ holes on the Bi $6s$ and O $2p_{\delta}$ orbitals, respectively. The unit cells are indexed by ${\bf r}=n_{x}{\bf a}+n_{y}{\bf b}$, where ${\bf a}=(a,0)$, ${\bf b}=(0,a)$ are the primitive lattice vectors along $x$- and $y$-directions, respectively, and $a$ is the Bi-Bi bond length (and our unit of length). To simplify the notation, we have introduced shorthand notation $p_{{\bf r},-x,\sigma}=p_{{\bf r}-{\bf a},x,\sigma}$ and $p_{{\bf r},-y,\sigma}=p_{{\bf r}-{\bf b},y,\sigma}$. The operators $\hat{n}^{s}_{{\bf r},\sigma}=s_{{\bf r},\sigma}^{\dagger}s_{{\bf r},\sigma}^{\phantom{\dagger}}$ and $\hat{n}^{p_{\alpha}}_{{\bf r},\sigma}=p_{{\bf r},\alpha,\sigma}^{\dagger}p_{{\bf r},\alpha,\sigma}^{\phantom{\dagger}}$ are the number operators for $s$ and $p_{\alpha}$ ($\alpha=x,y$) orbitals, respectively; $\epsilon_{s}$ and $\epsilon_{p}$ are the site energies; $\mu$ is the chemical potential; $t_{sp}$ and $t_{pp}$ are the Bi-O and O-O hopping integrals in the undistorted crystal; and $\alpha$ is the e-ph coupling constant. The phase factors are $P_{x(y)}=-P_{-x(-y)}=1$, and $P_{\pm x,\pm y}=P_{\pm y,\pm x}=-P_{\pm x,\mp y}=-P_{\mp y,\pm x}=1$. The motion of the O atoms 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 the neighboring Bi atoms. Thus, bare phonon frequency is $\Omega=\sqrt{2K/M}$. Finally, the atomic displacements modulate the hopping integral as $t_{sp}(P_{\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}$.

We study the model on a square lattice with $N=4\times 4$ Bi atoms (48 orbitals in total) using DQMC. We stress that the model considered here is free of the Fermion sign problem. The details are provided in the supplementary materials Sup , along with expressions for the standard quantities measured in this work, and supplementary exact diagonalization calculations. The details of all non-standard quantities are provided in the main text. Throughout, we adopt $t_{sp}=2.08$, $t_{pp}=0.056$, $\epsilon_{s}=6.42$, and $\epsilon_{p}=2.42$ (in units of eV), which are obtained from DFT calculations of BaBiO3 Khazraie et al. (2018). We adopt a phonon energy $\Omega=\sqrt{2}t_{sp}$ and e-ph coupling strength $\alpha=4a^{-1}$, which gives average displacement’s squared of $\frac{1}{N}\sum_{\bf r}\langle\hat{X}^{2}_{{\bf r},x}\rangle=\frac{1}{N}\sum_{\bf r}\langle\hat{X}^{2}_{{\bf r},y}\rangle=0.0356a^{2}$ at half-filling, indicating that the oxygen atoms do not cross the bismuth atoms during the sampling. (Here, we are limited to large $\Omega$ by long autocorrelation times.)

Results — Figures 1(b)-1(d) plots 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 inverse temperature $\beta=10/t_{sp}$, which provides evidence for a bond disproportionated insulating state at $\langle\hat{n}\rangle=1$. 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. This behavior reflects the breathing distortion sketched in Fig. 1(a), and is consistent with the bond disproportionation observed in the insulating phase of the BKBO Cox and Sleight (1979); Rice and Sneddon (1981); Kim et al. (2015).

Figure 2(a) plots the dc conductivity $\sigma_{\mathrm{dc}}$ and orbital-resolved spectral weight $\beta G_{\gamma\gamma}({\bf r}=0,\tau=\beta/2)$, where $\gamma$ is the orbital index, for $\langle\hat{n}\rangle=1$ Trivedi et al. (1996); Sup . The conductivity (black dots) initially increases as the temperature is lowered until reaching a maximum at $\beta\approx 5/t_{sp}$ then it is suppressed. All three orbital spectral weights follow a similar trend, indicating a concomitant removal of spectral weight at the Fermi level. The insulating phase is characterized by a ${\bf q}=(\pi,\pi)$ charge order, as evidenced by the charge susceptibility $\chi^{\mathrm{C}}_{\gamma\gamma}({\bf q})$ plotted in Fig. 2(b) as a function of temperature. Below $1/\beta t_{ps}=0.2$, the charge correlations rapidly increase on the $s$ orbital, while there is little change on the $p$ orbitals. This observation implies that the charge density on the O sublattice is uniform, even in the bond disproportionated structure, while a charge modulation forms on the Bi sites in the insulating state. An examination of the real-space charge density, as shown in the inset of Fig. 2(b), confirms this. We stress, however, that the charge transfer between the Bi sites is on the order of 0.1 holes/Bi.

From this analysis, it is clear that the model has a bond-disproportionated structure and a small charge modulation on the Bi atoms in the insulating state. This result supports the bond disproportionation scenario proposed for the bismuthates Foyevtsova et al. (2015). We now examine how this state evolves upon hole doping. Here, our focus is on the possible formation of lattice polarons, where holes are bound to local breathing distortions of the oxygen sublattice. These objects can be studied by considering the polaron number operator $\hat{p}({\bf r})=\hat{x}_{{\bf r},L_{s}}(\hat{n}_{{\bf r},s}+\hat{n}_{{\bf r},L_{s}})$, where $\hat{n}_{{\bf r},L_{s}}=\sum_{\sigma}L^{\dagger}_{{\bf r},s,\sigma}L^{\phantom{\dagger}}_{{\bf r},s,\sigma}$ is the number operator for the $A_{1g}$ combination of the ligand oxygen orbitals $L_{{\bf r},s,\sigma}=\frac{1}{2}(p_{{\bf r},x,\sigma}+p_{{\bf r},y,\sigma}-p_{{\bf r},-x,\sigma}-p_{{\bf r},-y,\sigma})$  Sup and $\hat{x}_{{\bf r},L_{s}}=(\hat{X}_{{\bf r},x}+\hat{X}_{{\bf r},y}-\hat{X}_{{\bf r},-x}-\hat{X}_{{\bf r},-y})$. This operator measures the combined presense of holes in the $A_{1g}$ molecular orbital surrounding a Bi site and a local contraction of those same orbitals, and can be used to trace the evolution of polarons with doping.

With increasing hole concentrations, we observe a MIT at $\beta=10/t_{sp}$. Figure 3(a) plots $\sigma_{\mathrm{dc}}$ as a function of filling, where it increases upon hole doping until saturating at $\langle\hat{n}\rangle\approx 1.4$, indicating metallic behavior. At the same time, the number of polarons $\frac{1}{N}\sum_{\bf r}\langle\hat{p}({\bf r})\rangle$ decreases as additional holes are introduced but remains nonzero even at the largest dopings [Fig.3(b)], indicating that the free carriers have polaronic character. We also study polaron correlations in real space using the staggered polaron correlation function $\langle P({\bf r})\rangle=(-1)^{n_{x}+n_{y}}\langle\hat{p}({\bf r})\hat{p}(0)\rangle$, which is plotted in Figs.3(c)-(f) for selected hole concentrations. At half filling, $\langle P({\bf r})\rangle$ is positive for all ${\bf r}$, indicating that the polarons are frozen into a long-range two-sublattice order, consistent with the patterns inferred from Figs. 1 and 2. With increasing hole concentrations, $\langle P({\bf r})\rangle$ decreases at the larger distances, signalling an overall relaxation of the bond disproportionated on long length scales but the persistence of short-range correlations. Such behavior could reflect nanoscale phase separation Naamneh et al. (2018); however, studies on large clusters are likely 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 at most one or two lattice constants.

We also examined the doping evolution of 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 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 enormous 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 3(b) 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 amount of bipolarons are present in the system. The staggered bipolaron correlation function is plotted in Figs. 3(g)-(j). 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 suppressed at all length scales, while a finite number of bipolarons are present, as indicated in Fig. 3(b). These results can again be easily understood if the metallic phase is a polaron liquid.

Given the presence of lattice polarons in the metallic phase, we computed the $s$-wave orbital-resolved pair field susceptibility $\chi_{\gamma}^{\mathrm{sc}}$ Sup . Figure 4(a) plots $\chi_{\gamma}^{\mathrm{sc}}$ as a function of temperature at $\langle\hat{n}\rangle=1.59$, and compares it against the dominant charge correlations $\chi^{C}_{ss}(\pi,\pi/2)$ at this doping. All three susceptibilities increase with decreasing temperature, but $\chi^{\mathrm{sc}}_{p_{x}}=\chi^{\mathrm{sc}}_{p_{y}}\equiv\chi^{\mathrm{sc}}_{p}$ dominates below $1/\beta t_{sp}\approx 0.04$. This observation implies pairing appears predominantly in the oxygen atoms. Extrapolating $1/\chi_{p}^{\mathrm{sc}}$ to zero (inset), yields an estimate for the superconducting $\beta_{c}\approx 63.29/t_{sp}$). This value is artificially high, due to the large value of $\Omega$ used in our calculations. Nevertheless, our results provide evidence that the bipolaronic rich metallic phase has a superconducting ground state. Fig. 4(b) plots $\chi^{\mathrm{sc}}_{p}$ as a function of doping at $1/\beta t_{sp}=0.03$, where we find that pairing susceptibility is suppressed in proximity to the insulating phase, suggesting the presence of a superconducting dome induced by competition with the insulating phase.

Summary — 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 for the first time. We obtained several results consistent with the observed properties of bismuthate high-Tc superconductors. At half filling, we find a bond disproportionated state that can be viewed as a lattice of localized bipolarons. 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 BKBO Bischofs et al. (2002); Giraldo-Gallo et al. (2015); Climent-Pascual et al. (2011); Tajima et al. (1985); Nagata et al. (1999); Naamneh et al. (2018). 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.

Acknowledgments — 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.