# Cooper Instability in the Occupation Dependent Hopping Hamiltonians

## I Formulation of the model

High temperature superconductivity in the lanthanum , yttrium  and related copper-oxide compounds remains a subject of intensive investigation and controversy. It was suggested that electron-phonon interaction mechanism, which is very successful in understanding of conventional (“low temperature”) superconductors within the Bardeen-Cooper-Schriffer scheme , may not be adequate for high-$`T_c`$ cuprates, and even the conventional Fermi liquid model of metallic state may require reconsideration. This opens an area for investigation of mechanisms of electron-electron interaction which can be relevant in understanding peculiarities of superconducting, as well as normal state, properties of cuprates. Specific to all of them is the existence of oxide orbitals. Band calculations  suggest that hopping between the oxygen $`p_x,p_y`$ orbitals and between the copper $`d_{x^2y^2}`$ orbitals may be of comparable magnitude. On the experimental side, spectroscopic studies  clearly show that the oxygen band appears in the same region of oxygen concentration in which superconductivity in cuprates is the strongest. Therefore there exists a possibility that specific features of oxide compounds may be related to oxygen-oxygen hopping, or to the interaction between the copper and the rotational $`p_xp_y`$ collective modes. If the oxygen hopping is significant then it immediately follows that intrinsic oxygen carriers ($`p_x,p_y`$ oxygen holes) should be different from the more familiar generic $`s`$-orbital derived itinerant carriers. The difference is related to low atomic number of oxygen such that removing or adding of one electron to atom induces a substantial change in the Coulomb field near the remaining ion and therefore results in the change of the effective radius of atomic orbitals near the ion. This will strongly influence the hopping amplitude between this atom and the atoms in its neighborhood. Such “orbital contraction” effect represents a source of strong interaction which does not simply reduce to the Coulomb (or phonon) repulsion (or attraction) between the charge carriers. It was suggested by Hirsch and coauthors , and by the present authors  that the occupation dependent hopping can have relevance to the appearance of superconductivity in high-temperature oxide compounds. In the present paper, we investigate the generic occupation-dependent hopping Hamiltonians with respect to peculiarities of the normal state, and to the range of existence of the superconducting state. Theoretical investigation of Cooper instability is supplemented by numeric study of pairing and diamagnetic currents in finite atomic clusters. We study the effect of Cooper pairing between the carriers and show that at certain values and magnitudes of the appropriate coupling parameters, the system is actually superconducting. The properties of such superconducting state are in fact only slightly different from the properties of conventional (low-$`T_c`$) superconductors. Among those we so far can only mention the change in the fluctuation conductivity above or near the critical temperature $`T_c`$. Relaxation of the pairing parameter to equilibrium acquires a small real part due to the asymmetry of contraction-derived interaction between the quasi-particles above and below the Fermi energy.

Oxygen atoms in the copper-oxygen layers of the cuprates (Figure 1) have simple quadratic lattice. We assume that $`p_z`$ orbitals of oxygen ($`z`$ is the direction perpendicular to the cuprate plane) are bound to the near cuprate layers whereas carriers at the $`p_x,p_y`$ orbitals may hop between the oxygen ions in the plane.

Let $`t_1`$ be the hopping amplitude of $`p_x(p_y)`$ and $`t_2`$ the hopping amplitude of $`p_y(p_x)`$ oxygen orbitals between the nearest lattice sites in the $`x(y)`$ direction in a square lattice with a lattice parameter $`a`$. Then the non-interacting Hamiltonian is

$$H_0=t_1\underset{<ij>_x}{}a_i^{}a_jt_2\underset{<ij>_y}{}a_i^{}a_jt_1\underset{<ij>_y}{}b_i^{}b_jt_2\underset{<ij>_x}{}b_i^{}b_j$$
(1)   
where $`a_i^{}(a_i)`$ is the creation (annihilation) operator for $`p_x`$ and correspondingly $`b_i^{}(b_i)`$ for $`p_y`$ orbitals. The interaction Hamiltonian includes the terms

$$H_1=\underset{<ij>}{}a_i^{}a_j\left[Vm_im_j+W(m_i+m_j)\right]+\underset{<ij>}{}b_i^{}b_j\left[Vn_in_j+W(n_i+n_j)\right]$$
(2)   
where $`n_i=a_i^{}a_i`$, $`m_i=b_i^{}b_i`$. This corresponds to the dependence of the hopping amplitude on the occupation numbers $`n_i,m_i`$ of the form

$$\left(\widehat{t}_{ij}\right)_{a_ia_j}=\tau _0(1m_i)(1m_j)+\tau _1\left[(1m_i)m_j+m_i(1m_j)\right]+\tau _2m_im_j$$
(3)   
and correspondingly $`(\widehat{t}_{ij})_{b_ib_j}`$ of the same form with $`m_i`$ replaced with $`n_i`$. The amplitudes $`\tau _0,\tau _1,\tau _2`$ correspond to the transitions between the ionic configurations of oxygen:

$`\tau _0:`$ $`O_i^{}+O_j^2`$ $`O_i^2+O_j^{}`$ (4)   
$`\tau _1:`$ $`O_i+O_j^2`$ $`O_i^2+O_j`$ (5)   
$`\tau _2:`$ $`O_i+O_j^{}`$ $`O_i^{}+O_j`$ (6)   
$`O`$ corresponds to the neutral oxygen ion whereas $`O^{}`$ to the single charged and $`O^2`$ to the double charged negative ions. Since oxygen atom has $`1s^22p^42s^2`$ configuration in its ground state, filling of the $`p`$ shell to the full occupied configuration $`2p^6`$ is the most favorable. Amplitudes $`V`$ and $`W`$ relate to the parameter $`\tau _0,\tau _1,\tau _2`$ according to

$$V=\tau _02\tau _1+\tau _2,W=\tau _1\tau _2.$$
(7)   
Assuming $`t_1=t_2`$ and replacing $`a_i,b_i`$ with $`a_i`$ with the pseudo-spin indices $`\sigma =,`$ we write the Hamiltonian Eq.(1) in the form

$$H=t\underset{<ij>\sigma }{}a_{i\sigma }^{}a_{j\sigma }+H_U+H_V+H_W$$
(8)   
where

$`H_U`$ $`=`$ $`U{\displaystyle \underset{i}{}}n_in_i`$ (9)   
$`H_V`$ $`=`$ $`V{\displaystyle \underset{<ij>\sigma }{}}a_{i\sigma }^{}a_{j\sigma }n_{i,\overline{\sigma }}n_{j,\overline{\sigma }}`$ (10)   
$`H_W`$ $`=`$ $`W{\displaystyle \underset{<ij>\sigma }{}}a_{i\sigma }^{}a_{j\sigma }(n_{i,\overline{\sigma }}+n_{j,\overline{\sigma }})`$ (11)   
where we also included the in-site Coulomb interaction (U) between the dissimilar orbitals at the same site. $`\sigma `$ can also be considered as a real spin projection of electrons at the site. In that case, the pairing will originate between the spin-up and spin-down orbitals, rather than between $`p_x`$ and $`p_y`$ orbitals. More complex mixed spin -and orbital- pairing configurations can also be possible within the same idea of orbital contraction (or expansion) at hole localization but are not considered in this paper. The following discussion does not distinguish between the real spin and the pseudo-spin pairing. The Hamiltonian, Eq.(6), is a model one which can not refer to the reliable values of the parameters appropriate to the oxide materials. The purpose of our study is rather to investigate the properties of superconducting transition specific to the model chosen and to find the range of the $`U,V,W`$ values which may correspond to superconductivity. This will be done along the lines of the standard BCS model  in the weak coupling limit, $`U,V,W\mathrm{\hspace{0.17em}0}`$, and by an exact diagonalization of the Hamiltonian for a finite atomic cluster at large and intermediate coupling.

In the momentum representation, the Hamiltonian becomes $`H=H_0+H_1+H_2`$ with

$`H_0={\displaystyle \underset{𝐩\sigma }{}}\xi _𝐩a_{𝐩\sigma }^{}a_{𝐩\sigma }`$ (12)   
$`H_1={\displaystyle \frac{1}{4}}{\displaystyle \underset{p_1p_2p_3p_4,\alpha \beta \gamma \delta }{}}a_{𝐩_1\alpha }^{}a_{𝐩_2\beta }^{}\mathrm{\Gamma }_{\alpha \beta \gamma \delta }^0(p_1,p_2,p_3,p_4)a_{𝐩_4\delta }a_{𝐩_3\gamma }`$ (13)   
where

$$\xi _𝐩=t\sigma _𝐩\mu ,\sigma _𝐩=2(\mathrm{cos}p_xa+\mathrm{cos}p_ya),$$
(14)   
and $`\mu `$ is the chemical potential. $`\mathrm{\Gamma }_{\alpha \beta \gamma \delta }^0`$ is the zero order vertex part defined as

$$\mathrm{\Gamma }_{\alpha \beta \gamma \delta }^0(p_1,p_2,p_3,p_4)=\left[U+(W+\frac{1}{2}\nu V)(\sigma _{𝐩_1}+\sigma _{𝐩_2}+\sigma _{𝐩_3}+\sigma _{𝐩_4})\right]\tau _{\alpha \beta }^x\tau _{\gamma \delta }^x(\delta _{\alpha \gamma }\delta _{\beta \delta }\delta _{\alpha \delta }\delta _{\beta \gamma })\delta _{𝐩_1+𝐩_2,𝐩_3+𝐩_4}$$
(15)   
where $`\tau _{\alpha \beta }^x`$ is a Pauli matrix

$`\left(\begin{array}{cc}0\hfill & \hfill 1\\ 1\hfill & \hfill 0\end{array}\right).`$   
For reasons which will be clear later, we separated $`H_V`$ and put some part of it into the $`H_1`$ term, while the remaining part is included in the $`H_2`$ term, thus giving

$$H_2=V\underset{<ij>\sigma }{}a_{i\sigma }^{}a_{j\sigma }(a_{i\overline{\sigma }}^{}a_{i\overline{\sigma }}\frac{\nu }{2})(a_{j\overline{\sigma }}^{}a_{j\overline{\sigma }}\frac{\nu }{2})$$
(17)   
with $`\nu =<n_i>`$ being the average occupation of the site.

## II The Cooper instability in the occupation-dependent hopping Hamiltonians

The Cooper instability realizes at certain temperature $`T=T_c`$ as a singularity in a two-particle scattering amplitude at zero total momentum. Let’s introduce a function

$$\mathrm{\Gamma }(p_1p_2,\tau \tau ^{})=<T_\tau a_{𝐩_1}(\tau )a_{𝐩_1}(\tau )\overline{a}_{𝐩_2}(\tau ^{})\overline{a}_{𝐩_2}(\tau ^{})>$$
(18)   
where $`\overline{a}_{𝐩\alpha }(\tau )=\mathrm{exp}(H\tau )a_{𝐩\alpha }^{}\mathrm{exp}(H\tau )`$, $`a_{𝐩\alpha }=\mathrm{exp}(H\tau )a_{𝐩\alpha }\mathrm{exp}(H\tau )`$ are the imaginary time $`(\tau )`$ creation and annihilation operators. At $`𝐩_1=𝐩_2`$, $`𝐩_3=𝐩_4`$, the kernel of $`\mathrm{\Gamma }_{\alpha \beta \gamma \delta }`$ is proportional to $`G_{\alpha \beta }^x`$ $`G_{\gamma \delta }^y`$ ($`G`$ is one-electron Green function). We keep notation $`\mathrm{\Gamma }(\mathrm{𝐩𝐩}^{})`$ for such a reduced Green function specifying only momenta $`𝐩=𝐩_1=𝐩_2`$ and $`𝐩^{}=𝐩_3=𝐩_4`$. By assuming temporarily $`V=0`$, this Hamiltonian results in an equation for the Fourier transform $`\mathrm{\Gamma }(p,p^{},\mathrm{\Omega })`$

$$\mathrm{\Gamma }(𝐩,𝐩^{},\mathrm{\Omega })=\mathrm{\Gamma }^0(𝐩,𝐩^{})T\underset{\omega }{}\underset{𝐤}{}\mathrm{\Gamma }^0(𝐩,𝐤)G_\omega (𝐤)G_{\omega +\mathrm{\Omega }}(𝐤)\mathrm{\Gamma }(𝐤,𝐩^{},\mathrm{\Omega })$$
(19)   
corresponding to summation of Feynmann graphs shown in Figure 2. In the above formulas, $`\omega =(2n+1)\pi T`$ and $`\mathrm{\Omega }=2\pi mT`$ ($`n,m`$ integers) are the discrete odd and even frequencies of the thermodynamic perturbation theory . $`G(𝐤,\omega )`$ is a one-particle Green function in a Fourier representation

$$G(𝐤,\omega )=\frac{1}{\xi _ki\omega }.$$
(20)   
Diagrams of Figure 2 are singular since equal momenta of two parallel running lines bring together singularities of both Green functions $`G(𝐤,\omega )`$ and $`G(𝐤,\omega )`$.

6-vertex interaction, Eq.(8), is not generally considered in the theories of strongly-correlated fermionic systems. Such interaction also results in singular diagrams for $`𝐩𝐩`$ scattering shown in Figure 3. Since a closed loop in this figure does not carry any momentum to the vertex, it reduces to the average value of $`\overline{G}`$ which in turn is the average of the number operator, $`<a^{}a>`$. Taking into consideration of such diagrams is equivalent to replacing one of the $`n_i`$’s in Eq.(8) to its thermodynamical average $`\nu =<a_{i\sigma }^{}a_{i\sigma }>`$. Then the $`V`$ term can be added to the renormalized value of $`W`$,

$`WW+{\displaystyle \frac{1}{2}}\nu V`$   
We will check by numeric analysis in Sec. III to which extent such an approximation may be justified.

Solution to Eq.(16) can be received by putting

$$\mathrm{\Gamma }(𝐩,𝐩^{},\mathrm{\Omega })=A(\mathrm{\Omega })+B_1(\mathrm{\Omega })\sigma _𝐩+B_2(\mathrm{\Omega })\sigma _𝐩^{}+C(\mathrm{\Omega })\sigma _𝐩\sigma _𝐩^{}.$$
(21)   
Substituting this expression into Eq.(16) and introducing the quantities

$$S_n(\mathrm{\Omega })=T\underset{\omega }{}\underset{𝐤}{}\sigma _𝐤^nG_\omega (𝐤)G_{\omega +\mathrm{\Omega }}(𝐤)$$
(22)   
we receive a system of coupled equations for $`A,B_1,B_2,C`$

$`\left(\begin{array}{cccc}1+US_0+\stackrel{~}{W}S_1& US_1+\stackrel{~}{W}S_2& 0& 0\\ \stackrel{~}{W}S_0& 1+\stackrel{~}{W}S_1& 0& 0\\ 0& 0& 1+US_0+\stackrel{~}{W}S_1& US_1+\stackrel{~}{W}S_2\\ 0& 0& \stackrel{~}{W}S_0& 1+\stackrel{~}{W}S_1\end{array}\right)`$ $`\left(\begin{array}{c}A\\ B_1\\ B_2\\ C\end{array}\right)=\left(\begin{array}{c}U\\ \stackrel{~}{W}\\ \stackrel{~}{W}\\ 0\end{array}\right)`$ (35)   
where $`\stackrel{~}{W}=W+\frac{1}{2}\nu V`$, which are solved to give

$$A=\frac{U\stackrel{~}{W}^2S_2}{D},B_1=B_2=\frac{\stackrel{~}{W}(1+\stackrel{~}{W}S_1)}{D},C=\frac{\stackrel{~}{W}^2S_0}{D}$$
(36)   
where $`D`$ is a determinant

$$D=\left|\begin{array}{cc}1+US_0+\stackrel{~}{W}S_1& US_1+\stackrel{~}{W}S_2\\ WS_0& 1+WS_1\end{array}\right|.$$
(37)   
The determinant becomes zero at some temperature which means an instability in the two-particle scattering amplitude ($`\mathrm{\Gamma }\mathrm{}`$). This temperature is the superconducting transition temperature $`T_c`$. At $`T_c`$, Eq.(16) is singular, which means that two-particle scattering amplitude gets infinite. Below $`T_c`$, the finite value of $`\mathrm{\Gamma }`$ is established by including the non-zero thermal averages (the order parameters), $`<a_𝐩^{}a_𝐩^{}>`$, $`<a_𝐩a_𝐩>`$. We first analyze the case of non-retarded, non-contraction interaction $`U`$, and after that will consider the effect of the occupation-dependent hopping terms, $`V`$ and $`W`$.

### A Direct non-retarded interaction

Neglecting contraction parameters $`V,W`$, solution to Eq.(16) reduces to

$$\frac{1}{U}=T\underset{\omega }{}\underset{𝐤}{}\frac{1}{\xi _𝐤^2+\omega ^2}$$
(38)   
which after the summation over the discrete frequencies reduces to the conventional BCS equation (at negative $`U`$)

$$\frac{1}{|U|}=\underset{𝐤}{}\frac{12n_𝐤}{2\xi _𝐤},$$
(39)   
with $`n_𝐤=(\mathrm{exp}(\beta \xi _𝐤)+1)^1`$. At finite frequency $`\mathrm{\Omega }`$, Eq.(23) reduces to

$$\mathrm{ln}\frac{T}{T_c}=T\underset{\omega }{}_{E_1}^{E_2}𝑑\xi \frac{i\mathrm{\Omega }}{(\xi ^2+\omega ^2)(\xi +i\omega +i\mathrm{\Omega })}$$
(40)   
where we replaced for simplicity an integration over the Brillouin, zone $`d^3k`$, by the integration over the energy assuming that the density of states near the Fermi energy $`\mu `$ is flat. $`E_1`$ and $`E_2`$ are the lower and upper limits of integration equal to $`4t\mu `$ and $`4t\mu `$, respectively. Such an approximation is not very bad since most singular contribution to integral comes from the point $`\xi _p=0`$ where the integrand is the largest.

Above $`T_c`$, Eq.(25) determines the frequency of the order parameter relaxation . There is a small change in this frequency compared to the BCS model in which limits of the integration $`(E_1,E_2)`$ are symmetric with respect to the Fermi energy, and small in comparison to $`\epsilon _F`$, therefore we briefly discuss it now.

To receive a real-time relaxation frequency, Eq.(25) needs to be analytically continued to a real frequency domain from the discrete imaginary frequencies $`i\omega _n=(2n+1)\pi iT`$ . Using the identity

$$T\underset{\omega }{}\frac{1}{(\omega +i\xi _1)(\omega +i\xi _2)\mathrm{}(\omega +i\xi _n)}=(i)^n\underset{i=1}{\overset{n}{}}\underset{ij}{}\frac{n(\xi _i)}{\xi _i\xi _j}$$
(41)   
where $`n(\xi )`$ is a Fermi function $`n(\xi )=(\mathrm{exp}(\beta \xi )+1)^1`$ gives

$$\mathrm{ln}\frac{T}{T_c}=\frac{i\mathrm{\Omega }}{2}_{E_1}^{E_2}\frac{\mathrm{tanh}\frac{\xi }{2T}}{\xi (2\xi +i\mathrm{\Omega })}𝑑\xi $$
(42)   
where

$$T_c=\frac{2\gamma }{\pi }\sqrt{E_1E_2}\mathrm{exp}\left(\frac{1}{N(\epsilon _F)|U|}\right),\mathrm{ln}\gamma =C=0.577.$$
(43)   
$`C`$ is Euler constant. Analytic continuation is now simple: we change $`\mathrm{\Omega }`$ to $`i(\omega i\delta )`$, $`\delta =+0`$, to receive a function which will be analytic in the upper half plane of complex $`\omega `$, $`Im\omega >0`$. The order parameter relaxation equation becomes

$$\left(\mathrm{ln}\frac{T}{T_c}\frac{\omega }{4}_{E_1}^{E_2}\frac{\mathrm{tanh}\frac{\xi }{2T}}{\xi (\xi \frac{\omega }{2}+i\delta )}𝑑\xi \right)\mathrm{\Delta }=0.$$
(44)   
At $`\omega T_c`$ and $`TT_cT_c`$, the real and imaginary parts of Eq.(29) are easily evaluated to give

$$\left(TT_c\frac{\pi i\omega }{8T_c}+\omega \frac{E_1E_2}{4E_1E_2}\right)\mathrm{\Delta }=0.$$
(45)   
Thus, the order parameter relaxation equation at $`T>T_c`$ becomes

$$(1+i\lambda )\frac{\mathrm{\Delta }}{t}+\mathrm{\Gamma }\mathrm{\Delta }=0$$
(46)   
where

$$\mathrm{\Gamma }=\frac{8}{\pi }(TT_c),\lambda =\frac{2(E_1E_2)}{\pi E_1E_2}T_c.$$
(47)   
In comparison to the BCS theory in which $`E_1=E_2=\omega _D`$ ($`\omega _D`$ is the Debye frequency) and therefore $`\lambda =0`$, we receive the relaxation which has a non-zero “inductive” component, $`i\lambda \mathrm{\Gamma }`$. Typically, $`E_1E_2\epsilon _F`$ and therefore $`|\lambda |`$ is a small quantity. It increases however near the low ($`\nu 1`$) or near the maximal ($`\nu 2`$) occupation where $`E_1`$ or $`E_2`$ become small. Such mode of relaxation is specific to a non-retarded (non-phonon) interaction which is not symmetric near $`\epsilon _F`$ and spans over the large volume of the $`𝐤`$-space rather than is restricted to a narrow energy $`\omega _D\epsilon _F`$ near the Fermi energy.

### B Occupation-dependent hopping instability and relaxation

Neglecting direct interaction, we put $`U=0`$ in Eq.(22) and receive

$$\frac{1}{\stackrel{~}{W}}=S_1(\omega )\pm \sqrt{S_0(\omega )S_2(\omega )}$$
(48)   
where at finite frequency $`\omega `$

$$S_n(\omega )=N(\epsilon _F)T\underset{\omega }{}_{E_1}^{E_2}\left(\frac{\xi +\mu }{t}\right)^n\frac{\mathrm{tanh}\frac{\xi }{2T}}{2\xi \omega +i\delta }𝑑\xi .$$
(49)   
Putting $`\omega =0`$ we receive from Eq.(33) a transition temperature $`T_c`$. The equation has a solution at $`\stackrel{~}{W}<0`$, $`\mu <0`$, or at $`\stackrel{~}{W}>0`$, $`\mu >0`$ (we assume that $`t>0`$). Plus or minus sign is chosen to receive the maximal value of $`T_c`$ (the second solution corresponding to smaller $`T`$, then, has to be disregarded since at $`T<T_c`$ the order parameter will be finite and therefore Eqs.(20)-(22) do not apply). This gives an expression for $`T_c`$

$$T_c=\frac{2\gamma }{\pi }\sqrt{E_1E_2}\mathrm{exp}\left[\frac{E_1E_2}{2|\mu |t}(t|\mu |)+\frac{E_2^2E_1^2}{8\mu ^2}\right]\mathrm{exp}\left(\frac{t}{2|\stackrel{~}{W}|N(\epsilon _F)}\right)$$
(50)   
where $`\mu <0`$, $`\stackrel{~}{W}<0`$ (second exponent is dominating the first one in the weak coupling limit $`\stackrel{~}{W}0`$). Real and imaginary parts of $`S_n(\omega )`$ are calculated at $`\omega T_c`$

$`ImS_n(\omega ){\displaystyle \frac{\pi \omega }{8T_c}}\left({\displaystyle \frac{\mu }{t}}\right)^nN(\epsilon _F)`$ (51)   
$`ReS_n(\omega )={\displaystyle \frac{\omega }{4}}N(\epsilon _F)\left({\displaystyle \frac{\mu }{t}}\right)^n\times \{\begin{array}{cc}\frac{E_2E_1}{E_1E_2},\hfill & n=0\hfill \\ \frac{E_2E_1}{E_1E_2}+\frac{2}{\mu }\mathrm{ln}\frac{\gamma \sqrt{E_1E_2}}{T_c},\hfill & n=1\hfill \\ \frac{E_2E_1}{E_1E_2}+\frac{E_1E_2}{\mu ^2}+\frac{2}{\mu }\mathrm{ln}\frac{\gamma \sqrt{E_1E_2}}{T_c},\hfill & n=2\hfill \end{array}`$ (55)   
Equation for $`\lambda `$ is received with a value larger than the previous one (Eq.(32))

$$\lambda \frac{T_c}{\mu }\left(3\mathrm{ln}\frac{2\gamma \sqrt{E_1E_2}}{\pi T_c}+\frac{2\mu (E_2E_1)}{E_1E_2}+\frac{E_1E_2}{2\mu }\right).$$
(56)   
Eigenvalue equation gives the $`𝐩`$-dependence of the two particle correlator $`\mathrm{\Gamma }(𝐩,𝐩^{})=<a_𝐩^{}a_𝐩^{}a_𝐩^{}a_𝐩^{}>`$ near $`T_c`$

$$\mathrm{\Gamma }(𝐩,𝐩^{})=C\left[S_2S_1(\sigma _𝐩+\sigma _𝐩^{})+S_0\sigma _𝐩\sigma _𝐩^{}\right].$$
(57)   
Since $`C`$ diverges at $`T_c`$, this determines that order parameter becomes macroscopic at $`T<T_c`$. Then, the pair creation operator, $`a_𝐩^{}a_𝐩^{}`$, will almost be a number, i.e., we may decompose Eq.(39) into a product

$$\mathrm{\Delta }_𝐩^{}\mathrm{\Delta }_𝐩=<a_𝐩^{}a_𝐩^{}><a_𝐩^{}a_𝐩^{}>$$
(58)   
and, to be consistent with the $`𝐩`$, $`𝐩^{}`$ dependences, by putting $`\xi _𝐩=\xi _𝐩^{}`$ we receive

$$\mathrm{\Delta }_𝐩=C_1\left(\mathrm{exp}(i\theta /2)\sqrt{S_2(0)}+\mathrm{exp}(i\theta /2)\sqrt{S_0(0)}\right)\mathrm{exp}(i\phi )$$
(59)   
where

$$\mathrm{cos}\theta =S_1(0)/\sqrt{S_0(0)S_2(0)}$$
(60)   
and $`\phi `$ is an overall phase which is irrelevant for a single superconductor but is important for calculating currents in multiple or weakly coupled superconductors. Therefore, system undergoes a pairing transition at temperature found from the Eq.(35). Since the pairs are charged, the state below $`T_c`$ can not be non-superconducting.

We have not calculated the Meissner response but in the following section we present numerical calculation of flux quantization which supports the above statement.

## III Exact diagonalization of the occupation-dependent hopping Hamiltonians in finite cluster

We calculate the ground state energy of a cubic system as shown in Figure 4. A magnetic flux $`\mathrm{\Phi }`$ is produced by a solenoid passing through the cube. Corners of the cube are the lattice sites, which can be occupied by electrons. With the inclusion of the magnetic flux, model Hamiltonian, Eq. 6, becomes

$`H=t{\displaystyle \underset{<ij>\sigma }{}}a_{i\sigma }^{}a_{j\sigma }\mathrm{exp}(i\alpha _{ij})+h.c.+U{\displaystyle \underset{i}{}}n_in_i+`$ (61)   
$`+{\displaystyle \underset{<ij>\sigma }{}}a_{i\sigma }^{}a_{j\sigma }\left[Vn_{i\overline{\sigma }}n_{j\overline{\sigma }}+W(n_{i\overline{\sigma }}+n_{j\overline{\sigma }})\right]\mathrm{exp}(i\alpha _{ij})+h.c.`$ (62)   
where

$$\alpha _{ij}=(2\pi /\mathrm{\Phi }_0)_{𝐫_i}^{𝐫_j}𝐀𝑑𝐥$$
(63)   
and $`\mathrm{\Phi }_0=hc/e`$ is the magnetic flux quantum. Throughout the calculations we take $`t=1`$.

We start with constructing the model Hamiltonian. In a Hilbert space of one electron

$`a=\left(\begin{array}{cc}0\hfill & \hfill 1\\ 0\hfill & \hfill 0\end{array}\right),a^{}=\left(\begin{array}{cc}0\hfill & \hfill 0\\ 1\hfill & \hfill 0\end{array}\right).`$ (68)   
with a basis specified as $`\psi _0=(0,1)`$ for the ground state ($`n=0`$) and $`\psi _1=(1,0)`$ for the excited state ($`n=1`$). In case of $`N`$ states operator of annihilation $`a_n`$ takes the form

$$a_n=s^{n1}ae^{Nn}$$
(69)   
where $`e`$ is the unit matrix and $`s`$ is unitary matrix

$`e=\left(\begin{array}{cc}1\hfill & \hfill 0\\ 0\hfill & \hfill 1\end{array}\right),s=\left(\begin{array}{cc}1\hfill & \hfill 0\\ 0\hfill & \hfill 1\end{array}\right)`$ (74)   
and $``$ stands for the Kronecker matrix multiplication. Explicitly, we have

$`a_1`$ $`=`$ $`aeee\mathrm{}e`$   
$`a_2`$ $`=`$ $`saee\mathrm{}e`$   
$`\mathrm{}`$   
$`a_N`$ $`=`$ $`sss\mathrm{}sa`$   
Thus, for example, for two states

$`a_1=\left(\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 0& 0\end{array}\right),a_2=\left(\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right).`$ (83)   
These matrices, which are annihilation operators, and corresponding Hermitian conjugate matrices, which are the creation operators, satisfy the Fermi anti-commutation relation. These operators are sparse matrices with only $`N/2`$ non-zero elements, which are equal to $`\pm 1`$. Next we solve the Schrödinger equation $`H\psi =E\psi `$. We implemented a novel algorithm for solving such sparse systems, which will be described elsewhere.

The cubic cluster within the Hubbard Hamiltonian and no external flux applied to the system was studied previously by Callaway et. al. . Quantum Monte Carlo methods applicable to large systems within the Hubbard model (both attractive and repulsive), but not the occupation-dependent hopping Hamiltonians, are reviewed in a paper of Dagotto .

### A The number parity effect

Superconductivity reveals itself in the lowering of the ground state energy as electrons get paired. Therefore the energy needs to be minimal for even number of electrons $`n`$ and will attain a larger value when $`n`$ is odd. We consider a “gap” parameter 

$$\mathrm{\Delta }_l=E_{2l+1}\frac{1}{2}\left(E_{2l}+E_{2l+2}\right)$$
(84)   
as a possible “signature” of superconductivity (where $`E_m`$ corresponds to the ground state energy for $`m`$ fermions). For all interaction parameters set to zero ($`U=V=W=0`$), no sign of pairing is observed. To check our analytic results of Sec. IIb and the argument following Eq.(34), we calculated $`\mathrm{\Delta }`$ above and below the half-filling ($`n=8`$ in case of cubic cluster). Below the half-filling chemical potential is negative ($`\mu <0`$) and above the half-filling it is positive ($`\mu >0`$). We first checked that the $`W0^+`$, $`W0^{}`$ and $`V0^+`$, $`V0^{}`$ calculation is consistent with an exact solution available for a non-interacting system of $`n`$ electrons.

We then test our program for the case of negative-$`U`$ Hubbard Hamiltonian ($`U<0,V=0,W=0`$) which is known to be superconducting (e.g. Refs. 22,23). Positive-$`U`$ Hubbard model does not show any sign of superconductivity, in disagreement with some statements in the literature . Our calculations can not disprove the (possible) non-pairing mechanisms of superconductivity but these seem to be unlikely models for the problem of superconductivity in oxides which clearly shows pairing of electrons (holes) in the Josephson effect and in the Abrikosov vortices. The relation $`2eV=\mathrm{}\omega `$ is justified in the first case  and flux quantum of a vortex is $`hc/2e`$ in the second , both with the value of the charge equal to twice the electronic charge, $`e`$.

Figure 5 shows the dependence of the ground state energy upon the number of particles in case of negative-$`U`$ and positive-$`U`$ Hubbard models assuming $`V=0`$ and $`W=0`$. Such dependences are typical for any value of $`|U|`$. There clearly is the pairing effect when $`U<0`$ and there is no sign of pairing at $`U>0`$.

Tests for pairing in the contraction $`V`$, $`W`$-models ($`V0,U=W=0`$ and $`W0,U=V=0`$, respectively) are shown in Figs. 6,7. The results are in agreement with our perturbative calculation of Sec. II and with its extension for the intermediate and strong coupling limits $`|V|t`$, $`|W|t`$. Since chemical potential is negative below the half-filling and positive above the half-filling, there is no pairing in the former case ($`\stackrel{~}{W}0^+`$) and there is a sign of pairing in the latter case ($`\stackrel{~}{W}>0`$), in accord with the value of the effective coupling constant $`\stackrel{~}{W}=W+\frac{1}{2}\nu V`$. Similarly, for $`\stackrel{~}{W}0^{}`$ below the half-filling there is a sign of pairing ($`\mathrm{\Delta }0`$) while above the half-filling there is no pairing. These results are summarized in Table 1.

For larger values of the interaction parameters, the perturbative results do not remain applicable anymore. Figure 8b shows the dependence of the parity gap $`\mathrm{\Delta }`$ on the strength of the interaction. From Figure 8 it is understood that the $`W`$ interaction introduces a “signature” of pairing in a similar way as the negative-$`U`$ interaction does. The possibility of the “contraction” pairing has been investigated formerly in the papers .

### B Flux quantization

Flux quantization is another signature of superconductivity which is a consequence of the Meissner effect. We also tested for the periodicity of the energy versus flux dependence with the period $`\mathrm{\Phi }_1=hc/2e`$ as compared to the period $`\mathrm{\Phi }_0=hc/e`$ in the non-interacting system . Unfortunately, the even harmonics of $`\mathrm{\Phi }_0`$-periodic dependence of the ground state energy (and related to it, the harmonics of the persistent current $`J=E/\mathrm{\Phi }`$ ) may simulate the pairing in a non-superconductive system. Small-size (mesoscopic) system can mask the superconducting behavior . Flux quantization in Hubbard Hamiltonians was studied formerly in Refs. 29-31.

We first demonstrate the behavior of the ground state energy with respect to flux, Figure 9. A characteristic feature of mesoscopic system suggests that addition of one extra particle to the system changes the sign of the derivative of the ground state energy with respect to magnetic flux at $`\mathrm{\Phi }=0`$. That is, depending on the parity of the number of particles and on the number of sites, system can change from paramagnetic to diamagnetic state or vice versa. But this behavior is not always observed for the cubic geometry studied. Except the sign change from $`n=2`$ to $`n=3`$ and from $`n=7`$ to $`n=8`$, no such behavior is seen. As mentioned above, however, the $`\mathrm{\Phi }_1`$-periodic component of the $`E(\mathrm{\Phi })`$ dependence begins to appear at the higher value of $`n`$ (Figure 9,c). For both contraction parameters equal to zero, i.e. $`W=V=0`$, we observe appearance of the $`hc/2e`$-periodic component for some values of $`U`$ (Figure 10). Even for positive (repulsive) values of $`U`$, it is possible to see a local minimum appearing at $`\mathrm{\Phi }=hc/2e`$ (Figure 10b). This is in agreement with the authors’ previous works . But this minimum, which does not lead to an exact periodicity of the ground state energy with a period $`\mathrm{\Phi }_0/2`$, should not be attributed to superconductivity, this is rather a characteristic behavior in mesoscopic systems.

For $`U<0`$ (while $`W=V=0`$), the expected mesoscopic behavior, that is the change of the sign of the slope of ground state energy at $`\mathrm{\Phi }=0`$, starts to demonstrate itself (Figure 11). But this happens at sufficiently large absolute values of (negative) $`U`$. For other values of $`U`$, however, there is no such change.

More pronounced $`hc/2e`$-periodic components are observed with the introduction of non-zero interaction parameters. The role of $`W`$ on the ground state energy, when both $`U`$ and $`V`$ are zero, is shown in Figure 12. Meanwhile setting both $`U`$ and $`W`$ to zero and observing the effect of the non-zero $`V`$ shows that $`V`$ does not play a role as significant as the other two interaction parameters do. There is not much difference in the behavior of the ground state energy upon magnetic flux between the zero and non-zero $`V`$ (for example $`V=1`$) cases.

## IV Conclusions

We studied the peculiarity of electron conduction in systems in which conduction band is derived from the atomic shells with a small number of electrons ($`N_e`$) in an atom. Such materials may include oxygen ($`N_e=8`$) in the oxides, carbon ($`N_e=6`$) in borocarbides (e.g. $`LuNi_2B_2C`$), hydrogen ($`N_e=1`$) in some metals (e.g., $`PdH`$). Some materials of this kind are superconductors. It was argued that the Coulomb effects within the atoms strongly influence the inter-atom wave function overlap between the atomic sites and therefore the electron hopping amplitude between the sites. The phenomenology of such conduction mechanism results in a novel, to the conventional solid state theory, Hamiltonians called the occupation-dependent-hopping (or contraction) Hamiltonians, specified with the two coupling parameters $`V`$, $`W`$. We then attempted a study of superconductivity in such systems within the BCS-type approach assuming the Cooper pairing of electrons. The weak-coupling limit allows determination of the range of parameters $`V`$, $`W`$ values and also of the in-site Coulomb interaction $`U`$ value which show the Cooper instability. The strong-coupling limit was addressed by a numeric calculation on finite clusters using the novel algorithm (of non-Lanczos type) for eigenvalues of large sparse matrices. One of the results of this numeric calculation was that the positive-$`U`$ Hubbard model, sometimes believed to be a candidate for high-$`T_c`$ superconductivity, does not comply with the goal.

This work was partially supported by the Scientific and Technical Research Council of Turkey (TÜBİTAK) through the BDP program.