in case of the in-plane bipolaron [@alekor2], Fig.2.

-0.5mm

The Bloch bands of these bipolarons are obtained using the canonical transformation, described above, projecting the transformed Hamiltonian, Eq.(\[trans\]), onto a reduced Hilbert space containing only empty or doubly occupied elementary cells [@alebook]. The wave function of the apex bipolaron localized, say in cell ${\bf m}$ is written as $$|{\bf m}\rangle =\sum_{i=1}^{4}A_{i}c_{i}^{\dagger
}c_{apex}^{\dagger }|0\rangle ,$$ where $i$ denotes the $p_{x,y}$ orbitals and spins of the four plane oxygen ions in the cell, Fig.1, and $c_{apex}^{\dagger }$ is the creation operator for the hole in one of the three apex oxygen orbitals with the spin, which is the same or opposite to the spin of the in-plane hole depending on the total spin of the bipolaron. The probability amplitudes $A_{i}$ are normalized by the condition $|A_{i}|=1/2,$ if four plane orbitals $%
p_{x1},p_{y2},p_{x3}$ and $p_{y4}$ are involved, or by $|A_{i}|=1/\sqrt{2}$ if only two of them are relevant. Then a matrix element of the Hamiltonian Eq.(\[trans\]) describing the *bipolaron* tunnelling to the nearest neighbor cell ${\bf
m+a}$ is found as $$t_b=\langle {\bf m}|\tilde{H}|{\bf m}+{\bf a}\rangle
=|A_{i}|^{2}T_{pp^{\prime }}^{apex}e^{-g^{2}},$$ where $T_{pp^{\prime }}^{apex}e^{-g^{2}}$ is a *single polaron* hopping integral between two apex ions. The inter-site bipolaron tunnelling 