Isotope effects in high-Tc cuprate superconductors: Ultimate proof for bipolaron theory of superconductivity

Developing a theory of high-temperature superconductivity in copper oxides is one of the outstanding problems in physics. Twenty-five years after its discovery, no consensus on the microscopic theory has been reached despite tremendous theoretical and experimental efforts. Attempts to understand this problem are hindered by the subtle interplay among a few mechanisms and the presence of several nearly degenerate and competing phases in these systems. Here we provide unified parameter-free explanation of the observed oxygen-isotope effects on the critical temperature, the magnetic-field penetration depth, and on the normal-state pseudogap for underdoped cuprate superconductors within the framework of the bipolaron theory compatible with the strong Coulomb and Froehlich interactions, and with many other independent observations in these highly polarizable doped insulators. Remarkably, we also quantitatively explain measured critical temperatures and magnitudes of the magnetic-field penetration depth. The present work thus represents an ultimate proof of the bipolaron theory of high-temperature superconductivity, which takes into account essential Coulomb and electron-phonon interactions.

In 1911 curiosity concerning the electrical properties of metals at low temperatures led the Dutch physicists, Kamerling Onnes and his assistant G. Holst to discover superconductivity at 4.2 K in mercury [2]. This discovery was one of the most important experimental findings in low temperature physics. On the long way towards a microscopic understanding of superconductivity, the observation of an isotope effect on the critical temperature, T c in 1950 [3,4] gave an important clue to the microscopic mechanism of superconductivity. The presence of an isotope effect thus implies that superconductivity is not of purely electronic origin. In the same year Fröhlich [5] pointed out that the electron-phonon interaction gave rise to an attractive interaction between electrons, which might be responsible for superconductivity. Fröhlich's theory played a decisive role in establishing the correct mechanism. Finally, in 1957, Bardeen, Cooper and Schrieffer [7] (BCS) developed the BCS theory that was the first successful microscopic theory of superconductivity. The BCS theory implies an isotopemass dependence of T c , with an isotope-effect exponent α = −d ln T c /d ln M = 1/2, in excellent agreement with the reported isotope exponents in simple metallic superconductors like Hg, Sn and Pb.
The doping dependent oxygen-isotope effect (OIE) on the critical temperature T c , α O = −d ln T c /d ln M O (where M O is the oxygen-isotope mass) [8] and the substantial OIE on the in-plane supercarrier mass m * * ab , α O m * = dm * * ab /d ln M O (ref. [9][10][11][12][13][14]), provide direct evidence for a significant electron-phonon interaction (EPI) also in high-temperature cuprate superconductors. High resolution angle-resolved photoemission spectroscopy (ARPES) [15] provides further evidence for the strong EPI apparently with c-axis-polarised optical phonons. These results along with optical [16], neutron scattering [17,18], and tunneling data [19][20][21] unambiguously show that lattice vibrations play a significant but unconventional role in high-temperature superconductivity. The interpretation of the optical spectra of high-T c materials as due to multi-polaron absorption [22] strengthens the view [23] that the Fröhlich EPI is important in those structures. Operating together with a shorter-range deformation potential and molecular-type (e.g., Jahn-Teller [24]) EPIs, the Fröhlich EPI can readily overcome the Coulomb repulsion at a short distance of about the lattice constant for electrons to form realspace intersite bipolarons (see [25] and Supplementary  Information [26]).
Despite all these remarkable and well-done experiments that lead to the consistent conclusion about the important role of EPI in high-temperature superconductors, there is no consensus on the microscopic origin of the observed unconventional isotope effects on the inplane magnetic-field penetration depth and the normalstate pseudogap. The doping dependent α O has been explained as due to the doping independent oxygenisotope effect on the in-plane carrier concentration n, that is, α O n = −d ln n/d ln M O = 0.146 (ref. [27]). This interpretation contradicts other independent experiments [11,12,14] which consistently show that the carrier concentrations of the two oxygen-isotope samples are the same within 0.0004 per Cu. This is also in sharp contrast to the observed very large oxygen-isotope effect on the low-temperature magnetic-field penetration depth, λ ab ∝ (m * * ab /n) 1/2 , in both La to α O n ≃ 2 if one would assume that the supercarrier mass is independent of the oxygen mass. Another model based on the pair-breaking effects due to impurities, disorder, and/or pseudogap can also explain the observed oxygenisotope effects on the penetration depth and the critical temperature in deeply underdoped samples [29]. But this model cannot consistently explain the negligibly small α O but the large OIE on the penetration depth in optimally doped samples [13].
Alternatively, the bipolaron theory of superconductivity [25] can naturally account for the substantial α O m * and large α O in deeply underdoped cuprates [30]. There is a qualitative difference between ordinary metals and polaronic conductors. The renormalized effective mass of electrons is independent of the ion mass M in ordinary metals (where the Migdal adiabatic approximation is believed to be valid), because the EPI coupling constant λ does not depend on the isotope mass (see, for instance ref. [26]). However, when electrons form polarons (new quasiparticles dressed by lattice distortions), their effective mass m * depends on M [30].
Although the bipolaron theory can qualitatively explain both α O m * and α O in deeply underdoped samples, some important issues have not been well addressed by the theory. The first issue is why α O m * is not equal to α O even in deeply underdoped cuprates. The second issue is why α O is much smaller than α O m * for slightly underdoped samples. The third issue is why there is a giant negative oxygen-isotope effect on the pesudogap formation temperature T * in HoBa 2 Cu 4 O 8 (ref. [31]). Here we provide parameter-free explanations to the observed oxygen-isotope effects on the critical temperature, the inplane supercarrier mass, and on the normal-state pseudogap in HoBa 2 Cu 4 O 8 within the framework of the bipolaron theory. The present work thus represents an ultimate proof of the bipolaron theory of high-temperature superconductivity.
We adopt the bipolaronic low-energy excitation structure of cuprates, Fig. 1, derived from the microscopic Hamiltonian with the strong Coulomb and Fröhlich interactions [25,26]. The critical temperature of quasi-twodimensional (2D) bipolarons, which are hard-core bosons, depends on their density n b (T c ) at the critical temperature and the in-plane bipolaron mass m * * ab : The bipolaron density slightly depends on temperature due to bipolaron depletion into unbound single polarons, n b (T ) = [x − n p (T )]/2, with the density n p given by where W is the polaron half-bandwidth, x is the in-plane doping level, and ∆ is the bipolaron binding energy, Fig. 1. This expression is obtained by integrating the Fermi-Dirac distribution function with a constant (2D) density of states in the polaron band, N (E) = 1/2W and assuming that the polaron half-bandwidth is large Here we take into account that the chemical potential is zero in the superconducting state at and below T c if all energies are taking with respect to the bipolaron ground state. The polaron in-plane effective mass, m * ab ∝ exp(A √ M ) (A is a constant), the polaron inverse bandwidth, 1/W , and the inter-site bipolaron mass m * * ab have the same isotope exponent d ln m * * ab / ln M ≡ α m * = (1/2) ln(m * ab /m) [30], where m is the band mass in a rigid lattice. The isotope effect on the pseudogap is given by: and The above expressions are obtained by taking into account that the pseudogap in Fig. 1 is given by ∆/2 = J p /2−(W −t/2) = J p /2−3W/4, where J p is the phononinduced intersite attraction, which is independent of the ionic mass [26], and t ≈ W/2 in the intermediate coupling regime [13]. Then using Eqs. (1,2) and neglecting the terms on the order of k B T c /W ≪ 1, one readily obtains the ratio: Eq. 6 can naturally explain why α O m * is always larger than α O [11,12]. It is worth noting that Eq. 6 is valid only if d ln ∆/d ln M is small. When d ln ∆/d ln M is large, we need to use Eqs. 1 and 2 to directly calculate the T c and n p (T c ) changes upon the isotope exchange, that is, Using the above equations we can quantitatively explain the oxygen-isotope effects on the pesudogap and the critical temperature in slightly underdoped HoBa 2 Cu 4 O 8 (ref. [31]). The oxygen isotope effect on the relaxation rate of crystal-field excitations in this compound was investigated by means of inelastic neutron scattering [31]. The relaxation rate, which is related to the free-carrier spin density, clearly shows a large oxygen-isotope effect (see Fig. 2). For the 16 O sample there is evidence for the opening of an electronic gap in the normal state at T * ≃ 170 K while for the 18 O sample T * is shifted to about 220 K. In contrast, the T c is shifted from 79.0 to 78.5 K upon replacing 16 In a normal state with no pseudogap, the relaxation is the exchange integral between the 4f electrons of the Ho 3+ ions and the charge carriers and N (E F ) is the electronic density of states at the Fermi energy [31]. Within the polaron/bipolaron framework, N (E F ) = 1/2W ∝ m * ab ∝ m * * ab , so the oxygen-isotope effect on N (E F ) is the same as the oxygen isotope effect on m * * ab . For a slightly underdoped YBa 2 Cu 3 O 7−y film, δm * * ab /m * * ab was found to be 5.5% upon replacing 16 18 O concentration is about 95%) [14]. If we assume that the oxygen-isotope effect on m * * ab for HoBa 2 Cu 4 O 8 is similar to that for the slightly underdoped YBa 2 Cu 3 O 7−y film, we expect that δm * * ab /m * * ab = 4.3% in HoBa 2 Cu 4 O 8 . In the superconducting state or below T * , the relaxation rate is suppressed due to openning of the gap. Then the relaxation rate is given by [31]: The solid lines in Fig. 2 represent the best fits of Eq. 8 to the data below T * . The best fits yield ∆/k B = 94.6±8.8 K for the 16 O sample and ∆/k B = 286±19 K for the 18 O sample. Therefore, there is a giant oxygenisotope effect on the pseudogap, in agreement with Eq. 4. Substituting δ∆ = 16.5 meV and δm * * ab /m * * ab = 4.3% into Eq. 4, we find W = 0.256 eV and t = 0.128 eV. Using m * * ab = 2h 2 /ta 2 , we calculate m * * ab = 8.1 m e , which is very close to that (8.3 m e ) inferred for the slightly underdoped YBa 2 Cu 3 O 6.88 with T c = 87.9 K (see ref. [33]).  [39].
In the absence of charge localization, as the case of the stoichiometric HoBa 2 Cu 4 O 8 , the Bose-Einstein condensation temperature T c is given by [16] where t c is related to the out-of-plane bipolaron mass m * * c as t c =h 2 /2m * * c d 2 . The out-of-plane bipolaron mass m * * c is deduced to be 518 m e using m * * c /m * * ab = 64 at T c (ref. [36]) and m * * ab = 8.1 m e . Then we calculate t c = 0.160 meV. Substituting t = 0.128 eV, t c = 0.160 meV, and n b (T c ) = 0.1077 into the above equation, we calculate T c = 78.9 K, in quantitative agreement with the measured T c (79.0 K). We further show that (see Supplementary Information [26]) Eq. 20 can also quantitatively explain the underdoped La 1.90 Sr 0.10 CuO 4 with T c = 29 K but overestimates the T c 's of some optimally doped samples. This implies that optimally doped cuprates are in the crossover regime from bipolaronic real-space pairing to the Cooper pairing of polarons [25].
Apart from the striking isotope effects explained quantitatively here there is abundant independent evidence in favor of bipolarons and the Bose-Einstein condensation in underdoped cuprate superconductors. In particular, the parameter-free estimates of the Fermi energy using the magnetic-field penetration depth [41] and the magnetic quantum oscillations [42] yielded a very low value (below 50 meV) supporting the real-space pairing in underdoped cuprate superconductors. Magnetotransport and thermal magnetotransport data strongly support preformed bosons in cuprates. In particular, many high-magneticfield studies revealed the non-BCS upward curvature of the upper critical field H c2 (T ) as a function of temperature [43] as predicted for the Bose-Einstein condensation of charged bosons in the magnetic field [44]. The Lorenz number differs significantly from the conventional Sommerfeld value of the standard Fermi-liquid theory because the carriers are double charged bosons [45]. Direct measurements of the Lorenz number using the thermal Hall effect just above T c [46] produce its value, which is about the same as that predicted by the bipolaron model. The unusual normal-state diamagnetism uncovered by torque magnetometery has been convincingly explained as the normal state (Landau) diamagnetism of charged bosons [47]. Single polarons, localised within an impurity bandtail, coexist with bipolarons in the charge-transfer doped Mott-Hubbard insulator, where the chemical potential is pinned within the charge-transfer gap due to bipolaron formation. This band-tail model accounts for two energy scales in ARPES and in the extrinsic and intrinsic tunnelling, their temperature and doping dependence, and for the asymmetry and inhomogeneity of extrinsic tunnelling spectra of cuprates [48].

Microscopic derivation of the low-energy bandstructure
In highly polarizable ionic lattices like cuprate superconductors both the Coulomb repulsion and the Fröhlich electron-phonon interaction (EPI) are quite strong (of the order of 1 eV) compared with the low Fermi energy of doped carriers because of a poor screening by non-or near-adiabatic carriers [1]. In those conditions the BCS-Eliashberg theory [2] breaks down because of the polaronic collapse of the electron bandwidth [3] so that one has to apply a non-adiabatic small polaron theory [4].
Here we sketch the microscopic derivation of the lowenergy band structure, Fig. 1, using the analytical multipolaron theory in the strong-coupling regime for highly polarisable lattices (more details are found in Refs. [4]).
Quantitative calculations of the interaction matrix elements can be performed from pseudopotentials using the density functional theory (DFT) [5]. On the other hand, one can express the bare Coulomb repulsion and EPI through material parameters rather than computing them from first principles in many physically important cases [6]. In particular, for a polar coupling to longitudinal optical phonons (the Fröhlich EPI), which is the major EPI in polar crystals, both the momentum depen-dence of the matrix element, M (q), and its magnitude are well known, |M (q)| = γ(q))hω 0 / √ 2N with a dimensionless γ(q) = 4πe 2 /κΩhω 0 q 2 , where Ω is a unit cell volume, N is the number of unit cells in a crystal, ω 0 is the optical phonon frequency, and κ = ǫ ∞ ǫ 0 /(ǫ 0 − ǫ ∞ ). The high-frequency, ǫ ∞ and the static, ǫ 0 dielectric constants are both measurable in a parent polar insulator.
The dielectric response function of strongly correlated electrons is apriori unknown. Hence one has to start with a generic Hamiltonian including unscreened Coulomb and Fröhlich interactions operating on the same scale since any ad-hoc assumption on their range and relative magnitude might fail, Here T ij ≡ T (m − n) is the bare hopping integral, µ is the chemical potential, i = m, s and j = n, s ′ include both site (m, n) and spin (s, s ′ ) states, u(m, q) = (2N ) −1/2 γ(q) exp(iq · m), c i , d q are electron and phonon operators, respectively,n i = c † i c i is a site occupation operator, and H ph = qh ω 0 (d † q d q + 1/2) is the polar vibration energy.
In highly polarisable lattices with ǫ 0 → ∞ the familiar Lang-Firsov (LF) [7] whereσ ij = T (m − n)X † iX j is the renormalised hopping integral involving the multi-phonon transitions described withX i = exp q u(m, q)d q − H.c. , andμ = µ+E p is the chemical potential shifted by the polaron level shift, Here, the integration goes over the Brillouin zone (BZ) and E p = 0.647 eV in La 2 CuO 4 [1]. The electron-phonon coupling constant is defined as λ = 2E p N (0). In the case of 2D carriers with a constant bare density of states, N (0) = ma 2 /2πh 2 per spin, Eq.(12) places cuprates in the intermediate to strong-coupling regime, λ > ∼ 0.5, if the bare band mass m > m e (here a is the in-plane lattice constant).
The number of virtual phonons in the polaron cloud is large in oxides and some other polar lattices, E p /hω 0 > 1 with the characteristic (oxygen) optical phonon frequencyhω 0 < ∼ 80 meV, so that multi-phonon vertexes are essential in the expansion of the hopping operator σ ij . To deal with this challenging problem let us single out the coherent hopping in Eq.(11) averagingσ ij with respect to the phonon vacuum, and consider the remaining terms as perturbation,H = H 0 + H p−ph . Here describes free phonons and polarons coherently propagating in a narrow band with the exponentially diminished hopping integral, and is the residual polaron-multiphonon interaction, which is a perturbation at large λ. In the diagrammatic technique the corresponding vertexes have any number of phonon lines. The second-order in H p−ph polaron selfenergy (Σ p ≈ −E p /2zλ 2 ) and the phonon self-energy (Σ ph ≈ −xhω 0 /zλ 2 ) are small, if λ ≫ 1/ √ 2z [8](here z is the lattice coordination number and x is the atomic density of carriers). Hence the perturbation expansion in 1/λ is applied. Importantly there is no structural instability in the strong coupling regime since |Σ ph | ≪hω 0 [8].
The LF transformation, Eq. (11) is exact for any adiabatic ratiohω 0 /T (a). However, if the perturbation expansion in 1/λ is restricted by lowest orders, then it significantly overestimates polaron masses in the adiabatic regime,hω 0 /T (a) < 1, for the case of the short-range (Holstein) EPI (here T (a) is the nearest-neighbor bare hopping integral). The polaronic band narrowing factor, exp(−g 2 ) becomes very small for this EPI in the strongcoupling regime, which would eliminate any possibility of high temperature superconductivity and even metallicity of the small Hosltein polarons.
However in the case of the long-range (Fröhlich) EPI, Quantum Monte-Carlo simulations [9] show that the LF transformation provides numerically accurate polaron masses already in the zero order of the inverse-coupling expansion both in the adiabatic regime as well as in the non-adiabatic one for any strength of the Fröhlich EPI. Moreover, such small polarons [9] and small bipolarons [10] are perfectly mobile in the relevant range of the coupling and the adiabatic ratio.
The perturbation H p−ph has no diagonal matrix elements with respect to phonon occupation numbers. Hence it can be removed from the Hamiltonian in the first order using a second canonical transforma- where E n , E n ′ and |n , |n ′ are the energy levels and the eigenstates of H 0 , respectively. Taking into account that the polaron Fermi energy is small compared with the phonon energy at strong coupling and/or sufficiently low doping [1], one can neglect the polaron contribution to E n ′ − E n ≈hω 0 q n ′ q and project the second-order in 1/λ Hamiltonian H onto the phonon vacuum |0 with the following result where is the Heisenberg multi-phonon operator obtained by replacing d q inX † i with d q exp(iω 0 t). Calculating the integral, Eq.(17) with δ → +0 yields All matrix elements, Eq. (18), of the polaron-polaron interaction are small compared with the polaron kinetic energy except the exchange interaction, J p (m − n) ≡ V nm mn such that f (mn, m ′ n ′ ) = 2g 2 (m − n). Using ∞ k=1 y k /k!k = −C − ln(y) + Ei ⋆ (y) with C ≈ 0.577 and Ei ⋆ (y) ≈ e y /y (for large y) one obtains a substantial J p (m) = T 2 (m)/2g 2 (m)hω 0 , which is larger than the nearest-neighbour polaron hopping integral, t(a)/J p ∝ 2hω 0 g 2 e −g 2 /T (a) < 1. Keeping only this exchange we finally arrive with the polaronic "t-J p " Hamiltonian [4], where S m = (1/2) s,s ′ c † ms τ ss ′ c ms ′ is the spin 1/2 operator ( τ are the Pauli matrices),n m = sn i , and µ =μ + m J p (m) is the chemical potential further renormalized by H p−ph .
There is a striking difference between this polaronic t-J p Hamiltonian and the familiar t-J model derived from the repulsive Hubbard U Hamiltonian in the limit U ≫ t omitting the so-called three-site hoppings and EPI [11]. The latter model acts in a projected Hilbert space constrained to no double occupancy. On the contrary in the polaronic t-J p Hamiltonian, Eq. (19) there is no constraint on the double on-site occupancy since the Coulomb repulsion is negated by the Fröhlich EPI. The polaronic hopping integral t(a) leads to the coherent (bi)polaron band and the antiferromagnetic exchange of purely phononic origin J p bounds polarons into small superlight inter-site bipolarons. Last but not least the difference is in the "+" sign in the last term of Eq. (19) proportional ton mnn , which protects the ground superconducting state from the bipolaron clustering, in contrast with the "-" sign in the similar term of the standard t-J model, where the phase separation is expected at sufficiently large J [12].
The polaronic t-J p Hamiltonian, Eq. (19) is analytically solvable in the limit of sufficiently low atomic density of carriers [4]. Neglecting the first term in H, which is the polaron kinetic energy proportional to t(a) < J p , one can readily diagonalise the remaining spin-exchange part of the Hamiltonian. Its ground state is an ensemble of inter-site singlet bipolarons with the binding energy ∆ b = J p localised on nearest neighbor sites. Such small bipolarons repel each other and single polarons via a short-range repulsion of about J p .
The kinetic energy operator in Eq. (19) connects singlet configurations in the first and higher orders with respect to the polaronic hopping integrals. Taking into account only the lowest-energy degenerate singlet configurations and discarding all other configurations one can project the t-J p Hamiltonian onto the inter-site bipolaronic Hamiltonian using the bipolaron annihilation operators B m = 2 −1/2 (c m↑ c m+a↓ − c m↓ c m+a↑ ), where a connects nearest neighbors [10]. Such inter-site bipolarons are perfectly mobile since they tunnel via single-polaron transitions [10,13]. At finite temperatures single polarons, thermally excited above the pseudogap, coexists with these bipolarons as shown in Fig. 1 of the main text.
Small bipolarons are hard-core bosons with the shortrange repulsion and a huge anisotropy of their effective mass since their inter-plane hopping is possible only in the second order of the polaron hopping integral [14]. The occurrence of superconductivity in bipolaronic systems is not controlled by the pairing strength, but by the phase coherence among the electron pairs below the Bose-Einstein condensation temperature [15]. 2 Parameter-free calculations of the Bose-Einstein consensation temperatures of some cuprate superconductors In the absence of charge localization, the Bose-Einstein condensation temperature T c is given by [16] where t is the bipolaron half-bandwidth and t c is related to the out-of-plane bipolaron mass m * * c as t c = h 2 /2m * * c d 2 (d is the inter-plane distance). The above equation can be written in terms of the measurable parameters such as the inplane penetration depth and the supercarrier mass anisotropy constant γ 2 = m * * c /m * * ab ,