{'fraction_non_alphanumeric': 0.03271822938816264, 'fraction_numerical': 0.010095593809515005, 'mean_word_length': 4.559284576245053, 'pattern_counts': {'":': 0, '<': 40, '<?xml version=': 0, '>': 58, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Small Bipolarons in the 2-dimensional Holstein-Hubbard Model. II Quantum Bipolarons

We study the effective mass of the bipolarons and essentially the possibility to get both light and strongly bound bipolarons in the Holstein-Hubbard model and some variations in the vicinity of the adiabatic limit. Several approaches to investigate the quantum mobility of polarons and bipolarons are proposed for this model. It is found that the bipolaron mass generally remains very large except in the vicinity of the triple point of the phase diagram, where the bipolarons have several degenerate configurations at the adiabatic limit (single site (S0), two sites (S1) and quadrisinglet (QS)), while the polarons are much lighter. This degeneracy reduces the bipolaron mass significantly. The triple point of the phase diagram is washed out by the lattice quantum fluctuations which thus suppress the light bipolarons. We show that some model variations, for example a phonon dispersion may increase the stability of the (QS) bipolaron against the quantum lattice fluctuations. The triple point of the phase diagram may be stable to quantum lattice fluctuations and a very sharp mass reduction may occur, leading to bipolaron masses of the order of 100 bare electronic mass for realistic parameters. Thus such very light bipolarons could condense as a superconducting state at relatively high temperature when their interactions are not too large, that is, their density is small enough. This effect might be relevant for understanding the origin of the high Tc superconductivity of doped cuprates far enough from half filling.


Specific Problem for high TC Superconductivity
Superconducting materials at temperatures significantly higher than the maximum Tc predicted by MacMillan [6] for the standard BCS superconductivity are exceptional. Up to now, there is a wide variety of such materials which are all cuprates built with CuO 2 planes and with many kinds of interlayer dopants [2].
When the electron-phonon coupling increases too much, it is known for several decades (Migdal [7]) that the BCS theory should break down because of lattice instabilities. These instabilities are associated with the formation of polarons and bipolarons. Alexandrov et al. [17] developed later a theory of bipolaronic superconductivity where the electrons form strongly bound on-site bipolarons. They are described by a hard core boson model which could become superfluid. Unfortunately, their calculations also show that when the electron phonon coupling increases beyond the Migdal instability, the effective mass of these bosons grows exponentially fast and becomes so huge that it seems hopeless to get superconductivity in this model, at least at non-negligible temperatures.
One of us(SA) conjectured in [19,20] that the interplay of the electron-phonon coupling with a direct electronelectron repulsion could reduce significantly the bipolaron effective mass and thus favor high Tc superconductivity.
In the present paper we support this conjecture by calculating the effective mass of a single bipolaron in a 2d model which involves both an electron-phonon coupling and a direct electron-electron repulsion. For this purpose we choose to investigate first this effect in the Holstein-Hubbard model because of its simplicity.
In the absence of Hubbard repulsion, we confirm that the effective mass of the bipolaron is indeed very large [13], which is incompatible with a high Tc superconducting phase. When the Hubbard term is increased new bipolaronic states become stable. They are 2-site bipolarons which consist of two neighboring polarons bound by their magnetic interaction in a singlet state and also a bipolaron called a "quadrisinglet" (QS) which consists of the combination of four singlets sharing one common site. In the parameter region where these bipolaronic states have nearly degenerate energies, the effective bipolaron mass is sharply reduced. The drastic mass reduction is due to resonance between the different bipolarons. Certain variations of the model, such as a phonon dispersion, might even increase the binding energy of the bipolaron while allowing it to keep a very light effective mass for realistic parameters.
We first discuss some early known results about the effective masses of polarons and bipolarons. Adiabatic results (section 2) described in [25] are briefly recalled. Next we treat the quantum lattice fluctuations as a perturbation of the adiabatic limit. The main effect is to lift the bipolaron degeneracy both due to translational invariance and to the possible existence of several kinds of adiabatic bipolarons with almost the same energy (section 3). This correction is only valid for a very small quantum lattice parameter.
To extend the field of validity of our calculations we next propose to use the Toyozawa variational form. Quantum polarons and bipolarons are approximated by a selfconsistent Bloch wave that is exact at the adiabatic limit. At this limit it is demonstrated in [25] that there is not a great loss of accuracy if the shapes of polaron and bipolaron are exponentials. Thanks to this approximation we gain much simplicity for the variational form. We first apply this method to polarons (section 4) and next to bipolarons (section 5). For very small quantum lattice fluctuations the results obtained by perturbation of the adiabatic limit are practically recorvered but there are significant deviations when the fluctuations become larger. First order transitions that cannot exist physically are washed out by hybridizing several Toyozawa variational forms. Actually, the regime where the bipolaron mass is sharply reduced is swallowed up in the domain where the ground state is unbound polarons. We demonstrate in section (6) that this undesirable phenomena can be avoided by variations ofthe model which increase the stability of the (QS) bipolaron (e.g. a phonon dispersion with the appropriate sign).

The Holstein-Hubbard Model
Let us first recall our notations for the model we study here. Its Hamiltonian is +gn j (a + j + a j ) + υn j,↑ n j,↓ where j and k represent lattice sites, T is the transfer integral between nearest neighbor sites < j, k >. The electrons are Fermions represented by the standard anticommuting operators C + j,σ and C j,σ at site j with spin σ =↑ or ↓. n i,σ = C + j,σ C j,σ and n i = n i↑ + n i↓ . a + j and a j are standard creation and annihilation boson operators of phonons andhω 0 is the phonon energy of a dispersion-less optical phonon branch. g is the on-site electron phonon coupling and υ the on-site electron-electron repulsion (Hubbard interaction).
We choose E 0 = 8g 2 /hω 0 as energy unit as in [25]. Defining the position and momentum operators as with the commutation relation [u j , p j ] = i (4) we obtain the dimensionless Hamiltonian Our reduced dimensionless parameters are Despite the primitive nature of our model, it may catch important aspects of reality. However, we shall also demonstrate at the end of this paper that certain model variations could be favorable for bipolaron mass reduction. For example, we may introduce a coupling between nearest neighbor atoms so that the phonon branch is no longer dispersionless. Then the new Hamiltonian is the sum of (5) and the extra energy term −C <i,j> u i u j where < i, j > represents all the pairs of nearest neighbor sites i and j. When C > 0, the bipolaron mass reduction is enhanced while it remains strongly bound. Thus we demonstrate that relatively minor changes in the model may favor (or disfavor) superconductivity at relatively high temperature.

Polaron and Bipolaron Effective Mass
Let us first briefly recall some standard results about the effective masses of polarons and bipolarons. The Lang-Firsov unitary transformation [1] yields a new Hamiltonian where H LF = e −iSLF He iSLF with S LF = 1 2 j p j n j After this transformation, the creation operatorC + j,σ at site j acts on the vacuum by creating both an electron and a lattice distortioñ that is, a polaron. A standard but rough mean-field approximation consists in taking the average of the transfer integral for unperturbed phonons. We obtain an approximate formula for the transfer integral of this polaron When there is a single electron in the system the electronelectron interaction does not play any role. Then the effective mass of a single polaron is defined as the inverse of the second order derivative versus the wave-vector q of the polaronic energy E(q), that is, T LF . The effective mass of the polaron is that of the bare electron multiplied by exp [1/(8α)]. It becomes exponentially large when α is small.
This approximation tends to become right only when the operator corresponding to the transfer integral in eq.7 has small fluctuations. This condition is fulfilled when the pre-factor t is small. It is also fulfilled when α is large, that is, for a weak electron-phonon coupling g compared to the phonon energyhω 0 . In the antiadiabatic limit (α large), the model becomes a Hubbard model with an onsite effective electron-electron interactionŨ = (U − 1/4) which is attractive for U < 1/4 and repulsive for U > 1/4 and where the transfer integral has been renormalized. The negative U model is expected to have superconducting phases [14] for non-vanishing band filling. However, we treat here the opposite case α small which is close to the adiabatic limit.
Actually when t is small, our numerical results agree with formula (9). For larger t, the effective mass of the polaron given by (9) becomes larger than the mass we compute. Note that our result should be more reliable because it yields a lower variational energy for a single polaron.
The effective mass of the bipolaron has been calculated by Alexandrov et al in the same limit (t small) [17] for strongly bound bipolarons (that is, for U small) and considering the kinetic energy term in eq.7 as a perturbation. In our dimensionless units they found the transfer integral t b for a bipolaron If one extrapolates naively this formula for larger U , one would find that t b becomes infinite. Of course, the associated effective mass of the bipolaron cannot vanish, but our results nevertheless demonstrate that it is sharply depressed not far from the region U ≈ 0.25. Comparison of formula (9) and 10 shows that in the region where both α and t are small the effective mass, of the bipolaron is much larger by many order of magnitude than the polaron mass which is itself much larger than that of the bare electron. In most physically realistic situations, the bipolaron masses are so huge that it is unreasonable to consider that they could exhibit a Bose condensation [13].
We perform here a numerical calculation of the effective mass of the bipolaron (and also the polaron) in order to show that in some specific regions of the parameter space, when the Hubbard term becomes comparable with the electron-phonon binding energy, these effective masses can be drastically reduced so that Bose condensations of bipolarons become plausible.

The Mean-field Holstein-Hubbard Model
We calculate first the adiabatic bipolarons which are groundstate of a mean-field Hamiltonian. They are the exact solutions in the adiabatic limit when γ = α 2 is zero (that is, when the atomic kinetic energy is negligible). These spatially localized solutions are degenerate under lattice translation. This degeneracy is lifted when the atomic kinetic energy is taken into account. Within a perturbative treatment, this explicitly gives bands of extended bipolarons characterized by a wave-vector. The inverse curvature of the lowest band at zero wave-vector yields the bipolaron effective mass. This calculation become exact in principle in the limit of small γ. Note that similar methods were already developed in [22] to calculate the effective masses of discommensurations in Charge Density Waves.

The adiabatic regime
The Hamiltonian eq.5 can be written as the sum of three terms H = H el + H ph + H f (11) where H el and H ph are decoupled electron and phonon Hamiltonians respectively and H f is a fluctuation term.
n i andū i are variational parameters which are determined by minimizing the ground-state energy of the effective Hamiltonian H ad = H el + H ph . It comes out that u i = u i andn i = n i are the average of the corresponding operators. The standard mean-field approximation for the electron phonon coupling consists in neglecting the fluctuation energy H f .
Minimizing the ground-state energy of Hamiltonian eq.13 also yields Then, the ground-state of the mean field Hamiltonian H ad has the form A pair of electrons with the electronic wave function ψ i,j is created as well as the corresponding lattice distortionū i . The electronic wave function is a singlet state, that is, a symmetric function of (i, j): ψ i,j = ψ j,i . It fulfills an extended nonlinear Schroedinger equation which is exactly the same as in the adiabatic case at α = 0 [25].
The square root of the mean square lattice fluctuation (u i −ū i ) 2 1/2 = α = γ 1/2 is small of order α. Thus, the mean-field approximation obviously becomes exact in the adiabatic limit α → 0 when there are no lattice fluctuations.

Bipolarons from Anti-integrable limit and variational approximations
For an easy understanding the reader should refer to our early paper [25] were the adiabatic (or mean field) bipolarons were investigated in detail in the two-dimensional model by continuation from the anti-integrable limit (t = 0).
The main result of [25] is that we found a quite rich phase diagram with first order transition lines in the parameter space (U, t). For large t the electrons remain extended and do not self localize as bipolarons. For small t there are several kinds of structures that compete to be the bipolaron ground-state. These bipolarons were denoted (S0), (S1) and (QS). Bipolaron (S0) is mostly localized at a single site and has square symmetry. Bipolaron (S1) consists into a bound pair of polarons in a magnetic singlet state localized on two neighboring sites. It breaks the square symmetry and is oriented either in the x direction (S1) x or the y direction (S1) y . The quadrisinglet bipolaron (QS) is a combination of four singlet states with a common central site and has square symmetry.
Interesting properties are obtained at a triple point corresponding to the intersection of three first-order transition lines. At that point, and apart from the translational degeneracy, there are four different degenerate bipolarons (S0), (QS), (S1) x and (S1) y . We shall see that the quantum lattice perturbations hybridize these degenerate states and hence drastically enhance the bandwidth of the bipolaron or, equivalently, reduce its effective mass. Within a classical picture we already noticed that the energy barrier ( Peierls-Nabarro barrier) which has to be overcome to move the bipolaron through the lattice was drastically reduced.
We also investigated in [25] some approximations with exponential variational forms for the bipolarons that allow analytical calculations. The exact phase diagram calculated numerically was reproduced with a quite good accuracy with the following forms for bipolarons (S0), (S1) and (QS) respectively. A, B, and C are normalization factors and the parameters λ, λ 1 and λ 2 are optimized for energy minimization. We shall develop here a quantum analogous version of these approximations to improve our methods.

Quantum Lattice Corrections
We now treat the mean-field fluctuation H f = 1/2 i (u i − u i )(n i −n i ) as a perturbation that lifts the translational degeneracy of the mean field bipolarons (16), whose wave functions are denoted |Ω S (j) > where S represents bipolarons (S0), (S1) x , (S1) y or (QS). The index j is the site where the bipolaron (S) is located. For bipolarons (S1) x and (S1) y which occupy two adjacent sites (j x , j y ) and (j x + 1, j y ) or (j x , j y ) and (j x , j y + 1), respectively, we choose by convention j = (j x , j y ). To treat the mean-field fluctuation in lowest order, the initial Hamiltonian (11) should be projected and diagonalized in the subspace generated by all these translated wave functions. We already noticed that the bipolaron energies might be degenerate or almost degenerate so that we should take into account their possible hybridization. The eigenstates should have the general form where a S,j are coefficients to be determined by extremalization of < Ω|H|Ω > with the normalization constraint < Ω|Ω >= 1. Both < Ω|Ω > and < Ω|H|Ω > are quadratic functions of a S,j where the coefficients of matrices P and M are defined as It is important to take into account the fact that the eigenstates of the self-consistent mean-field Hamiltonian H ad are not orthogonal one with each other, since the matrix of scalar products is not diagonal. For two normalized wavefunctions |Ω > and |Ω ′ > with the form (16) and with electronic wave functions {ψ i,j } and {ψ ′ i,j }, electronic densitiesn i = −2ū i andn ′ i = −2ū ′ i respectively the scalar products eqs. (26), (25) can be calculated explicitly for the Hamiltonian (5).
where ∆ψ i,j = ψ k,l is the discrete Laplacian on a 4d lattice.
The extremalization equation of < Ω|H|Ω > with respect to A = {a S,j } with the normalization condition A * .P.A = 1, is M.A = E.P.A where E is the Lagrange parameter, which is also the eigenenergy. It can be written as an eigenvalue problem for the normalized vector Note that the extensive term N α 2 < Ω|Ω ′ > in the second term of eq.(28) does not disturb the calculations. It yields a constant term N α 2 in the effective Hamiltonian P −1/2 .M.P −1/2 which is nothing but the zero point phonon energy of the system with size N (without electrons).
Because of the translation invariance of the model, M (S,i),(S ′ ,j) and P (S,i),(S ′ ,j) only depends on j − i = n. As a result, eq.29 can be diagonalized as combinations of plane waves with the form B S,j (K) = S B S (K)e iKj with wave vector K which fulfills the eigenequation P −1/2 (K).M(K).P −1/2 (K).B = E ν (K).B(K) (30) with the Fourier coefficients Then the diagonalization of the 4 × 4 matrix P −1/2 (K).M(K).P −1/2 (K) yields the eigenenergies E ν (K). Figure (1) shows an example of calculation of these four bands in the vicinity of the triple point. Thus when there are four bipolarons that are metastable (e.g. in the vicinity of the triple point of the phase diagram) one obtains four bipolaron bands E ν (K). Within our approach the number of bipolaron bands is equal to the number of metastable states for the adiabatic bipolaron which provides the base about which we expand the eigenstates. In other regions of the phase diagram the number of metastable bipolarons changes, which induces (unrealistic) discontinuities for the number of bands. For example, when U = 0 only the bipolaron (S0) is metastable, and there is only one bipolaron band. However the lowest bipolaron band does not exhibit very sharp changes despite a small discontinuous variation. The reason that the upper bands are not reliable is that the energies of these states might be also degenerate with phonon excitations of the bipolaronic states of the lower band. The real excited states involve complex hybridization between these states.
Conversely, the bipolaronic states with the lowest energies should not hybridize significantly with the higher energy states involving phonon excitations. Thus we consider that the lowest-energy bipolaron band provides a reliable description of the bipolaron excitations close to its ground state. We use it to measure the bipolaron effective mass, that is, the inverse of the curvature T b at zero wave vector K = 0. T b is constant in all directions because of the square symmetry ( fig.1). It can be viewed as the effective hopping coefficient for the bipolaron tunnelling through the lattice and can be compared with the prediction of [17] given by formula (10), which is valid at both U and t small. Fig.2 exhibits the ratio T b /t b as a function of the effective transfer integral t for U = 0 for several values of the quantum parameter α. For U small this ratio goes to 1 when t goes to 0, which confirms the validity of formula (10) predicted by Alexandrov et al. [17] in that regime. We also note that beyond this regime when the parameters (U, t) are larger than 0, T b becomes significantly larger than t b , or equivalently the bipolaron effective mass calculated numerically drops faster than predicted by (10).
The insert of fig.3 shows the bipolaron energy gain compared to a pair of free electrons 1 at U = 0 and the fig.3 shows the effective transfer integral compared to the bare electronic transfer integral, which is negligible at the scale of the electronic energy. The bipolaron effective mass appears to be much beyond than 10 10 electronic masses even when its binding energy vanishes. It is clear that this regime U = 0 is not favorable at all for the Bose condensation of such bipolarons that should occur below a critical temperature inversely proportional to the effective mass of the quasi-particle.  When the Hubbard term increases for relatively small t, fig.4 (t = 0.04) shows that a sharp discontinuity occurs when the ground-state bipolaron becomes (S1) U > 0.25. There is a sharp increase of the tunnelling energy by five orders of magnitude for this bipolaron at α = 0.01. In that case the Peierls-Nabarro barrier calculated in the previous paper [25] is still very high and consequently there is almost no hybridization between (S0) and (S1). The smoothing of the discontinuity of the tunnelling energy is thus hardly visible.
When t is larger, the bipolaron (QS) becomes stable for U ≈ 0.23 and hybridization between (S0), (S1) and (QS) becomes significant. Actually the most important contribution to the tunnelling energy of the bipolaron comes from the hybridization between (QS) and (S1). It is responsible for the sharp increase of the tunnelling energy or equivalently the sharp drop of the bipolaron effective mass. This quantum mobility is favored when (QS) and (S1) are degenerate in energy and separate by a weak Peierls-Nabarro barrier. Then (QS) may tunnel to one of the four neighboring (S1) and the latter tunnels to its neighboring (QS) that corresponds to the initial one translated by one lattice spacing in the direction of (S1) and so on. The bipolaron tunnelling energy could reach 10 −3 the bare electronic energy which is not negligible anymore.
Let us point out that such a high mobility cannot be obtained within the approximations used in [17] which do not consider the possible degeneracies of several bipolarons. The conclusions of [13] about the physical impossibility of bipolaronic superconductivity are irrelevant for that situation.

Variational Calculation of Quantum Polarons
In principle the above approach is valid only for very small α: that is, when the quantum lattice fluctuations are small. However these fluctuations may increase drastically, especially close to the first-order transitions when there are several degenerate bipolarons that we are especially interested in. Thus it is worthwhile to improve our previous calculations by a variational approach which should be equivalent to the mean field perturbation for small quantum lattice fluctuations.
Our purpose is now to develop a quantum version of the variational forms [26] used and tested in the adiabatic case but which could hold for larger values of α. Our approach is a simplified version of those of Toyozawa (see [34] and [35] for a recent application to the polaron in 1D). We first test this method for the single polaron and will extend it in the next section for the bipolarons of the Holstein-Hubbard model.
Because of the invariance of the system under translations the wave function of a quantum single polaron is written as a Bloch wave: where Λ is a normalization factor and |Ψ P (n) > is obtained from a unique wave function |Ψ P (0) > changing all the indices i of its electronic and atomic variables into i + n. This transformation is nothing but a shift of the wave function from site 0 to n.

Toyozawa approximation
A simple variational approximation proposed by Toyozawa is to assume that the local wave function is similar to the mean-field polaron: To simplify the spin of the electron is omitted. For each wave vector K the variational energy is a function of the scalar products which does not depend on j and has to be extremalized with respect to the 2N parameters corresponding to the electronic wave function {ψ P j } and the lattice distortion {v l }. This form becomes exact in the adiabatic limit and should improve the previous perturbation theory as it is self-consistent.
A relation between the electronic density and the average of the atomic displacement can be easily taken into account in this variational form. First, let us recall that the true eigenfunctions of H are extrema of < Ψ |H|Ψ > in the full space of normalized functions Ψ . For a given normalized eigenfunction Ω of H we can consider the one parameter family of normalized functions Ψ (δ) = exp (iδ.p j )|Ω > where the coordinate u j of the atom j is changed into u j + δ. The variational energy of this wave function is should be extremal for δ = 0, which implies This result is nothing but an extension to the non-adiabatic case of the standard relation between the average atomic positions and the electronic densities.

Toyozawa Exponential Approximation: TEA
Minimizing the variational form (35) for the whole set of 2N − 1 parameters {ψ i } and {v i } with condition (40) is a complex numerical task which moreover will become even more complex when extended to the bipolaron problem. However, we can expect that the behavior of the variational parameters {φ n } and {v n } will not be far from exponential at infinity. Thus assuming simple exponentials for {φ n } and {v n } should not be a bad approximation as proposed in [26] at the adiabatic limit. Taking into account the normalization and condition (40) we postulate that the electronic wave function and the atomic modulation have the form: for each wave vector K there are only two variational parameters λ(K) and µ(K) instead of 2N for the original Toyozawa ansatz which allows much simpler calculations although we still need a numerical minimization of (35). To that aim we use a simplex method [37], which is the most efficient algorithm we tested because it avoids any precision problem due to the numerical computation of the derivatives.
The Toyozawa Exponential Ansatz (TEA) turns out to be almost as good as the full ansatz when the polarons are small, since in that case the exponential approximates quite well its shape. When the size of the polaron becomes larger, the TEA (as well as the original Toyozawa ansatz) yields a first order transition. This first order transition is well-known to exist at the adiabatic limit α = 0 where the ground state of a single electron undergoes a first order transition from a small polaron to a free electron [4,26] at t = t p ≈ 0.07486.
We define the binding energy of the quantum polaron as the difference between the energy of the extended electron at zero wave vector K = 0 and that of the bottom of the polaron band. Fig.5 shows the variation of the binding energy versus t for the quantum polaron calculated in several different approximations including the assumption that: 1. The polaron band is calculated as for the bipolaron bands (section 3) from perturbation of the mean field polaron (thin dashed line) ; 2. The polaron band is hybridized with the free electron band (thin full line) ; 3. The polaron band is calculated with the TEA approximation (thick dashed line) ; and 4. The polaron band is calculated with the HTEA approximation where small and large polarons are hybridized (thick full line) (see next section for details) .
When the quantum lattice fluctuations are small (which occurs either at the adiabatic limit α = 0 or when t is small), these approximations yield practically the same result. When α = 0 the best variational form is that which gives the lowest energy for the ground-state (that is, the largest binding energy). The results of these approximations become significantly different when t approaches the critical value t p at which the adiabatic first order transition occurs. Each of these approximations improves the previous one, since the polaron energy becomes lower at each step.
It is clear that approximations (1) and (2), which keep the polaron shape rigidly fixed to that at the adiabatic limit, are not appropriate to remove the first-order transition (see 5). The TEA approximation (3) also yields firstorder transitions, but there are two distinct transitions occurring at t = t 1 p (α) < t p and t = t 2 p (α) > t p and the amplitudes of energy discontinuities are weak because the polaron shape is determined self consistently. Fig.6 shows λ and µ values that minimize (35) at K = 0. The first of the TEA transitions (t 1 p (α) < t p ) occurs between a small and a large polaron (see fig.6) and t 1 p (α) decreases when α increases before that transition disappears for α > 0.03. The second TEA transition (t 2 p (α) > t p ) persists for large (α > 0.7) but it is hardly distinguishable on the binding energy plot ( fig.7). The transition occurs between a large polaron and a quasi-free electron with an extended phonon part (that is, µ(K = 0) tends to 1 when t is large). Note also that λ(K = 0) may become negative in the regime of large t and small α but then the polaron binding energy becomes negligible so that it is meaningless to use a polaron picture for a regime that is better described as a Fermi liquid.

Smoothing the First Order transition: HTEA
Actually, any first-order transition for the polaron groundstate (or the bipolaron) which would be obtained by any variational method cannot exist physically. The reason is that at the transition point there are two approximate wavefunctions with the same variational energy which are supposed to approximate the ground-state. It is possible to hybridize these two degenerate states to obtain a new state with a lower energy. The same arguments hold for the exact ground-state, which cannot exhibit any firstorder transition.
On the basis of these arguments, we demonstrate numerically that the two first-order transitions obtained with the TEA of polaron can be smoothed using a variational form for the wave function Ψ P (0) in eq.33 which hybridizes three wave functions, Each of these wave functions has the TEA form (41) and (42) with parameters λ 1 , µ 1 , λ 2 , µ 2 and λ 3 , µ 3 respectively. Hybridizing three wave functions instead of two has the advantage of sweeping out simultaneously the two successive first order transitions. The variational energy (35) now depends on 9 parameters λ S , µ S , β S with S ∈ {1, 2, 3}. Let us note and where (S, S ′ ) ∈ {1, 2, 3} 2 . We point out that because of the central symmetry of the TEA, the 3 × 3 matrices M and P are real. Then the energy of the ground-state E(K) has the following variational form: The extremalization of E(K) eq.46 with respect to β * 1 β * 2 and β * 3 yields the set of three equations and therefore we have to solve eigenvalue problem M β = E(K)P β: that is, E(K) is the lowest eigenvalue of the matrix P −1/2 (K)M (K)P −1/2 . That calculation is very similar to the perturbative method of the mean-field described previously in the case of the bipolaron but here the lowest eigenvalue E(K) has still to be minimized with respect to the set of six parameters (λ S , µ S ). For small t we recover the TEA results (that is, only one β is nonnegligible, fig.5). Close to the TEA first-order transitions the variational groundstate appears as the hybridization of either a small polaron and a large polaron or a large polaron and a quasi-free electron (very large polaron). A significant increase of the binding energy of polaron results from this hybridization in these crossover regions where the first order transitions are smoothed and thus removed fig.7. Furthermore our calculations show that the energy gain due to hybridization persists for large t values. In that regime the fluctuations of the quantum lattice are strong enough to hybridize two TEA states, the large polaron and the quasi-free electron, whose energies differ only slightly.
A consequence of the hybridization can be also observed on the shapes of the polaron bands. In the adiabatic limit (α = 0), the small polaron is degenerate under arbitrary lattice translations, which means that the polaron band is perfectly flat, as shown on fig.8. In the regime where the polaron is metastable for t > t p , the flat polaron band intersects the free electron band so that there is a line of wave vectors where the small polaron state and the free electron state are degenerate (see fig.8). With nonvanishing quantum lattice fluctuations ( α = 0), the degeneracies are lifted along the intersection line. Approximation 2 provides an important correction in the vicinity of t p where the adiabatic polaron becomes extended. Around the degenerate line there is a cross-over region in wave vector where the component of the free electron to the groundstate varies from almost 1 to almost 0 when K goes from 0 to π and the opposite for the component of the small polaron. Thus there is a smooth exchange of the quantum state from the large to the small polaron ( fig.10). This exchange also occurs for the upper band but in reverse order from the small polaron to the large polaron.
In the band of the TEA, for t(α) 1 < t < t(α) 2 only one first-order transition is observed in K space between a small polaron and a large one. For t > t(α) 2 two first order transition might be observed in K space at different   K values and they occur between the three kinds of polaron describded previously. They are smoothed with the HTEA.
The effective polaron mass is the inverse curvature 1/T p at the bottom of the lowest polaron band at K = 0. It can be calculated as a function of t and α and compared with the value 1/T LF obtained from the Lang-Firsov transformation (9). The variation versus t of the ratio T p /T LF is shown fig.10 for different α. For α small, this ratio T p /T LF is almost one, which confirms that the mean-field approximation used to establish formula (9) is valid for both quan- tum lattice fluctuations and t small. When t increases from zero the ratio T p /T LF increases from unity, which means that eq.9 overestimates the polaron effective mass. We already observed this effect for the bipolaron case in the formula eq.10. However this effect does not imply sharp variations.
To compare the polaron mass and the bare electronic mass, the ratio T p /t is plotted fig.11 for different values of α. For large t as well as for large α, the polaron effective mass reduces to the bare electron mass. In other words the electron becomes practically free.
When α is small, there is a sharp drop in the inverse effective polaron mass, which is reminiscent of the first order transition at α = 0 between the localized small polaron and the extended electron. This sharp variation becomes smoother and smoother as α increases. Fig.11 compares the binding energy of the polaron and its tunnelling energy T p . For small t, the binding energy is much larger than T p , while for large t it becomes much smaller. There is a value t = t c (α) where both energies are equal. In the vicinity of this region the polaron has maximum mobility while it remains reasonably well-bound (compared to this tunnelling energy!).

Variational Calculation of Quantum Bipolarons
The variational methods (34) we used for the single polaron can be extended to bipolarons with variational forms (S0), (S1) and (QS). For this purpose we write the bipolaron wave function as a Bloch wave: and we postulate an extended Toyozawa form for the local wave function

TEA quantum Bipolarons
The simple TEA approximation for the bipolaron consists in choosing ψ B j,k with the form (19) for B=(S0), (20) for B=(S1) or (21) for B=(QS) and v l with exponential forms which depend on the type of bipolaron as follows: The same arguments used to prove equation (40) imply which determines the parameters C S0 , C S1 and C QS . Using the scalar product formula (27) and (28), the variational energy (35) is calculated numerically and minimized with respect to both λ and µ parameters for each value of the wave vector K and for each bipolaron (S0),(S1) or (QS). This variational form still has a small number of parameters which allows a fast numerical minimization. The minimum energy is always found to be at the bottom of the lowest band at K = 0. The quantum corrections to the energies of bipolarons (S0), (S1) and (QS) are compared with the energy of two polarons far apart. We use the HTEA result described in the previous section ( fig.12), since we know that it yields the lowest and thus the most accurate energy for the quantum polaronic ground-state.
As for the TEA polaron each TEA bipolarons (S0), (S1) or (QS) exhibits a first-order transition when t increases between a small and a large bipolaron with the same symmetry. Actually if one compares the energies of all the possible solutions these large bipolarons are found never to be the groundstate whatever α is, because a pair of single quantum polarons has always less energy. As a result, these bipolarons always gain energy by breaking up into two polarons ( fig.12) even for large α.
In the adiabatic limit (α = 0), these TEA calculations become identical to the variational calculation which was described in [25] (see fig.14). Comparing the energies of these TEA bipolarons (without any hybridization) we construct a new phase diagram for α non zero with first-order transition lines and test how it changes when the quantum lattice parameter α increases.
The approximate calculations of the tunnelling energy for the polaron (9) and for the bipolaron (10) suggests that for U sufficiently different from 1/4 the tunnelling of a single polaron with one electron is much easier than that of a bipolaron, which contains two electrons and moreover involves a bigger lattice distortion. Thus, one should expect more generally that the energy gain generated by the quantum lattice fluctuations for the single polaron is systematically much larger than that for the bipolarons. As a result, the domain of parameters where the groundstate consists of an unbound pair of large polarons should extend at the expense of the domains of the bipolarons when α increases.  Indeed fig.14 confirms that the first-order transition lines which exist in the adiabatic limit shift to lower values of t when α increases. As a consequence the domain of existence of the (QS) bipolaronic ground-state shrinks to zero for a rather small value (approximately 0.002) of α and completely disappears for larger values.
The disappearance of the triple point of the phase diagram between bipolarons (S0),(S1) and (QS) for relatively small values of α seems to rule out our suggestion that bipolarons could become very light. We show how to recover this possibility in the last section by minor changes in the model which may restore this triple point for relatively large values of α.

HTEA quantum Bipolarons
We carefully examined whether the HTEA calculation of bipolarons, which in principle should be more accurate, could change this conclusion. Actually, it will not change it, and to not confuse or bother the reader, all details of our unsuccessful (but useful) numerical investigations are not presented.
As said previously, in principle no first-order transitions could exist for the ground-state of a pair of electrons interacting with the lattice. They are removed by hybridization of all (or only those which are relevant), degenerate bipolaron solutions (S0),(S1) or (QS) both small and large, which necessarily generates some energy gain.
The HTEA calculation for the bipolaron is similar to that for the polaron except that it may involve more bipolaronic states. We assume generally that the wave function Ψ B (0) (47) is a normalized combination of n wave functions (n depends on the number of TEA which hybridize) which may have different bipolaronic forms S=(S0), (S1 x ), (S1 y ), (QS) each of which can be small and large, so that in principle there are 8 different states. However, we have not use simultaneously all these states since there are no situations in the parameter space t, U where all their energies are simultaneously degenerate but only a relevant subset 2 . Then the energy of the ground-state E(K) has the following variational form: where and The extremalization of (54) with respect to β S is done by a diagonalization of the matrix P −1/2 (K)M (K)P −1/2 of size n × n. The variational energy is minimized with respect to parameters of eq.(53).
In all regions of the phase diagram, for small α the HTEA energy corrections for the bipolarons (S0),(S1) or (QS) are systematically much smaller than those involved by the polarons. The tunnelling energy of bipolarons is much smaller than those of the polaron.
The hybridization cross-overs which are found at each smoothed first-order transition of the TEA phase diagram remain very narrow and the hybridization energy gain is negligible. One needs to have a high bipolaronic degeneracy such as the triple point or a relatively large value of α (α ≥ 0.05) to observe non negligible crossovers. Even in that case the energy gains remain small compared to those of an unbound pair of the HTEA polarons.
If the HTEA bipolarons keep almost the same energy as the TEA bipolarons, the phase diagram fig.14 is practically unchanged. Of course, the first-order transition lines which appear in this phase diagram should now be viewed as sharp crossover lines. The crossover between the bound bipolarons and the unbound pair of polarons has been investigated with a general HTEA bipolaron form (including the latest) but no significant hybridization has been found between these two kinds of states so that we can not draw a conclusion about the nature of this transition.
The triple point is a special point of the phase diagram where the bipolarons (S0), (S1) and (QS) are degenerate and where we should expect a higher energy gain by hybridization when α is not too small. Unfortunately, this triple point disappears when α increases beyond approximately 0.002. When it just disappears the TEA bipolaron binding energy referred to two unbound HTEA polaron is just zero but then its tunnelling energy T b is maximum (but still only 10 −7 × t : that is the bipolaron effective mass is seven order of magnitude larger than those of the bare electron).
The negative conclusion of this section is that more sophisticated variational calculations does not confirm the conclusion of section (3) which was based on the assumption α small extrapolated to larger α.
The present study also shows that in the domain of small U one may have a quantum bipolaron ground-state with a large tunneling energy occurring at very large α > 0.1. This result is simply obtained with only the TEA of the small bipolaron (S0) that is proved to have a nonnegligible binding energy for both t and U small enough. Nevertheless, this result is not relevant for such large α, our approach based on a perturbative theory of the adiabatic limit fails because of too large quantum lattice fluctuations.

Phonon Dispersion Effect
We intend to show that highly degenerate point that could persist under large quantum lattice fluctuations implies very light bipolarons. To achieve that goald, a simple procedure consists in changing the model so as to favor the bipolaron (QS). If we could make it more robust to quantum lattice fluctuations it should become very light for reasonably large α by hybridization with the other degenerate bipolarons at the triple point.
We choose to introduce a phonon dispersion, but this might not be the unique way. When an electron is present at a given site it will also distort the lattice at the neighboring sites. If the sign of the dispersion is appropriate, the lattice potential at the neighboring sites is lower which favors its occupancy by electrons and thus the spatial extension of the bipolaron. The bipolaron (QS) which is more extended than the bipolaron (S0) should be favored. 3 We consider the new Hamiltonian where H is the Holstein-Hubbard Hamiltonian (1) and its reduced Hamiltonian corresponds to H (5) that gives When the coupling C is positive the dispersive term generates an effective attractive interaction between polarons. This coupling cannot exceed the value 1/2 beyond which the low wavevector phonons becomes unstable.

Adiabatic Limit
At the adiabatic limit the equation eq.15 becomes where D is the matrix : Bipolarons (S0),(S1), (QS)... which were found at the anti-integrable limit of the Holstein-Hubbard model at t = 0 persist as ground-states in this model with nonzero coupling C [25] (see the diagrams figs.15 and 17). The domain where bipolaron (QS) is the ground-state enlarges when C increases up to its maximum value 1/2. As expected the existence of a positive dispersion favors the quadrisinglet ground-state.
The first-order transition between bipolaron (S0) and (QS) becomes almost second-order and difficult to distinguish numerically since there is no symmetry breaking 3 Phonon dispersion may induce other important effects in the bipolaron structure as shown in [38] for CDW's. between these two bipolaronic states. Then as expected, there is a soft internal mode which almost vanishes at the transition on both side of the transition which corresponds to a breathing mode of the bipolaron with the same symmetry. Simultaneously the Peierls-Nabarro barrier almost vanishes. This soft mode which does not break the bipolaron symmetry is not a pinning mode and does not favor the classical mobility of this bipolaron. To that purpose the pinning mode which also softens at the first order transition between (QS) and (S1) is the most appropriate ( see fig.16 and refs. [24], [25]). When the coupling C is too large 0.2 < C < 0.5 our accuracy is limited in practice because of the bipolaron (QS) extension, which requires large system sizes we cannot afford. This problem occurs especially close to the firstorder transition between bipolaron (QS) and the extended state (see diagram fig.17). However the exponential ansatz eq.19,20,21 still fits quite well the bipolaron ground-state as shown on the diag.15. These variational forms allow us to compute quickly the bipolarons even for large C values and to determine approximately the groundstate with a reasonable accuracy (see diagrams figs.15 and 17).
As we already know the flaw of this approximate method is that spurious first-order transitions may occur. This situation happens nearby the first-order transition between (QS) and the extended state as seen on diag.17. It is due to the exponential ansatz which does not provide a good fit of the bipolaron when it becomes more extended.
However, at the triple point the bipolaron ground-state is still localized on very few sites ( fig.18) for C = 0.3 and the exponential ansatz remains sufficiently accurate.

Quantum Corrections
Same methods, as those used above for the original Holstein-Hubbard model are applied to deal with the quantum lattice fluctuations of the modified model. The degeneracy due to the translation invariance of the model is lifted according to standard perturbation theory. One gets a tight binding model as in section 3 which yields both binding and tunnelling energies of the quantum groundstate. Figs.19) shows these quantities for a strong coupling (C = 0.3).
The binding energy of the bipolaron refers to two noninteracting polarons calculated with the HTEA method, which is the most accurate. For a single polaron condition (40) becomes and we choose to write the displacement as For a large enough phonon coupling C > 0.2, in the region we investigate t < 0.1 the HTEA method only requires the hybridization between a small polaron and a large polaron. The almost second-order transition displayed by the TEA at C = 0 occurs now at a larger t P (α) 2 .
The binding energy of the quantum bipolaron is still large in that region and one notices the optimal regime where both tunnelling and binding energies have the same value.
To obtain the optimal region, a fine tuning of the parameters is required because changing them slightly can either reduce the binding energy so that the bipolaron becomes fragile against temperature or sharply increase its effective mass, killing its quantum mobility.
Phonon dispersion favors the mobility of the bipolaron because it extends the lattice distortion around the bipolaron (see fig.18) as well as the electronic wave function. Classically, this effect is manifested by internal mode softening and by the depression of the Peierls-Nabarro energy barrier (not calculated here see paper I [25]) between the different bipolarons. As a result, when the lattice is quantum the hybridization between the different bipolarons is increased, which increases the band width and decreases the effective mass.
In the vicinity of the (QS) region (see figs. 19) the effective mass of the HTEA for the bipolaron is about five times larger than the effective mass computed with the perturbative method, but the bipolaron mass is still very small. The comparison of the binding energy calculated with the two methods shows that the variational HTEA method is not accurate in the area of the QS region. Indeed the perturbative method gives a stronger binding energy and thus it is variationally better. This is likely due to the fact that when the bipolaron extends too much the TEA is not accurate because the bipolaron shape is not well approximated by the exponential. Fig.19 shows for α = 0.017 the effective mass of the bipolaron in the optimal regime that ranges not far from 100 bare electronic mass. We choose as an example the realistic optical phonon frequencyhω 0 = 1.10 −1 eV and to be in the optimal regime C = 0.3, α = 0.017 U = 0.25 t = 0.1 the initial parameter of Hamiltonian (1) must be g = 3.10 −2 eV E 0 = 6eV υ = 1.5eV t = 0.6eV c = 0.3eV . The tunnelling energy as well as the bipolaron binding energy are about 6.10 −3 eV . With such characteristic values and a bipolaron concentration not too large, a superfluid state could be expected at relatively high temperatures, perhaps few hundred degrees K. This estimate neglects the bipolaron interactions, but when their concentration becomes large these interactions cannot be neglected, especially at half filling when there is one polaron per site. Close to this close packing regime the bipolaronic structure cannot exist anymore for sure. Instead, a magnetic spatially ordered polaronic structures could occur. Further studies should investigate the situation with large electron densities.

Concluding Remarks
In some circumstances the bipolaron might become unusually light, which allows in principle the formation of superconducting states at rather high temperature with physically realistic parameters. This effect is due to the degeneracy of several bipolaronic states in the adiabatic limit for some specific regions of the phase diagram. In this situation there are small Peierls Nabarro barriers and phonon softening for the different bipolaronic states. Then the quantum lattice fluctuations lift the degeneracy between the degenerate states and may yield very light hybridized bipolarons, which however are well-bound.
We realized this situation in a modified Holstein-Hubbard model, which involves both an electron-phonon interaction and a direct repulsive electron-electron interaction.
The superconducting state of such very light bipolarons occurs for weak concentrations. When the concentration becomes larger there are strong interactions between the bipolarons, which may both break them into polarons and organize different structures (for example, magnetic).
This situation may happen in superconducting cuprates. In the undoped regime where the band of electrons is half filled, the structure can be viewed as close-packed polarons with an antiferromagnetic ordering. This polaron structure should persist for low doping till a certain electron concentration where the holes are polaron vacancies. For a sufficiently large doping the electron concentration may become low enough in order that a (first order) transition toward a superfluid of light quantum bipolarons takes place. The real phenomenology should be more complex because one should expect that the model parameters depend on the doping and thus that the system does not remain always close from the optimal regime with strongly bonded light bipolarons but move around this point. Otherwise, we suggested in [25] that in some appropriate models the (QS) bipolaron could have a d-symmetry. We have not yet realized an explicit model where such an effect occurs, but we hope to.
The numerical techniques we used (Toyozawa Exponential Ansatz) and its improvement (HTEA) where the hybridization between different states is taken into account, turned out to be very efficient to study the bipolaron mass. It should be developed to consider models with many electrons. In [28] it was proven that at adiabatic limit, the ground-state at large electron-phonon coupling was bipolaronic. Variation of the exponential ansatz may provide strong simplifications for these case and a qualitative understanding of the many-polaron problem first in the adiabatic limit, next with quantum lattice fluctuations. Finally, the problem of quantization of discrete breathers can be approached with similar techniques [36]. Abstract. We study the effective mass of the bipolarons and essentially the possibility to get both light and strongly bound bipolarons in the Holstein-Hubbard model and some variations in the vicinity of the adiabatic limit. Several approaches to investigate the quantum mobility of polarons and bipolarons are proposed for this model. First, the quantum fluctuations are treated as perturbations of the mean-field (or adiabatic) approximation of the electron-phonon coupling in order to calculate the bipolaron bands. It is found that the bipolaron mass generally remains very large except in the vicinity of the triple point of the phase diagram (see [1]), where the bipolarons have several degenerate configurations at the adiabatic limit (single site (S0), two sites (S1) and quadrisinglet (QS)), while the polarons are much lighter. This degeneracy reduces the bipolaron mass significantly. Next we improve this result by variational methods (modified Toyozawa Exponential Ansatz or TEA) valid for larger quantum perturbations away from the adiabatic limit. We first test this new method for the single polaron. We find that the triple point of the phase diagram is washed out by the lattice quantum fluctuations which thus suppress the light bipolarons. Further improvements of the method by hybridization of several TEA states do not change this conclusion. Next we show that some model variations , for example a phonon dispersion may increase the stability of the (QS) bipolaron against the quantum lattice fluctuations. We show that the triple point of the phase diagram may be stable to quantum lattice fluctuations and a very sharp mass reduction may occur, leading to bipolaron masses of the order of 100 bare electronic mass for realistic parameters. Thus we argue that such very light bipolarons could condense as a superconducting state at relatively high temperature when their interactions are not too large, that is, their density is small enough. This effect might be relevant for understanding the origin of the high Tc superconductivity of doped cuprates far enough from half filling.

Specific Problem for high TC Superconductivity
Superconducting materials at temperatures significantly higher than the maximum Tc predicted by MacMillan [2] for the standard BCS superconductivity [3] are exceptional [4]. Up to now, there is a wide variety of such materials which are all cuprates built with CuO 2 planes and with many kinds of interlayer dopants [5]. When the electron-phonon coupling increases too much, it is known for several decades (Migdal [6]) that the BCS theory should break down because of lattice instabilities. These instabilities are associated with the formation of polarons and bipolarons. Alexandrov et al. [7] developed later a theory of bipolaronic superconductivity where the electrons form strongly bound on-site bipolarons. They are described by a hard core boson model which could become superfluid. Unfortunately, their calculations also show that when the electron phonon coupling increases beyond the Migdal instability, the effective mass of these bosons grows exponentially fast and becomes so huge that it seems hopeless to get superconductivity in this model, at least at non-negligible temperatures.
One of us(SA) conjectured in [8,9] that the interplay of the electron-phonon coupling with a direct electronelectron repulsion could reduce significantly the bipolaron effective mass and thus favor high Tc superconductivity.
In the present paper we support this conjecture by calculating the effective mass of a single bipolaron in a 2d model which involves both an electron-phonon coupling and a direct electron-electron repulsion. For this purpose we choose to investigate first this effect in the Holstein-Hubbard model because of its simplicity.
In the absence of Hubbard repulsion, we confirm that the effective mass of the bipolaron is indeed very large [10], which is incompatible with a high Tc superconducting phase. When the Hubbard term is increased new bipolaronic states become stable. They are 2-site bipolarons which consist of two neighboring polarons bound by their magnetic interaction in a singlet state and also a bipolaron called a "quadrisinglet" (QS) which consists of the combination of four singlets sharing one common site. In the parameter region where these bipolaronic states have nearly degenerate energies, the effective bipolaron mass is sharply reduced. The drastic mass reduction is due to resonance between the different bipolarons. Certain variations of the model, such as a phonon dispersion, might even increase the binding energy of the bipolaron while allowing it to keep a very light effective mass for realistic parameters.
We first discuss some early known results about the effective masses of polarons and bipolarons. Adiabatic results (section 2) described in [1] are briefly recalled. Next we treat the quantum lattice fluctuations as a perturbation of the adiabatic limit. The main effect is to lift the bipolaron degeneracy both due to translational invariance and to the possible existence of several kinds of adiabatic bipolarons with almost the same energy (section 3). This correction is only valid for a very small quantum lattice parameter.
To extend the field of validity of our calculations we next propose to use the Toyozawa variational form. Quantum polarons and bipolarons are approximated by a selfconsistent Bloch wave that is exact at the adiabatic limit. At this limit it is demonstrated in [1] that there is not a great loss of accuracy if the shapes of polaron and bipolaron are exponentials. Thanks to this approximation we gain much simplicity for the variational form. We first apply this method to polarons (section 4) and next to bipolarons (section 5). For very small quantum lattice fluctuations the results obtained by perturbation of the adiabatic limit are practically recorvered but there are significant deviations when the fluctuations become larger. First order transitions that cannot exist physically are washed out by hybridizing several Toyozawa variational forms. Actually, the regime where the bipolaron mass is sharply reduced is swallowed up in the domain where the groundstate is unbound polarons. We demonstrate in section (6) that this undesirable phenomena can be avoided by variations ofthe model which increase the stability of the (QS) bipolaron (e.g. a phonon dispersion with the appropriate sign).

The Holstein-Hubbard Model
Let us first recall our notations for the model we study here. Its Hamiltonian is where j and k represent lattice sites, T is the transfer integral between nearest neighbor sites < j, k >. The electrons are Fermions represented by the standard anticommuting operators C + j,σ and C j,σ at site j with spin σ =↑ or ↓. n i,σ = C + j,σ C j,σ and n i = n i↑ + n i↓ . a + j and a j are standard creation and annihilation boson operators of phonons andhω 0 is the phonon energy of a dispersion-less optical phonon branch. g is the on-site electron phonon coupling and υ the on-site electron-electron repulsion (Hubbard interaction).
We choose E 0 = 8g 2 /hω 0 as energy unit as in [1]. Defining the position and momentum operators as with the commutation relation [u j , p j ] = i (4) we obtain the dimensionless Hamiltonian Our reduced dimensionless parameters are Despite the primitive nature of our model, it may catch important aspects of reality. However, we shall also demonstrate at the end of this paper that certain model variations could be favorable for bipolaron mass reduction. For example, we may introduce a coupling between nearest neighbor atoms so that the phonon branch is no longer dispersionless. Then the new Hamiltonian is the sum of (5) and the extra energy term −C <i,j> u i u j where < i, j > represents all the pairs of nearest neighbor sites i and j. When C > 0, the bipolaron mass reduction is enhanced while it remains strongly bound. Thus we demonstrate that relatively minor changes in the model may favor (or disfavor) superconductivity at relatively high temperature.

Polaron and Bipolaron Effective Mass
Let us first briefly recall some standard results about the effective masses of polarons and bipolarons. The Lang-Firsov unitary transformation [11] yields a new Hamiltonian where H LF = e −iSLF He iSLF with S LF = 1 2 j p j n j After this transformation, the creation operatorC + j,σ at site j acts on the vacuum by creating both an electron and a lattice distortioñ that is, a polaron. A standard but rough mean-field approximation consists in taking the average of the transfer integral for unperturbed phonons. We obtain an approximate formula for the transfer integral of this polaron When there is a single electron in the system the electronelectron interaction does not play any role. Then the effective mass of a single polaron is defined as the inverse of the second order derivative versus the wave-vector q of the polaronic energy E(q), that is, T LF . The effective mass of the polaron is that of the bare electron multiplied by exp [1/(8α)]. It becomes exponentially large when α is small.
This approximation tends to become right only when the operator corresponding to the transfer integral in eq.7 has small fluctuations. This condition is fulfilled when the pre-factor t is small. It is also fulfilled when α is large, that is, for a weak electron-phonon coupling g compared to the phonon energyhω 0 . In the antiadiabatic limit (α large), the model becomes a Hubbard model with an onsite effective electron-electron interactionŨ = (U − 1/4) which is attractive for U < 1/4 and repulsive for U > 1/4 and where the transfer integral has been renormalized. The negative U model is expected to have superconducting phases [12] for non-vanishing band filling. However, we treat here the opposite case α small which is close to the adiabatic limit.
Actually when t is small, our numerical results agree with formula (9). For larger t, the effective mass of the polaron given by (9) becomes larger than the mass we compute. Note that our result should be more reliable because it yields a lower variational energy for a single polaron.
The effective mass of the bipolaron has been calculated by Alexandrov et al in the same limit (t small) [7] for strongly bound bipolarons (that is, for U small) and considering the kinetic energy term in eq.7 as a perturbation. In our dimensionless units they found the transfer integral t b for a bipolaron If one extrapolates naively this formula for larger U , one would find that t b becomes infinite. Of course, the associated effective mass of the bipolaron cannot vanish, but our results nevertheless demonstrate that it is sharply depressed not far from the region U ≈ 0.25. Comparison of formula (9) and 10 shows that in the region where both α and t are small the effective mass, of the bipolaron is much larger by many order of magnitude than the polaron mass which is itself much larger than that of the bare electron. In most physically realistic situations, the bipolaron masses are so huge that it is unreasonable to consider that they could exhibit a Bose condensation [10].
We perform here a numerical calculation of the effective mass of the bipolaron (and also the polaron) in order to show that in some specific regions of the parameter space, when the Hubbard term becomes comparable with the electron-phonon binding energy, these effective masses can be drastically reduced so that Bose condensations of bipolarons become plausible.

The Mean-field Holstein-Hubbard Model
We calculate first the adiabatic bipolarons which are groundstate of a mean-field Hamiltonian. They are the exact solutions in the adiabatic limit when γ = α 2 is zero (that is, when the atomic kinetic energy is negligible). These spatially localized solutions are degenerate under lattice translation. This degeneracy is lifted when the atomic kinetic energy is taken into account. Within a perturbative treatment, this explicitly gives bands of extended bipolarons characterized by a wave-vector. The inverse curvature of the lowest band at zero wave-vector yields the bipolaron effective mass. This calculation become exact in principle in the limit of small γ. Note that similar methods were already developed in [13] to calculate the effective masses of discommensurations in Charge Density Waves.

The adiabatic regime
The Hamiltonian eq.5 can be written as the sum of three terms where H el and H ph are decoupled electron and phonon Hamiltonians respectively and H f is a fluctuation term.
n i andū i are variational parameters which are determined by minimizing the ground-state energy of the effective Hamiltonian H ad = H el + H ph . It comes out that u i = u i andn i = n i are the average of the corresponding operators. The standard mean-field approximation for the electron phonon coupling consists in neglecting the fluctuation energy H f .
Minimizing the ground-state energy of Hamiltonian eq.13 also yields Then, the ground-state of the mean field Hamiltonian H ad has the form A pair of electrons with the electronic wave function ψ i,j is created as well as the corresponding lattice distortionū i . The electronic wave function is a singlet state, that is, a symmetric function of (i, j): ψ i,j = ψ j,i . It fulfills an extended nonlinear Schroedinger equation which is exactly the same as in the adiabatic case at α = 0 [1].
The square root of the mean square lattice fluctuation (u i −ū i ) 2 1/2 = α = γ 1/2 is small of order α. Thus, the mean-field approximation obviously becomes exact in the adiabatic limit α → 0 when there are no lattice fluctuations.

Bipolarons from Anti-integrable limit and variational approximations
For an easy understanding the reader should refer to our early paper [1] were the adiabatic (or mean field) bipolarons were investigated in detail in the two-dimensional model by continuation from the anti-integrable limit (t = 0).
The main result of [1] is that we found a quite rich phase diagram with first order transition lines in the parameter space (U, t). For large t the electrons remain extended and do not self localize as bipolarons. For small t there are several kinds of structures that compete to be the bipolaron ground-state. These bipolarons were denoted (S0), (S1) and (QS). Bipolaron (S0) is mostly localized at a single site and has square symmetry. Bipolaron (S1) consists into a bound pair of polarons in a magnetic singlet state localized on two neighboring sites. It breaks the square symmetry and is oriented either in the x direction (S1) x or the y direction (S1) y . The quadrisinglet bipolaron (QS) is a combination of four singlet states with a common central site and has square symmetry.
Interesting properties are obtained at a triple point corresponding to the intersection of three first-order transition lines. At that point, and apart from the translational degeneracy, there are four different degenerate bipolarons (S0), (QS), (S1) x and (S1) y . We shall see that the quantum lattice perturbations hybridize these degenerate states and hence drastically enhance the bandwidth of the bipolaron or, equivalently, reduce its effective mass. Within a classical picture we already noticed that the energy barrier ( Peierls-Nabarro barrier) which has to be overcome to move the bipolaron through the lattice was drastically reduced.
We also investigated in [1] some approximations with exponential variational forms for the bipolarons that allow analytical calculations. The exact phase diagram calculated numerically was reproduced with a quite good accuracy with the following forms for bipolarons (S0), (S1) and (QS) respectively. A, B, and C are normalization factors and the parameters λ, λ 1 and λ 2 are optimized for energy minimization. We shall develop here a quantum analogous version of these approximations to improve our methods.

Quantum Lattice Corrections
We now treat the mean-field fluctuation H f = 1/2 i (u i − u i )(n i −n i ) as a perturbation that lifts the translational degeneracy of the mean field bipolarons (16), whose wave functions are denoted |Ω S (j) > where S represents bipolarons (S0), (S1) x , (S1) y or (QS). The index j is the site where the bipolaron (S) is located. For bipolarons (S1) x and (S1) y which occupy two adjacent sites (j x , j y ) and (j x + 1, j y ) or (j x , j y ) and (j x , j y + 1), respectively, we choose by convention j = (j x , j y ). To treat the mean-field fluctuation in lowest order, the initial Hamiltonian (11) should be projected and diagonalized in the subspace generated by all these translated wave functions. We already noticed that the bipolaron energies might be degenerate or almost degenerate so that we should take into account their possible hybridization. The eigenstates should have the general form where a S,j are coefficients to be determined by extremalization of < Ω|H|Ω > with the normalization constraint < Ω|Ω >= 1. Both < Ω|Ω > and < Ω|H|Ω > are quadratic functions of a S,j < Ω|H|Ω > = where the coefficients of matrices P and M are defined as It is important to take into account the fact that the eigenstates of the self-consistent mean-field Hamiltonian H ad are not orthogonal one with each other, since the matrix of scalar products is not diagonal. For two normalized wavefunctions |Ω > and |Ω ′ > with the form (16) and with electronic wave functions {ψ i,j } and {ψ ′ i,j }, electronic densitiesn i = −2ū i andn ′ i = −2ū ′ i respectively the scalar products eqs. (26), (25) can be calculated explicitly for the Hamiltonian (5).
where ∆ψ i,j = ψ k,l is the discrete Laplacian on a 4d lattice.
The extremalization equation of < Ω|H|Ω > with respect to A = {a S,j } with the normalization condition A * .P.A = 1, is M.A = E.P.A where E is the Lagrange parameter, which is also the eigenenergy. It can be written as an eigenvalue problem for the normalized vector Note that the extensive term N α 2 < Ω|Ω ′ > in the second term of eq.(28) does not disturb the calculations. It yields a constant term N α 2 in the effective Hamiltonian P −1/2 .M.P −1/2 which is nothing but the zero point phonon energy of the system with size N (without electrons).
Because of the translation invariance of the model, M (S,i),(S ′ ,j) and P (S,i),(S ′ ,j) only depends on j − i = n. As a result, eq.29 can be diagonalized as combinations of plane waves with the form B S,j (K) = S B S (K)e iKj with wave vector K which fulfills the eigenequation P −1/2 (K).M(K).P −1/2 (K).B = E ν (K).B(K) (30) with the Fourier coefficients P S,S ′ (K) = n P (S,j),(S ′ ,j+n) e iKn (31) Then the diagonalization of the 4 × 4 matrix P −1/2 (K).M(K).P −1/2 (K) yields the eigenenergies E ν (K). Figure (1) shows an example of calculation of these four bands in the vicinity of the triple point. Thus when there are four bipolarons that are metastable (e.g. in the vicinity of the triple point of the phase diagram) one obtains four bipolaron bands E ν (K). Within our approach the number of bipolaron bands is equal to the number of metastable states for the adiabatic bipolaron which provides the base about which we expand the eigenstates. In other regions of the phase diagram the number of metastable bipolarons changes, which induces (unrealistic) discontinuities for the number of bands. For example, when U = 0 only the bipolaron (S0) is metastable, and there is only one bipolaron band. However the lowest bipolaron band does not exhibit very sharp changes despite a small discontinuous variation. The reason that the upper bands are not reliable is that the energies of these states might be also degenerate with phonon excitations of the bipolaronic states of the lower band. The real excited states involve complex hybridization between these states.
Conversely, the bipolaronic states with the lowest energies should not hybridize significantly with the higher energy states involving phonon excitations. Thus we consider that the lowest-energy bipolaron band provides a reliable description of the bipolaron excitations close to its ground-state. We use it to measure the bipolaron effective mass, that is, the inverse of the curvature T b at zero wave vector K = 0. T b is constant in all directions because of the square symmetry ( fig.1). It can be viewed as the effective hopping coefficient for the bipolaron tunnelling through the lattice and can be compared with the prediction of [7] given by formula (10), which is valid at both U and t small. Fig.2 exhibits the ratio T b /t b as a function of the effective transfer integral t for U = 0 for several values of the quantum parameter α. For U small this ratio goes to 1 when t goes to 0, which confirms the validity of formula (10) predicted by Alexandrov et al. [7] in that regime. We also note that beyond this regime when the parameters (U, t) are larger than 0, T b becomes significantly larger than t b , or equivalently the bipolaron effective mass calculated numerically drops faster than predicted by (10).
The insert of fig.3 shows the bipolaron energy gain compared to a pair of free electrons 1 at U = 0 and the fig.3 shows the effective transfer integral compared to the bare electronic transfer integral, which is negligible at the scale of the electronic energy. The bipolaron effective mass appears to be much beyond than 10 10 electronic masses even when its binding energy vanishes. It is clear that this regime U = 0 is not favorable at all for the Bose condensation of such bipolarons that should occur below a critical temperature inversely proportional to the effective mass of the quasi-particle. When the Hubbard term increases for relatively small t, fig.4 (t = 0.04) shows that a sharp discontinuity occurs when the ground-state bipolaron becomes (S1) U > 0.25. There is a sharp increase of the tunnelling energy by five orders of magnitude for this bipolaron at α = 0.01. In that case the Peierls-Nabarro barrier calculated in the previous paper [1] is still very high and consequently there is almost no hybridization between (S0) and (S1). The smoothing of the discontinuity of the tunnelling energy is thus hardly visible.
When t is larger, the bipolaron (QS) becomes stable for U ≈ 0.23 and hybridization between (S0), (S1) and (QS) becomes significant. Actually the most important contribution to the tunnelling energy of the bipolaron comes from the hybridization between (QS) and (S1). It is responsible for the sharp increase of the tunnelling energy or equivalently the sharp drop of the bipolaron effective mass. This quantum mobility is favored when (QS) and (S1) are degenerate in energy and separate by a weak Peierls-Nabarro barrier. Then (QS) may tunnel to one of the four neighboring (S1) and the latter tunnels to its neighboring (QS) that corresponds to the initial one translated by one lattice spacing in the direction of (S1) and so on. The bipolaron tunnelling energy could reach 10 −3 the bare electronic energy which is not negligible anymore.
Let us point out that such a high mobility cannot be obtained within the approximations used in [7] which do not consider the possible degeneracies of several bipolarons. The conclusions of [10] about the physical impossibility of bipolaronic superconductivity are irrelevant for that situation.

Variational Calculation of Quantum Polarons
In principle the above approach is valid only for very small α: that is, when the quantum lattice fluctuations are small. However these fluctuations may increase drastically, especially close to the first-order transitions when there are several degenerate bipolarons that we are especially interested in. Thus it is worthwhile to improve our previous calculations by a variational approach which should be equivalent to the mean field perturbation for small quantum lattice fluctuations.
Our purpose is now to develop a quantum version of the variational forms [14] used and tested in the adiabatic case but which could hold for larger values of α. Our approach is a simplified version of those of Toyozawa (see [15] and [16] for a recent application to the polaron in 1D). We first test this method for the single polaron and will extend it in the next section for the bipolarons of the Holstein-Hubbard model.
Because of the invariance of the system under translations the wave function of a quantum single polaron is written as a Bloch wave: where Λ is a normalization factor and |Ψ P (n) > is obtained from a unique wave function |Ψ P (0) > changing all the indices i of its electronic and atomic variables into i + n. This transformation is nothing but a shift of the wave function from site 0 to n.

Toyozawa approximation
A simple variational approximation proposed by Toyozawa is to assume that the local wave function is similar to the mean-field polaron: To simplify the spin of the electron is omitted. For each wave vector K the variational energy < Ω P (K) |H|Ω P (K) >= p e iKp < Ψ P (j)|H|Ψ P (j + p) > p e iKp < Ψ P (j)|Ψ P (j + p) > is a function of the scalar products which does not depend on j and has to be extremalized with respect to the 2N parameters corresponding to the electronic wave function {ψ P j } and the lattice distortion {v l }. This form becomes exact in the adiabatic limit and should improve the previous perturbation theory as it is self-consistent.
A relation between the electronic density and the average of the atomic displacement can be easily taken into account in this variational form. First, let us recall that the true eigenfunctions of H are extrema of < Ψ |H|Ψ > in the full space of normalized functions Ψ . For a given normalized eigenfunction Ω of H we can consider the one parameter family of normalized functions Ψ (δ) = exp (iδ.p j )|Ω > where the coordinate u j of the atom j is changed into u j + δ. The variational energy of this wave function is should be extremal for δ = 0, which implies This result is nothing but an extension to the non-adiabatic case of the standard relation between the average atomic positions and the electronic densities.

Toyozawa Exponential Approximation: TEA
Minimizing the variational form (35) for the whole set of 2N − 1 parameters {ψ i } and {v i } with condition (40) is a complex numerical task which moreover will become even more complex when extended to the bipolaron problem. However, we can expect that the behavior of the variational parameters {φ n } and {v n } will not be far from exponential at infinity. Thus assuming simple exponentials for {φ n } and {v n } should not be a bad approximation as proposed in [14] at the adiabatic limit. Taking into account the normalization and condition (40) we postulate that the electronic wave function and the atomic modulation have the form: for each wave vector K there are only two variational parameters λ(K) and µ(K) instead of 2N for the original Toyozawa ansatz which allows much simpler calculations although we still need a numerical minimization of (35). To that aim we use a simplex method [17], which is the most efficient algorithm we tested because it avoids any precision problem due to the numerical computation of the derivatives.
The Toyozawa Exponential Ansatz (TEA) turns out to be almost as good as the full ansatz when the polarons are small, since in that case the exponential approximates quite well its shape. When the size of the polaron becomes larger, the TEA (as well as the original Toyozawa ansatz) yields a first order transition. This first order transition is well-known to exist at the adiabatic limit α = 0 where the ground state of a single electron undergoes a first order transition from a small polaron to a free electron [18,14] at t = t p ≈ 0.07486.
We define the binding energy of the quantum polaron as the difference between the energy of the extended electron at zero wave vector K = 0 and that of the bottom of the polaron band. Fig.5 shows the variation of the binding energy versus t for the quantum polaron calculated in several different approximations including the assumption that: 1. The polaron band is calculated as for the bipolaron bands (section 3) from perturbation of the mean field polaron (thin dashed line) ; 2. The polaron band is hybridized with the free electron band (thin full line) ; 3. The polaron band is calculated with the TEA approximation (thick dashed line) ; and 4. The polaron band is calculated with the HTEA approximation where small and large polarons are hybridized (thick full line) (see next section for details) .
When the quantum lattice fluctuations are small (which occurs either at the adiabatic limit α = 0 or when t is small), these approximations yield practically the same result. When α = 0 the best variational form is that which gives the lowest energy for the ground-state (that is, the largest binding energy). The results of these approximations become significantly different when t approaches the critical value t p at which the adiabatic first order transition occurs. Each of these approximations improves the previous one, since the polaron energy becomes lower at each step.
It is clear that approximations (1) and (2), which keep the polaron shape rigidly fixed to that at the adiabatic limit, are not appropriate to remove the first-order transition (see 5). The TEA approximation (3) also yields firstorder transitions, but there are two distinct transitions occurring at t = t 1 p (α) < t p and t = t 2 p (α) > t p and the amplitudes of energy discontinuities are weak because the polaron shape is determined self consistently. Fig.6 shows λ and µ values that minimize (35) at K = 0. The first of the TEA transitions (t 1 p (α) < t p ) occurs between a small and a large polaron (see fig.6) and t 1 p (α) decreases when α increases before that transition disappears for α > 0.03. The second TEA transition (t 2 p (α) > t p ) persists for large (α > 0.7) but it is hardly distinguishable on the binding energy plot (fig.7). The transition occurs between a large polaron and a quasi-free electron with an extended phonon part (that is, µ(K = 0) tends to 1 when t is large). Note also that λ(K = 0) may become negative in the regime of large t and small α but then the polaron binding energy becomes negligible so that it is meaningless to use a polaron picture for a regime that is better described as a Fermi liquid.

Smoothing the First Order transition: HTEA
Actually, any first-order transition for the polaron groundstate (or the bipolaron) which would be obtained by any variational method cannot exist physically. The reason is that at the transition point there are two approximate wavefunctions with the same variational energy which are supposed to approximate the ground-state. It is possible to hybridize these two degenerate states to obtain a new state with a lower energy. The same arguments hold for the exact ground-state, which cannot exhibit any firstorder transition.
On the basis of these arguments, we demonstrate numerically that the two first-order transitions obtained with the TEA of polaron can be smoothed using a variational form for the wave function Ψ P (0) in eq.33 which hybridizes three wave functions, Each of these wave functions has the TEA form (41) and (42) with parameters λ 1 , µ 1 , λ 2 , µ 2 and λ 3 , µ 3 respectively. Hybridizing three wave functions instead of two has the advantage of sweeping out simultaneously the two successive first order transitions. The variational energy (35) now depends on 9 parameters λ S , µ S , β S with S ∈ {1, 2, 3}. Let us note and where (S, S ′ ) ∈ {1, 2, 3} 2 . We point out that because of the central symmetry of the TEA, the 3 × 3 matrices M and P are real. Then the energy of the ground-state E(K) has the following variational form: The extremalization of E(K) eq.46 with respect to β * 1 β * 2 and β * 3 yields the set of three equations and therefore we have to solve eigenvalue problem M β = E(K)P β: that is, E(K) is the lowest eigenvalue of the matrix P −1/2 (K)M (K)P −1/2 . That calculation is very similar to the perturbative method of the mean-field described previously in the case of the bipolaron but here the lowest eigenvalue E(K) has still to be minimized with respect to the set of six parameters (λ S , µ S ). For small t we recover the TEA results (that is, only one β is nonnegligible, fig.5). Close to the TEA first-order transitions the variational ground-state appears as the hybridization of either a small polaron and a large polaron or a large polaron and a quasi-free electron (very large polaron). A significant increase of the binding energy of polaron results from this hybridization in these crossover regions where the first order transitions are smoothed and thus removed fig.7. Furthermore our calculations show that the energy gain due to hybridization persists for large t values. In that regime the fluctuations of the quantum lattice are strong enough to hybridize two TEA states, the large polaron and the quasi-free electron, whose energies differ only slightly.
A consequence of the hybridization can be also observed on the shapes of the polaron bands. In the adiabatic limit (α = 0), the small polaron is degenerate under arbitrary lattice translations, which means that the polaron band is perfectly flat, as shown on fig.8. In the regime where the polaron is metastable for t > t p , the flat polaron band intersects the free electron band so that there is a line of wave vectors where the small polaron state and the free electron state are degenerate (see fig.8). With nonvanishing quantum lattice fluctuations ( α = 0), the degeneracies are lifted along the intersection line. Approximation 2 provides an important correction in the vicinity of t p where the adiabatic polaron becomes extended. Around the degenerate line there is a cross-over region in wave vector where the component of the free electron to the ground-state varies from almost 1 to almost 0 when K goes from 0 to π and the opposite for the component of the small polaron. Thus there is a smooth exchange of the quantum state from the large to the small polaron ( fig.9). This exchange also occurs for the upper band but in reverse order from the small polaron to the large polaron.
In the band of the TEA, for t 1 (α) < t < t 2 (α) only one first-order transition is observed in K space between a small polaron and a large one. For t > t 2 (α) two first order transition might be observed in K space at different   K values and they occur between the three kinds of polaron describded previously. They are smoothed with the HTEA.
The effective polaron mass is the inverse curvature 1/T p at the bottom of the lowest polaron band at K = 0. It can be calculated as a function of t and α and compared with the value 1/T LF obtained from the Lang-Firsov transformation (9). The variation versus t of the ratio T p /T LF is shown fig.10 for different α. For α small, this ratio T p /T LF is almost one, which confirms that the mean-field approximation used to establish formula (9) is valid for both quan- tum lattice fluctuations and t small. When t increases from zero the ratio T p /T LF increases from unity, which means that eq.9 overestimates the polaron effective mass. We already observed this effect for the bipolaron case in the formula eq.10. However this effect does not imply sharp variations.
To compare the polaron mass and the bare electronic mass, the ratio T p /t is plotted fig.11 for different values of α. For large t as well as for large α, the polaron effective mass reduces to the bare electron mass. In other words the electron becomes practically free.
When α is small, there is a sharp drop in the inverse effective polaron mass, which is reminiscent of the first order transition at α = 0 between the localized small polaron and the extended electron. This sharp variation becomes smoother and smoother as α increases. Fig.11 compares the binding energy of the polaron and its tunnelling energy T p . For small t, the binding energy is much larger than T p , while for large t it becomes much smaller. There is a value t = t c (α) where both energies are equal. In the vicinity of this region the polaron has maximum mobility while it remains reasonably well-bound (compared to this tunnelling energy!).

Variational Calculation of Quantum Bipolarons
The variational methods (34) we used for the single polaron can be extended to bipolarons with variational forms (S0), (S1) and (QS). For this purpose we write the bipolaron wave function as a Bloch wave: and we postulate an extended Toyozawa form for the local wave function

TEA quantum Bipolarons
The simple TEA approximation for the bipolaron consists in choosing ψ B j,k with the form (19) for B=(S0), (20) for B=(S1) or (21) for B=(QS) and v l with exponential forms which depend on the type of bipolaron as follows: The same arguments used to prove equation (40) imply which determines the parameters C S0 , C S1 and C QS . Using the scalar product formula (27) and (28), the variational energy (35) is calculated numerically and minimized with respect to both λ and µ parameters for each value of the wave vector K and for each bipolaron (S0),(S1) or (QS) (see fig.13). This variational form still has a small number of parameters which allows a fast numerical minimization. The minimum energy is always found to be at the bottom of the lowest band at K = 0. The quantum corrections to the energies of bipolarons (S0), (S1) and (QS) are compared with the energy of two polarons far apart. We use the HTEA result described in the previous section ( fig.12), since we know that it yields the lowest and thus the most accurate energy for the quantum polaronic ground-state.
As for the TEA polaron each TEA bipolarons (S0), (S1) or (QS) exhibits a first-order transition when t increases between a small and a large bipolaron with the same symmetry. Actually if one compares the energies of all the possible solutions these large bipolarons are found never to be the ground-state whatever α is, because a pair of single quantum polarons has always less energy. As a result, these bipolarons always gain energy by breaking up into two polarons ( fig.12) even for large α.
In the adiabatic limit (α = 0), these TEA calculations become identical to the variational calculation which was described in [1] (see fig.14). Comparing the energies of these TEA bipolarons (without any hybridization) we construct a new phase diagram for α non zero with firstorder transition lines and test how it changes when the quantum lattice parameter α increases.
The approximate calculations of the tunnelling energy for the polaron (9) and for the bipolaron (10) suggests that for U sufficiently different from 1/4 the tunnelling of a single polaron with one electron is much easier than that of a bipolaron, which contains two electrons and moreover involves a bigger lattice distortion. Thus, one should expect more generally that the energy gain generated by the quantum lattice fluctuations for the single polaron is systematically much larger than that for the bipolarons. As a result, the domain of parameters where the ground-state consists of an unbound pair of large polarons should extend at the expense of the domains of the bipolarons when α increases.  Indeed fig.14 confirms that the first-order transition lines which exist in the adiabatic limit shift to lower values of t when α increases. As a consequence the domain of existence of the (QS) bipolaronic ground-state shrinks to zero for a rather small value (approximately 0.002) of α and completely disappears for larger values.
The disappearance of the triple point of the phase diagram between bipolarons (S0),(S1) and (QS) for relatively small values of α seems to rule out our suggestion that bipolarons could become very light. We show how to recover this possibility in the last section by minor changes in the model which may restore this triple point for relatively large values of α.

HTEA quantum Bipolarons
We carefully examined whether the HTEA calculation of bipolarons, which in principle should be more accurate, could change this conclusion. Actually, it will not change it, and to not confuse or bother the reader, all details of our unsuccessful (but useful) numerical investigations are not presented.
As said previously, in principle no first-order transitions could exist for the ground-state of a pair of electrons interacting with the lattice. They are removed by hybridization of all (or only those which are relevant), degenerate bipolaron solutions (S0),(S1) or (QS) both small and large, which necessarily generates some energy gain.
The HTEA calculation for the bipolaron is similar to that for the polaron except that it may involve more bipolaronic states. We assume generally that the wave function Ψ B (0) (47) is a normalized combination of n wave functions (n depends on the number of TEA which hybridize) which may have different bipolaronic forms S=(S0), (S1 x ), (S1 y ), (QS) each of which can be small and large, so that in principle there are 8 different states. However, we have not use simultaneously all these states since there are no situations in the parameter space t, U where all their energies are simultaneously degenerate but only a relevant subset 2 . Then the energy of the ground-state E(K) has the following variational form: where and The extremalization of (54) with respect to β S is done by a diagonalization of the matrix P −1/2 (K)M (K)P −1/2 of size n × n. The variational energy is minimized with respect to parameters of eq.(53).
In all regions of the phase diagram, for small α the HTEA energy corrections for the bipolarons (S0),(S1) or (QS) are systematically much smaller than those involved by the polarons. The tunnelling energy of bipolarons is much smaller than those of the polaron.
The hybridization cross-overs which are found at each smoothed first-order transition of the TEA phase diagram remain very narrow and the hybridization energy gain is negligible. One needs to have a high bipolaronic degeneracy such as the triple point or a relatively large value of α (α ≥ 0.05) to observe non negligible crossovers. Even in that case the energy gains remain small compared to those of an unbound pair of the HTEA polarons.
If the HTEA bipolarons keep almost the same energy as the TEA bipolarons, the phase diagram fig.14 is practically unchanged. Of course, the first-order transition lines which appear in this phase diagram should now be viewed as sharp crossover lines. The crossover between the bound bipolarons and the unbound pair of polarons has been investigated with a general HTEA bipolaron form (including the latest) but no significant hybridization has been found between these two kinds of states so that we can not draw a conclusion about the nature of this transition.
The triple point is a special point of the phase diagram where the bipolarons (S0), (S1) and (QS) are degenerate and where we should expect a higher energy gain by hybridization when α is not too small. Unfortunately, this triple point disappears when α increases beyond approximately 0.002. When it just disappears the TEA bipolaron binding energy referred to two unbound HTEA polaron is just zero but then its tunnelling energy T b is maximum (but still only 10 −7 × t : that is the bipolaron effective mass is seven order of magnitude larger than those of the bare electron).
The negative conclusion of this section is that more sophisticated variational calculations does not confirm the conclusion of section (3) which was based on the assumption α small extrapolated to larger α.
The present study also shows that in the domain of small U one may have a quantum bipolaron ground-state with a large tunneling energy occurring at very large α > 0.1. This result is simply obtained with only the TEA of the small bipolaron (S0) that is proved to have a nonnegligible binding energy for both t and U small enough. Nevertheless, this result is not relevant for such large α, our approach based on a perturbative theory of the adiabatic limit fails because of too large quantum lattice fluctuations.

Phonon Dispersion Effect
We intend to show that highly degenerate point that could persist under large quantum lattice fluctuations implies very light bipolarons. To achieve that goald, a simple procedure consists in changing the model so as to favor the bipolaron (QS). If we could make it more robust to quantum lattice fluctuations it should become very light for reasonably large α by hybridization with the other degenerate bipolarons at the triple point.
We choose to introduce a phonon dispersion, but this might not be the unique way. When an electron is present at a given site it will also distort the lattice at the neighboring sites. If the sign of the dispersion is appropriate, the lattice potential at the neighboring sites is lower which favors its occupancy by electrons and thus the spatial extension of the bipolaron. The bipolaron (QS) which is more extended than the bipolaron (S0) should be favored. 3 We consider the new Hamiltonian where H is the Holstein-Hubbard Hamiltonian (1) and its reduced Hamiltonian corresponds to H (5) that gives When the coupling C is positive the dispersive term generates an effective attractive interaction between polarons. This coupling cannot exceed the value 1/2 beyond which the low wavevector phonons becomes unstable.

Adiabatic Limit
At the adiabatic limit the equation eq.15 becomes where D is the matrix : Bipolarons (S0),(S1), (QS)... which were found at the anti-integrable limit of the Holstein-Hubbard model at t = 0 persist as ground-states in this model with nonzero coupling C [1] (see the diagrams figs.15 and 17). The domain where bipolaron (QS) is the ground-state enlarges when C increases up to its maximum value 1/2. As expected the existence of a positive dispersion favors the quadrisinglet ground-state.
The first-order transition between bipolaron (S0) and (QS) becomes almost second-order and difficult to distinguish numerically since there is no symmetry breaking 3 Phonon dispersion may induce other important effects in the bipolaron structure as shown in [19] for CDW's. between these two bipolaronic states. Then as expected, there is a soft internal mode which almost vanishes at the transition on both side of the transition which corresponds to a breathing mode of the bipolaron with the same symmetry. Simultaneously the Peierls-Nabarro barrier almost vanishes. This soft mode which does not break the bipolaron symmetry is not a pinning mode and does not favor the classical mobility of this bipolaron. To that purpose the pinning mode which also softens at the first order transition between (QS) and (S1) is the most appropriate ( see fig.16 and refs. [20], [1]).  When the coupling C is too large 0.2 < C < 0.5 our accuracy is limited in practice because of the bipolaron (QS) extension, which requires large system sizes we cannot afford. This problem occurs especially close to the firstorder transition between bipolaron (QS) and the extended state (see diagram fig.17). However the exponential ansatz eq.19,20,21 still fits quite well the bipolaron ground-state as shown on the diag.15. These variational forms allow us to compute quickly the bipolarons even for large C values and to determine approximately the ground-state with a reasonable accuracy (see diagrams figs.15 and 17).
As we already know the flaw of this approximate method is that spurious first-order transitions may occur. This situation happens nearby the first-order transition between (QS) and the extended state as seen on diag.17. It is due to the exponential ansatz which does not provide a good fit of the bipolaron when it becomes more extended.
However, at the triple point the bipolaron ground-state is still localized on very few sites ( fig.18) for C = 0.3 and the exponential ansatz remains sufficiently accurate.

Quantum Corrections
Same methods, as those used above for the original Holstein-Hubbard model are applied to deal with the quantum lattice fluctuations of the modified model. The degeneracy due to the translation invariance of the model is lifted according to standard perturbation theory. One gets a tight binding model as in section 3 which yields both binding and tunnelling energies of the quantum groundstate. Figs.19) shows these quantities for a strong coupling (C = 0.3).
The binding energy of the bipolaron refers to two noninteracting polarons calculated with the HTEA method, which is the most accurate. For a single polaron condition (40) becomes and we choose to write the displacement as v P = D −1 v where v is given by For a large enough phonon coupling C > 0.2, in the region we investigate t < 0.1 the HTEA method only requires the hybridization between a small polaron and a large polaron. The almost second-order transition displayed by the TEA at C = 0 occurs now at a larger t 2 p (α). The binding energy of the quantum bipolaron is still large in that region and one notices the optimal regime where both tunnelling and binding energies have the same value.
To obtain the optimal region, a fine tuning of the parameters is required because changing them slightly can either reduce the binding energy so that the bipolaron becomes fragile against temperature or sharply increase its effective mass, killing its quantum mobility.
Phonon dispersion favors the mobility of the bipolaron because it extends the lattice distortion around the bipolaron (see fig.18) as well as the electronic wave function. Classically, this effect is manifested by internal mode softening and by the depression of the Peierls-Nabarro energy barrier (not calculated here see paper I [1]) between the different bipolarons. As a result, when the lattice is quantum the hybridization between the different bipolarons is increased, which increases the band width and decreases the effective mass.
In the vicinity of the (QS) region (see figs. 19) the effective mass of the HTEA for the bipolaron is about five times larger than the effective mass computed with the perturbative method, but the bipolaron mass is still very small. The comparison of the binding energy calculated with the two methods shows that the variational HTEA method is not accurate in the area of the QS region. Indeed the perturbative method gives a stronger binding energy and thus it is variationally better. This is likely due to the fact that when the bipolaron extends too much the TEA is not accurate because the bipolaron shape is not well approximated by the exponential. Fig.19 shows for α = 0.017 the effective mass of the bipolaron in the optimal regime that ranges not far from 100 bare electronic mass. We choose as an example the realistic optical phonon frequencyhω 0 = 1.10 −1 eV and to be in the optimal regime C = 0.3, α = 0.017 U = 0.25 t = 0.1 the initial parameter of Hamiltonian (1) must be g = 3.10 −2 eV E 0 = 6eV υ = 1.5eV t = 0.6eV c = 0.3eV . The tunnelling energy as well as the bipolaron binding energy are about 6.10 −3 eV . With such characteristic values and a bipolaron concentration not too large, a superfluid state could be expected at relatively high temperatures, perhaps few hundred degrees K. This estimate neglects the bipolaron interactions, but when their concentration becomes large these interactions cannot be neglected, especially at half filling when there is one polaron per site. Close to this close packing regime the bipolaronic structure cannot exist anymore for sure. Instead, a magnetic spatially ordered polaronic structures could occur. Further studies should investigate the situation with large electron densities.

Concluding Remarks
In some circumstances the bipolaron might become unusually light, which allows in principle the formation of superconducting states at rather high temperature with physically realistic parameters. This effect is due to the degeneracy of several bipolaronic states in the adiabatic limit for some specific regions of the phase diagram. In this situation there are small Peierls Nabarro barriers and phonon softening for the different bipolaronic states. Then the quantum lattice fluctuations lift the degeneracy between the degenerate states and may yield very light hybridized bipolarons, which however are well-bound.
We realized this situation in a modified Holstein-Hubbard model, which involves both an electron-phonon interaction and a direct repulsive electron-electron interaction.
The superconducting state of such very light bipolarons occurs for weak concentrations. When the concentration becomes larger there are strong interactions between the bipolarons, which may both break them into polarons and organize different structures (for example, magnetic).
This situation may happen in superconducting cuprates. In the undoped regime where the band of electrons is half filled, the structure can be viewed as close-packed polarons with an antiferromagnetic ordering. This polaron structure should persist for low doping till a certain electron concentration where the holes are polaron vacancies. For a sufficiently large doping the electron concentration may become low enough in order that a (first order) transition toward a superfluid of light quantum bipolarons takes place. The real phenomenology should be more complex because one should expect that the model parameters depend on the doping and thus that the system does not remain always close from the optimal regime with strongly bonded light bipolarons but move around this point. Otherwise, we suggested in [1] that in some appropriate models the (QS) bipolaron could have a d-symmetry. We have not yet realized an explicit model where such an effect occurs, but we hope to.
The numerical techniques we used (Toyozawa Exponential Ansatz) and its improvement (HTEA) where the hybridization between different states is taken into account, turned out to be very efficient to study the bipolaron mass. It should be developed to consider models with many electrons. In [21] it was proven that at adiabatic limit, the ground-state at large electron-phonon coupling was bipolaronic. Variation of the exponential ansatz may provide strong simplifications for these case and a qualitative understanding of the many-polaron problem first in the adiabatic limit, next with quantum lattice fluctuations. Finally, the problem of quantization of discrete breathers can be approached with similar techniques [22].