Bipolaronic high-temperature superconductivity

,


I. INTRODUCTION
The quest for new pathways to high-temperature superconductivity from physically simple ideas has captivated researchers for decades.Conventional superconductivity is well understood within the standard framework of the Bardeen-Cooper-Schrieffer (BCS) theory, in which the exchange of phonons between electrons acts as a pairing glue to produce superconductivity with a transition-temperature T c that is a small fraction of the phonon frequency Ω ( = 1).In this framework, high T c can arise in systems with very large phonon frequencies, and indeed the pressure-stabilized hydrides with very large phonon frequencies exhibit superconductivity at remarkably high temperatures [1,2].A great deal of experimental and theoretical work has also focused on routes to high T c based on unconventional superconductivity with non-s-wave pairing symmetry arising from electronic correlations [3].
The pursuit of phonon-mediated high-temperature superconductivity from strong electron-phonon coupling at fixed phonon frequency Ω is, however, challenged by several important constraints.In a conventional superconducting material, such as Al, phonon-induced attraction between electrons induces a Cooper instability, giving rise to superconductivity with a T c /Ω that is small and vanishes as the strength of the electron-phonon coupling λ → 0. As λ increases, T c /Ω increases.However, the standard theoretical treatment of the large-λ problem is based on the Migdal-Eliashberg approximation which breaks down for λ larger than a critical value of order 1 [4][5][6][7][8][9][10] because of lattice reconstruction or formation of heavy bipolarons [4,7,8,10].In high-electron-density materials, increasing λ beyond λ ≈ 1 typically leads to lattice reconstruction into a new structure with a reduced λ, leading to a maximum value of T c at λ ≈ 1.At lower carrier density, the electron-lattice interaction may not induce lattice reconstruction; instead, bipolarons emerge and these bipolarons either form a charge-localized nonsuperconducting state [4][5][6] or undergo a superfluid transition at a T c determined by the inverse of their effective mass.However, the effective mass of strongly-bound bipolarons has generically been believed to be large and to increase rapidly with λ [11][12][13], implying generically low values of T c from bipolaronic superconductivity [11].
Here, we challenge the widely held view that bipo-arXiv:2203.07380v5[cond-mat.supr-con] 2 Feb 2023 laron formation is not favorable for high transitiontemperature superconductivity by providing a concrete, experimentally relevant model for phonon-mediated bipolaronic high-T c superconductivity with a T c that can be significantly higher than previously established upper bounds, see Fig. 1.Our work is based on the observation that in the dilute limit bipolarons are in effect interacting bosons, with a transition temperature that depends both on the mass and the density.At fixed mass, the transition temperature increases as the density increases-until either the transition temperature becomes of the order of the bipolaron binding energy or the density becomes large enough that the bipolarons significantly overlap, at which point the theory breaks down and, we suspect, the superconducting transition temperature saturates.Thus the maximum transition temperature is set by a combination of binding strength, inverse mass and inverse size, with small-size, light-mass, strongly bound bipolarons optimizing the maximum transition temperature.In the extreme strong-coupling regime, the size saturates to a value of the order of the lattice constant, while it appears that in all models the polaronic mass enhancement grows exponentially in λ.While these qualitative considerations are generic, the specifics and hence the maximum value of T c will depend on the specifics of the underlying microscopic model studied.General understanding of this physics has been obtained from studies of the Holstein model, in which lattice distortions couple to the electron density (potential energy), and the focus has been on the bipolaron mass, with the size receiving less attention.In these models the mass increases very rapidly as λ is increased.Fröhlich or extended-Holstein models in which the coupling to the electron density is longer ranged have also been studied; light masses have been found in circumstances involving an interplay of competing forces [14][15][16], but recent studies indicate that this light mass occurs in a very limited region of parameter space to be relevant in realistic systems [17].However, alternative forms of electron-phonon coupling are also important.In particular, any material with a unit cell consisting of more than a single atom will experience distortions of its atomic bonds that locally modulate the electronic hopping, as in the Peierls [18] or Su-Schrieffer-Heeger (SSH) models [19].Models of this type exhibit significant differences in polaron formation relative to Holstein-type models [20][21][22][23][24][25] including, importantly, lighter polarons [22][23][24] and, as shown in Ref. [25], lighter bipolarons.This previous analysis of bipolarons [25], however, did not address superconductivity explicitly, did not examine the bipolaron size systematically and was limited to the one-dimensional case and constrained to the regime t ∼ Ω (t is the amplitude of the electronic hopping) leaving the possibility that mass enhancement would become more severe in higher dimensions in the adiabatic limit t Ω relevant to most materials, thus destroying any hope of phononinduced high-T c behavior.
In this work, we use a recently developed numerically Here EBP(K) is the bipolaron dispersion, ΨBP is the bipolaron ground state wavefunction, and K is the bipolaron momentum.We contrast the behavior of the bipolaronic superconductivity in the bond-Peierls model against superconductivity of Holstein (H) bipolarons (filled circles, dotted orange line) for t/Ω = 2 and U = 8t and against the prediction of Migdal-Eliashberg (ME) theory of strong-coupling superconductivity in the bond-Peierls (bP; empty squares) and Holstein (H; empty circles) models.Here we use λ = α 2 /2Ωt for bond-Peierls bipolarons, λ = α 2 H /8Ωt for Holstein bipolarons, λ = 8α 2 /2πΩt in Migdal-Eliashberg theory of the bond-Peierls model, and λ = α 2 H /2πΩt in Migdal-Eliashberg theory of the Holstein model, where α and αH are the electron-phonon coupling constants of the bond-Peierls and Holstein models, respectively, see Appendix A for more details about the conventions used.Error bars represent statistical errors in QMC simulations corresponding to one standard deviation.This comparison illustrates that Tc of the bond-Peierls bipolaronic superconductor can exceed previously expected upper bounds for phonon-mediated superconductivity in a wide swath of parameter space.
exact, sign-problem-free quantum Monte Carlo (QMC) approach [26] to study bipolaronic superconductivity in a minimal model of Peierls electron-lattice coupling in two dimensions (2D).We compute bipolaron properties in various regimes of coupling and adiabaticity (including deep in the adiabatic limit) and combine these results with analytical understanding of the bipolaron superfluid phase diagram to determine T c at which a liquid of bipolarons undergoes a Berezinskii-Kosterlitz-Thouless transition into a superfluid.Our main results are as follows.First, any non-zero λ mediates an attractive interaction between electrons, giving rise to an s-wave bipolaronic superconductor with a T c /Ω that can become significantly larger than the upper bound predicted from Migdal-Eliashberg theory of strong-coupling superconductivity or from Holstein bipolaron superconductivity nearly everywhere in parameter space, see Fig. 1.Second, we find that Coulomb repulsion, modeled phenomenologically as a Hubbard U term, enhances the magnitude of T c of the s-wave bipolaronic superconductor up to a critical value of U/t, beyond which T c becomes suppressed (at intermediate to strong coupling, this behavior extends to large values of U/t), see Fig. 2. Finally, in our theory T c is largest when t/Ω ∼ 1 -2, implying that manipulating the stiffness of a crystal via structural or moiré engineering or fabricating crystals with light atoms offers a path towards realizing high-temperature superconductors.Some aspects of the physics we discuss here may be operative in known materials, possibly including the iron-based pnictide superconductors.

A. Model of bond-phonon-coupled electrons
We consider a minimal model of Peierls electronphonon coupling, the bond-Peierls model (also referred to as the bond-SSH model) on a 2D square lattice.In this model the electronic hopping between two sites is modulated by an oscillator associated with the bond connecting the two sites.This coupling can arise if transverse oscillations of out-of-plane atoms modulate the hopping of electrons between atoms in the plane.This is believed to occur in the superconducting iron pnictide materials, as discussed in Section IV and Appendix H.The Hamiltonian is Ĥ = Ĥe + Ĥph + Ve-ph . ( Here electrons with spin σ ∈ {↑, ↓} are governed by a Hubbard model Ĥe = −t i,j ,σ ĉ † i,σ ĉj,σ + h.c.+ U i ni,↑ ni,↓ with onsite repulsion U and ni,σ = ĉ † i,σ ĉi,σ at site i.The notation i, j refers to nearest-neighbor sites.We set the lattice constant a = 1 in what follows.We model distortions of the bonds connecting sites i and j as Einstein oscillators Ĥph = i,j b † i,j bi,j with frequency ( = 1) Ω = K/M (note Xi,j is a single oscillator associated with the bond, not a composite variable representing a difference of displacements of the atoms at the two ends of the bond).We take the interaction between electrons and phonons to be the simplest coupling term within the family of Peierls models describing the modulation of electron hopping by an oscillator Xi,j associated with the bond connecting sites i and j with coupling coefficient α = α/ √ 2M Ω.We henceforth set M = 1.The relevant parameters are a dimensionless coupling constant λ = (α 2 /K)/4t = α 2 /(2Ωt), the ratio of the typical polaronic energy scale to the free electron energy scale, and an adiabaticity parameter t/Ω.As discussed in Appendix B, this model does not have a sign problem in the singlet two-electron sector (unlike other, related models, e.g. the site-Peierls model [25]) and therefore can be studied to great accuracy using quantum Monte Carlo methods.

B. Method
Using a QMC approach based on a path-integral formulation of the electronic sector combined with either a real-space diagrammatic or a Fock-space path-integral representation of the phononic sector [26] allows us to study pairing and singlet bipolaron formation in the twoelectron sector of the model.The absence of a sign problem, specific to this particular microscopic formulation and sector of the model, enables numerically exact results with small statistical errors on large lattices even in the challenging regime of t Ω all the way to the heavy mass limit, see Appendix B for details.While the coupling of electronic hopping to phonons can take various forms, our approach here allows us to draw generic conclusions about high-T c bipolaronic superconductivity due to the modulation of electronic hopping by lattice distortions in the entire parameter space.1) at adiabaticity parameter t/Ω = 10/3 as a function of the electron-phonon coupling λ = α 2 /(2Ωt) for different onsite Hubbard U (in units of the electron hopping t): a. Bipolaron binding energy ∆BP in units of the electron hopping t, b.Bipolaron radial size probability density distribution (absolute value squared of the bipolaron wavefunction) PBP(R) for λ = 0.5, and c.Bipolaron effective mass m BP in units of the mass of two free electrons m0 = 2me = 1/t and its mean squared-radius R 2 BP .Error bars represent statistical errors in QMC simulations corresponding to one standard deviation.Error bars in PBP(R) correspond to statistical errors smaller than the symbol size and therefore are not shown.

III. SUPERFLUID OF BIPOLARONS
In 2D, bipolarons undergo a superfluid transition at a temperature T ≤ T c , where T c is determined by the bipolaron density and effective mass and depends only double-logarithmically weakly on the effective bipolaronbipolaron interactions [27][28][29].We can thus safely ignore bipolaron-bipolaron interactions, barring any competing instability, e.g.phase separation or Wigner crystallization.Based on prior work [30] we argue that these instabilities are unlikely.These considerations reduce our problem to that of the superfluidity of a gas of hardcore bipolarons in 2D, for which T c ≈ 1.84ρ BP /m BP [29], where ρ BP is the density of bipolarons and m BP is the bipolaron effective mass [31].This formula for T c remains valid in a broad density range so long as bipolarons do not overlap.The largest T c from this mechanism thus arises for a ρ BP corresponding to a liquid of bosons with an inter-particle separation that is on the order of the bipolaron radial size R BP , which, after lattice regularization, must be at least unity, i.e. for ρ BP = min{(1/(πR 2 BP ), 1/π)}.From this we obtain an estimate for the maximum T c of the Berezinskii-Kosterlitz-Thouless transition of the bipolaronic liquid that depends only on the bipolaron properties given by a. Bipolaronic high-T c superconductivity.BP obtained from QMC simulations (see Appendix B2) of Eq. (1) as a function of λ for different t/Ω at U = 8t.Our results shown in Fig. 1 prove that high-T c bipolaronic superconductivity is not only possible, but is robust even in the presence of a large Coulomb repulsion U = 8t.To appreciate this result, we contrast our computed T c /Ω against upper bounds based on superfluidity of Holstein bipolarons computed using the same methodology as for the bond-Peierls bipolarons or Migdal-Eliashberg theory of strong-coupling superconductivity out of a Fermi liquid applied to the bond-Peierls and Holstein models.We find that T c of the bond-Peierls bipolaronic superconductor generically exceeds these bounds in a large swath of parameter space.T c of the Holstein bipolaronic superconductor, for t/Ω = 2 at U = 0, rapidly drops with λ from a maximum of ∼ 0.05Ω at λ ∼ 0.25 (Appendix D).As U/t increases, the bipolaron mass decreases [12,32], but the binding energy drops very rapidly and the radius correspondingly increases [32], limiting T c to even smaller values, see Appendix D for more details.In contrast, the bond-Peierls bipolaron becomes strongly bound but remains relatively light as λ increases and can resist a large U/t, consistent with the predictions of Ref. [25], see Fig. 3.The comparison of our calculations to that of Migdal-Eliashberg theory (Appendix E) reveals that T c of the bond-Peierls bipolaronic superconductor also exceeds the maximum T c corresponding to strong-coupling superconductivity out of a Fermi liquid.This maximum T c is in qualitative agreement with a typical upper bound based on McMillan's approach to conventional superconducting materials [33].McMillan's approach is valid only in the regimes of validity of Migdal-Eliashberg theory, i.e. up to λ ≈ 1 [7,8,10,34].Thus, a typical upper bound from the McMillan approach can be estimated at λ = 1 to give a maximum T c /Ω ∼ 0.05 for a Coulomb pseudopotential µ = 0.12 (see Appendix F), in qualitative agreement with our results of Migdal-Eliashberg theory applied to the two models presented in Fig. 1.These comparisons illustrate that T c obtained in our model of bipolaronic superconductivity significantly and generically exceed previously held expectations.
Figure 1 demonstrates a remarkable phenomenology of the high-T c bipolaronic superconductivity.In particular, T c /Ω exhibits a non-monotonic, dome-like dependence on λ with a peak that shifts to smaller values with larger t/Ω.We can understand this behavior as follows.In the anti-adiabatic limit t/Ω 1, phonon exchange mediates an instantaneous pair-hopping interaction between electrons ĉ † i,↓ ĉj,↓ ĉj,↑ + h.c., which induces formation of light-mass, strongly bound bipolarons with an s-wave wavefunction (see Appendix G and Ref. [25]).By continuity, we expect this behavior to qualitatively persist and develop a frequency dependence (retardation) as t/Ω increases, accompanied by a proclivity for the bipolaron mass to increase as the number of phonons in the bipolaronic cloud grows.This competition between phonon-mediated kinetic energy-enhancing electron pairhopping interactions and a tendency to gain mass determines the properties of bipolaronic superconductivity.Our numerics reveal a parameterically large, physically relevant regime in which the bipolaron mass m BP exhibits weak to moderate enhancement whilst retaining a relatively small radial size R 2 BP and a large binding energy ∆ BP [25], see Fig. 3.This behavior is completely absent in the standard Holstein model in which bipolarons rapidly become heavy in a manner that depends exponentially on the electron-phonon coupling strength [12,13].The characteristic behavior of our model, which we believe to hold generically for Peierls-coupled systems, explains the increase in T c /Ω with λ up to an optimal λ op beyond which bipolarons enter a regime of exponential mass enhancement that becomes prominent for λ > λ op and larger t/Ω.Nonetheless, over a broad swath of parameter space we find the simulated T c /Ω curves to surpass all previously held expectations.
b. Coulomb repulsion enhancement of T c of bipolaronic superconductivity.Most importantly, T c /Ω in Fig. 1 exceeds these bounds even for large values of U/t, demonstrating the robustness of the bipolaronic mechanism against Coulomb repulsion.Figure 2 shows that T c /Ω exhibits a dome-like dependence also on U/t.The value of U that maximizes T c depends on λ (not shown).At intermediate coupling λ ∼ 0.5 -0.6, T c /Ω peaks around a value of U = 8t, which varies little with t/Ω.This unconventional enhancement of the value of T c of our s-wave bipolaronic superconductor up to such large values of U/t can be understood from Fig. 3, which reveals that the bipolaron's effective mass m BP and its squared-radius R 2 BP depend on U/t, but in opposite ways.For larger λ, R 2 BP depends very weakly on U/t, while m BP decreases with increasing U/t, explaining the growth in T c .This implies that in this limit the bipolaron size is already sufficiently large to enable the bound pair to avoid the Hubbard repulsion so that increasing U bears no effect on the symmetry of the pairing wavefunction.However, the bipolaron binding energy ∆ BP also decreases with increasing U/t, and this decrease ultimately limits T c , which cannot be greater than ∆ BP .In the Holstein model, the U term directly competes with an onsite phonon-mediated attraction [12,13] and while Fate of bipolaronic high-Tc superconductivity in the limit of large electronic densities.Schematic diagram illustrating the expected dependence of the superconducting transition temperature Tc in the model of Eqs. ( 1) -( 2) on the average electronic density ni .For ni min{2/πR 2 BP , 2/π}, bipolarons form a dilute superconductor.For larger densities bipolarons overlap and a new strongly correlated state may emerge.We envision a few possibilities: either pairing correlations continue to dominate and Tc saturates or grows (dashed line) or competing effects suppress Tc (dotted line).Ultimately, at or near half filling, barring superconductivity at weak λ in absence of nesting, Tc vanishes and, depending on the value of the electron-phonon coupling λ, a valence bond solid (VBS) or an antiferromagnet (AFM) emerges [35][36][37][38] (gray region).
this results in a decrease in m BP , it also causes a rapid increase in R 2 BP accompanied by a fast decrease in ∆ BP , and the latter two factors are more important and mean that ultimately U does not significantly enhance T c , see Appendix D. This analysis reveals that the Peierls bipolarons are generally much less sensitive to Coulomb repulsion than their Holstein counterparts [25].
c. Evolution with density of bipolaronic superconductivity.Our estimates of T c for bipolaronic superconductivity correspond to the largest electronic density for which bipolarons do not overlap.At higher densities, studies of the emerging strongly correlated state require new techniques capable of discerning between competing phases.We imagine at least two possibilities, depicted in Fig. 4. Either strong pairing correlations between the electrons [41] result in saturation of T c or even higher values of T c (dashed line), or competing effects resulting from the bipolaron overlap combined with the Hubbard repulsion and lattice reconstruction suppress T c (dotted line).In either scenario, the system eventually becomes non-superconducting at sufficiently large densities near or at half filling (unless unnested) and, depending on the value of λ, antiferromagnetism or valence-bond-solid charge order develops [35][36][37][38].

IV. OUTLOOK
Obtaining high-transition-temperature superconductivity from Bose condensation of bipolarons has for many Here there are two hopping processes: one resulting from direct overlap between dxy orbitals on neighboring iron atoms with amplitude t < 0 (dashed double-headed arrow) and another involving a second-order process in which the dxy orbitals on iron atoms overlap with the px orbital of the pnictogen atom (solid double-headed arrows), resulting in a hopping with amplitude t > 0. Destructive interference between these two processes results in a reduced net hopping along this pathway, e.g. in FeSe, see Appendix H for more details.
years been thought to be very difficult.Recent work [25] revealed conditions under which light bipolarons can arise, opening a door for phonon-mediated high-T c superconductivity.This paper uses an exact quantum Monte Carlo treatment of a precisely defined 2D square lattice model to demonstrate that bond electron-phonon coupling in fact gives rise to bipolaronic superconductivity with a transition temperature that is much higher than that obtained from the more extensively studied Holstein (density-coupled) bipolarons or from Migdal-Eliashberg theory of superconductivity out of a Fermi liquid.Perhaps more importantly, we find that T c of the bondcoupled bipolaronic superconductor is enhanced by local Coulomb repulsion.The key ingredient in this bipolaronic mechanism that gives rise to high T c is the combination of light mass and relatively small size of bipolarons.While the calculations reported here employ the bond-Peierls coupling in which the amplitude for an electron to hop between two sites is modulated by a phonon defined on the bond connecting the two sites, bipolaronic high-T c superconductivity may also arise in the site-coupled Peierls (SSH) model [25] where the hopping is modulated by the relative distance between the atoms at the two ends of the bond, but an analysis including anharmonic couplings is required.These remarkable properties call for consideration of the physics of Peierls electron-phonon coupling in quantum materials, and motivate further theoretical study of other physical situations that may give rise to small-size, light-mass bipolarons needed in order to support a state with high-T c superconductivity.Models of potential interest include ones in which a phonon on a bond couples to two sites or a plaquette with either Holstein (e.g.Ref. [42]) or Peierls coupling or both, or a situation in which N phonons, each couples to N electronic sites, see e.g.Ref. [43].
The bond-Peierls coupling studied here arises generically in materials where the orbitals of out-of-plane atoms mix with the bonding orbitals of in-plane atoms so that transverse fluctuations of the distance of the out-of-plane atoms give rise to modulation of the hopping.This physics may be operative in several families of materials, including the 90 • -bonded [44,45] and corner-sharing [46] perovskites.One particularly intriguing example occurs in the iron pnictides where electron transfer between two adjacent Fe ions can occur via a state on an intermediate pnictogen ion so that modulation of the pnictogen height strongly modulates particular hopping pathways [39,40], and where for certain geometries different pathways destructively interfere, producing a bond-Peierls coupling with a large coupling constant.Figure 5 demonstrates this scenario.See Appendix H for more details and discussion.Although the model studied here is not directly applicable to the pnictides, which are complex materials involving multi-orbital "Hund's metal" physics, it is intriguing to note that an "extended s-wave" state with some similarities to our bipolaron state has been proposed for its superconductivity.
From a materials science perspective, it is worth highlighting that the high-T c bipolaronic superconductivity becomes notably pronounced in the most "quantal" regime of t/Ω ∼ 1 (Fig. 1) in which phonons and electrons are competitive energetically.Here, T c ∼ 0.2Ω, which could easily give rise to a T c ∼ 70K for a typical value of phonon frequency ∼ 0.03 eV if and only if the unusual limit of t ∼ Ω can be achieved.For the material FeSe we expect the ratio of the nearest-neighbor hopping amplitude in the x -y plane to the relevant outof-plane phonon frequency to be t/Ω ∼ 2 -3 [39,47] (see Appendix H), which is close to this ideal regime.In addition, structural engineering of a crystal's electronic stiffness [48], for example, by strain or moiré engineering [49] as found in twisted bilayer graphene [50], can produce a reduction in t without significantly changing Ω, perhaps realizing t ∼ Ω.Alternatively, functional superatomic crystal engineering [51] may enable synthesis of crystals with light atoms on the bonds of the lattice, providing a path to high-Ω phonons, achieving t Ω and an even higher value of T c .
Important future research directions raised by our work include full theoretical characterization of the phenomenology of bipolaronic high-T c superconductivity, extension of our results to the richer models describing the Hund's and multiorbital physics of materials such as the pnictides, and search for other bond-Peierls coupled compounds, thus opening a door to a new route to design principles of novel high-temperature superconductors.In this appendix we detail the conventions used to define the dimensionless electron-phonon coupling constant λ for the different models considered in the main text.
The electron-phonon coupling terms in the bond-Peierls (BP) and Holstein (H) models are where α = α √ 2M Ω.Our notation here differs slightly from that of the main text to clearly distinguish the BP and H models and to make explicit that there are two BP phonons in each unit cell of the square lattice, one associated with a bond in the x direction and one with a bond in the y direction.
In the polaronic limit, λ is normally defined as the ratio of the typical polaron energy scale to the free electron energy scale, which for the BP model is λ = α 2 BP /(2Ωt) and for the H model is conventionally λ = α 2 H /(4Ωt) [13].However, in order to contrast T c /Ω for the bipolaronic superfluid in the two models on the same scale, we rescale λ → λ/2 in the case of the H model so that it becomes λ = α 2 H /(8Ωt).In Migdal-Eliashberg (ME) theory of strong-coupling superconductivity out of a Fermi liquid, λ, which we denote as λ ME , is conventionally defined via the Fermi surface mass enhancement: Thus, λ ME acquires a dependence on the density of states (see details in Appendix E) and in the low electron density limit becomes λ ME = 8α We employ a recently developed quantum Monte Carlo (QMC) approach based on a path-integral representation of the electronic sector combined with either a path-integral or a diagrammatic representation of the phononic sector.The method and its details can be found in Ref. [26].Here we provide a brief overview of the approach, and details of the numerical simulations.

Methodology
We use Monte Carlo (MC) to stochastically sample the imaginary-time propagator G ba (τ ) ≡ b| e −τ Ĥ |a , where τ is imaginary time, and |a and |b are any two-electron states in the singlet sector on a two-dimensional (2D) square lattice.We can evaluate the ground-state (GS) where O ba (τ ) = b| e −(τ /2) Ĥ Ôe −(τ /2) Ĥ |a and W ab are MC weights.
For the model in Eqs. ( 1), ( 2), the basis B corresponds to site Fock states for the electrons and bond Fock states for the phonons.Using the interaction representation for the evolution operator in imaginary time: e −τ Ĥ = e −τ ĥ σ(τ ), and expanding σ(τ ) one finds [52,53]: where E βα = E β − E α .In this representation, there are three types of the elementary processes: 1. bare electron hopping, 2. electron hopping assisted by phonon creation, and 3. electron hopping assisted by phonon annihilation.Phonons can also be treated diagrammatically within the same expansion whilst maintaining a pathintegral representation only for the electronic sector.The MC scheme employed is based on the statistical interpretation of the right-hand side of Eq. (B2) as an average over an ensemble of graphs, in which each graph represents a string in space-time coordinates characterized by the number and types of kinks V γi+1γi .For the model in Eqs. ( 1), (2), the choice of basis B ensures that −V βα are non-negative, real numbers, and thus graphs are sampled according to non-negative weights given by the values of the corresponding integrands in Eq. (B2), rendering this a sign-problem-free MC approach.The sign problem is present for the site-Peierls model because the coupling involves the difference between phonon displacement operators on different sites, which means that the sign of a subset of the −V βα factors in Eq. (B2) will be negative.In contrast, in the bond-Peierls model, the sign of the coupling can be always gauged away and therefore the sign problem is absent.We use the following updates in the MC scheme: we stochastically 1. add and remove the last bare-hopping kink using a pair of complementary updates [52,53], 2. switch between the three types of the hopping terms, and 3. sample the length of the last time interval separating the last kink from the state b|.This scheme yields states b| that admit any allowed configuration of excited phonon modes, and, as a result, ensures ergodicity.

Computation of bipolaron properties
First, we note that the imaginary-time dependence of G ba (τ ) contains direct information about the groundstate energy E GS within a given momentum symmetry sector of the Hilbert space: and thus we can use this relation to extract E GS in the τ → ∞ limit for a given momentum symmetry sector, see Fig. 6.Furthermore, in the center-of-mass coordinate representation of the two-electron states |a and |b , the propagator G ba (τ, R) which depends on the relative center-of-mass displacement R = R b − R a , assumes a universal form given by the propagator of a free particle whose effective mass is given by the bipolaron mass m BP (see, for example, Ref. [54]), which in 2D takes the form: where A ba is a non-universal coefficient which depends on the choice of states |a and |b .We can thus extract m BP , by averaging over the states |a and |b , from the meansquare fluctuations of the center-of-mass displacement in FIG. 6. Computation of the bipolaron mass using a QMC approach based on a path-integral representation of the electrons combined with a diagrammatic treatment of the phonons.a. Green function G ba (τ ) as a function of imaginary time τ in different bipolaron momentum (K) sectors (solid lines) and fits to the long-τ asymptotic behavior of Eq. (B3) (dashed lines).b.The bipolaron dispersion EBP(K) constructed from EGS in the different K sectors obtained in a from the fits of G ba (τ ) in the large-τ limit to the asymptotic form in Eq. (B3).We compute the bipolaron mass by fitting the dispersion.For this data set, we find the following fitting function: EBP(K) = −9.38883+ 0.00188511K + 0.0791534K 2 , which yields m BP /m0 = 6. the large-τ limit: R 2 (τ ) ultimately saturates to a straight line at sufficiently large τ leading to an accurate estimate of the effective mass, see Fig. 7.
We can thus compute the bipolaron mass either by constructing the bipolaron dispersion E BP (K) as a function of the bipolaron momentum K from the asymptotic behavior of Eq. (B3) as shown in Fig. 6 or directly from the asymptotic behavior of R 2 (τ ) in Eq. (B5) as shown in Fig. 7.We use a QMC approach based on a pathintegral representation of the electrons combined with a diagrammatic treatment of the phonons to simulate the asymptotic behavior of Eq. (B3), and a QMC approach based on a path-integral representation of both the electrons and the phonons to simulate asymptotic behavior of Eq. (B5).Estimates of the mass obtained using these two approaches agree within the error bars, see Fig. 8. Importantly, the dispersion we obtain in our simulations retains a parabolic form for all energies lower than T c even when K becomes on the order of the inverse size of the bipolaron, justifying the use of a continuous space description.
To compute the bipolaron mean squared-radius R 2

BP
we use Monte Carlo estimators (Eq.(B1)) to evaluate the probability distribution P (R 12 ) of finding the two electrons at a distance R 12 from their center-of-mass po- sition within a bound bipolaron, from which we compute ). (Note that R 2 BP and R 2 (τ ) of the two-electron state are unrelated; the former is the square of the relative distance measured from the center-of-mass coordinate, i.e. half the distance between the two electrons, while the latter is the square of the center-of-mass displacement-during imaginary-time evolution-averaged over worldline configurations.) We use the two QMC approaches to simulate the behavior of two electrons in the singlet sector of the model (Eqs.( 1), ( 2)) on a 2D square lattice with linear size L = 128 sites.The accuracy of the QMC results is controlled by the maximum value of τ , τ max which determines the accuracy of the projection of the propagation onto the ground state, with numerically exact results recovered in the limit τ max → ∞.All results presented in this work are converged with respect to τ max .Our calculations of the mass and size of bipolarons in the ground state allow us to estimate T c reliably at temperatures below the binding gap, i.e., so long as T c ≤ ∆ BP .Error bars shown in the figures represent statistical errors in QMC simulations corresponding to one standard deviation, and, when applicable, account for errors in fitting.

Appendix C: Supplementary results on superfluidity of bond-Peierls bipolarons
In the main text, in Fig. 1, we show T c of a superfluid of bond-Peierls bipolarons at U/t = 8 at various values of the adiabaticity parameter t/Ω.The behavior of the mass and size of the bipolaron determines the value of T c as can be seen from Eq. (3). Figure 3 provides informa-tion about the bipolaron mass and size at t/Ω = 10/3.Here we present in Figs. 9, 10 additional results on the properties of the bipolaron for various values of t/Ω.
Figure 9 shows that the overall trend of ∆ BP , m * BP and R 2 BP with λ shifts to larger values of λ as t/Ω decreases, but attains a qualitatively similar profile.Similarly, the spatial structure of the bipolaron P BP (R) at the optimal λ that maximizes T c appears to not depend strongly on the value of t/Ω, see Fig. 10.From Fig. 9 we see that upon increasing λ, both m * BP increases and R 2 BP decreases in a fashion in which there exists an optimal λ (e.g. for t/Ω = 10/3, the optimal λ is ∼ 0.5) for which m * BP is not too large yet R 2 BP is relatively small, leading to a maximum in the T c curve, see Fig. 1.

Appendix D: Supplementary results on superfluidity of of Holstein bipolarons
In the main text, in Fig. 1, we show T c of a superfluid of Holstein bipolarons at U/t = 8.Here we present additional results on the behavior of T c of a superfluid of Holstein bipolarons computed from Eq. ( 3) using the same methodology as for the bond-Peierls bipolarons, and provide details on the behavior of the mass and radial size of the bipolaron.
In Fig. 11, we show the behavior of T c /Ω for Holstein bipolarons at t/Ω = 2 as a function of λ for U/t = 0 and 8 (Fig. 11a), and as a function of U/t for λ = 0.5 (Fig. 11b).T c /Ω never exceeds ∼ 0.05.Increasing λ past an optimal but small value leads to rapid bipolaron mass enhancement (see below) and suppression of T c .Increasing U/t leads to a decrease in the mass, but the binding energy drops very rapidly and the radius correspondingly increases, limiting T c to even smaller values.This behavior can be seen clearly in Fig. 12 which shows that T c is much smaller in the Holstein model than in the bond-Peierls model because the Holstein bipolarons are always very heavy and small when strongly bound.They can become lighter only when weakly bound which also results in large R 2 BP , once again producing a small T c .maximum physically attainable electron-phonon coupling strength in realistic models in which the electron-phonon interaction is computed from microscopics in a theory of electrons and ions coupled by the physical Coulomb interactions.Here we focus on the properties of the model systems discussed in the main text, with the Hubbard U set to zero.The results for the H model presented here are consistent with recent work of Esterlis and collabo-rators [34].
Contact with conventional ME theory is more clearly made in momentum (k) space, where the H (Eq. (A2)) and BP (Eq.(A1)) couplings are We have rewritten the BP coupling in terms of the phonon operators: 1. Phonon stiffness and stability limits The ME theory of these models proceeds by first computing a physical phonon stiffness given by the difference of the bare phonon stiffness denoted by K and α2 multiplied by the zero frequency limit of an appropriate electron correlation function.(In the adiabatic limit Ω 0 E F , the frequency dependence of the correlator is negligible; in other words, the phonon mass is not renormalized, and the correlator may be computed us- ing the bare electron Green functions.Here Ω 0 is the bare phonon frequency, which is the same as Ω in the main text.) In the H model the relevant correlator is the electron density-density correlation function and we have Since a positive phonon stiffness is required for stability of the oscillator, the maximum coupling is bounded (in the ME approximation) by α2 < max q K χρρ(q,0) .Standard density functional theory (DFT) computations of phonon frequencies include (within the approximations of DFT and within the adiabatic limit) renormalization of the phonon frequency, i.e. a material crystal structure is by construction stable (within the DFT approximation).Esterlis and collaborators [10,34] studied the stability issue using numerical methods which allowed them to go beyond the ME approximation; the main focus of their work was electron densities near half filling where the susceptibility has a substantial density dependence and density waves provide a competing ground state.A qualitative result of their work is that while the ME approximation is not accurate near the stability limit, the estimate α2 < max q K χρρ(q,0) still provides a reasonable bound.
Our interest here is specifically in the ME approximation and in the low density limit.For the electronic model dispersion used here ε k = −2t (cos k x + cos k y ), at low electronic densities (say, n ≡ ni 0.25 electron per FIG. 13.Static density-density correlation function χρρ(q) ≡ χρρ(q, 0) of non-interacting electrons on a 2D square lattice tight-binding model with dispersion ε k = −2t(cos kx + cos ky) at carrier concentration n = 0.25 carrier per site as function of momentum (q) along direction (1, 0) (solid line) and (1, 1) (dashed line), with n → 0 value of density of states 1/(2πt) (dotted line).
cal phonon frequency Ω phys is, to a good approximation, constant and given by where the bare ME coupling constant is and the second, approximate, equality becomes exact in the low density limit.Thus, the stability limit of the Holstein model in the ME approximation (see Eq. (E5)) at low densities is λ 0 H;ME = 1/2.Note that the λ defined for Holstein bipolarons differs from this definition, see Table I.
In the BP model there are two phonon modes per unit cell of the square lattice and the couplings are momentum-dependent.In the ± basis of Eq. (E2) Eq. (E5) becomes χ ++ (q, 0) χ +− (q, 0) χ +− (q, 0) χ −− (q, 0) , (E8) where χ ±± is the Λ ± Λ ± correlator and χ +− is the Λ + Λ − correlator (see Eq. (E4)).Figure 14 shows numerical calculation of these correlators.The cross correlator χ +− (q, 0) is in general very small; we neglect it here.In the very low density limit (n 0.05/site), χ ++ (q, 0) = 8χ ρ,ρ = 4 πt and χ −− (q, 0) ≈ 0. For moderately low densities (e.g.n = 0.25/site) χ −− (q, 0) remains small relative to χ ++ (q, 0) but χ ++ (q, 0) acquires nonnegligible momentum dependence with the largest value being at q = 2k F .For larger n, χ −− can become larger than χ ++ and in fact sets the limit of stability.The momentum dependence of the electron-phonon coupling leads to some ambiguity in the definition of λ 0 BP;ME .Here we adopt a definition appropriate to the very low density limit, writing where χ ± (q, 0) ≡ χ ±± (q, 0).The λ defined for BP bipolarons in the main text is a factor of π/8 smaller than λ 0 BP;ME , see also Table I.The factor πt 4 χ + (q, 0) becomes 1, independent of q in the very low density limit, and the − mode decouples, so the physics becomes identical to the H model apart from the relation between α and λ.At larger n ∼ 0.25 the − mode still decouples but the factor is less than 1 at all q and is q dependent, being largest at q = 2k F .At larger densities the − mode does not decouple.

Electron self energy
The ME approximation to the electron self energy Σ (making use of the adiabatic limit, which allows us to average the self energy over all wave vectors on the Fermi surface) gives where the phonon propagator for mode a is and the Fermi surface-projected electron Green function is  In the H model the phonon propagator and coupling coefficient are independent of momentum and its angle on the Fermi surface, and it is convenient to define the coupling constant as λ H;ME = λ 0 H;ME Ω0 Ω phys using the physical, renormalized phonon frequency.After linearizing in the anomalous self energy we have, explicitly These are the forms in which the Migdal-Eliashberg equations for the electron self energy are traditionally presented and solved.If the temperature is low compared to the phonon frequency one finds from Eq. (E15) that the low frequency mass enhancement is 1 + λ H;ME .The critical transition temperature may easily be determined following Bergmann and Rainer [55] by recasting Eq. (E16) as an eigenvalue equation for the vector W (ω n )/|ω n Z(ω n )| and defining the transition temperature T c as the temperature at which the leading eigen-value vanishes.The result is a T c which, as a fraction of the renormalized phonon frequency, evolves from the small λ Bardeen-Cooper-Schrieffer (BCS) form of e −1/λ to the Allen-Dynes form of √ λ as λ is increased from small to large values.Figure 15 shows the ratio of T c to Ω phys as a function of λ H;ME calculated from Eq. (E16).
In the present context it is of greater relevance to present the results as the ratio of T c to the bare oscillator frequency Ω 0 , as a function of the bare coupling λ 0 H;ME .These are obtained by a simple scaling of the results in Fig. 15 and are presented in the main text in Fig. 1.
In the BP model the presence of two phonon modes in the unit cell and the momentum dependence of their coupling coefficients and phonon frequencies make the analysis more involved.In order to evaluate Eq. (E11) for the BP model we need the coupling function Λ(k x , k y ) (Eq. (E4)) at the momentum (k + k )/2.For simplicity, we make the circular Fermi surface approximation and define ψ = (θ + θ )/2 and φ = θ − θ , then Eq. (E14) can be written as Figure 16a shows the dependence of the coupling constant for the + mode as a function of phonon momentum (parametrized by φ) calculated within the circular Fermi surface approximation for densities n = 0.25, 0.1 and 0.05.Also shown is the coupling for the − mode at n = 0.25 (for the lower densities the coupling function is essentially indistinguishable from zero).We see that at all of the relevant densities the − mode decouples.We solve the linearized gap equations for the BP model in the explicit form with We calculate the transition temperature by proceeding as in the H case. Note however that we have formulated the equations form the outset in terms of the bare λ, λ 0 , not the conventional ME λ defined in terms of the renormalized phonon frequency Ω phys .Results for n = 0.25 are shown in the right panel of Fig. 16b and in the main text in Fig. 1.

Appendix F: McMillan's phenomenological approach to phonon-mediated strong-coupling superconductivity
McMillan's approach to strong-coupling superconductivity is based on a phenomenological treatment of experimental data on conventional superconducting materials within the framework of Migdal-Eliashberg theory.This approach makes use of a coupling constant λ estimated directly from experiment by considering an electron-phonon interaction averaged over the Fermi surface, and supplements this treatment by empirical parameters used to mimic the effect of Coulomb interaction in order to better fit experimental data.The Migdal-Eliashberg theory itself is valid in the adiabatic limit t Ω and for moderate values of λ, because at λ 1, apart from structural instability [7], the Fermi liquid becomes a metastable state, higher in energy than a state formed of bipolarons [4-6, 8-10, 34, 56].Within its domain of applicability, McMillan's formula finds that strong electron-phonon coupling induces a superconducting instability out of a Fermi liquid at [33] T where µ = 0.12 is the value of the Coulomb pseudopotential found in many materials [33].itative agreement with the results of Migdal-Eliashberg theory applied to the two models presented in Fig. 1 (see also Appendix E).These comparisons illustrate that T c obtained in our model of bipolaronic superconductivity generically exceeds previously held expectations.
Appendix G: Effective electronic Hamiltonian in the antiadiabatic limit To unravel the pairing mechanism responsible for the formation of bipolarons, we derive an effective electronic Hamiltonian in the antiadiabatic limit t Ω by projecting out high-energy subspaces with one or more phonons [25,57].This procedure is valid at strong coupling λ 1 within the antiadiabatic regime if t α Ω such that α 2 Ωt.To second order, we find an effective two-electron Hamiltonian given by: Ĥeff.= ĥe + Ûe-e + Ve-e , Ω is the polaron formation energy, Ũ = U − 8α 2 Ω−U + 8α 2 Ω is the amplitude of a phonon-renormalized effective onsite density-density interaction, T = 2α 2 Ω−U is the amplitude of a phonon-mediated onsite electron pair hopping interaction [25], Ṽ = 2α 2 Ω − α 2 Ω+U is the amplitude of a phonon-mediated nearest-neighbor density-density interaction [58], and J = 4α 2 Ω+U is the amplitude of a phononmediated SU(2)-preserving nearest-neighbor spin-spin interaction [25,58].
Our numerical results away from the antiadiabatic limit indicate that the salient features embodied by Ĥeff., specifically the phonon-mediated kinetic energy enhancing pair-hopping interaction [25], continue to hold qualitatively as t/Ω increases, as evidenced by the light bipolaron masses.In contrast, in the adiabatic t/Ω 1, phonons should behave classically and have no dynamics, thus a pair of electrons experience a retarded phononmediated attraction and form a singlet bipolaron localized on a lattice bond in order to minimize the total energy.Our numerical results indicate that away from these asymptotic limits a competition between the phononmediated kinetic energy enhancing electron pair-hopping interaction and the tendency to localize electron pairs is at play and determines the fate of bipolarons, which, nonetheless, appear to be light in a large region of parameter space.
Finally, we see that, at least in the antiadiabatic limit, the Ve-e contains a phonon-mediated repulsive part [58], which disfavors pairing in all but the rotationally symmetric s-channel, and a large U/t, e.g.U = 8t discussed in the main text, further enhances this tendency.ment of the out-of-plane atoms give rise to modulation of the barrier for electron tunneling across bonds, precisely as embodied by the model in Eqs.(1), (2).Interestingly, Refs.[39,40] show that this physics is operative in the iron pnictides wherein modulation of the pnictogen height can strongly modulate particular hopping pathways.Figure 5 demonstrates this scenario.Here a pnictogen atom sits at the apex of an octahedron with four iron atoms residing in the x-y plane in the middle of the octahedron.The transverse motion of the pnictogen atom out of the x-y plane in the z direction causes fluctuations in the barrier for electronic tunneling between the iron atoms within the x-y plane.[39,40] (Fig. 5a).A dominant pathway for direct electronic hopping between iron atoms in this class of materials involves overlaps between lopes of d xy orbitals of opposite sign on neighboring iron atoms in the x-y plane, resulting in a negative hopping t < 0 [39].An indirect electronic pathway, resulting from a second order, superexchange-like process, involves the overlap of the apex atom's p x orbital with each lobe of the two d xy orbitals, resulting in a net positive hopping t > 0 [39].These two pathways (Fig. 5b) nearly cancel in FeSe because the magnitudes of t and t are nearly equal, and more generally the ratio of t to t vary in other pnictide materials resulting in an overall reduction in the magnitude of the net electronic hopping between the iron atoms [39].This interference effect combined with the large modulation of the tunneling barrier by the displacement of the pnictogen atom along this particular hopping pathway means that the value of the dimensionless electron-bond-phonon coupling strength λ relevant to this mode in this class of materials can be large, see Refs.[39,47].As an example, Ref. [60] suggests a value of λ ∼ 0.5 in one member of this family of materials, but more work is needed to accurately determine the strength of electron-phonon coupling in specific compounds.The net overall hopping in the x-y plane is roughly ∼ 50 meV [39] and the transverse phonon frequency is estimated to be ∼ 5.3 THz ≈ 22 meV in FeSe, and ∼ 17 meV in FeTe [47], implying a ratio of the relevant phonon frequency to the magnitude of the relevant (net) electron hopping of ∼ 2 -3 in these materials.We also note that while the pairing symmetry in these materials has not yet been fully determined, an "extended s-wave" state with some similarities to our bond bipolaron state is a leading candidate.
This analysis reveals that the bond-Peierls coupling may be operative in the pnictides.However, additional electron-phonon interaction terms may exist in these materials.For example, since a bond between two iron atoms connects two octahedra, the motion of a single pnictogen atom out of plane within one octahedron is correlated with that of another pnictogen atom in the neighboring octahedron, and this correlated pnictogenpnictogen motion can simultaneously modulate the hopping across two iron-iron bonds, giving rise to yet another electron-phonon coupling term in these materials.The true extent to which the bond-Peierls coupling is important in determining the behavior of these materials requires more work to understand the interplay between the electron-phonon interaction terms and other features such as those presented by the multiple electronic bands, Hund's coupling, and the form and range of the effective electron-electron interactions near the Fermi surface.

20 T 5 FIG. 2 .
FIG. 2.Coulomb repulsion-mediated enhancement of bipolaronic high-Tc superconductivity.Tc of the bond-Peierls bipolaronic superconductor in units of the phonon frequency Ω for different adiabaticity ratios of the electron hopping t to Ω at intermediate coupling λ ∼ 0.5 -0.6 as a function of the onsite Hubbard repulsion U in units of t computed according to Eq. (3) from QMC simulations of the bipolaron effective mass m BP and mean squared-radius R 2 BP .Error bars represent statistical errors in QMC simulations corresponding to one standard deviation.Tc of the bond-Peierls bipolaronic superconductor exceeds the largest value of Tc ∼ 0.05Ω predicted by strong-coupling Migdal-Eliashberg theory (see Fig.1) even for large values of U/t.

FIG. 3 .
FIG.3.Bipolaron properties in the bond-Peierls model.Bipolaron properties computed from QMC calculations performed on Eq. (1) at adiabaticity parameter t/Ω = 10/3 as a function of the electron-phonon coupling λ = α 2 /(2Ωt) for different onsite Hubbard U (in units of the electron hopping t): a. Bipolaron binding energy ∆BP in units of the electron hopping t, b.Bipolaron radial size probability density distribution (absolute value squared of the bipolaron wavefunction) PBP(R) for λ = 0.5, and c.Bipolaron effective mass m BP in units of the mass of two free electrons m0 = 2me = 1/t and its mean squared-radius R 2 BP .Error bars represent statistical errors in QMC simulations corresponding to one standard deviation.Error bars in PBP(R) correspond to statistical errors smaller than the symbol size and therefore are not shown.
Figure 1 presents T c /Ω computed from Eq. (3) using m BP and R 2

FIG. 5 .
FIG.5.Bond-Peierls coupling in the pnictides.A pnictogen atom sits at the apex of an octahedron with four iron atoms residing in the x-y plane in the middle of the octahedron.a.The phonon associated with the transverse motion of the pnictogen atom out of the x-y plane in the z direction causes fluctuations in the barrier for electronic tunneling between the iron atoms within the x-y plane[39,40].b.Interference pattern along the bond connecting the iron atoms.Here there are two hopping processes: one resulting from direct overlap between dxy orbitals on neighboring iron atoms with amplitude t < 0 (dashed double-headed arrow) and another involving a second-order process in which the dxy orbitals on iron atoms overlap with the px orbital of the pnictogen atom (solid double-headed arrows), resulting in a hopping with amplitude t > 0. Destructive interference between these two processes results in a reduced net hopping along this pathway, e.g. in FeSe, see Appendix H for more details.
Appendix A: Conventions used for the definition of λ in the different models

4 FIG. 7 .
FIG.6.Computation of the bipolaron mass using a QMC approach based on a path-integral representation of the electrons combined with a diagrammatic treatment of the phonons.a. Green function G ba (τ ) as a function of imaginary time τ in different bipolaron momentum (K) sectors (solid lines) and fits to the long-τ asymptotic behavior of Eq. (B3) (dashed lines).b.The bipolaron dispersion EBP(K) constructed from EGS in the different K sectors obtained in a from the fits of G ba (τ ) in the large-τ limit to the asymptotic form in Eq. (B3).We compute the bipolaron mass by fitting the dispersion.For this data set, we find the following fitting function: EBP(K) = −9.38883+ 0.00188511K + 0.0791534K 2 , which yields m BP /m0 = 6.33 ± 1.0.Error bars represent statistical errors in QMC simulations corresponding to one standard deviation.Results shown in this figure are for the BP bipolarons at λ = 0.5, t/Ω = 10/3 and U/t = 8.

FIG. 8 .
FIG.8.Bipolaron mass m BP in units of the mass of two free electrons m0 = 2me = 1/t obtained using a QMC approach based on a path-integral representation of the electrons combined with a diagrammatic treatment of the phonons (solid line) versus that obtained using a QMC approach based on a path-integral representation of both the electrons and the phonons (dashed line).Results shown in this figure are for the BP bipolarons at t/Ω = 10/3 and U/t = 8.

FIG. 9 .
FIG. 9. Bipolaron properties in the bond-Peierls model.Bipolaron properties computed from QMC calculations performed on Eq. (1) at various values of the adiabaticity parameter t/Ω (different rows) as a function of the electron-phonon coupling λ = α 2 /(2Ωt) for different onsite Hubbard U (in units of the electron hopping t): a1,a2, a3.Bipolaron binding energy ∆BP in units of the electron hopping t, b1,b2, b3.Bipolaron effective mass m BP in units of the mass of two free electrons m0 = 2me = 1/t, and c1,c2, c3.Bipolaron mean squared-radius R 2 BP .Error bars represent statistical errors in QMC simulations corresponding to one standard deviation.

FIG. 10 .FIG. 11 .
FIG.10.Bipolaron radial size probability density distribution PBP(R) in the bond-Peierls model.PBP(R) computed from QMC calculations performed on Eq. (1) at various values of the adiabaticity parameter t/Ω for the value of λ = α 2 /(2Ωt) which maximizes Tc at onsite Hubbard repulsion U/t = 0, 8, 12; see Fig.1.Error bars in PBP(R) correspond to statistical errors smaller than the symbol size and therefore are not shown.

FIG. 12 .
FIG.12.Bipolaron properties in the Holstein model.Bipolaron properties computed from QMC calculations performed on Eq. (A2) at adiabaticity parameter t/Ω = 2 as a function of the electron-phonon coupling λ = α 2 /(8Ωt) for an onsite Hubbard U/t = 8: a. Bipolaron effective mass m BP in units of the mass of two free electrons m0 = 2me = 1/t, and b.Bipolaron mean squared-radius R 2 BP .Error bars represent statistical errors in QMC simulations corresponding to one standard deviation.

5 T
FIG. 15.Transition temperature Tc in Migdal-Eliashberg theory calculations of the Holstein model divided by the physical (renormalized) phonon frequency Ω phys calculated as a function of the conventional Migdal-Eliashberg coupling constant λH;ME from Eq. (E16).
FIG. 16. a. Momentum dependence of the electron-phonon coupling parameter in Eq. (E18) for the + mode of the bond-Peierls model calculated for densities n = 0.05, 0.1, and 0.25 carrier per site (solid lines; highest to lowest) and for the − mode of the model for density n = 0.25 carrier per site (dashed line).b.Transition temperature Tc in Migdal-Eliashberg theory calculations of the bond-Peierls model divided by bare phonon frequency Ω0 as function of the bare coupling λ 0 BP;ME for density n = 0.25 carrier per site.
σ α,β ĉi,β and σ is a vector of SU(2) Pauli operators.Here = 4α 2 2 BP /(2πΩt) in the BP model and λ ME = α 2 H /(2πΩt) in the H model. Table I summarizes the conventions we use for the definition of λ in the different models considered in this work.

TABLE I .
Definition of λ for bipolarons and in Migdal-Eliashberg (ME) theory in the bond-Peierls (BP) and Holstein (H) models as employed in this work.
)where ν and ω n are Matsubara frequencies for bosons and fermions, respectively, and τ 1 in the first of the SU(2) Pauli matrices.Here the normal component of the self energy is Σ n = iω n (1 − Z(ω n )) (where Z(ω n ) is known as the Z factor) and the anomalous component is W (ω n ).In the circular Fermi surface approximation which is reasonably accurate for n ≤ 0.25 we have ε k ≈ −4t + k 2 /(2m)