Electron pairing: from metastable electron pair to bipolaron

Starting from the shell structure in atoms and the significant correlation within electron pairs, we distinguish the exchange-correlation effects between two electrons of opposite spins occupying the same orbital from the average correlation among many electrons in a crystal. In the periodic potential of the crystal with lattice constant larger than the effective Bohr radius of the valence electrons, these correlated electron pairs can form a metastable energy band above the corresponding single-electron band separated by an energy gap. In order to determine if these metastable electron pairs can be stabilized, we calculate the many-electron exchange-correlation renormalization and the polaron correction to the two-band system with single electrons and electron pairs. We find that the electron-phonon interaction is essential to counterbalance the Coulomb repulsion and to stabilize the electron pairs. The interplay of the electron-electron and electron-phonon interactions, manifested in the exchange-correlation energies, polaron effects, and screening, is responsible for the formation of electron pairs (bipolarons) that are located on the Fermi surface of the single-electron band.


I. INTRODUCTION
Since the discovery of high T c superconductivity by Bednorz and Müller [1] in 1986, great progresses have been made in the experimental and theoretical investigation of unconventional superconductivity. However, the mechanism of electron pairing in unconventional superconductors remains one of the most challenging and unresolved problems in condensed matter physics. [2][3][4] Vast experimental evidences have shown that electron pairing and unconventional superconductivity occur in many different materials, such as cuprates, [1][2][3][4] iron-based superconductors, [5][6][7][8] and carbon-based superconductors, [9,10] etc. Although there are many different theories for unconventional superconductivity, almost all theories follow the basic idea of the BCS theory [11]. They presume that there is some effective attraction between electrons leading to Cooper pairing which spontaneously condense into a collective non-Fermi liquid state.
We would like to mention some very recent experimental results related to the electron-pairing mechanism in unconventional superconductors. Božović et al. reported very impressive and accurate results on superconductivity in high-T c cuprates. [12] They synthesized atomically perfect thin films and multilayers of cuprates La 2x Sr x CuO 4 (LSCO) and measured the absolute value of the magnetic penetration depth and the phase stiffness with high accuracy in thousands of samples. The large statistics revealed clear trends in the intrinsic proper- * Corresponding e-mail: hai@ifsc.usp.br ties of the cuprate superconductors. They found that the obtained results disagree with the BCS theory in any variant, i.e. clean or dirty, including the Migdal-Eliashberg theory. Rather, the experimental data indicated small (local) and very light electron pairs with mass on the order of an electron mass. These pairs are preformed well above T c and at T c undergo Bose-Einstein condensation. [12,13] Investigations performed by Zhong et al. [14] and by Ren et al. [15] challenged the d-wave pairing mechanism in cuprates. With scanning tunneling spectroscopy, they revealed anisotropic and nodeless superconducting gaps in the cuprate superconductors Bi 2 Sr 2 CaCu 2 O 8+δ (Bi2212) and YBa 2 Cu 3 O 7−x (YBCO). In their paper, Zhong et al. [14] affirmed that "this is contradictory to the nodal d-wave pairing scenario that is often thought to be the most important result in the 30-year study of the HTS mechanism of cuprates". Important progresses in recent investigations on the electron pairing mechanism in iron-based superconductors indicate small and preformed Cooper pairs. [6][7][8] For instance, using Bogoliubov quasiparticle interference imaging, Sprau et al. [6] found that the superconducting energy gap in FeSe is extremely anisotropic and nodeless. Their investigation discovered the existence of orbital-selective Cooper pairing in FeSe. Gerber et al. [8] combined two time-domain experiments into a coherent lock-in measurement in the terahertz regime and was able to quantify the electron-phonon coupling strength in FeSe. Their study revealed a strong enhancement of the electron-phonon coupling strength in FeSe owing to electron correlations and highlighted the importance of the cooperative interplay between electron-electron and electron-phonon interactions In this paper, we present a theory for electron pairing in crystals where we consider the electron-electron correlation, the periodic potential of the crystal lattice, and the electron-phonon interaction. Our theory is different from all previous ones. We obtain preformed small electron pairs in the periodic potential of the crystal. This is essentially an orbital dependent electron-pairing theory in which the exchange-correlation between two electrons occupying the same orbital is decisive for pair formation. However such pairs are metastable unless the electronphonon interaction is included. The calculations show that the electron-phonon coupling and the polaron effects are responsible for the stabilization of the electron pairs. Our basic idea for the electron pairing process is that the electron-electron correlation is orbital dependent. Within the jellium model for the electron gas in solids, [16] the atomic nuclei that form the periodic lattice are smeared out into a uniform positive charge distribution. Each electron is totally delocalized. Therefore, many electrons "see" each other with their fluctuation potential at the same time and thus correlate all at once, giving rise to collective screening and oscillation effects. However, in atoms and molecules, significant correlations occur within electron pairs. [17][18][19] Strong exchange-correlation interaction between two electrons in the same orbital manifests in the shell structure of atoms and also in the covalent and ionic bonding in molecules. Our starting point in this study is to distinguish the exchange-correlation effects between two electrons of opposite spins occupying the same atomic orbital from the average correlation among many electrons in a crystal. This may happen in a crystal but the electrons have to "feel" the nuclei potential well. This leads to a preliminary condition that the effective Bohr radius of the valence electrons in the crystal has to be comparable or smaller than the lattice constant. For instance, for a cuprate crystal with effective electron mass m * ≃ 5m 0 and static dielectric constant ǫ 0 ≃ 30, the effective Bohr radius a B ≃ 3.2Å is smaller than the lattice constant of about 3.8Å. Because we want to show that the electronpair correlation in atoms can manifest themselves in electron transport in crystals, our calculations have to start first with the formation of energy bands. [20] In order to find out the electron-pair states in the crystal, we will first establish a simple crystal model to discuss the physical process. We consider a "hydrogen solid" model with single-electron state of H atom and electron-pair state of the H − ion. We will show that, besides the energy bands from the single-electron energy levels of individual atoms in the crystal, there can exist a metastable electron-pair energy band from the correlated electron pairs of the H − state for the lattice constant λ being larger than the effective Bohr radius a B . [21] The electron pairs are metastable because the Coulomb repulsion is strong overwhelming the exchange-correlation. In order to stabilize them we have to include the electronphonon interaction to counterbalance the Coulomb repulsion. Therefore, the electron-phonon interaction is necessitated in a natural way in the electron pairing process.
This paper is organized as follows. In Sec. II we calculate and discuss the metastable electron-pair energy band in two-dimensional square-lattice crystals. In Sec. III we present the many-particle Hamiltonian consisting of electrons in both the single-electron and electron-pair bands being coupled to the LO-phonons. In Sec. IV and Sec. V, we calculate the many-particle exchangecorrelation (xc) energies due to electron-electron interactions and polaron energies due to electron-phonon interactions. Then, we show in Sec. VI the conditions under which the metastable electron pairs can be stabilized in the ground state including many-body effects. Finally, we summarize our work in Sec. VII.
In the calculations, we will use effective Bohr radius a B = ǫ 0h 2 /m * e 2 and effective Rydberg R y =h 2 /2m * a 2 B as the units for length and energy, respectively.

II. METASTABLE ELECTRON PAIRS IN CRYSTAL
Electronic band structure is of fundamental importance for our understanding of many physical properties of solids. Within the independent electron approximation, the electron states are of Bloch form in the periodic potential of the crystal lattice. The effects of electronelectron interactions are accounted for by an effective potential which repeats this periodicity. [20] In this section, we will show that, besides the energy bands from the single-electron states there can exist a metastable electron-pair energy band depending upon the crystal structure and potential. This metastable electron-pair band originates from two correlated electrons of opposite spins occupying the same atomic orbital.

A. Two-electron atoms
Our study will start with the simplest electron-pair system, i.e., two-electron atoms. These helium-like atoms with two electrons of opposite spins occupying the same orbital, e.g. helium atom He and negatively charged hydrogen anion H − have played an important role in the development of theoretical physics in the last century. [22] It is a challenge to determine accurately the correlation energy, even in simple systems such as He atom and H − ion. [22][23][24][25] Hylleraas' result for the ground-state energy of He atom obtained in 1929 was -5.80648 R y [26]. After generations of calculation, very accurate (nonrelativistic) ground-state energies [27][28][29][30][31][32] of two electron atoms have been obtained: -5.807448754068· · · R y for He and -1.055502033088· · ·R y for H − . Recently, using high-precision variational calculations Estienne et al. [33] determined the critical nuclear charge Z c =0.911 028 224 077 255 73(4) which is the minimum charge required to bind two electrons in a helium-like atom.
On the other hand, the famous experiment on twoelectron atoms by Madden and Codling [34] revealed that the simple model based on independent particle picture is inappropriate to characterize a series of doubly excited states because of strong electron-electron correlation. [22,34] In comparison with the single-electron states of the H atom, the H − ion is a closed-shell system with two strongly correlated electrons. Such an electron-pair state is different in its nature from the single-electron states because of the strong correlation. It should be recognized as a new strongly correlated electronic state.
Negative hydrogen ion H − in two-dimensional (2D) system has also been investigated in the last decades mostly because of the discovery of its counterpart D − center in 2D semiconductor quantum wells [35]. The D − center is a negatively charged shallow donor impurity center in semiconductors, such as a negatively charged Si impurity in a GaAs quantum well. It is an H − -like state in solid-state environment but with very different energy and length scales (e.g., in GaAs, the effective Rydberg R y =5.9 meV and effective Bohr radius a B =98Å). Therefore, the D − center in semiconductors is considered as an ideal "laboratory" to study the H − properties, for instance, in high magnetic fields. [36] Earlier variational calculation by Phelps and Bajaj found the energy of the H − in 2D is -4.48 R y . [37] Further numerical calculations obtained -4.48054 R y by Ivanov and Schmelcher [38] and -4.4804798 R y by Ruan et al. [39].
In Table I we compare the ground-state energies of the 2D and 3D H − states. The binding energy E b is defined as the difference between the energies E H of the neutral H atom and E H − of the H − ion. This is the energy required to remove one of the two electrons from the H − ion to infinity. It is also called electron affinity of the hydrogen atom. One sees that the binding energy E b of the H − in 2D is almost 9 times larger than that in 3D because electron correlation in 2D is much stronger. In the last column we also give the energy per electron ε p in the H − state. In order to explain the so-called Mott insulator and metal-insulator transition [40,41], Mott considered a hydrogen solid model with the single-electron energy band only, i.e., a simple cubic lattice crystal of one-electron atoms and made the following discussion [40]. For small values of the lattice constant λ, there is a half-filled band in such a crystal and thus it is metallic. If one varies the lattice constant to large values (but not so large as to prevent tunnelling), the Coulomb interaction U for two electrons occupying the same atomic site overcomes the kinetic energy (characterized by the band width W ). In this case, each electron should be assigned to its parent atom. The crystal must be nearly the same as a collection of isolated neutral hydrogen atoms and thus it is an insulator. This reveals a competition between potential and kinetic effects. At large λ (small W ) the Coulomb repulsion U dominates, the electrons are localized, and the system is insulating. [40,41]. The idea of Mott led to the theoretical model introduced by Hubbard. [42] The Hubbard model traces the insulating behavior to strong Coulomb repulsion between electrons occupying the same orbital. The competition between the kinetic and Coulomb energies gives rise to strong electron-electron correlations. The Hubbard model was proposed originally to describe the transition between conducting and insulating systems. It has also been widely used to study materials with strongly correlated electrons and high-temperature superconductivity. [43] In this paper we will consider both the single-electron and electron-pair states in the "hydrogen solid" model. Our calculation will be performed for 2D systems because we can obtain more accurate numerical results in 2D. Another reason why 2D systems are more interesting is that electron-electron correlations are stronger. Many unconventional superconductor materials are found to be essentially two dimensional. Fig. 1 shows diagrams representing (a) a 2D H atom and (b) a 2D H − ion with their respective energy levels. The nuclear potential V a (r − R m ) of the atom is represented by the black curves, where R m is the position of the nucleus. The single-electron levels of a 2D H atom are given by ε i = −R y /(i + 1/2) 2 (for i=0,1,2,...). The energy level of a correlated electron pair in H − ion, i.e., the energy per electron in the ground state, is given by ε p = E H − /2 = −2.24 R y . We remind that a single H − ion is stable.
We now consider the following "hydrogen solid": N atoms are arranged into a simple-lattice crystal at positions R m (for m = 0, 1, 2, · · · , N − 1) with lattice constant λ of the order of the Bohr radius a B as indicated in Fig. 1(c). The crystal potential for an electron at r j is given by It is known that the single-electron levels ε i of individual atoms form energy bands E i (k) in such a crystal [20] as indicated by the horizontal thick-blue lines in Fig. 1(c). In principle, there is also the possibility that two electrons of opposite spins occupy the same atomic orbital forming an H − -like state in the crystal, but it becomes unstable due to the presence of the neighbor atoms. Therefore, the counterpart of the H − state in a crystal has never been investigated. In this section we will show that, though such an electron pair is unstable in a crystal due to the Coulomb repulsion, they may form a metastable energy band (indicated by the horizontal thick-green line in Fig.  1(c)) depending on the crystal structure and potential.
In such a crystal the lattice constant λ should not be so small as to prevent individual atoms to bind two electrons, but not so large as to prevent co-tunnelling of an electron pair between the neighbor unit cells.
In the center of mass and relative coordinates (R,r) of the two electrons at r 1 and r 2 , defined by R = 1 2 (r 1 + r 2 ) and r = r 1 − r 2 , the wavefunction of an individual electron pair with energy 2ε p bound to the atom at R m is given by φ(R − R m , r). The Schrödinger equation for two electrons in the crystal potential given by Eq. (1) can be written as, where r 1 = R+ 1 2 r and r 2 = R− 1 2 r. The term 2/|r| is the Coulomb repulsion potential between the two electrons. We should bear in mind that, due to electron-electron repulsion, the ground state of this system can be found for |r 1 − r 2 | = |r| → ∞. In other words, the ground state of this two-electron system corresponds to two noninteracting single electrons separated by an infinitely long distance. But we are looking for the quantum states of two-correlated electrons occupying the same orbital in the same unit cell in the crystal with the average separation r being less than the lattice constant λ. Therefore, the electron-pair states in the crystal are metastable. The calculations in Sec. II-C will confirm that such a metastable state does exist in 2D periodic potentials.
If the electron-pair states of two correlated electrons in the crystal can be approximated by a linear combination of the electron-pair wavefunctions φ(R − R m , r) of single atoms, written as, for r < λ, we can obtain the following homogeneous linear equations, with α p (R l = 0) = 1, and We observe that Eq. (5) for c m depends only on R l = R m − R n . This is an eigenvalue problem of a block circulant matrix. [44] The solution has the following form where the vector k should be a reduced wavevector in the first Brillouin zone and C is the normalization constant.
We finally obtain the electron-pair wavefunction in the crystal given by The above wavefunction is a Bloch wavefunction in the center of mass coordinates R because it can be written as where the part in the square brackets is a periodic function in the coordinates R with period of the crystal lattice. The dispersion relation of the electron-pair band is given by Because the considered electron-pair wavefunction φ(R−R m , r) of a two-electron atom has essentially the ssymmetry [24,37], the dispersion relation of the electronpair band in the 2D square-lattice crystal with lattice constant λ can be approximated as (13) where only the nearest-neighbor tunneling term J p (R 1 ) is considered. The value of J p (R 1 ) determines the bandwidth of the electron-pair band. Because this is a cotunneling process of two paired electrons between the neighbor sites and the effective mass of the electron pair is twice of a single electron, the electron-pair bandwidth should be much smaller than that of the single-electron band. This dispersion relation will be confirmed in the next section by making numerical calculations of the metastable electron-pair band in 2D periodic potential of a square lattice.

C. Metastable electron-pair band in 2D square lattices
For a quantitative demonstration of the metastable electron-pair band in a crystal and its renormalization due to many-body effects, we will use the following 2D periodic potential. For an electron at r j = (x j , y j ) in a 2D square lattice with the lattice constant λ, the considered potential is given by where q = 2π/λ and V 0 is the amplitude of the crystal potential. Notice that, V 0 is not a measurable quantity, e.g., the amplitude of the crystal potential in Fig. 1(c) should be infinity. The potential defined in Eq. (14) with two parameters λ and V 0 will simplify our numerical calculations without losing any essential features of the theory. In this 2D periodic potential, the energy has a continuous spectrum for E ≥ 0. Therefore, two electrons can possibly bind into a pair for E < 0 only. We have calculated the single-electron and metastable electron-pair states in this periodic potential. For the calculation details we refer to Ref. 21.
The single-electron states are well known in this potential. The Schrödinger equation for a single electron is given by with where k is the wavevector in the first Brillouin zone, l is the band index, and G l = l x qi + l y qj (with l x , l y = 0, ±1, ±2, ...) is the reciprocal-lattice vector; E l (k) and ψ k+G l (r) are the eigenvalue and eigenfunction, respectively.
When we consider two electrons in this periodic potential, their Hamiltonian is given by where the last term is the Coulomb repulsion potential between the two electrons. In the center of mass R = (X, Y ) and relative r = (x, y) coordinates defined in Eq. (2), the two-electron Hamiltonian becomes This Hamiltonian is periodic in X and Y with period λ. We can choose a Bloch wavefunction in the center-ofmass coordinates for our basis. As to the function in the relative coordinates r = (r, θ), we have to consider the symmetry of the electron-electron Coulomb potential and the periodic potential representing a 2D square lattice. We use the following basis for our wavefunction, and where n = 0, 1, 2, · · ·, m = 0, 1, 2, · · · , n, ξ n = 2/(2n + 1), is the generalized Laguerre polynomial. The function R n,m (r) is taken from the wavefunction of a 2D hydrogen atom [45,46] with a modification introduced by a dimensionless scaling parameter β. The two-electron wavefunction can be written as Ψ k (R, r) = lx,ly n,m a lx,ly;n,m (k)ψ lx,ly;n,m (R, r). (22) Considering the antisymmetry of the electron wavefunctions with spin states, the two-electron wavefunction of the spin singlet state is given by the above expression with the sum over even m only.
Solving the corresponding eigenvalue equation of the two-electron Hamiltonian given by Eq. (18) with the above basis, we find a metastable electron-pair state of spin-singlet in the 2D square lattice potential. As shown in Ref. [21], for fixed period λ, a metastable electronpair state can be found when V 0 is larger than a certain value. A minimum potential V 0 required for a metastable electron-pair state in the periodic potential corresponds to the critical nuclear charge Z c for a helium-like twoelectron atom. The metastable electron-pair state exist for E < 0 only. This indicates that co-tunneling of the paired electrons occurs in the formation of the electronpair band in this 2D periodic potential. In the calculations we found that the average separation r between two electrons in a metastable pair state is always smaller than half the lattice constant, r < λ/2. The global minimum of the eigenenergy of the two-electron system occurs at r → ∞ corresponding to two non-interacting single electrons. We also want to emphasize that the parameter β in the wavefunction in Eq. (20) plays the role of variational parameter to improve the correlation energy of the electron pair. Because this parameter is directly related to the average distance between the two electrons in a pair, it helps us to understand better the metastable electron-pair state. However, the parameter β does not determine the existence of the electron-pair state in the periodic potential.
In the following, we will present numerical results for the band structure and show that the dispersion relation of the electron-pair band in this potential fits very well the expression given by Eq. (13). Fig. 2(a) shows the dispersion relations of the electron states in the 2D crystal potential with λ = 1.3 a B and V 0 = 15 R y . The dispersion relation E p (k) of the spin-singlet metastable electron-pair band is given together with that of the lowest single-electron band E 0 (k). The electron-pair band remains above the corresponding single-electron band because the Coulomb repulsion between the two electrons is stronger than their correlation. However, the shape of the dispersion relations of the two bands are very similar. This confirms the dispersion relation of the electron-pair band discussed in previous section within the framework of the tight-binding approach. The similarity is due to the fact that two paired electrons are closely bound in the relative coordinates in real space with an average separation r less than half of the lattice period. Furthermore, the single-electron and electron-pair states in the relative coordinates are of the same symmetry. The average distance between the two electrons in the case of Fig. 2(a) is r = 0.44λ. The energy gap between the electron-pair band and the single-electron band is E (0) g . The energy difference between the bottoms of the two bands at Γ point is defined as ν 0 . The electron pair behaves as a larger particle with both the mass and charge twice of a single electron. Consequently, the tunneling probability of an electron pair to its neighbor site is much smaller than that of a single electron leading to a much narrower electron-pair band.
We also find that the dispersion relation of the electron-pair band can indeed be described by Eq. (13). For instance, the dispersion of the electron-pair band in Fig. 2 , where E p,0 = −1.35401 ± 0.00002 R y and J p,1 = 0.03368 ± 0.00002 R y . The fitting gives an extremely small fitting parameter χ 2 = 3×10 −8 .
In Fig. 2(b) we plot the electron-pair band as a function of the lattice constant λ together with the two lowest single-electron bands for fixed potential amplitude V 0 =15 R y . The two curves for each band indicate the minimum and maximum energies of the band. We see that the electron-pair band appears for λ > ∼ 1.22 a B with a bandwidth < ∼ 0.2 R y . It stays above the lowest single-electron band with a gap E (0) g of about 3 to 5 R y . The bandwidth of the electron-pair band is at least one order of magnitude smaller than that of the single-electron band. Fig. 2(c) shows the electron-pair bands as a function of λ for different V 0 . We see that the electron-pair band appears for E < 0 because the energy spectrum is continuous for positive energy in the crystal potential given by Eq. (14). It means that co-tunneling of the paired electrons is required to form the energy band. Therefore, their bandwidth is much smaller than that of the single-electron band. The existence of the metastable electron-pair band is a result of the local confinement in each unit cell, the electron-electron correlation, and the co-tunneling of the electron pair in the crystal. In Fig. 2(d) we plot the energy gap E (0) g together with ν 0 . They are important quantities for the renormalization of the electron-pair states.
In the rest of the paper, we will demonstrate that the metastable electron pairs can be stabilized at certain electron densities by including electron-electron and electron-phonon interactions in the crystal. Since the electron pairs are spin singlet, they will be considered as bosonic quasiparticles and mostly distributed at the bottom of the electron-pair band at low temperature. The many-body effects in the crystal renormalize the band structure. If the band renormalization can bring down the bottom of the electron-pair band to the Fermi surface of the single-electron band, the electron pairs at the bottom of the band cannot decay into two single electrons because of the Pauli exclusion principle. In such a case, the electron pairs can be stabilized.

III. HAMILTONIAN OF THE TWO-BAND MANY-ELECTRON SYSTEM INTERACTING WITH LO-PHONONS
We consider the 2D electron system with two energy bands: a lower single-electron band E s (k) and a higher metastable electron-pair band ν 0 + E p (k). The bottom of the single-electron band is taken as reference for energy E = 0 and the bottom of the electron-pair band is at E = ν 0 (see Fig. 2). Assuming that there are N t electrons in the system consisting of N s single electrons and N p electron pairs with N t = N s + 2N p , the many-particle exchange-correlation (xc) interactions will renormalize the energy bands reducing the energies of both the single electrons and electron pairs. In real materials, a 2D electron system can be found in the interface and surface of bulk materials or in a 2D layer of layered crystal structure such as superconducting cuprate, therefore interaction between the electron system and crystal lattice vibration and polarization affects the electron states. Considering the ionic and polar-covalent characteristics of many su-perconducting materials, the electron-LO-phonon interactions can be significant and their contribution to the energy-band renormalization is important. [16,47] In this section, we will present the many-particle Hamiltonian including electron-electron (e-e) and electron-phonon (eph) interactions assuming that the 2D electron layer is immersed in a 3D phonon field. The Hamiltonian of the considered system is giving by where H el is for the electronic part with single electrons and electron pairs, H ph for 3D LO-phonons, and H el−ph for e-ph interaction. The Hamiltonian of the electronic part H el was derived in Ref. [48], and is given by where H single is for the single-electron (se) band, H pair for the electron-pair (ep) band, H s−p for single-electronelectron-pair (se-ep) interaction. For the N s electrons in the single-electron band, their Hamiltonian is given by where v q = v ss (q) = 2 2π q is the single-electron-singleelectron (se-se) Coulomb potential, the operators c † k,σ and c k,σ are creation and annihilation operators, respectively, for a single electron of momentumhk and spin σ. They obey the fermion anti-commutation rela- For N p electron pairs in the electronpair band, the Hamiltonian is given by where v pp (q) = 4v q f pp (q) is the electron-pair-electronpair (ep-ep) interaction potential with the form factor f pp (q) given in Ref. [48], the operators b † k and b k are creation and annihilation operators, respectively, for a spinsinglet electron pair of momentumhk. They obey the bo-

The se-ep interband interaction is give by
with and where v sp (q) = 2v q f sp (q) is the se-ep interband scattering potential (without breaking the electron pair) with form factor f sp (q) and v t (q) is the se-ep interaction potential for interband transition (breaking or forming electron pairs). They are given in Ref. 48 and f 2 sp (q) = f pp (q). Notice that in the above Hamiltonian, the summation over q does not include q=0 because it is cancelled out with the background ion-ion interaction and the system is neutral.
As we discussed in the previous section, the single electrons and the paired electrons share the same space in the crystal. The maximum total electron density is two electrons per unit cell (per atom), paired or not. However, when the potential amplitude V 0 of the crystal is larger than a certain value, two electrons in the same unit cell will occupy the same atomic orbital in the relative coordinates forming an electron pair. It means that in this case, each unit cell can be occupied by a single electron or by an electron pair, but not both at the same time. Therefore, for a certain single-electron density n s = N s /A (where A is the area of the sample), the maximum electron-pair density n p = N p /A in the square-lattice crystal is given by where λ −2 is the density of the unit cell of the crystal. Notice that n s = λ −2 corresponds to half filling of the single-electron band. The above condition indicates that the single-electron band should be less than half filled if there are any electron pairs in the crystal. The Hamiltonian of the optical-phonon modes in bulk materials with energyhω LO and 3D wavevector Q =(q, q z ) is given by where a † Q (a Q ) is the creation (annihilation) operator of the LO-phonons. The interaction Hamiltonian of a many-electron system with the LO-phonons is given by, [16] with the Fourier coefficient of the e-ph interaction potential where r j is the position of the electron j with band mass m * . The Fröhlich electron-phonon coupling constant α is defined by where ǫ 0 and ǫ ∞ are the static and high-frequency dielectric constant, respectively. In the considered two-band system, the gap between the two bands is E (0) g in the zero density limit. Although the e-e and e-ph interactions in the crystal reduce the energies of both the single-electron and electron-pair bands, their effects on the electrons in the paired states are larger than those in the unpaired single-electron band. Therefore, we expect that many-body effects will reduce the energy gap between the two bands. We will calculate the exchange-correlation corrections as well as the polaron energies for both single electrons and electron pairs including the screening effects in order to find out whether the electron pairs can be stabilized or not.
In order to obtain the ground-state energy of the present many-particle system with the single electrons, electron pairs and LO-phonons, we will employ the Lee-Low-Pines transformation [49] in dealing with the manypolaron system [49][50][51][52]. For weak e-ph coupling, we can assume the ground state of the electron-phonon system as |GS = |GS el |VAC ph , where |GS el is the ground state of the electronic part and |VAC ph is the phonon vacuum state with zero real phonons. [49,50] The above approximation is valid for weak and intermediate e-ph coupling strength [47,49,50,53,54] and it allows us to calculate separately the electron exchange-correlation and polaron contributions to the ground-state energy of the system, given by where E GS is the ground-state energy of the electronic part without interaction with phonons and E (tot) pol is the total polaron correction due to electron-phonon interaction.
In the next two sections, we will calculate the contributions of the electronic exchange-correlation interaction and electron-phonon interaction to the renormalization of the single-electron end electron-pair energies.

IV. ELECTRON-ELECTRON INTERACTIONS WITH SINGLE ELECTRONS AND ELECTRON PAIRS
It is known that in both 3D and 2D systems the exchange-correlation energy for band-gap renormaliza-tion (BGR) is almost independent of the band characteristics. The many-particle exchange-correlation energy depends only on the inter-particle distance r s (determined by the particle density) in appropriate rescaled natural units in a universal manner [55][56][57]. The contribution of the electron-electron interaction to the BGR can be obtained by calculating the average exchangecorrelation energies per particle or by calculating the selfenergies of the particles involved. The kinetic energy is usually assumed to be unchanged in a renormalization process.
In this section, we will calculate the ground-state energy E (el) GS of the electronic part consisting of the single electrons and electron pairs. The corresponding Hamiltonian is given by Eqs. (24)(25)(26)(27)(28)(29). Although the electron pairs are metastable, we will treat them in the calculations as if they were stable particles. Only the final results including full many-body corrections will tell us whether they can really be stabilized or not. In the calculations, we will first take the single-electron density n s and electron-pair density n p as inputs. The single electrons are considered as fermions and electron pairs as bosons. Therefore, we are dealing with a many-particle system consisting of a boson-fermion mixture. [58] The exchange-correlation energies are obtained as a function of n s and n p . We then determine their contributions to the ground-state energy and the band renormalization. Within such a scheme, the se-ep interaction Hamiltonian H t int given in Eq. (29) will not be invoked explicitly in the calculations. However, transitions between the singleelectron and electron-pair bands are permitted because of this term. The many-particle interaction energy in such a two-component system of boson-fermion mixture can be obtained by [59] where the inter-particle interaction potential E int ij (ξ) is given by The above potential depends on the "bare" inter-particle potential v ij (q) and static structure factor S ij (q), for i, j = s (single electron) and p (electron pair). The static structure factor can be calculated by, Within the linear response theory, the density-density response function χ ij (q, ω) of this two-component plasma is given by [60,61], where χ ii (q, ω) is the non-interacting polarizability [60][61][62][63] of the ith component and ϕ ij (q) is the static effective interaction potential. The function χ (0) ss (q, ω) is for a noninteracting 2D electron gas given in Ref. [62]. The function χ (0) pp (q, ω) is the polarizability of a non-interacting 2D boson gas of electron pairs with density n p . The effective potential ϕ ij (q) defines a local field correction [60,64] in terms of the "bare" potential v ij (q). Within the randomphase approximation (RPA), the local field correction on the static effective interaction is neglected and therefore ϕ ij (q) = v ij (q). The potential v ij (q) has been determined in the previous section given by v ss (q) = v q , v pp (q) = 4v q f pp (q), and v sp (q) = v ps (q) = 2v q f sp (q).
For a non-interacting 2D boson (i.e., electron-pair) gas with density n p , we can assume that all the electron pairs are in the same state at the bottom of the electronpair band at zero temperature, i.e., in the condensate phase [61,63]. The polarizability of the non-interacting boson gas is given by, where ε q,p = q 2 /2. The contribution of the e-e interaction to the singleelectron band renormalization is given by ∆E s = E ss + E sp and to the electron-pair band given by ∆E p = (E pp + E sp )/2. The ground-state energy of the two-band system is given by where E kin is the kinetic energy of the many-electron system. We calculate these energies within the RPA. Although the RPA overestimates the exchange-correlation energies ∆E s and ∆E p , only the difference between them contributes to the band-gap renormalization. The errors resulting from the RPA should be partially cancelled in the process when determining the condition of stability of the electron pairs. Therefore, we consider the RPA a reasonable approximation for our purpose. As mentioned above, the kinetic energy E kin will be assumed unchanged in the renormalization. It is given by E kin = N s ε F /2, where ε F is the Fermi energy of the single-electron band in relation to its band bottom. In a two-dimensional system, the average kinetic energy of a single electron is ε F /2. The electron pairs have no kinetic energy because they are assumed to be at the bottom of the electron-pair band in the condensate phase.
If we further assume that v sp (q) = 0, the single electrons and electron pairs become independent in different energy bands. The single-electron system is an usual onecomponent 2D electron gas which has been widely investigated within different methods and approximations. [62,65,66] The system of electron pairs is new. Although the electron pairs are metastable states, we will treat them as stable particles when searching for their manyparticle ground state. The static structure factor of a one-component boson system of electron pairs within the RPA is given by, Consequently, we obtained the many-electron-pair exchange-correlation energy, where r p s = (πn p ) − 1 2 , and The calculations in Ref. [48] showed that f pp (q) is a monotone decreasing function with f pp (0) = 1 and f pp (∞) = 0. If we take f pp (q) ≡ 1 (consequently r = 0), i.e., assuming the electron pair as an ideal boson with charge -2e and mass 2m * , the above integral becomes I = I 0 = −1.29. This is the well-known RPA result for an ideal charged 2D boson gas [67]. Notice that the factor 8 in Eq. (43) is due to the units used here. In Fig. 3 we show the many-particle interaction energies E ss , E sp , and E pp within the RPA in the system keeping the same number of electrons in the singleelectron and electron-pair bands, i.e., n s = 2n p . The energy E sp = E ps is positive because of the se-ep repulsion. Fig. 3(a) is for the 2D crystal with λ = 1.5 a B and V 0 = 15 R y . It shows the density dependence of the interaction energies and the effects of seep interaction. The blue-dash and red-dash curves are the exchange-correlation energies E (1) pp and E (1) ss , respectively, without the se-ep interaction. In this case, the energy E (1) pp is obtained from Eqs. (43) and (44). We see that the nonzero distance between the two electrons in the pair (i.e., r > 0) affects the energy E pp . In the calculations we found that λ/2 > r > λ/3. If we assume r = 0 for the electron pairs, we obtain for an ideal charged boson system in 2D indicated by the black-dotted curve.
We observe that the se-ep interband interaction not only introduces the energy E sp but also reduces the energies E pp and E ss being evident in the difference between the solid and dashed curves in Fig. 3(a). The density dependence of the ep-ep interaction energy E pp is different from that of the se-se interaction E ss . For example, E pp is about 7 times larger than E ss at lower density n s = 2n p = 0.001 a −2 B and this ratio is reduced to 3 times at higher density n s = 2n p = 1.0 a −2 B . In Fig. 3(b) we show the many-particle interaction energies for different potential V 0 . With increasing V 0 for the same lattice constant λ, the average distance r between the two electrons in the same pair decreases and, consequently, the form factor f pp (q) increases (as shown in the inset) and the energy E pp becomes larger.

V. POLARON EFFECTS IN THE TWO-BAND SYSTEM WITH SINGLE ELECTRONS AND ELECTRON PAIRS
In this section, we will study the polaron effects on the N t electrons in the two-band system interacting with LO-phonons assuming N s electrons in the single-electron band and N p electron pairs in the electron-pair band, N s + 2N p = N t . The Hamiltonian of the system is given in Eq. (23) and the electron-phonon interaction given by Eq. (32). Considering a boson-fermion mixture (namely, the electron pairs and single electrons) interacting with the phonons, the electron-phonon interaction Hamiltonian in Eq. (32) can be separated into two parts, with for N s single electrons, and for N p electron pairs, where V p Q (r) = 2V Q cos( Q·r 2 ) and R j and r are the center-of-mass and relative coordinates of the electron pair, respectively.
In the calculations of the polaron energies, we will ignore the direct participation of the single-electronelectron-pair interaction potential v sp (q). The potential v sp (q) has a minor effect on the e-ph interaction. Its influence on the e-ph interaction is mostly indirect through the electronic screening and it is taken into account in the static structure-factor. Within such an approximationignoring the potential v sp (q) in dealing with the e-ph interactions, the Hamiltonian of the whole system defined in Eq. (23) can be separated into two subsystems. They are the one consisting of single electrons interacting with phonons and the other of electron pairs interacting with phonons. In this way, the contribution of the electronphonon interaction to the ground-state energy E tot pol in Eq. (35) can be calculated as where E single pol and E pair pol are the polaron energies of the single electron and electron pair, respectively.
The subsystem composed of single electrons interacting with LO-phonons is where H single is given by Eq. (25), H ph by Eq. (31), and H single e−ph by Eq. (46). The e-ph coupling leads to a polaron consisting of an electron and a surrounding phonon cloud. When the e-ph interaction is not too strong, the polaron correction to the ground-state energy of an electron gas can be calculated within the Lee-Low-Pines (LLP) unitary transformation method. [49,50] This method has been used to study polaron gases in bulk materials and also in low-dimensional systems. [51,52] It is known that the polaron energy obtained from the LLP transformation is exact for α → 0. In the low electron density limit, the polaron energy obtained from the LLP method for e-ph coupling constant α = 6 is 90% of the exact value. [47,54] Here we are dealing with a polaron gas in which the screening reduces the electron-phonon interaction strength. Therefore, the LLP method should yield a reasonable polaron energy for α < 6.
We calculate the polaron energy within the LLP method for the above Hamiltonian in Eq. (49), given by, where S ss (q) is the static structure factor of the singleelectron gas. In the low electron-density limit, S ss (q) = 1. This leads to the well known perturbation result E single pol = −(π/2)αhω LO for the polaron energy of an electron in 2D coupled with 3D-phonons. [51,68] Fig. 4(a) shows the polaron energy as a function of the single-electron density n s in the 2D square-lattice crystal with λ = 1.5 a B coupled with the 3D LO-phonons. The polaron energy in the low-density limit without screening is indicated by the thin-dotted line. It is seen that the screening of the electron gas considered in the structure factor S ss (q) in Eq. (50) reduces the polaron effect. At higher electron density n s = 1.0 a −2 B , the polaron energy is only about 20% of its low-density value. When there are also electron pairs in the crystal and the se-ep interaction is included in the structure factor S ss (q) in the screening, the calculation results show that the se-ep interaction reduces the screening and consequently, enhances the polaron energy. The dashed curves in Fig. 4(a) are obtained with electron-pair density n p = n s /2 interacting with single-electrons. The red, green and blue dashed curves are for V 0 = 12, 15, and 18 R y , respectively. The potential V 0 affects the potentials v sp and v pp . Its influences on the polaron energy is indirectly through the structure factor S ss (q). Therefore, we obtained almost the same value of E single pol for different V 0 .
The subsystem consisting of electron pairs and LOphonons is given by the following Hamiltonian, where H pair is given by Eq. (26), H ph by Eq. (31), and H pair e−ph by Eq. (47). Following a similar procedure as before and applying the LLP transformation to the electron-pair-phonon interactions, we obtain the polaron contribution to the energy of the electron pairs, given by where S pp (q) is the static structure factor of the electron pairs. The matrix element is given by where |Ψ k (R, r); Q = |Ψ k (R, r) |Q , with |Ψ k (R, r) for the electron-pair state and |Q for the phonon state. Using the metastable electron-pair wavefunction in Eq. (22), we obtain where V Q is the e-ph interaction potential given by Eq. (33), f pp (q) is the form-factor that appears in the pair-pair interaction potential in Eq. (26).
We now obtain a Fröhlich bipolaron formed by an electron pair coupled with LO-phonons. Fig. 4(b) shows the bipolaron energy in the crystals with λ = 1.5 a B and V 0 = 12, 15, and 18 R y . The solid curves give the bipolaron energies for different V 0 with only screening of the electron pairs. For larger V 0 , the average distance r between the two electrons in a pair is smaller, i.e., the size of the electron pair is smaller. Therefore, the bipolaron energy is larger. But the screening reduces significantly the bipolaron energy. We also see that, when se-ep interaction is included in the structure factor S pp (q), the screening of the electron-pair-phonon coupling is reduced enhancing the bipolaron energy as shown by the dashed curves for n s = 2n p in the figure.
If we assume S pp (q) = 1, we obtain the low-density limit of the bipolaron energy without screening given by the horizontal dotted-dash lines. This is the case of a single bipolaron without screening. If we further assume f pp (q) = 1, we obtain the upper bound of the bipolaron energy, E pair pol = −4 √ 2(π/2)αhω LO . [69] This is equivalent to assuming the electron pair as a "larger single electron" with charge -2e and mass 2m * coupled to the LOphonons.
As a matter of fact, the effects of electron-phonon interaction and possible bipolaron formation have been extensively studied as the electron pairing mechanism for unconventional superconductivity. [47,[70][71][72] In order to form a bipolaron in the crystal, a crucial point is that the electron-phonon coupling induced attraction has to overcome the electron-electron Coulomb repulsion. [47] This requires not only a small ratio ǫ ∞ /ǫ 0 of dielectric constants but also a large enough electron-phonon coupling constant α leading to a critical e-ph coupling constant α c = 2.9 in 2D and α c = 6.8 in 3D for bipolaron formation.
However, the bipolaron formation mechanism in the present theory is distinct from the traditional bipolaron theory in the literature. In this paper, we show a preformed metastable electron pair due to strong correlation of two electrons occupying the same orbital. The electron-phonon interaction dresses the electron pairs up with a phonon cloud forming bipolarons. Therefore, in the present context we introduce an internal interaction due to orbital dependent electron correlation of the electron pair. The electron-phonon coupling involves in the renormalization of the preformed electron pairs. Notice that polaron effects can be much larger on the electrons in the correlated pairs than on the single electrons, and therefore the bipolaron contribution to the stabilization of the electron pairs overcoming the Coulomb repulsion becomes essential. We will show in the next section that the certain electron pairs can indeed be stabilized as bipolarons.

VI. STABILIZATION OF THE ELECTRON PAIRS AND THE GROUND STATE OF THE MANY-PARTICLE SYSTEM
The condition for stabilization of the electron pairs when including renormalization is that the bottom of the electron-pair band occurs at the Fermi energy of the single-electron band. Using the same energy reference defined in Sec. III, i.e., taking the bottom of the singleelectron band at Γ point before the renormalization as E = 0, the many-particle interactions lower the singleelectron band bottom to where E ss and E sp are determined by Eq. (36) and E single pol is given by Eq. (50). The Fermi energy of the singleelectron band can be obtained by E F = E Γ s + ε F where ε F is the energy difference between the Fermi energy and the bottom of the band. It is determined by the singleelectron density and the density of states of the band. We will not consider the many-body effects on ε F . The average kinetic energy per single electron in this band is given by ε F /2 in a two-dimensional system.
Including the band renormalization, the bottom of the electron-pair band is given by where E pp and E sp are given by Eq. (36) and E pair pol by Eq. (52). Notice that the energy of an electron pair is given by 2E Γ p . Fig. 5 shows the energies for the renormalized singleelectron and electron-pair bands in the 2D crystal with λ = 1.3 a B and V 0 = 15 R y coupled with the 3D LO-phonons. The band structure without many-body corrections was given in Fig. 2(a). The many-particle interaction energies are obtained for a constant total electron density n t = n s + 2n p = 0.5 a −2 B with different e-ph coupling constant α. Because the polaron energy is given by the coupling constant α and the LOphonon energyhω LO , the ratio R y /hω LO is required in the calculations. [73] This ratio is a material dependent parameter and has been used in the bound-polaron and bound-bipolaron problems [73,74]. Its value is in the range from 0.5 to 2 for different materials. [74,75] For the present discussion, we take R y /hω LO =1. Fig. 5(a) shows for fixed total electron density n t = 0.5 a −2 B and e-ph coupling constant α = 4, the dependence of the single-electron band bottom E Γ s , its Fermi energy E F , and the electron-pair band bottom E Γ p on the electron distribution in the two bands n s and n p . The total energy per electron E tot is also given in the figure. It is the average ground-state energy per electron in the system calculated by Eqs. (35), (41), and (48), and is given by GS + E where E Γ s and E Γ p are given by Eqs. (55) and (56), respectively. It is seen that, though the total electron density n t is a constant, the above obtained energies depend on the distribution of the electrons between the two bands. Especially the energy of the electron-pair band E Γ p depends strongly on n p . In Fig. 5(a) we see that with increasing the single-electron density n s , the single-electron band E Γ s and the Fermi energy E F vary slowly. Because n t is a constant, increasing n s means decreasing n p . At low density n p , the weak screening in the electron-pair band enhances greatly the bipolaron energy resulting in a lowering of the bottom of the electron-pair band E Γ p which reaches the Fermi surface. We find that under the considered condition, E Γ p touches the Fermi surface at n s = 0.4774 a −2 B and n p = 0.0113 a −2 B as indicated by the black dot. For these densities, the electron pairs (or the bipolarons) are stabilized on the Fermi surface. It means that only 4.5% of the electrons in the system form stable electron pairs in this case. If we look at the total energy E tot , it tends to reduce the energy of the system at the higher single-electron density side. But the electron pairs on the Fermi surface cannot decay into two single electrons due to the Pauli exclusion principle. Moreover, breaking an electron pair needs a cost to overcome their correlation energy. Therefore, the electron-pair density is stabilized at n p = 0.0113 a −2 B in the ground state of the system with α = 4 and n t = 0.5 a −2 B . Fig. 5(b) shows the electron-pair band energy E Γ p and the single-electron band Fermi energy E F in the same system with n t = 0.5 a −2 B but for different e-ph coupling constant α. We see that for α = 2 no stable electron pairs are found. The calculations indicate that for α > ∼ 3, part of the electrons can form electron pairs on the Fermi surface at low electron-pair density. The stabilized electron-pair densities are n p = 0.001, 0.011, and 0.026 a −2 B for α = 3, 4 and 5, respectively.
In order to determine the electron-pair density in the ground-state of the system, we solve the following equation E Γ p (n s , n p ) = E F (n s , n p ), as a function of n s and n p for fixed λ, V 0 and α. Figs. 6(a) and 6(b) show the density of stabilized electron-pairs or bipolarons in linear and logarithmic scale, respectively, as a function of the single-electron density in the system with λ = 1.3 a B and V 0 =15 R y for α = 2.8, 3.0, 3.5, 4.0, 4.5, and 5.0. We see that for α = 2.8 a very low density of electron pairs become stabilized. With increasing α, more electron pairs appear in the system. For α = 5, their density may reach n p =0.032 a −2 B . This corresponds to about 11% of the electrons in the system are in the electron-pair band. The thick-red curve in the figure is the possible maximum electron-pair density n max p given by Eq. (30). It is important to notice that the relation between n max p and n s not only restrict the electron-pair density being less than n max p . It also means that if n s is larger than λ −2 , i.e., the single-electron band filling is more than half, there are no stable electron pairs in the system. For the crystal with λ = 1.3 a B , the half filling of the single-electron band occurs at n s = 0.5917 a −2 B . In Fig. 6(c) we show the electron-pair density in a different 2D potential with λ = 2.8 a B and V 0 =5 R y . In this case, the band-gap E (0) g and also ν 0 are small, as shown in Fig. 2. The half filling of the single-electron band is at n s = 0.1276 a −2 B . We find that for α ≃ 2, some electron pairs can already be stabilized with small density. For α = 5, the electron-pair density may reach 0.02 a −2 B corresponding to 31% of the total electrons being in the paired state.

VII. SUMMARY AND OUTLOOK
The starting point of our work is that the electronelectron correlation is orbital dependent. The exchangecorrelation energy of an electron pair occupying the same orbital is larger than the average of the exchangecorrelation energy of many electrons.
Depending upon the crystal structure and potential, such electron pairs can form a metastable electron-pair band. The metastable electron-pair band is obtained in two different ways. One follows the tight-binding method and the other is similar to the nearly-free electron model. They give consistent results for the dispersion relation of the electron-pair states in the period potential.
Combining the metastable electron-pair band together with the corresponding single-electron band, we constructed a many-particle system consisting of single electrons and metastable electron pairs. When further considering many-body interaction renormalization including single electrons, electron pairs, and optical phonons, we found that the metastable electron pairs can be stabilized. The calculations show that the polaron effects play an essential role in counterbalancing the Coulomb repulsion in the stabilization of the electron pairs. The electron-phonon coupling is enhanced in the strongly correlated electron pairs and leads to the formation of bipolarons. On the other hand, screening affects significantly both the polaron and bipolaron energies manifesting a cooperative interplay of electron-electron and electronphonon interactions. The obtained ground state of the system consists of polarons in the single-electron band and bipolarons in the electron-pair band sitting on the Fermi surface of the single-electron band.
The numerical calculations presented in this paper are performed for simple potentials and with the s-orbital in 2D square-lattice crystals. However, the physical processes and numerical calculations for electron pairing can be extended to a quasi-2D crystal of a single atomic layer with finite thickness or to 3D systems. In principle, it can also be extended to study the electron pairing in pand d-orbitals. The obtained results within the present simple model predict light and small electron pairs (bipolarons) in the crystal. In the center-of-mass coordinates, the electron pair is a Bloch wavefunction with an effective mass of twice a single electron. The pair is small and local because the average separation between the two electrons is less than half of the lattice constant and they are localized in the same unit cell when expressed in the relative coordinates. Furthermore, only a fraction of the electrons in the system form pairs.
We have obtained a many-particle ground state with spin-singlet electron pairs in the condensate phase. We expect that these preformed electron pairs at certain densities in coherent state will contribute to superconductivity. Finally, we want to comment on the possible "superconducting energy gap". From our calculations we naturally infer such a gap to the transition energy required to break the stabilized electron pairs sitting on the Fermi surface. This transition is determined by the Hamiltonian H t int in Eq. (29) where the potential v t (q) was given by Eq. (11) in Ref. [48]. Because the electron pairs are stabilized on the Fermi surface in the center of the Brillouin zone at k = 0, an external excitation has to overcome the internal correlation energy of an electron pair to bring two electrons to the single-electron states above the Fermi surface at ±k = 0. The minimum energy required (or the gap) is primarily determined by the potential v t (q) and the final momenta ±hk of the single-electron states. It is band structure dependent, anisotropic, and nodeless.