# Two-particle bound states on a lattice

Pavel E. Kornilovitch<EMAIL_ADDRESS>
###### Abstract

Two-particle lattice states are important for physics of magnetism,
superconducting oxides, and cold quantum gases. The quantum-mechanical lattice
problem is exactly solvable for finite-range interaction potentials. A two-
body Schrödinder equation can be reduced to a system of linear equations whose
numbers scale with the number of interacting sites. For the simplest cases
such as on-site or nearest-neighbor attractions, many pair properties can be
derived analytically, although final expressions can be quite complicated. In
this work, we systematically investigate bound pairs in one-, two-, and three-
dimensional lattices. We derive pairing conditions, plot phase diagrams, and
compute energies, effective masses and radii. Along the way, we analyze
nontrivial physical effects such as light pairs and the dependence of binding
thresholds on pair momenta. At the end, we discuss the preformed-pair
mechanism of superconductivity and stability of many-pair systems against
phase separation. The paper is a combination of original work and pedagogical
tutorial.

###### keywords:

Two-particle states , bound states , lattice , preformed pairs ,
superconductivity

††journal: Annals of Physics

Department of Physics, Oregon State University, Corvallis, Oregon 97331, USA

###### Contents

  1. 1  Introduction
  2. 2  General theory
    1. 2.1  Model
    2. 2.2  Unsymmetrized solution
    3. 2.3  Symmetrized solution. Spin singlets
    4. 2.4  Anti-symmetrized solution. Spin triplets
    5. 2.5  Pair size
  3. 3  Negative-$U$ Hubbard model
    1. 3.1  General expressions
    2. 3.2  1D. One dimensional chain
    3. 3.3  2D. Square lattice
    4. 3.4  2D. Triangular lattice
    5. 3.5  3D. Simple cubic lattice
    6. 3.6  3D. Tetragonal lattice
    7. 3.7  3D. Body-centered cubic (BCC) lattice
    8. 3.8  3D. Face-centered cubic (FCC) lattice
  4. 4  1D. $UV$ model on the one-dimensional chain
    1. 4.1  Singlet states
    2. 4.2  Light bound pairs
    3. 4.3  Triplet states
  5. 5  2D. $UV$ model on the square lattice
    1. 5.1  Singlet states. $\Gamma$-point
    2. 5.2  Singlet states. Arbitrary momentum ${\bf P}$
    3. 5.3  Triplet states
    4. 5.4  Longer-range attractions
  6. 6  2D. $UV$ model on the rectangular lattice
    1. 6.1  Singlet states. $\Gamma$-point
    2. 6.2  Triplet states. $\Gamma$-point
  7. 7  2D. $UV$ model on the triangular lattice
    1. 7.1  General dispersion relations
    2. 7.2  Pairing thresholds at the $\Gamma$-point
  8. 8  3D. $UV$ model on the simple cubic lattice
    1. 8.1  Singlet states
    2. 8.2  Triplet states
  9. 9  3D. Tetragonal $UV$ model
    1. 9.1  Singlet states. Pairing thresholds at the $\Gamma$-point
    2. 9.2  Triplet states. Pairing thresholds at the $\Gamma$-point
    3. 9.3  Bose-Einstein condensation of real-space pairs
  10. 10  Light pairs in the strong coupling limit
    1. 10.1  Staggered square planes
    2. 10.2  Square $UV$ model with resonant attraction, $U=V$
    3. 10.3  Generalization to other lattices with resonant attraction, $U=V$
  11. 11  More than two particles
  12. 12  Summary and discussion
    1. 12.1  Common properties of lattice bound pairs
      1. 12.1.1  Mathematically exact solutions
      2. 12.1.2  Pair binding energy depends on its momentum
      3. 12.1.3  Light pairs
      4. 12.1.4  Anisotropic 3D pairs
      5. 12.1.5  Stability of fermion pairs against phase separation
    2. 12.2  Relevance to cuprate superconductivity
    3. 12.3  Future outlook
  13. A  Omissions and extensions
    1. A.1  Miscellaneous two particle problems
    2. A.2  Extension to non-Bravais lattices and multi-orbital models
  14. B  Green’s functions of the 2D square and rectangular lattices
    1. B.1  Definitions
    2. B.2  Basic integrals for the square case, $\alpha=\beta$
    3. B.3  Recurrence relations for all $M_{nm}$ in the square case, $\alpha=\beta$
    4. B.4  $M^{\rm sq}_{nm}-M^{\rm sq}_{00}$ at threshold. Square case, $\alpha=\beta$
    5. B.5  Basic integrals in the rectangular case, $\alpha\neq\beta$
    6. B.6  Recurrence relations for the rectangular case, $\alpha\neq\beta$
    7. B.7  Table of $M_{nm}-M_{00}$ at threshold in the rectangular case, $\alpha\neq\beta$
  15. C  $UV$ model on the 2D triangular lattice
    1. C.1  Evaluation of basic integral $M^{\rm tr}_{00}$, Eq. (15), for the triangular lattice
    2. C.2  Numerical evaluation of $M^{\pm}_{\bf nm}$ for the triangular lattice
    3. C.3  Table of $M^{\rm tr}_{\bf nm}-M^{\rm tr}_{\bf 00}$ at threshold at the $\Gamma$-point
    4. C.4  Application of group theory
  16. D  $UV$ model on the simple cubic and tetragonal lattices
    1. D.1  Green’s functions, Eq. (219), in the simple cubic case, $\alpha=\beta=\gamma$
    2. D.2  Green’s functions, Eq. (219), in the tetragonal case, $\alpha=\beta\neq\gamma$
    3. D.3  Green’s functions, Eq. (219), in the orthorhombic case, $\alpha\neq\beta\neq\gamma$
    4. D.4  Determinant of Eq. (234)
    5. D.5  Derivation of Eq. (237)
  17. E  $UV$ model on the BCC lattice
    1. E.1  Derivation of Eq. (102)

##  1  Introduction

In 1986, Daniel Mattis published a review on few-body lattice problems [1].
The motivation for that work was deeper understanding of magnetism. Many
successful models of magnetism are formulated in terms of localized spins
arranged in regular lattices [2]. Exact results obtained for two and three
interacting magnons provided valuable physical insights. Mattis also mentioned
other areas that could benefit from similar analysis, specifically excitons
and superconductivity (‘$\ldots$ do Cooper pairs have bound states of “Cooper
molecules”?’, [1], p. 362). At about the same time, high-temperature
superconductivity (HTSC) was discovered by Bednorz and Müller [3], which
generated enormous interest in unconventional theories of superconductivity.
The unusual HTSC properties such as a low carrier density, a short coherence
length [4], and a Bose-gas-like scaling of the magnetic penetration depth with
$T_{c}$ [5, 6] revived the pre-BCS proposal [7, 8, 9, 10] that
superconductivity could be a Bose-Einstein condensation (BEC) of charged
bosons formed by binding of electrons or holes into real-space pairs. Many
properties of superconducting oxides have been shown to be well described by
charged Bose-gas phenomenology [11, 12, 13, 14, 15, 16]. More recently,
normal-state pairs were observed in iron-based superconductors [17, 18] and in
the shot noise in copper oxide junctions [19]. The debates about the preformed
pairs mechanism of superconductivity and its relevance to HTSC and pseudogap
physics continue today [20].

Physics of real-space pairs is tightly linked with the pairing mechanism. In
superconducting oxides, main candidates are: strong electron-phonon
interaction [21], the Jahn–Teller effect [3, 22, 23, 24], and spin
fluctuations [25], although other mechanisms and combinations have been
proposed [11, 26, 27, 28, 29]. However, the complexity of the problem suggests
splitting one big puzzle into smaller ones. By postulating a simple
phenomenological attraction of some kind, one can investigate the BEC-BCS
crossover [16, 30, 31], the physics of pseudogap, the role of anisotropy [32],
phase separation [33, 34, 35], electrodynamics, and other nontrivial topics,
all without arguing about the specific nature of pairing interaction.

An added benefit of this approach is that the two-particle lattice problem is
exactly solvable for a wide class of non-retarded interaction potentials,
which provides a welcome rigor and often an analytical formula. The binding
threshold, pair energy, dispersion, effective mass $m^{\ast}_{p}$, effective
radius $r^{\ast}_{p}$, and wave function: all can be determined without
approximations. These properties provide important physical insights. For
example, mass anisotropy of pairs differs from the bare anisotropy of the
member particles, see, e.g., Sections 3.6 and 9.3. This relationship may be
helpful in understanding the anisotropy of transport and electromagnetic
properties of HTSC. Additionally, pair wave function is directly proportional
to a macroscopic superconducting order parameter [36], In particular, both
share the same orbital symmetry. Knowledge of the exact pair wave function
helps elucidate the relationship between the order parameter and other
elements of the system. For example, it was proposed [37] that correlation-
induced diagonal hopping may produce a $d$-symmetric ground-state pair.

Figure 1: Evolution of BEC of real-space pairs with density. $n_{\rm
cp}=\Omega^{-1}_{p}$ is the close-packing density, $T_{c}$ is a “critical
temperature of phase transition”. $T^{\ast}_{\rm BEC}$ is determined by the
mass and size of a single pair, see Eq. (2). The symmetry of a macroscopic
order parameter is defined by the pair wave function [36].

A striking application of two-particle properties to unconventional
superconductivity comes from analyzing the BEC temperature of pairs:

$T_{\rm BEC}=C\,\frac{\hbar^{2}}{m^{\ast}_{p}}\,n^{2/3}_{p}\>.$ (1)

Here, $C$ is a numerical coefficient and $n_{p}$ is the pair number density
assumed to be known from chemical composition. $m^{\ast}_{p}$ is the pair mass
provided by an exact solution. Consider system’s evolution with $n_{p}$ shown
in Fig. 1. (A) At low $n_{p}$, the average distance between pairs is larger
than their size, and Eq. (1) applies. (B) With $n_{p}$ increasing, the system
reaches “close-packing” when the pairs begin to overlap. The corresponding
density is approximately given by an inverse pair volume, $n_{\rm
cp}=\Omega^{-1}_{p}=(r^{\ast}_{p})^{-3}$. The pairs interact strongly but the
phase transition is still of BEC type with a transition temperature
approximately given by

$T_{\rm BEC}(n_{\rm cp})=T^{\ast}_{\rm
BEC}=C\,\frac{\hbar^{2}}{m^{\ast}_{p}}\frac{1}{r^{\ast
2}_{p}}\propto\frac{1}{m^{\ast}_{p}\,r^{\ast 2}_{p}}\>.$ (2)

(C) Upon further increase of density, the pairs overlap much more, with the
average distance smaller than the pair size. The constituent fermions begin to
form a Fermi sea [38], and the phase transition gradually shifts to the BCS
type. Equation (1) no longer applies and $T_{c}$ begins to fall. Thus, the
“maximal attainable” transition temperature is given by Eq. (2) [32, 39, 40].
Remarkably, the latter contains only the properties of a single pair, with
both $m^{\ast}_{p}$ and $r^{\ast}_{p}$ supplied by the exact solution. Then,
various effects on $T^{\ast}_{\rm BEC}$ may be studied rigorously. This
methodology was applied, for example, to understanding the effects of
interlayer hopping in Ref. [32]. This topic is the subject of Section 9.3.

Early developments of short-range attractive models of superconductivity [41,
42, 43] were summarized by Micnas, Ranninger, and Robaszkiewicz [11]. Those
authors considered two groups of models. The first group included on-site
attraction, and was essentially derived from the attractive Hubbard model.
Those models were useful to follow the BCS-BEC crossover, to understand
system’s thermodynamics, electrodynamics and other properties. At the same
time, those models were too simplistic to describe real superconducting
materials. It is hard to come up with a physical mechanism strong enough to
overcome Coulomb repulsion between charge carriers and still keep the pairs
mobile. The second group of models included intersite attraction and were more
realistic. Moving attraction to finite distances allowed keeping a strong
repulsive core (Hubbard repulsion) that could model a screened Coulomb
repulsion. In addition, intersite models add the possibility of
antiferromagnetic ordering, of $p$\- and $d$-pairing, and make a better
connection to lattice models of magnetism. Such models with on-site repulsion
and intersite attraction will be the main focus of the present work. They will
be referred hereafter as “$UV$ models”. These models continue to be actively
investigated today [44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56].

Even nearest-neighbor attraction is an approximation to inter-electron or
inter-hole attractions of real materials. The latter results from mediation by
an intermediary subsystem such as phonons, and typically extend beyond nearest
neighbors [57, 58]. Such two-particle models are still exactly solvable but
the complexity of solution increases sharply with the range. One example is
analyzed in Section 5.4. What is important, however, is that $UV$ models
comprise the two most essential elements: some repulsion $U$ representing
strong correlations and some attraction $V$ representing a mediating
subsystem. As such, the $UV$ potential is the simplest representative of an
entire class of potentials that possess common properties. For example, all
such models have a threshold function $V_{\rm cr}(U)$ that separates bound and
unbound states and saturates in the $U\rightarrow\infty$ limit. Effects of
other factors, such as anisotropy or pair motion, on the delicate balance
between the two forces can be studied within the $UV$ model, at least
qualitatively.

Intersite attraction also appear within another popular theory of
unconventional superconductivity based on the t-J model [59, 33, 60, 61, 34,
62, 63]. In the dilute limit (many holes and few electrons), the exchange
interaction is equivalent to a nearest-neighbor attraction and the system
reduces to a gas of bound pairs. As the t-J model has mostly been studied on
square lattices, this case will be covered in Section 5.

In the last two decades, optical lattices and cold atoms emerged as another
physical realization of models with short-range attraction [64, 65, 66, 67].
Whereas in the cuprates the $UV$ model is an approximation to real inter-
particle potentials, in optical lattices it can be precisely engineered and
studied in pure form. The onsite interaction is controlled via Feshbach
resonances [68] and can be made either repulsive [69] or attractive [70]. The
intersite interaction can be controlled either by exciting dressed Rydberg
atoms to large quantum numbers [71] or via proper alignment of dipolar quantum
gases [72]. Precise manipulation of few particles in optical traps [73, 74,
75] and BEC of molecules in an attractive Fermi gas [76] have been
demonstrated. Local pairing can now be measured directly using gas microscopy
[77, 78].

Bound states of two spin waves also appear in models of quantum magnetism [2,
79, 80, 81, 82, 83] including frustrated magnetics [84, 85, 86].

Basic properties of two-body states in $UV$ models were discussed by Micnas,
Ranninger, and Robaszkiewicz [11]. In particular, conditions for pair
formation and binding energies were found for 1D, 2D square, and 3D simple
cubic $UV$ models. Since then, the field has seen several developments. For
example, it was realized that threshold of pair formation depends on the pair
momentum [87]. $UV$ models have been solved for triangular [88], tetragonal
[32], BCC [89], and FCC [90] lattices. New analytical results obtained for
cubic Watson integrals by Joyce and coworkers [91, 92, 93, 94, 95, 96, 97, 98]
suggest revisiting the cubic models for deeper analysis. Given also that the
HTSC puzzle remains largely unsolved and the recent observation of real-space
pairs in the normal state [17, 18, 19], a fresh review of two-body states on a
lattice seems to be worthwhile.

The purpose of this work is twofold. First, we collect results obtained for
different lattices [87, 32, 99, 100, 88, 101, 102, 89, 90] and develop them
systematically in one place and with unified formalism. Second, we present a
large body of new results that were developed over the last 30 years but
remained unpublished until now. Most of the material in the main text is
original, including the general theory of Section 2. If a particular result or
a formula is not accompanied by a specific reference to prior work, one should
assume the result is published for the first time. The paper has a pedagogical
side as well as we included in Appendixes some textbook material to make
exposition self-contained. In particular, we collected available information
on lattice Green’s functions in different geometries. We also included in C.4
an example of utilizing group theory to separate states of different orbital
symmetries. In general, we find the two-particle problem to be an excellent
primer on non-relativistic quantum mechanics that teaches wave-function
symmetry, emergence of complex band dispersions, the Galilean invariance and
lack thereof, multi-component wave functions, scattering states, and other
topics.

The paper is organized as follows. After formulating the general theory
(Section 2), we analyze bound states in the attractive Hubbard model in
different dimensions (Section 3), and then the $UV$ model in 1D (Section 4),
2D (Sections 5, 6, and 7), and 3D (Sections 8 and 9). The next two sections
are devoted to two special topics: light pairs (Section 10), and the stability
of pairs against phase separation (Section 11). Section 12.1 summarizes the
most important common physical properties of lattice bound states. Sections
9.3 and 12.2 discuss relevance to HTSC. Extensive Appendixes contain
additional information intended for subject matter expects. A.1 briefly
summarizes two-particle problems not covered in this work in detail. A.2
explains how the present formalism can be extended to multi-orbital models
[39, 103, 104]. The latter are excluded from this review because they are not
yet ready for systematic exposition. The other Appendixes provide explicit
expressions for analytically known lattice Green’s functions, recipes for
efficient numerical integration when Green’s functions are not known
analytically, a group theory example, and details of several most cumbersome
derivations.

Figure 2: Illustration of a $UV$ model on the square lattice with nearest-
neighbor attraction and hopping. The cartoon in the center illustrates an
$s$-symmetric wave function in which the on-site amplitude is suppressed
relative to the inter-site ones.

##  2  General theory

###  2.1  Model

We begin with formulating an underlying model. We restrict consideration to
Bravais lattices only. (Two-body problems in more complex lattices can also be
solved [103, 104, 39] but the resulting algebraic expressions are much more
cumbersome, see A.2.) The inversion symmetry of Bravais lattices allows for
easier separation of singlet and triplet pair states, which greatly reduces
the complexity of equations. Second, we will consider only finite-range
potentials. The two-particle Schrödinger equation reduces to an algebraic
matrix equation with a size equal to the number of nonzero elements in the
inter-particle potential. An infinite-range potential will lead to an
infinitely large matrix and the advantages of an exact solution will be lost.
Even for finite ranges, the complexity of the matrix solution grows rapidly
with the interaction range. In practice, only solutions with short-range
potentials are simple enough to produce analytical results. Most space in this
paper will be devoted to models with on-site and nearest-neighbor
interactions, although one case of a longer range interaction will be
discussed in Section 5.4. A typical $UV$ model is illustrated in Fig. 2.
Third, we will consider only spin-$\frac{1}{2}$ fermions. Their coordinate
wave function can be either symmetric or antisymmetric; so the solutions will
also cover spin-0 bosons and spinless fermions. Finally, bound states will be
the main focus, and scattering states will not be discussed.

A second-quantized model Hamiltonian reads

$\hat{H}=-\sum_{{\bf m},\bar{\bf b},\sigma}t_{\bar{\bf b}}\,c^{\dagger}_{{\bf
m}+\bar{\bf b},\sigma}c_{{\bf m}\sigma}+\frac{U}{2}\sum_{\bf m}\hat{n}_{\bf
m}\left(\hat{n}_{\bf m}-1\right)+\frac{1}{2}\sum_{{\bf m},{\bf b}}V_{\bf
b}\,\hat{n}_{{\bf m}+{\bf b}}\,\hat{n}_{\bf m}\>.$ (3)

Here $\hat{n}_{\bf m}=\sum_{\sigma}\hat{n}_{{\bf m}\sigma}=c^{\dagger}_{{\bf
m}\uparrow}c_{{\bf m}\uparrow}+c^{\dagger}_{{\bf m}\downarrow}c_{{\bf
m}\downarrow}$ is the fermion number operator on site ${\bf m}$. $\bar{\bf b}$
and ${\bf b}$ are hopping and interaction neighbor vectors, respectively. The
first term in Eq. (3) is kinetic energy of free fermions defined by spin-
independent hopping integrals $t_{\bar{\bf b}}$. We will invariably consider
only negative hopping integrals, which is reflected by explicitly writing a
negative sign in front of the sum. Therefore, $t_{\bar{\bf b}}>0$ and
$t_{-\bar{\bf b}}=t_{\bar{\bf b}}$ for all $\bar{\bf b}$. The energy of atomic
orbitals is used as a zero energy, and the corresponding term is not written.
The second term is the on-site (Hubbard) repulsion with amplitude $U$. Because
of the property $\hat{n}^{2}_{\sigma}=\hat{n}_{\sigma}$, it is equivalent to
the usual form $U\sum_{\bf m}\hat{n}_{{\bf m}\uparrow}\hat{n}_{{\bf
m}\downarrow}$. The last term in Eq. (3) represents nearest-neighbor
interaction, where $V_{-{\bf b}}=V_{\bf b}$. The prefactor $\frac{1}{2}$ is
included to compensate double-counting. Equation (3) is written in the $UV$
form to emphasize the special role of on-site interaction $U$. For most of the
paper, we will set the lattice spacing to one, $a=1$, only restoring $a$ in
places where it provides additional physical insights.

###  2.2  Unsymmetrized solution

In first quantization, two-body wave function $\Psi({\bf m}_{1},{\bf m}_{2})$
must satisfy the Schrödinger equation:

$-\sum_{\bar{\bf b}}t_{\bar{\bf b}}\left[\Psi({\bf m}_{1}+\bar{\bf b},{\bf
m}_{2})+\Psi({\bf m}_{1},{\bf m}_{2}+\bar{\bf b})\right]+U\,\delta_{{\bf
m}_{1},{\bf m}_{2}}\Psi({\bf m}_{1},{\bf m}_{2})+\sum_{\bf b}V_{\bf
b}\,\delta_{{\bf m}_{1}-{\bf m}_{2},{\bf b}}\Psi({\bf m}_{1},{\bf
m}_{2})=E\,\Psi({\bf m}_{1},{\bf m}_{2})\>,$ (4)

where $E$ is the total energy. Equation (4) is converted in momentum space
using the transformation

$\displaystyle\Psi({\bf m}_{1},{\bf m}_{2})$
$\displaystyle=\frac{1}{N}\sum_{{\bf k}_{1}{\bf k}_{2}}\psi_{{\bf k}_{1}{\bf
k}_{2}}\,e^{i{\bf k}_{1}{\bf m}_{1}+i{\bf k}_{2}{\bf m}_{2}}\,,$ (5)
$\displaystyle\psi_{{\bf k}_{1}{\bf k}_{2}}$
$\displaystyle=\frac{1}{N}\\!\sum_{{\bf m}_{1}{\bf m}_{2}}\\!\\!\Psi({\bf
m}_{1},{\bf m}_{2})\,e^{-i{\bf k}_{1}{\bf m}_{1}-i{\bf k}_{2}{\bf m}_{2}}\,,$
(6)

where $N$ is the total number of lattice unit cells. A transformed equation
reads

$(E-\varepsilon_{{\bf k}_{1}}-\varepsilon_{{\bf k}_{2}})\,\psi_{{\bf
k}_{1}{\bf k}_{2}}=U\frac{1}{N}\sum_{\bf q}\psi_{{\bf q},{\bf k}_{1}+{\bf
k}_{2}-{\bf q}}+\frac{1}{N}\sum_{{\bf b}{\bf q}}V_{\bf b}\,e^{i({\bf q}-{{\bf
k}_{1}}){\bf b}}\,\psi_{{\bf q},{\bf k}_{1}+{\bf k}_{2}-{\bf q}}\>.$ (7)

Here,

$\varepsilon_{\bf k}=-\sum_{\bar{\bf b}}t_{\bar{\bf b}}\,e^{i{\bf k}\bar{\bf
b}}$ (8)

is the one-particle dispersion of the model. The right-hand-side of Eq. (7) is
a linear combination of quantities

$\displaystyle\Phi_{\bf 0}({\bf k}_{1}+{\bf k}_{2})$ $\displaystyle=\Phi_{\bf
0}({\bf P})\equiv\frac{1}{N}\sum_{\bf q}\psi_{{\bf q},{\bf k}_{1}+{\bf
k}_{2}-{\bf q}}=\frac{1}{N}\sum_{\bf q}\psi_{{\bf q},{\bf P}-{\bf q}}\>,$ (9)
$\displaystyle\Phi_{\bf b}({\bf k}_{1}+{\bf k}_{2})$ $\displaystyle=\Phi_{\bf
b}({\bf P})\equiv\frac{1}{N}\sum_{\bf q}e^{i{\bf q}{\bf b}}\,\psi_{{\bf
q},{\bf k}_{1}+{\bf k}_{2}-{\bf q}}=\frac{1}{N}\sum_{\bf q}e^{i{\bf q}{\bf
b}}\,\psi_{{\bf q},{\bf P}-{\bf q}}\>,$ (10)

where

${\bf P}={\bf k}_{1}+{\bf k}_{2}\>,$ (11)

is the total lattice momentum of the two fermions. It is critically important
for the existence of an exact solution that $\Phi_{\bf 0}$ and $\Phi_{\bf b}$
are functions of only one argument ${\bf P}$ rather than two separate
arguments ${\bf k}_{1}$ and ${\bf k}_{2}$. Utilizing the definitions, Eqs. (9)
and (10), the wave function is expressed from Eq. (7) as follows

$\psi_{{\bf k}_{1}{\bf k}_{2}}=\frac{U}{E-\varepsilon_{{\bf
k}_{1}}-\varepsilon_{{\bf k}_{2}}}\,\Phi_{\bf 0}({\bf P})+\sum_{\bf b}V_{\bf
b}\,\frac{e^{-i{\bf k}_{1}{\bf b}}}{E-\varepsilon_{{\bf
k}_{1}}-\varepsilon_{{\bf k}_{2}}}\,\Phi_{\bf b}({\bf P})\>.$ (12)

Substitution of Eq. (12) back in Eq. (9) yields a system of linear algebraic
equations for $\Phi$:

$\displaystyle\Phi_{\bf 0}({\bf P})$ $\displaystyle=-UM_{\bf 00}\,\Phi_{\bf
0}({\bf P})-\sum_{{\bf b}^{\prime}}V_{{\bf b}^{\prime}}\,M_{{\bf 0}{\bf
b}^{\prime}}(E,{\bf P})\,\Phi_{{\bf b}^{\prime}}({\bf P})\>,$ (13)
$\displaystyle\Phi_{\bf b}({\bf P})$ $\displaystyle=-UM_{\bf b0}\,\Phi_{\bf
0}({\bf P})-\sum_{{\bf b}^{\prime}}V_{{\bf b}^{\prime}}\,M_{{\bf b}{\bf
b}^{\prime}}(E,{\bf P})\,\Phi_{{\bf b}^{\prime}}({\bf P})\>,$ (14)

where

$\displaystyle M_{\bf 00}(E,{\bf P})$ $\displaystyle=\frac{1}{N}\sum_{\bf
q}\frac{1}{-E+\varepsilon_{\bf q}+\varepsilon_{{\bf P}-{\bf q}}}\>,\hskip
28.45274ptM_{{\bf 0}{\bf b}^{\prime}}(E,{\bf P})=\frac{1}{N}\sum_{\bf
q}\frac{e^{-i{\bf q}{\bf b}^{\prime}}}{-E+\varepsilon_{\bf
q}+\varepsilon_{{\bf P}-{\bf q}}}\>,$ (15) $\displaystyle M_{\bf b0}(E,{\bf
P})$ $\displaystyle=\frac{1}{N}\sum_{\bf q}\frac{e^{i{\bf q}{\bf
b}}}{-E+\varepsilon_{\bf q}+\varepsilon_{{\bf P}-{\bf q}}}\>,\hskip
28.45274ptM_{{\bf b}{\bf b}^{\prime}}(E,{\bf P})=\frac{1}{N}\sum_{\bf
q}\frac{e^{i{\bf q}({\bf b}-{\bf b}^{\prime})}}{-E+\varepsilon_{\bf
q}+\varepsilon_{{\bf P}-{\bf q}}}\>.$ (16)

Notice how in the process of substitution, ${\bf k}_{1}$ is replaced by ${\bf
q}$, and ${\bf k}_{2}$ by ${\bf P}-{\bf q}$. But the argument of $\Phi$
remains unchanged: ${\bf P}={\bf k}_{1}+{\bf k}_{2}\rightarrow{\bf q}+{\bf
P}-{\bf q}={\bf P}$. As a result, $\Phi({\bf P})$ can be moved outside the
${\bf q}$ sums, and the equations reduce from integral to algebraic. Thus the
entire method is predicated on conservation of total momentum.

Quantities $M_{{\bf b}{\bf b}^{\prime}}$ in Eqs. (15) and (16) are two-body
Green’s functions of underlying lattices. In some lattices, $M_{{\bf b}{\bf
b}^{\prime}}$ can be reduced to one-body Green’s functions by an appropriate
transformation. In solving a typical two-body problem, most effort is spent on
calculating and analyzing $M_{{\bf b}{\bf b}^{\prime}}$, and the existence of
close-form final formulas depends on whether $M_{{\bf b}{\bf b}^{\prime}}$ can
be evaluated analytically. In 1D, integration is elementary leading to
algebraic expressions. In 2D, $M_{{\bf b}{\bf b}^{\prime}}$ can be usually
reduced to the complete elliptic integrals of different kinds. And in 3D,
$M_{{\bf b}{\bf b}^{\prime}}$ are generalized Watson integrals for which many
new analytical results have been obtained in the last 30 years [91, 92, 93,
94, 95, 96, 97, 98]. (As mentioned in the Introduction, the latest advances in
mathematical physics is one of the reasons that prompted this paper.) Much of
the Appendixes are devoted to deriving and listing available results for
$M_{{\bf b}{\bf b}^{\prime}}$ on different lattices. Note that we have flipped
the sign of the energy denominator in Eq. (15) for future convenience. In this
work, we will be interested only in bounds states with $E<0$. The definitions,
Eqs. (15) and (16), render most of $M_{{\bf b}{\bf b}^{\prime}}$ positive. In
the following, we will often write $|E|$ instead of $-E$ to avoid any
confusion.

The consistency condition of Eqs. (13) and (14),

${\rm det}\left|\\!\\!\begin{array}[]{cc}UM_{\bf 00}(E,{\bf P})+1&V_{{\bf
b}^{\prime}}\,M_{{\bf 0}{\bf b}^{\prime}}(E,{\bf P})\\\ UM_{\bf b0}(E,{\bf
P})&V_{{\bf b}^{\prime}}\,M_{{\bf b}{\bf b}^{\prime}}(E,{\bf P})+\delta_{{\bf
b}{\bf b}^{\prime}}\end{array}\\!\\!\right|=0\>,$ (17)

determines system’s energy $E$ as a function of total momentum ${\bf P}$ and,
consequently, the pair’s energy, dispersion, and effective mass. If $n_{\bf
b}$ is the number of vectors ${\bf b}$ with nonzero $V_{\bf b}$, then in Eq.
(17), $M_{\bf b0}$ is an $(n_{\bf b}\times 1)$ column, $V_{{\bf
b}^{\prime}}\,M_{{\bf 0}{\bf b}^{\prime}}$ is a $(1\times n_{\bf b})$ row, and
$V_{{\bf b}^{\prime}}\,M_{{\bf b}{\bf b}^{\prime}}$ is an $(n_{\bf b}\times
n_{\bf b})$ matrix where each column is multiplied by its respective $V_{{\bf
b}^{\prime}}$. The eigenvector of Eqs. (13) and (14) determines a pair wave
function via Eq. (12). Equations (12)-(17) constitute the general solution of
the two-body lattice problem.

###  2.3  Symmetrized solution. Spin singlets

The size and complexity of the main system, Eq. (17), rapidly grows with the
radius of interaction. Consider for example the square lattice. For contact,
zero-range interaction, Eq. (17) is a $(1\times 1)$ matrix, for nearest-
neighbor $UV$ model it is a $(5\times 5)$ matrix (one zero-range potential
plus four nearest-neighbor potentials), and for next-nearest-neighbor $UV$
model it is already a $(9\times 9)$ matrix. The ability to perform analytical
calculations diminishes rapidly, especially for ${\bf P}$’s away from special
symmetry points. In this situation, permutation symmetry offers a way to
simplify the solution. The two-fermion wave function must be either symmetric
or antisymmetric with respect to argument exchange: $\Psi({\bf m}_{2},{\bf
m}_{1})=\pm\Psi({\bf m}_{1},{\bf m}_{2})$, corresponding to spin-singlet ($+$)
and spin-triplet ($-$) pair states. Including this symmetry from the beginning
reduces the final system’s size by about half. An additional benefit is the
$(+)$ solutions also describe bound pairs of spin-$0$ bosons and the $(-)$
solutions describe spinless fermions. In the rest of this section we derive
the $(+)$ solution. The $(-)$ solution is derived in Section 2.4.

In order to restrict the wave function symmetry, we use the following
transformation instead of Eq. (5):

$\displaystyle\Psi({\bf m}_{1},{\bf m}_{2})$
$\displaystyle=\frac{1}{2N}\sum_{{\bf k}_{1}{\bf k}_{2}}\left(e^{i{\bf
k}_{1}{\bf m}_{1}+i{\bf k}_{2}{\bf m}_{2}}+e^{i{\bf k}_{1}{\bf m}_{2}+i{\bf
k}_{2}{\bf m}_{1}}\right)\phi^{+}_{{\bf k}_{1}{\bf k}_{2}}\>,\hskip
14.22636pt\Psi({\bf m}_{2},{\bf m}_{1})=+\Psi({\bf m}_{1},{\bf m}_{2})\>,$
(18) $\displaystyle\phi^{+}_{{\bf k}_{1}{\bf k}_{2}}$
$\displaystyle=\frac{1}{2N}\sum_{{\bf m}_{1}{\bf m}_{2}}\left(e^{-i{\bf
k}_{1}{\bf m}_{1}-i{\bf k}_{2}{\bf m}_{2}}+e^{-i{\bf k}_{1}{\bf m}_{2}-i{\bf
k}_{2}{\bf m}_{1}}\right)\Psi({\bf m}_{1},{\bf m}_{2})\>,\hskip
14.22636pt\phi^{+}_{{\bf k}_{2}{\bf k}_{1}}=+\phi^{+}_{{\bf k}_{1}{\bf
k}_{2}}\>.$ (19)

Next, we multiply the Schrödinger equation, Eq. (4), by the expression in
parentheses in Eq. (19) and apply the operation $(2N)^{-1}\sum_{{\bf
m}_{1}{\bf m}_{2}}$. The result is

$\displaystyle(E-\varepsilon_{{\bf k}_{1}}-\varepsilon_{{\bf
k}_{2}})\,\phi^{+}_{{\bf k}_{1}{\bf k}_{2}}$
$\displaystyle=U\frac{1}{N}\sum_{\bf q}\phi^{+}_{{\bf q},{\bf k}_{1}+{\bf
k}_{2}-{\bf q}}$ $\displaystyle+\sum_{\bf b}V_{\bf b}\frac{1}{4N}\sum_{\bf
q}\\!\left[\left(e^{i{\bf k}_{1}{\bf b}}+e^{i{\bf k}_{2}{\bf
b}}\right)e^{-i{\bf q}{\bf b}}\\!+\left(e^{-i{\bf k}_{1}{\bf b}}+e^{-i{\bf
k}_{2}{\bf b}}\right)e^{i{\bf q}{\bf b}}\right]\phi^{+}_{{\bf q},{\bf
k}_{1}+{\bf k}_{2}-{\bf q}}\>.$ (20)

The $V$ term contains two groups of ${\bf q}$ integrals: one that contains
$\exp{(i{\bf qb})}$ and another $\exp{(-i{\bf qb})}$. One now observes that
the two groups transform into each other when ${\bf b}\rightarrow-{\bf b}$. To
make use of this symmetry, we arrange all vectors ${\bf b}$ into pairs $({\bf
b},-{\bf b})$, then select only one vector from each pair and collect them in
new group ${\bf b}_{+}$. Thus, the full group splits in two subgroups:

$\\{{\bf b}\\}\rightarrow\\{{\bf b}_{+}\\},-\\{{\bf b}_{+}\\}\>.$ (21)

For example, in the square lattice with nearest-neighbor interaction, $\\{{\bf
b}\\}=\\{+{\bf x},+{\bf y},-{\bf x},-{\bf y}\\}$ can be split in two pairs,
$\\{(+{\bf x},-{\bf x}),(+{\bf y},-{\bf y})\\}$. Then one can choose $\\{{\bf
b}_{+}\\}$ out of four equivalent possibilities: $\\{{\bf b}_{+}\\}=\\{+{\bf
x},+{\bf y}\\}$, $\\{{\bf b}_{+}\\}=\\{+{\bf x},-{\bf y}\\}$, $\\{{\bf
b}_{+}\\}=\\{-{\bf x},+{\bf y}\\}$, and $\\{{\bf b}_{+}\\}=\\{-{\bf x},-{\bf
y}\\}$. The ${\bf b}$ sum in the $V$ term in Eq. (20) splits in two:

$\displaystyle V\>{\rm term}=$ $\displaystyle\;\sum_{{\bf b}_{+}}V_{{\bf
b}_{+}}\frac{1}{4N}\sum_{\bf q}\left[\left(e^{i{\bf k}_{1}{\bf
b}_{+}}+e^{i{\bf k}_{2}{\bf b}_{+}}\right)e^{-i{\bf q}{\bf
b}_{+}}+\left(e^{-i{\bf k}_{1}{\bf b}_{+}}+e^{-i{\bf k}_{2}{\bf
b}_{+}}\right)e^{i{\bf q}{\bf b}_{+}}\right]\phi^{+}_{{\bf q},{\bf k}_{1}+{\bf
k}_{2}-{\bf q}}$ $\displaystyle+\sum_{{\bf b}_{+}}V_{-{\bf
b}_{+}}\frac{1}{4N}\sum_{\bf q}\left[\left(e^{-i{\bf k}_{1}{\bf
b}_{+}}+e^{-i{\bf k}_{2}{\bf b}_{+}}\right)e^{i{\bf q}{\bf
b}_{+}}+\left(e^{i{\bf k}_{1}{\bf b}_{+}}+e^{i{\bf k}_{2}{\bf
b}_{+}}\right)e^{-i{\bf q}{\bf b}_{+}}\right]\phi^{+}_{{\bf q},{\bf
k}_{1}+{\bf k}_{2}-{\bf q}}$ $\displaystyle=$ $\displaystyle\;\sum_{{\bf
b}_{+}}V_{{\bf b}_{+}}\frac{1}{2N}\sum_{\bf q}\left[\left(e^{i{\bf k}_{1}{\bf
b}_{+}}+e^{i{\bf k}_{2}{\bf b}_{+}}\right)e^{-i{\bf q}{\bf
b}_{+}}+\left(e^{-i{\bf k}_{1}{\bf b}_{+}}+e^{-i{\bf k}_{2}{\bf
b}_{+}}\right)e^{i{\bf q}{\bf b}_{+}}\right]\phi^{+}_{{\bf q},{\bf k}_{1}+{\bf
k}_{2}-{\bf q}}\>.$ (22)

The last equality is true because we consider only symmetric potentials,
$V_{{\bf b}_{+}}=V_{-{\bf b}_{+}}$. In the next crucial step, we apply a
variable change ${\bf q}^{\prime}={\bf k}_{1}+{\bf k}_{2}-{\bf q}$ to the
first half of Eq. (22). That renders the exponential terms equal to the
exponential terms of the second half but the wave function transforms to
$\phi^{+}_{{\bf k}_{1}+{\bf k}_{2}-{\bf q},{\bf q}}$. However due to the
permutation symmetry it is equal to $\phi^{+}_{{\bf q},{\bf k}_{1}+{\bf
k}_{2}-{\bf q}}$. This proves that the two halves of Eq. (22) are equal. The
entire $V$ term can be written as twice the second half (for example).
Returning to the full equation, Eq. (20), it reads

$(E-\varepsilon_{{\bf k}_{1}}-\varepsilon_{{\bf k}_{2}})\,\phi^{+}_{{\bf
k}_{1}{\bf k}_{2}}=U\frac{1}{N}\sum_{\bf q}\phi^{+}_{{\bf q},{\bf k}_{1}+{\bf
k}_{2}-{\bf q}}+\sum_{{\bf b}_{+}}V_{{\bf b}_{+}}\left(e^{-i{\bf k}_{1}{\bf
b}_{+}}+e^{-i{\bf k}_{2}{\bf b}_{+}}\right)\frac{1}{N}\sum_{\bf q}e^{i{\bf
q}{\bf b}_{+}}\phi^{+}_{{\bf q},{\bf k}_{1}+{\bf k}_{2}-{\bf q}}\>.$ (23)

Equation (23) has only half as many $V$ terms as its unsymmetrized
counterpart, Eq. (7). Accordingly, we introduce auxiliary functions

$\displaystyle\Phi^{+}_{\bf 0}({\bf k}_{1}+{\bf k}_{2})$
$\displaystyle=\Phi^{+}_{\bf 0}({\bf P})\equiv\frac{1}{N}\sum_{\bf
q}\phi^{+}_{{\bf q},{\bf k}_{1}+{\bf k}_{2}-{\bf q}}=\frac{1}{N}\sum_{\bf
q}\phi^{+}_{{\bf q},{\bf P}-{\bf q}}\>,$ (24) $\displaystyle\Phi^{+}_{\bf
b_{+}}({\bf k}_{1}+{\bf k}_{2})$ $\displaystyle=\Phi^{+}_{\bf b_{+}}({\bf
P})\equiv\frac{1}{N}\sum_{\bf q}e^{i{\bf q}{\bf b}_{+}}\,\phi^{+}_{{\bf
q},{\bf k}_{1}+{\bf k}_{2}-{\bf q}}=\frac{1}{N}\sum_{\bf q}e^{i{\bf q}{\bf
b}_{+}}\,\phi^{+}_{{\bf q},{\bf P}-{\bf q}}\>.$ (25)

The wave function is expressed from Eq. (23)

$\phi^{+}_{{\bf k}_{1}{\bf k}_{2}}=\frac{U}{E-\varepsilon_{{\bf
k}_{1}}-\varepsilon_{{\bf k}_{2}}}\,\Phi^{+}_{\bf 0}({\bf P})+\sum_{{\bf
b}_{+}}V_{{\bf b}_{+}}\,\frac{e^{-i{\bf k}_{1}{\bf b}_{+}}+e^{-i{\bf
k}_{2}{\bf b}_{+}}}{E-\varepsilon_{{\bf k}_{1}}-\varepsilon_{{\bf
k}_{2}}}\,\Phi^{+}_{{\bf b}_{+}}({\bf P})\>.$ (26)

Substituting $\phi^{+}$ back in the definitions of $\Phi^{+}$, one obtains:

$\displaystyle\Phi^{+}_{\bf 0}({\bf P})$ $\displaystyle=-UM^{+}_{\bf
00}\,\Phi^{+}_{\bf 0}({\bf P})-\sum_{{\bf b}^{\prime}_{+}}V_{{\bf
b}^{\prime}_{+}}\,M^{+}_{{\bf 0}{\bf b}^{\prime}_{+}}(E,{\bf
P})\,\Phi^{+}_{{\bf b}^{\prime}_{+}}({\bf P})\>,$ (27)
$\displaystyle\Phi^{+}_{{\bf b}_{+}}({\bf P})$ $\displaystyle=-UM^{+}_{{\bf
b}_{+}{\bf 0}}\,\Phi^{+}_{\bf 0}({\bf P})-\sum_{{\bf b}^{\prime}_{+}}V_{{\bf
b}^{\prime}_{+}}\,M^{+}_{{\bf b}_{+}{\bf b}^{\prime}_{+}}(E,{\bf
P})\,\Phi^{+}_{{\bf b}^{\prime}_{+}}({\bf P})\>,$ (28)

where

$\displaystyle M^{+}_{\bf 00}(E,{\bf P})$ $\displaystyle=\frac{1}{N}\sum_{\bf
q}\frac{1}{-E+\varepsilon_{\bf q}+\varepsilon_{{\bf P}-{\bf q}}}\>,\hskip
14.22636ptM^{+}_{{\bf 0}{\bf b}^{\prime}_{+}}(E,{\bf P})=\frac{1}{N}\sum_{\bf
q}\frac{e^{-i{\bf q}{\bf b}^{\prime}_{+}}+e^{-i({\bf P}-{\bf q}){\bf
b}^{\prime}_{+}}}{-E+\varepsilon_{\bf q}+\varepsilon_{{\bf P}-{\bf q}}}\>,$
(29) $\displaystyle M^{+}_{{\bf b}_{+}{\bf 0}}(E,{\bf P})$
$\displaystyle=\frac{1}{N}\sum_{\bf q}\frac{e^{i{\bf q}{\bf
b}_{+}}}{-E+\varepsilon_{\bf q}+\varepsilon_{{\bf P}-{\bf q}}}\>,\hskip
14.22636ptM^{+}_{{\bf b}_{+}{\bf b}^{\prime}_{+}}(E,{\bf
P})=\frac{1}{N}\sum_{\bf q}\frac{e^{i{\bf q}({\bf b}_{+}-{\bf
b}^{\prime}_{+})}+e^{i{\bf q}{\bf b}_{+}}e^{-i({\bf P}-{\bf q}){\bf
b}^{\prime}_{+}}}{-E+\varepsilon_{\bf q}+\varepsilon_{{\bf P}-{\bf q}}}\>.$
(30)

For nonzero ${\bf P}$, it may be more convenient to use an equivalent but a
more symmetric formulation. Let us perform a variable change in Eqs. (29) and
(30):

${\bf q}\to{\bf q}+\frac{\bf P}{2}\>,\hskip 28.45274pt{\bf P}-{\bf
q}\to\frac{\bf P}{2}-{\bf q}\>.$ (31)

That leads to the appearance of factors like $\exp{(\pm i{\bf P}{\bf
b}_{+}/2)}$ in front of the ${\bf q}$ sums. They can be absorbed into a new
definition of $\Phi^{+}$:

${\tilde{\Phi}}^{+}_{0}({\bf P})\equiv\Phi^{+}_{0}({\bf P})\>,\hskip
28.45274pt{\tilde{\Phi}}^{+}_{{\bf b}_{+}}({\bf P})\equiv e^{-i\frac{\bf
P}{2}{\bf b}_{+}}\Phi^{+}_{{\bf b}_{+}}({\bf P})\>.$ (32)

In terms of the new amplitudes, Eqs. (27) and (28) assume the form

$\displaystyle{\tilde{\Phi}}^{+}_{\bf 0}({\bf P})$
$\displaystyle=-U{\tilde{M}}^{+}_{\bf 00}\,{\tilde{\Phi}}^{+}_{\bf 0}({\bf
P})-\sum_{{\bf b}^{\prime}_{+}}V_{{\bf b}^{\prime}_{+}}\,{\tilde{M}}^{+}_{{\bf
0}{\bf b}^{\prime}_{+}}(E,{\bf P})\,{\tilde{\Phi}}^{+}_{{\bf
b}^{\prime}_{+}}({\bf P})\>,$ (33) $\displaystyle{\tilde{\Phi}}^{+}_{{\bf
b}_{+}}({\bf P})$ $\displaystyle=-U{\tilde{M}}^{+}_{{\bf b}_{+}{\bf
0}}\,{\tilde{\Phi}}^{+}_{\bf 0}({\bf P})-\sum_{{\bf b}^{\prime}_{+}}V_{{\bf
b}^{\prime}_{+}}\,{\tilde{M}}^{+}_{{\bf b}_{+}{\bf b}^{\prime}_{+}}(E,{\bf
P})\,{\tilde{\Phi}}^{+}_{{\bf b}^{\prime}_{+}}({\bf P})\>,$ (34)

where

$\displaystyle{\tilde{M}}^{+}_{\bf 00}(E,{\bf P})$
$\displaystyle=\frac{1}{N}\sum_{\bf q}\frac{1}{-E+\varepsilon_{{\bf
q}+\frac{\bf P}{2}}+\varepsilon_{{\bf q}-\frac{\bf P}{2}}}\>,\hskip
14.22636pt{\tilde{M}}^{+}_{{\bf 0}{\bf b}^{\prime}_{+}}(E,{\bf
P})=\frac{1}{N}\sum_{\bf q}\frac{2\cos{({\bf q}{\bf
b}_{+})}}{-E+\varepsilon_{{\bf q}+\frac{\bf P}{2}}+\varepsilon_{{\bf
q}-\frac{\bf P}{2}}}\>,$ (35) $\displaystyle{\tilde{M}}^{+}_{{\bf b}_{+}{\bf
0}}(E,{\bf P})$ $\displaystyle=\frac{1}{N}\sum_{\bf q}\frac{\cos{({\bf q}{\bf
b}_{+})}}{-E+\varepsilon_{{\bf q}+\frac{\bf P}{2}}+\varepsilon_{{\bf
q}-\frac{\bf P}{2}}}\>,\hskip 14.22636pt{\tilde{M}}^{+}_{{\bf b}_{+}{\bf
b}^{\prime}_{+}}(E,{\bf P})=\frac{1}{N}\sum_{\bf q}\frac{2\cos{({\bf q}{\bf
b}_{+})}\cos{({\bf q}{\bf b}^{\prime}_{+})}}{-E+\varepsilon_{{\bf q}+\frac{\bf
P}{2}}+\varepsilon_{{\bf q}-\frac{\bf P}{2}}}\>,$ (36)

and the condition $\varepsilon_{-{\bf q}}=\varepsilon_{\bf q}$ has been used.
There are several advantages of Eqs. (33)-(36) over Eqs. (27)-(30): (i) All
${\tilde{M}}$’s are manifestly real which simplifies analytics; (ii) Quite
generally, ${\tilde{M}}^{+}_{{\bf 0}{\bf
b}^{\prime}_{+}}=2{\tilde{M}}^{+}_{{\bf b}_{+}{\bf 0}}$; and (iii) Pair
momentum ${\bf P}$ only enters the denominators of ${\tilde{M}}$’s, which
simplifies analysis of pair effective mass and dispersion in some cases.
Otherwise, the two formulations are equivalent.

The consistency condition of the linear system, Eqs. (27)-(28) or Eqs.
(33)-(34), defines the pair’s energy, and its eigenvector defines the pair’s
wave function. The system comprises only $1+\frac{1}{2}n_{\bf b}$ linear
equations versus $1+n_{\bf b}$ equations in the unsymmetrized method. This is
significant simplification. For example, in the simple-cubic $UV$ model with
nearest neighbor interaction, symmetrization reduces the consistency condition
from a $(7\times 7)$ matrix to a $(4\times 4)$ matrix. The symmetrized
solution describes only spin-singlet pairs. Spin triplets are discussed next.

###  2.4  Anti-symmetrized solution. Spin triplets

In this section, we repeat the derivation of Section 2.3 for anti-symmetric
wave functions that describe spin-triplet pairs. The corresponding Fourier
transformation reads

$\displaystyle\Psi({\bf m}_{1},{\bf m}_{2})$
$\displaystyle=\frac{1}{2N}\sum_{{\bf k}_{1}{\bf k}_{2}}\left(e^{i{\bf
k}_{1}{\bf m}_{1}+i{\bf k}_{2}{\bf m}_{2}}-e^{i{\bf k}_{1}{\bf m}_{2}+i{\bf
k}_{2}{\bf m}_{1}}\right)\phi^{-}_{{\bf k}_{1}{\bf k}_{2}}\>,\hskip
14.22636pt\Psi({\bf m}_{2},{\bf m}_{1})=-\Psi({\bf m}_{1},{\bf m}_{2})\>,$
(37) $\displaystyle\phi^{-}_{{\bf k}_{1}{\bf k}_{2}}$
$\displaystyle=\frac{1}{2N}\sum_{{\bf m}_{1}{\bf m}_{2}}\left(e^{-i{\bf
k}_{1}{\bf m}_{1}-i{\bf k}_{2}{\bf m}_{2}}-e^{-i{\bf k}_{1}{\bf m}_{2}-i{\bf
k}_{2}{\bf m}_{1}}\right)\Psi({\bf m}_{1},{\bf m}_{2})\>,\hskip
14.22636pt\phi^{-}_{{\bf k}_{2}{\bf k}_{1}}=-\phi^{-}_{{\bf k}_{1}{\bf
k}_{2}}\>.$ (38)

We multiply the Schrödinger equation, Eq. (4), by the expression in
parentheses in Eq. (38) and apply operation $(2N)^{-1}\sum_{{\bf m}_{1}{\bf
m}_{2}}$. The result is

$(E-\varepsilon_{{\bf k}_{1}}-\varepsilon_{{\bf k}_{2}})\,\phi^{-}_{{\bf
k}_{1}{\bf k}_{2}}=\sum_{\bf b}V_{\bf b}\frac{1}{4N}\sum_{\bf
q}\\!\left[\left(e^{i{\bf k}_{1}{\bf b}}-e^{i{\bf k}_{2}{\bf
b}}\right)e^{-i{\bf q}{\bf b}}\\!+\left(e^{-i{\bf k}_{1}{\bf b}}-e^{-i{\bf
k}_{2}{\bf b}}\right)e^{i{\bf q}{\bf b}}\right]\phi^{-}_{{\bf q},{\bf
k}_{1}+{\bf k}_{2}-{\bf q}}\>.$ (39)

Of note here is the absence of a $U$ term that cancels out due to
antisymmetry. Next, we split the sum over ${\bf b}$ in two partial sums: over
${\bf b}_{+}$ and over $-{\bf b}_{+}$, and then change variables ${\bf
q}^{\prime}={\bf k}_{1}+{\bf k}_{2}-{\bf q}$ in the $\exp{(-i{\bf q}{\bf
b}_{+})}$ terms. The result is

$(E-\varepsilon_{{\bf k}_{1}}-\varepsilon_{{\bf k}_{2}})\,\phi^{-}_{{\bf
k}_{1}{\bf k}_{2}}=\sum_{{\bf b}_{+}}V_{{\bf b}_{+}}\left(e^{-i{\bf k}_{1}{\bf
b}_{+}}-e^{-i{\bf k}_{2}{\bf b}_{+}}\right)\frac{1}{N}\sum_{\bf q}e^{i{\bf
q}{\bf b}_{+}}\phi^{-}_{{\bf q},{\bf k}_{1}+{\bf k}_{2}-{\bf q}}\>.$ (40)

To convert this to linear equations, we introduce $\frac{1}{2}n_{\bf b}$
auxiliary functions

$\Phi^{-}_{\bf b_{+}}({\bf k}_{1}+{\bf k}_{2})=\Phi^{-}_{\bf b_{+}}({\bf
P})\equiv\frac{1}{N}\sum_{\bf q}e^{i{\bf q}{\bf b}_{+}}\,\phi^{-}_{{\bf
q},{\bf k}_{1}+{\bf k}_{2}-{\bf q}}=\frac{1}{N}\sum_{\bf q}e^{i{\bf q}{\bf
b}_{+}}\,\phi^{-}_{{\bf q},{\bf P}-{\bf q}}\>.$ (41)

Using these definitions, the pair wave function follows from Eq. (40):

$\phi^{-}_{{\bf k}_{1}{\bf k}_{2}}=\sum_{{\bf b}_{+}}V_{{\bf
b}_{+}}\,\frac{e^{-i{\bf k}_{1}{\bf b}_{+}}-e^{-i{\bf k}_{2}{\bf
b}_{+}}}{E-\varepsilon_{{\bf k}_{1}}-\varepsilon_{{\bf
k}_{2}}}\,\Phi^{-}_{{\bf b}_{+}}({\bf P})\>.$ (42)

Substituting $\phi^{-}_{{\bf k}_{1}{\bf k}_{2}}$ back in Eq. (41), one obtains
the final system

$\Phi^{-}_{{\bf b}_{+}}({\bf P})=-\sum_{{\bf b}^{\prime}_{+}}V_{{\bf
b}^{\prime}_{+}}\,M^{-}_{{\bf b}_{+}{\bf b}^{\prime}_{+}}(E,{\bf
P})\,\Phi^{-}_{{\bf b}^{\prime}_{+}}({\bf P})\>,$ (43) $M^{-}_{{\bf b}_{+}{\bf
b}^{\prime}_{+}}(E,{\bf P})=\frac{1}{N}\sum_{\bf q}\frac{e^{i{\bf q}({\bf
b}_{+}-{\bf b}^{\prime}_{+})}-e^{i{\bf q}{\bf b}_{+}}e^{-i({\bf P}-{\bf
q}){\bf b}^{\prime}_{+}}}{-E+\varepsilon_{\bf q}+\varepsilon_{{\bf P}-{\bf
q}}}\>.$ (44)

Applying the variable change defined in Eq. (31) and introducing new
amplitudes:

${\tilde{\Phi}}^{-}_{{\bf b}_{+}}({\bf P})\equiv e^{-i\frac{\bf P}{2}{\bf
b}_{+}}\Phi^{-}_{{\bf b}_{+}}({\bf P})\>,$ (45)

the system is cast in an alternative form:

${\tilde{\Phi}}^{-}_{{\bf b}_{+}}({\bf P})=-\sum_{{\bf b}^{\prime}_{+}}V_{{\bf
b}^{\prime}_{+}}\,{\tilde{M}}^{-}_{{\bf b}_{+}{\bf b}^{\prime}_{+}}(E,{\bf
P})\,{\tilde{\Phi}}^{-}_{{\bf b}^{\prime}_{+}}({\bf P})\>,$ (46)
${\tilde{M}}^{-}_{{\bf b}_{+}{\bf b}^{\prime}_{+}}(E,{\bf
P})=\frac{1}{N}\sum_{\bf q}\frac{2\sin{({\bf q}{\bf b}_{+})}\sin{({\bf q}{\bf
b}^{\prime}_{+})}}{-E+\varepsilon_{{\bf q}+\frac{\bf P}{2}}+\varepsilon_{{\bf
q}-\frac{\bf P}{2}}}\>.$ (47)

The size of the triplet system is one less that of the singlet system because
of the lack of the Hubbard term. Therefore, triplet pairs are usually easier
to deal with than singlets.

###  2.5  Pair size

In this section, we derive a general recipe of calculating pair’s effective
radius $r^{\ast}_{p}$. We consider only spin-singlet pairs with ${\bf P}=0$.
Other situations can be analyzed similarly. Radius components are defined as
follows

$\left(r^{\ast}_{pj}\right)^{2}=\langle m^{2}_{j}\rangle=\frac{\sum_{\bf
m}m^{2}_{j}\>\Psi^{\ast}({\bf m},0)\Psi({\bf m},0)}{\sum_{\bf
m}\Psi^{\ast}({\bf m},0)\Psi({\bf m},0)}\equiv\frac{J_{j}}{J_{0}}\>.$ (48)

The ${\bf k}_{1}+{\bf k}_{2}=0$ wave function follows from Eq. (18):

$\Psi({\bf m},0)=\frac{1}{2N}\sum_{\bf k}\left(e^{i{\bf k}{\bf m}}+e^{-i{\bf
k}{\bf m}}\right)\phi^{+}_{{\bf k},-{\bf k}}\>.$ (49)

Substitution in the denominator of Eq. (48) yields

$J_{0}=\sum_{\bf m}\Psi^{\ast}({\bf m},0)\Psi({\bf m},0)=\frac{1}{N}\sum_{\bf
k}\phi^{+}_{{\bf k},-{\bf k}}\phi^{+\ast}_{{\bf k},-{\bf
k}}=\frac{1}{N}\sum_{\bf k}\left|\phi^{+}_{{\bf k},-{\bf k}}\right|^{2}.$ (50)

Upon substitution in $J_{j}$, one notices that $m_{j}$ can be expresses as a
$k_{j}$ derivative of $e^{i{\bf k}{\bf m}}$. Taking into account the
permutation symmetry of $\phi^{+}_{{\bf k},-{\bf k}}$, one obtains

$J_{j}=\frac{1}{N^{2}}\sum_{\bf m}\left(\sum_{\bf k}\phi^{+}_{{\bf k},-{\bf
k}}\frac{\partial}{\partial k_{j}}e^{i{\bf k}{\bf m}}\right)\left(\sum_{{\bf
k}^{\prime}}\phi^{+\ast}_{{\bf k}^{\prime},-{\bf
k}^{\prime}}\frac{\partial}{\partial k^{\prime}_{j}}e^{-i{\bf k}^{\prime}{\bf
m}}\right).$ (51)

Both expressions in parentheses can be transformed utilizing the Green’s
theorem for periodic functions, see for example Ref. [105]. That transfers the
derivatives to $\phi^{+}$. The final result for the pair radius reads

$\left(r^{\ast}_{pj}\right)^{2}=\langle
m^{2}_{j}\rangle=\frac{J_{j}}{J_{0}}=\frac{\frac{1}{N}\sum_{\bf
k}\frac{\partial}{\partial k_{j}}(\phi^{+}_{{\bf k},-{\bf
k}})\frac{\partial}{\partial k_{j}}(\phi^{+\ast}_{{\bf k},-{\bf
k}})}{\frac{1}{N}\sum_{\bf k}(\phi^{+}_{{\bf k},-{\bf k}})(\phi^{+\ast}_{{\bf
k},-{\bf k}})}\>.$ (52)

The wave function is expressible via amplitudes $\Phi^{+}$ using Eq. (26)

$\phi^{+}_{{\bf k},-{\bf k}}=\frac{U}{E-2\varepsilon_{\bf k}}\,\Phi^{+}_{\bf
0}(0)+\sum_{{\bf b}_{+}}V_{{\bf b}_{+}}\,\frac{2\cos{({\bf k}{\bf
b}_{+})}}{E-2\varepsilon_{\bf k}}\,\Phi^{+}_{{\bf b}_{+}}(0)\>.$ (53)

From here, one derives

$\frac{\partial}{\partial k_{j}}(\phi^{+}_{{\bf k},-{\bf
k}})=\left[\frac{2U}{(E-2\varepsilon_{\bf
k})^{2}}\frac{\partial\varepsilon_{\bf k}}{\partial k_{j}}\right]\Phi^{+}_{\bf
0}(0)+\sum_{{\bf b}_{+}}V_{{\bf b}_{+}}\left[\frac{-2(b_{+j})\sin{({\bf k}{\bf
b}_{+})}}{E-2\varepsilon_{\bf k}}+\frac{4\cos{({\bf k}{\bf
b}_{+})}}{(E-2\varepsilon_{\bf k})^{2}}\frac{\partial\varepsilon_{\bf
k}}{\partial k_{j}}\right]\Phi^{+}_{{\bf b}_{+}}(0)\>.$ (54)

Thus, calculation of pair size comprises the following steps. (i) For given
$U$, $V_{{\bf b}_{+}}$, and $\varepsilon_{\bf k}$, one solves the eigenvalue
problem, Eqs. (27)-(28) or Eqs. (33)-(34), for ${\bf P}=0$. That provides pair
energy $E$ and eigenvector $\\{\Phi^{+}_{\bf 0},\Phi^{+}_{{\bf b}_{+}}\\}$.
The eigenvector does not have to be normalized since the normalization
constant cancels out in the ratio of two integrals. (ii) The wave function and
its derivatives are calculated according to Eqs. (53) and (54). (iii)
Everything is substituted in Eq. (52) and the two integrals are computed
numerically. This method is more efficient than direct calculation of
$\Psi({\bf m},0)$ as a Fourier transform of $\phi^{+}$ followed by a summation
over ${\bf m}$. In general, one expects $r^{\ast}_{p}$ to diverge near the
binding threshold and to be of order one lattice constant in the strong
coupling limit. In the simplest cases, $r^{\ast}_{p}$ can be calculated
analytically as shown in the next section.

##  3  Negative-$U$ Hubbard model

###  3.1  General expressions

Before getting to more complex $UV$ models, it is instructive to consider the
simpler case of zero-range interaction, that is the negative Hubbard model.
Several characteristic features of lattice bound states show up already at
this level. Additionally, due to the model’s relative simplicity, analytical
calculations can be carried out to the fullest extent. The model is defined by
the potential

$U=-|U|\>,\hskip 28.45274ptV_{\bf b}=0\>.$ (55)

One expects only one singlet bound state, so either unsymmetrized solution,
Eq. (13), or symmetrized one, Eq. (27), can be applied. In both cases, the
system reduces to a single equation for $\Phi_{\bf 0}$ with the consistency
condition

$|U|\,M_{\bf 00}(E,{\bf P})=1\>,$ (56)

which defines pair energy $E({\bf P})$. The pair wave function is

$\psi_{{\bf k}_{1}{\bf k}_{2}}=\frac{1}{E-\varepsilon_{{\bf
k}_{1}}-\varepsilon_{{\bf k}_{2}}}\>,$ (57)

up to a normalization constant. Since total momentum ${\bf P}$ is fixed,
$\psi$ is a function of only one argument:

$\psi_{{\bf P}}({\bf q})=\frac{1}{E-\varepsilon_{\bf q}-\varepsilon_{{\bf
P}-{\bf q}}}\>.$ (58)

The real-space wave function follows from Eq. (5)

$\Psi({\bf m}_{1},{\bf m}_{2})=\frac{1}{N}\sum_{\bf q}\frac{e^{i{\bf q}{\bf
m}_{1}+i({\bf P}-{\bf q}){\bf m}_{2}}}{E-\varepsilon_{\bf q}-\varepsilon_{{\bf
P}-{\bf q}}}=e^{i{\bf P}\frac{({\bf m}_{1}+{\bf
m}_{2})}{2}}\frac{1}{N}\sum_{\bf q}\frac{e^{i{\bf q}({\bf m}_{1}-{\bf
m}_{2})}}{E-\varepsilon_{\frac{{\bf P}}{2}+{\bf q}}-\varepsilon_{\frac{{\bf
P}}{2}-{\bf q}}}\>.$ (59)

The first factor describes center-of-mass motion while the integral over ${\bf
q}$ describes internal structure of the pair. Analysis will continue for
different lattices separately.

###  3.2  1D. One dimensional chain

The 1D attractive Hubbard model provides the simplest example of a lattice
bound state. Many pair properties can be derived analytically. The basic
integral, $M_{\bf 00}$ in Eq. (15), is

$M^{\rm 1D}_{\bf 00}=\int^{\pi}_{-\pi}\frac{{\rm
d}q}{2\pi}\frac{1}{|E|-4t\cos{\left(\frac{P}{2}\right)}\cos{q}}=\frac{1}{\sqrt{E^{2}-16\,t^{2}\cos^{2}{\left(\frac{P}{2}\right)}}}\>.$
(60)

Substitution in Eq. (56) yields pair energy

$E(P)=-\sqrt{|U|^{2}+16\,t^{2}\cos^{2}{\left(\frac{P}{2}\right)}}\>.$ (61)

This is a rare case when pair energy is known as an explicit formula. Based on
it, a number of interesting properties can be established. (i) The minimum
energy of two free particles with total momentum $P$ is
$E_{11}=-4t\cos{(P/2)}$. Comparing that with Eq. (61), one finds
$E(P)<E_{11}(P)$ for any $|U|>0$. In other words, $|U|=0$ is the threshold of
pair formation for any $P$. The same conclusion can be reached by noting that
$M_{\bf 00}$ diverges at $E\rightarrow E_{11}$, see Eq. (60). (ii) Energy
$E(P)$ is periodic with period $(2\pi)$, despite $P$ being a sum of two
single-particle momenta, each of which varying between $-\pi$ and $\pi$. Thus,
pairing leads to Brillouin zone (BZ) folding and the pair behaves as one
particle with $-\pi\leq P\leq\pi$. (iii) Pair energy in BZ corners is
$E(\pm\pi)=-|U|$. (iv) The pair binding energy is quadratic in coupling near
the threshold:

$E(|U|\ll
t)=-4t\cos{\left(\frac{P}{2}\right)}-\frac{|U|^{2}}{8t\cos{\left(\frac{P}{2}\right)}}\>.$
(62)

The first term here is the minimum energy of two free particles with total
momentum $P$. (v) Expansion of Eq. (61) for small $P$ yields the pair
effective mass [in units of the bare one-particle mass
$m_{0}=\hbar^{2}/(2ta^{2})$]:

$\frac{m^{\ast}_{p}}{m_{0}}=\frac{\sqrt{|U|^{2}+16\,t^{2}}}{2t}\>.$ (63)

The pair mass is not constant but increases with the binding energy. This is a
common property of bound states [1] related to the lack of Galilean invariance
on the lattice. Comparison between Eqs. (63) and (61) reveals a curious
relationship between the pair mass and its “rest energy” $E(0)$. Restoring for
a moment the intersite distance $a$, and using $m_{0}=\hbar^{2}/(2ta^{2})$,
one obtains

$|E(0)|=m^{\ast}_{p}\left(\frac{2ta}{\hbar}\right)^{2}.$ (64)

The expression in parentheses is recognized as the maximum group velocity on
the lattice. Thus, Eq. (64) has the form of $E=mc^{2}$ of relativistic
physics.

Transitioning to the wave function, the integral in Eq. (59) can be calculated
explicitly [106]:

$\int^{\pi}_{-\pi}\frac{{\rm
d}q}{2\pi}\,\frac{\cos{[q(m_{1}-m_{2})]}}{|E|-4t\cos{\frac{P}{2}}\,\cos{q}}=\frac{1}{\sqrt{|E|^{2}-\alpha^{2}}}\left[\frac{\alpha}{|E|+\sqrt{|E|^{2}-\alpha^{2}}}\right]^{|m_{1}-m_{2}|}=\frac{1}{|U|}\left[\frac{|U|+\sqrt{|U|^{2}+\alpha^{2}}}{\alpha}\right]^{-|m_{1}-m_{2}|}.$
(65)

where $\alpha\equiv 4t\cos(P/2)$ has been set for brevity. Thus, the un-
normalized wave function can be written as

$\Psi(m_{1},m_{2})=e^{i\frac{P}{2}(m_{1}+m_{2})}\cdot
e^{-\gamma|m_{1}-m_{2}|}\>,$ (66)

where

$\sinh{\gamma}=\frac{|U|}{\alpha}=\frac{|U|}{4t\cos{\frac{P}{2}}}\>,$ (67)

and

$E=-4t\cos{\frac{P}{2}}\cosh{\gamma}\>.$ (68)

The same expressions can be derived directly from a real-space Schrödinger
equation by means of two-particle Bethe ansatz. Using the explicit form of
$\Psi$, it is straightforward to compute potential energy, kinetic energy, and
effective radius of a moving pair:

$E_{\rm
pot}=\left\langle-|U|\,\delta_{m_{1},m_{2}}\right\rangle=-\frac{|U|^{2}}{\sqrt{|U|^{2}+16\,t^{2}\cos^{2}{\\!\frac{P}{2}}}}\>,$
(69) $E_{\rm kin}=E-E_{\rm
pot}=-\frac{16\,t^{2}\cos^{2}{\\!\frac{P}{2}}}{\sqrt{|U|^{2}+16\,t^{2}\cos^{2}{\\!\frac{P}{2}}}}\>,$
(70)
$r^{\ast}_{p}=\left\langle\left(m_{1}-m_{2}\right)^{2}\right\rangle^{1/2}=\frac{4t\cos{\\!\frac{P}{2}}}{\sqrt{2}\,|U|}\>.$
(71)

Interestingly, $r^{\ast}_{p}(P\rightarrow\pi)\rightarrow 0$, which means the
pair shrinks to a point. The same conclusion can also be derived directly from
Eq. (59). At $P=\pi$, the two kinetic terms in the denominator cancel out
which renders the internal wave function to be $\propto\delta_{{\bf m}_{1}{\bf
m}_{2}}$. The pair energy, mass, and radius of the 1D Hubbard pair are plotted
in Fig. 3.

Figure 3: Properties of ${\bf P}=0$ bound pairs in the attractive Hubbard
model on the 1D chain, 2D square lattice, and 3D simple cubic lattice. The
black circle marks the pair formation threshold in 3D, Eq. (87).

###  3.3  2D. Square lattice

For the square lattice, the basic integral $M_{\bf 00}$ in Eq. (15) can be
expressed via complete elliptic integral of the first kind ${\bf K}(z)$, see
B.5, Eq. (316), for details. Denoting $\alpha\equiv 4t\cos{\frac{P_{x}}{2}}$
and $\beta\equiv 4t\cos{\frac{P_{y}}{2}}$, one has

$M^{\rm sq}_{\bf
00}=\int\limits^{\pi}_{-\pi}\\!\\!\int\limits^{\pi}_{-\pi}\frac{{\rm
d}q_{x}\,{\rm
d}q_{y}}{(2\pi)^{2}}\frac{1}{|E|-\alpha\cos{q_{x}}-\beta\cos{q_{y}}}=\frac{2}{\pi\sqrt{|E|^{2}-(\alpha-\beta)^{2}}}\,{\bf
K}\\!\\!\left[\sqrt{\frac{4\alpha\beta}{|E|^{2}-(\alpha-\beta)^{2}}}\right]\\!.$
(72)

This result applies not only to the isotropic square model at arbitrary ${\bf
P}$, but also to the rectangular model with $t_{x}\neq t_{y}$. Inserting Eq.
(72) in Eq. (56) defines pair energy in the most general case. The minimum
energy of two free particles is $E_{11}=2\varepsilon_{{\bf
P}/2}=-(\alpha+\beta)$, at which the argument of ${\bf K}$ reaches 1 and
$M^{\rm sq}_{\bf 00}$ diverges logarithmically. Similar to 1D, the divergence
is interpreted as the existence of a bound state at any nonzero $|U|$. Thus,
$|U|=0$ is pair formation threshold at any ${\bf P}$. Let us determine pair
energy near threshold. Setting $E=-\alpha-\beta-\Delta$,
$\Delta\ll\alpha,\beta$, Eqs. (72) and (56) produce in the leading order

$\frac{|U|}{\pi\sqrt{\alpha\beta}}\,{\bf
K}\\!\left(1-\frac{\alpha+\beta}{4\alpha\beta}\,\Delta\right)=1\>.$ (73)

Using the asymptote ${\bf K}(1-z)\sim\frac{1}{2}\ln{\frac{8}{z}}$, at
$z\rightarrow+0$, one obtains

$\Delta=\frac{32\,\alpha\beta}{\alpha+\beta}\>{\rm
exp}\left(-\frac{2\pi\sqrt{\alpha\beta}}{|U|}\right).$ (74)

In the ground state, $\alpha=\beta=4t$, and the binding energy is

$\Delta_{0}=64\,t\>{\rm exp}\left(-\frac{8\pi t}{|U|}\right).$ (75)

Near BZ corners, pair energy is $E(\pm\pi,\pm\pi)=-|U|$, like in the 1D case.

Pair effective mass is derived next. It is most convenient to consider BZ
diagonal, $P_{x}=P_{y}=P$, where Eq. (72) simplifies considerably. Writing
$|E|=|E_{0}|-\frac{\hbar^{2}P^{2}}{m^{\ast}_{p}}$, expanding Eq. (72) for
$P\ll 1$, and applying the formula

$\frac{d{\bf K}(z)}{dz}=\frac{{\bf E}(z)}{z(1-z^{2})}-\frac{{\bf K}(z)}{z}\>,$
(76)

one obtains from Eq. (56):

$\frac{m_{0}}{m^{\ast}_{p}}=\frac{|E_{0}|}{16t}\left\\{1-\left[1-\frac{(8t)^{2}}{|E_{0}|^{2}}\right]\frac{{\bf
K}\left(\frac{8t}{|E_{0}|}\right)}{{\bf
E}\left(\frac{8t}{|E_{0}|}\right)}\right\\},$ (77)

where ${\bf E}(\kappa)$ is complete elliptic integral of the second kind.
Equation (77) has correct limits: $m^{\ast}_{p}(|U|\rightarrow 0)=2m_{0}$ and
$m^{\ast}_{p}(|U|\rightarrow\infty)=\frac{|U|}{2t}\,m_{0}$. The factor ${\bf
K}\left(\frac{8t}{|E_{0}|}\right)$ in Eq. (77) can be written as
$\frac{\pi|E_{0}|}{2|U|}$, which follows from Eq. (56).

Pair effective radius is discussed next. Unlike 1D, there is no explicit
formula for the pair wave function for all ${\bf m}$ in 2D. However, for each
given ${\bf m}$, the wave function can be derived from a few basic integrals
using recurrence relations, as explained in B. On BZ diagonals, $\Psi({\bf
m}_{1},{\bf m}_{2})$ is always a linear combination of ${\bf K}$ and ${\bf
E}$. At arbitrary ${\bf P}$, the wave function is a linear combination of all
three complete elliptic integrals ${\bf K}$, ${\bf E}$, and ${\bf\Pi}$. To
calculate effective radius, we substitute Eq. (59) in Eq. (48) and perform the
following transformation

$r^{\ast 2}_{pj}=\frac{\sum_{\bf m}\sum_{{\bf
q}_{1}}\frac{\frac{\partial}{\partial q_{1j}}(e^{-i{\bf q}_{1}{\bf
m}})}{E-\xi_{{\bf P},{\bf q}_{1}}}\sum_{{\bf
q}_{2}}\frac{\frac{\partial}{\partial q_{2\alpha}}(e^{i{\bf q}_{2}{\bf
m}})}{E-\xi_{{\bf P},{\bf q}_{2}}}}{\sum_{\bf m}\sum_{{\bf
q}_{1}}\frac{e^{-i{\bf q}_{1}{\bf m}}}{E-\xi_{{\bf P},{\bf q}_{1}}}\sum_{{\bf
q}_{2}}\frac{e^{i{\bf q}_{2}{\bf m}}}{E-\xi_{{\bf P},{\bf
q}_{2}}}}=\frac{\sum_{\bf q}\frac{\left(\frac{\partial\xi_{{\bf P},{\bf
q}}}{\partial q_{j}}\right)^{2}}{(E-\xi_{{\bf P},{\bf q}})^{4}}}{\sum_{\bf
q}\frac{1}{(E-\xi_{{\bf P},{\bf
q}})^{2}}}=\frac{\left(-\frac{1}{6}\right)\frac{\partial^{3}}{\partial|E|^{3}}\sum_{\bf
q}\frac{\left(\frac{\partial\xi_{{\bf P},{\bf q}}}{\partial
q_{j}}\right)^{2}}{|E|+\xi_{{\bf P},{\bf
q}}}}{(-1)\frac{\partial}{\partial|E|}\sum_{\bf q}\frac{1}{|E|+\xi_{{\bf
P},{\bf q}}}}\,,$ (78)

where $\xi_{{\bf P},{\bf q}}\equiv\varepsilon_{\frac{{\bf P}}{2}+{\bf
q}}+\varepsilon_{\frac{{\bf P}}{2}-{\bf q}}$. Both sums in Eq. (78) are
recognized as integrals $M^{\rm sq}_{nm}$ of the square lattice evaluated in
B. Let us limit consideration to BZ diagonals, $P_{x}=P_{x}\equiv P$. With the
notation of B.3, one writes

$r^{\ast 2}_{px}=r^{\ast
2}_{py}=\frac{(4t\cos{\frac{P}{2}})^{2}}{12}\frac{\frac{\partial^{3}}{\partial|E|^{3}}\left(M^{\rm
sq}_{00}-M^{\rm sq}_{20}\right)}{\frac{\partial}{\partial|E|}\left(M^{\rm
sq}_{00}\right)}\>.$ (79)

Using explicit expressions, Eqs. (295) and (298), one derives, after
transformations, a final formula

$r^{\ast 2}_{px}=r^{\ast
2}_{py}=\frac{1}{12}\left\\{\frac{1+\left(\frac{8t\cos{\frac{P}{2}}}{|E|}\right)^{2}}{1-\left(\frac{8t\cos{\frac{P}{2}}}{|E|}\right)^{2}}-\frac{{\bf
K}\left(\frac{8t\cos{\frac{P}{2}}}{|E|}\right)}{{\bf
E}\left(\frac{8t\cos{\frac{P}{2}}}{|E|}\right)}\right\\}.$ (80)

The pair shrinks to a point in the strong coupling limit,
$E\rightarrow-\infty$, and at $P=\pm\pi$ for any $|U|$. The radius diverges at
the threshold, $E\rightarrow-8t\cos{\frac{P}{2}}$, as expected on physical
reasoning. In the ground state, ${\bf P}=0$, the asymptotes are
$r^{\ast}_{px}(|U|\rightarrow\infty)\approx\frac{4t}{\sqrt{2}|U|}$ and
$r^{\ast}_{px}(|U|\rightarrow 0)\approx 96^{-\frac{1}{2}}\,\exp{\frac{4\pi
t}{|U|}}$. The energy, mass, and radius of the 2D Hubbard pair are plotted in
Fig. 3.

###  3.4  2D. Triangular lattice

In this section, we consider the two-dimensional triangular lattice with
nearest-neighbor isotropic hopping, $t_{\bf b}=t$. Single particle dispersion,
Eq. (8), is

$\varepsilon_{\bf
k}=-2t\cos{k_{x}}-4t\cos{\left(k_{x}/2\right)\cos{(\sqrt{3}k_{y}/2)}}.$ (81)

The double integral $M^{\rm tr}_{\bf 00}$, Eq. (15), corresponding to this
$\varepsilon_{\bf k}$ is evaluated in C.1. For the ground state,
$P_{x}=P_{y}=0$, the result is

$M^{\rm tr}_{\bf
00}=\frac{2}{\pi}\frac{1}{\sqrt{|E_{0}|^{2}-48t^{2}+16t\sqrt{2t|E_{0}|+12t^{2}}}}\>{\bf
K}\\!\left(\sqrt{\frac{32t\sqrt{2|E_{0}|t+12t^{2}}}{|E_{0}|^{2}-48t^{2}+16t\sqrt{2|E_{0}|t+12t^{2}}}}\right).$
(82)

The lowest energy of two free particles on the triangular lattice is
$E_{11}=-12t$. When $E\rightarrow E_{11}-0$ from below, the argument of ${\bf
K}$ approaches 1 and $M^{\rm tr}_{\bf 00}$ diverges logarithmically. Utilizing
Eq. (56), one concludes that a bound pair is formed for any attractive $U$.

Let us derive asymptotic behavior of $E$ at small couplings. Setting
$|E_{0}|=12t+\Delta$, $\Delta\ll t$ in Eq. (82), one obtains from Eq. (56)

$\frac{|U|}{4\pi\sqrt{3}t}\>{\bf K}\\!\left(1-\frac{\Delta}{18t}\right)=1\>,$
(83)

from where the binding energy is

$\Delta=144t\,\exp{\left(-\frac{8\pi\\!\sqrt{3}\,t}{|U|}\right)}.$ (84)

###  3.5  3D. Simple cubic lattice

The attractive Hubbard model in three dimensions possesses a new qualitative
feature: nonzero threshold of pair formation. In 3D, kinetic energy alone is
strong enough to counteract weak attraction. Mathematically, the triple
integral in $M_{\bf 00}$ converges when $E\rightarrow E_{11}$ for any ${\bf
P}$, rendering the critical potential depth $|U_{\rm cr}|$ finite. Calculation
of $|U_{\rm cr}|$ requires evaluation of the celebrated Watson integrals
[107]. Amazingly, on BZ diagonals they are now known analytically for a
general lattice point [91, 92, 93, 96, 98], see D. Setting
$P_{x}=P_{y}=P_{z}=P$ and $\alpha=4t\cos{(P/2)}$, one obtains

$M^{\rm sc}_{\bf
00}=\frac{1}{(2\pi)^{3}}\int\limits^{\pi}_{-\pi}\\!\\!\\!\int\limits^{\pi}_{-\pi}\\!\\!\\!\int\limits^{\pi}_{-\pi}\frac{{\rm
d}q_{x}\,{\rm d}q_{y}\,{\rm
d}q_{z}}{|E|-\alpha\,(\cos{q_{x}}+\cos{q_{y}}+\cos{q_{z}})}\>.$ (85)

The minimum energy of two free particles is $E_{11}=-3\alpha$ at which the
classic result of Watson’s reads [107]

$M^{\rm sc}_{\bf
00}(E=-3\alpha)=\frac{4(18+12\sqrt{2}-10\sqrt{3}-7\sqrt{6})}{\alpha\pi^{2}}\>{\bf
K}^{2}\left[(2-\sqrt{3})(\sqrt{3}-\sqrt{2})\right]=\frac{1}{\alpha}\,0.505462\ldots\>.$
(86)

The binding threshold then follows from Eq. (56):

$|U_{\rm cr}({\rm diag})|=(7.913552\ldots)\,t\cos{\left(P/2\right)}\>.$ (87)

The $P=0$ value of $7.914\ldots$ obtained by numerical integration was
reported in Ref. [108]. The threshold decreases along the diagonal and becomes
zero in the BZ corner. This also implies that for any $|U|<7.913552\ldots$,
there is a momentum $P_{0}$ such that the particles are not bound for
$P<P_{0}$ but bound at $P>P_{0}$. Thus, pairs become more stable at large
lattice momenta. It is a common property of lattice bound states, which will
be encountered many times later in this paper.

Figure 4: Binding threshold surface for the negative Hubbard model in the 3D
simple cubic lattice. What is shown is the surface’s intersections with plane
$P_{x}=P_{y}$ for several $|U|$. Filled circles indicate the threshold momenta
on the BZ diagonal, given by
$P_{z}=P_{x}=P_{y}=2\arccos{\left(|U|/7.913552\right)}$.

The binding threshold can also be computed on BZ planes that pass through four
BZ corners, for example on the plane $P_{x}=P_{y}$. This is possible thanks to
the extension of Watson’s result to the anisotropic case by Montroll [109,
98]:

$\displaystyle M^{\rm tg}_{\bf 00}(E=-8t\cos{(P_{x}/2)}-4t\cos{(P_{z}/2)})$
$\displaystyle=\frac{1}{\alpha\pi^{3}}\int\limits^{\pi}_{0}\\!\\!\\!\int\limits^{\pi}_{0}\\!\\!\\!\int\limits^{\pi}_{0}\frac{{\rm
d}q_{x}\,{\rm d}q_{y}\,{\rm
d}q_{z}}{2-\cos{q_{x}}-\cos{q_{y}}+\xi(1-\cos{q_{z}})}$
$\displaystyle=\frac{4}{\alpha\pi^{2}}\sqrt{\frac{2\kappa_{1}\kappa_{2}}{\xi}}\,{\bf
K}(\kappa_{1}){\bf K}(\kappa_{2})\>,$ (88)

where

$\kappa_{1,2}=\frac{1}{2\xi}\left(\sqrt{4+2\xi}\pm
2\right)\\!\left(2\sqrt{1+\xi}-\sqrt{4+2\xi}\right),$ (89) $\alpha\equiv
4t\cos{(P_{x}/2)}\>,\hskip
28.45274pt\xi\equiv\frac{\cos{(P_{z}/2)}}{\cos{(P_{x}/2)}}\>.$ (90)

The threshold value is an inverse of Eq. (88) by virtue of Eq. (56). Notice
how pair’s center-of-mass movement induces anisotropy. For fixed
$|U|<7.913552\ldots$ and $E=E_{11}$, Eq. (56) defines a surface in BZ that
separates bound and unbound states. Figure 4 shows the intersection of that
surface with the $P_{x}=P_{y}$ plane for several values of $|U|$.

The ground state energy, $E_{0}=E({\bf P}=0)$, is discussed next. A close-form
expression for $M^{\rm sc}_{\bf 00}$ in the isotropic simple cubic model was
found by Joyce [91, 92]

$M^{\rm sc}_{\bf
00}(E_{0}\leq-12t)=\frac{1}{|E_{0}|}\frac{(1-9\zeta^{4})}{(1-\zeta)^{3}(1+3\zeta)}\left[\frac{2}{\pi}{\bf
K}(\kappa)\right]^{2},$ (91)
$\kappa^{2}(\zeta)=\frac{16\,\zeta^{3}}{(1-\zeta)^{3}(1+3\zeta)}\>,$ (92)
$\zeta=\zeta(w)=\left[\frac{1-\sqrt{1-\frac{w^{2}}{9}}}{1+\sqrt{1-w^{2}}}\right]^{\frac{1}{2}},\hskip
14.22636ptw=\frac{12\,t}{|E_{0}|}\>.$ (93)

Using these formulas, Eq. (56) defines $E_{0}$ as a function of $|U|$. It is
plotted in Fig. 3(a).

To derive effective mass, we write
$|E|=|E_{0}|-\frac{3\hbar^{2}P^{2}}{2m^{\ast}_{p}}$ and expand Eq. (85) for
small $P$ to get

$\frac{m_{0}}{m^{\ast}_{p}}=\frac{1}{6}\frac{\int\limits^{\pi}_{-\pi}\\!\\!\\!\int\limits^{\pi}_{-\pi}\\!\\!\\!\int\limits^{\pi}_{-\pi}\frac{{\rm
d}q_{x}\,{\rm d}q_{y}\,{\rm
d}q_{z}\,(\cos{q_{x}}+\cos{q_{y}}+\cos{q_{z}})}{[|E_{0}|-4t\,(\cos{q_{x}}+\cos{q_{y}}+\cos{q_{z}})]^{2}}}{\int\limits^{\pi}_{-\pi}\\!\\!\\!\int\limits^{\pi}_{-\pi}\\!\\!\\!\int\limits^{\pi}_{-\pi}\frac{{\rm
d}q_{x}\,{\rm d}q_{y}\,{\rm
d}q_{z}}{[|E_{0}|-4t\,(\cos{q_{x}}+\cos{q_{y}}+\cos{q_{z}})]^{2}}}\>.$ (94)

Once $E_{0}$ is found from Eq. (56), the effective mass can be obtained by
numerical integration. Alternatively, the denominator in Eq. (94) is
$-\partial M^{\rm sc}_{\bf 00}/\partial|E_{0}|$, while the numerator can be
expressed using a similar derivative and a linear transformation to link it
also to $M^{\rm sc}_{\bf 00}$. The final formula is

$\frac{m_{0}}{m^{\ast}_{p}}=\frac{1}{24t}\left\\{|E_{0}|+\frac{1}{|U|(\partial
M^{\rm sc}_{\bf 00}/\partial|E_{0}|)}\right\\},$ (95)

where we have used $M^{\rm sc}_{\bf 00}=1/|U|$. Thus, the effective mass can
also be computed by numerical differentiation of $M^{\rm sc}_{\bf 00}$.
$m^{\ast}_{p}(U)$ is shown in Fig. 3(b).

To compute effective radius, we apply the 3D version of Eq. (78). For ${\bf
P}=0$, it yields

$r^{\ast
2}_{pj}=16t^{2}\>\frac{\int\limits^{\pi}_{-\pi}\\!\\!\\!\int\limits^{\pi}_{-\pi}\\!\\!\\!\int\limits^{\pi}_{-\pi}\frac{\sin^{2}{q_{j}}\>{\rm
d}q_{x}\,{\rm d}q_{y}\,{\rm
d}q_{z}}{[|E_{0}|-4t\,(\cos{q_{x}}+\cos{q_{y}}+\cos{q_{z}})]^{4}}}{\int\limits^{\pi}_{-\pi}\\!\\!\\!\int\limits^{\pi}_{-\pi}\\!\\!\\!\int\limits^{\pi}_{-\pi}\frac{{\rm
d}q_{x}\,{\rm d}q_{y}\,{\rm
d}q_{z}}{[|E_{0}|-4t\,(\cos{q_{x}}+\cos{q_{y}}+\cos{q_{z}})]^{2}}}\>.$ (96)

$r^{\ast}_{px}(U)$ is shown in Fig. 3(c).

Figure 5: Properties of bound pairs in the anisotropic tetragonal attractive
Hubbard model. (a) Binding threshold as a function of hopping anisotropy. (b)
Anisotropy of effective masses. (c) Anisotropy of effective radii. Circles
indicate threshold $t_{z}$ values for given $|U|$.

###  3.6  3D. Tetragonal lattice

The tetragonal lattice is of interest because of its relevance to high-$T_{c}$
superconductors. We explore several effects of hopping anisotropy in this
section. Equations (88) and (89) apply to the tetragonal attractive Hubbard
model, where the anisotropy is caused not by pair momentum but by the
difference in transfer integrals along the $z$ and $x,y$ directions. Limiting
consideration to the ground state, one sets $t_{x}=t_{y}=t$, $\alpha=4t$ and
$\xi=t_{z}/t$. Figure 5(a) shows the binding threshold as a function of $\xi$.
The function is sharp near $\xi=0$: already at $\xi=10^{-4}$, $|U_{\rm
cr}|=1.98\,t$. Even a tiny interlayer hopping suppresses long-range
logarithmic fluctuations and renders the system essentially three-dimensional.
An analytical expression for the Green function $M^{\rm tg}_{\bf 00}$ for
${\bf P}=0$ and arbitrary $E_{0}$ was derived by Delves, Joyce, and Zucker
[94, 95, 97]. The cumbersome expression is given in D.2.

Of physical interest is pair mass anisotropy $m^{\ast}_{px}/m^{\ast}_{pz}$.
Expanding $|U|M^{\rm tg}_{\bf 00}=1$ for small ${\bf P}$ one obtains

$\frac{m^{\ast}_{px}}{m^{\ast}_{pz}}=\frac{t_{z}}{t}\>\frac{\int\limits^{\pi}_{-\pi}\\!\\!\\!\int\limits^{\pi}_{-\pi}\\!\\!\\!\int\limits^{\pi}_{-\pi}\frac{\cos{q_{z}}\>{\rm
d}q_{x}\,{\rm d}q_{y}\,{\rm
d}q_{z}}{[|E_{0}|-4t\,(\cos{q_{x}}+\cos{q_{y}})-4t_{z}\cos{q_{z}}]^{2}}}{\int\limits^{\pi}_{-\pi}\\!\\!\\!\int\limits^{\pi}_{-\pi}\\!\\!\\!\int\limits^{\pi}_{-\pi}\frac{\cos{q_{x}}\>{\rm
d}q_{x}\,{\rm d}q_{y}\,{\rm
d}q_{z}}{[|E_{0}|-4t\,(\cos{q_{x}}+\cos{q_{y}})-4t_{z}\cos{q_{z}}]^{2}}}\>.$
(97)

The mass ratio is shown in Fig. 5(b). For a single particle,
$m_{x}/m_{z}=t_{z}/t$, and the graph would be a straight line. Pairing
enhances mass anisotropy which approaches $(t_{z}/t)^{2}$ in the strong
coupling limit. For intermediate $|U|$, the mass anisotropy lies between
$(t_{z}/t)$ and $(t_{z}/t)^{2}$.

Another interesting property is the ratio of effective radii. It follows from
Eq. (96) that

$\left(\frac{r^{\ast}_{pz}}{r^{\ast}_{px}}\right)^{2}=\frac{t^{2}_{z}}{t^{2}}\>\frac{\int\limits^{\pi}_{-\pi}\\!\\!\\!\int\limits^{\pi}_{-\pi}\\!\\!\\!\int\limits^{\pi}_{-\pi}\frac{\sin^{2}{q_{z}}\>{\rm
d}q_{x}\,{\rm d}q_{y}\,{\rm
d}q_{z}}{[|E_{0}|-4t\,(\cos{q_{x}}+\cos{q_{y}})-4t_{z}\cos{q_{z}}]^{4}}}{\int\limits^{\pi}_{-\pi}\\!\\!\\!\int\limits^{\pi}_{-\pi}\\!\\!\\!\int\limits^{\pi}_{-\pi}\frac{\sin^{2}{q_{x}}\>{\rm
d}q_{x}\,{\rm d}q_{y}\,{\rm
d}q_{z}}{[|E_{0}|-4t\,(\cos{q_{x}}+\cos{q_{y}})-4t_{z}\cos{q_{z}}]^{4}}}\>.$
(98)

The pair size anisotropy is shown in Fig. 5(c). Near binding threshold,
$r^{\ast}_{pz}/r^{\ast}_{px}=\sqrt{t_{z}/t}$ whereas in general it is confined
between $\sqrt{t_{z}/t}$ and $(t_{z}/t)$.

###  3.7  3D. Body-centered cubic (BCC) lattice

In a BCC lattice with nearest-neighbor hopping, $\varepsilon_{\bf
k}=-8t\cos{\frac{k_{x}}{2}}\cos{\frac{k_{y}}{2}}\cos{\frac{k_{z}}{2}}$. (In
this and the following sections, we set the cube length, $a=1$.) First, we
consider the ground state, ${\bf P}=0$. In calculating $M^{\rm bcc}_{\bf 00}$,
integration over the BCC Brillouin zone can be replaced with one fourth of the
integral over a cube with side length $4\pi$. Changing momentum variables
yields

$M^{\rm bcc}_{\bf
00}=\frac{1}{(2\pi)^{3}}\int\limits^{\pi}_{-\pi}\\!\\!\\!\int\limits^{\pi}_{-\pi}\\!\\!\\!\int\limits^{\pi}_{-\pi}\frac{{\rm
d}q_{x}\,{\rm d}q_{y}\,{\rm
d}q_{z}}{|E_{0}|-16t\,\cos{q_{x}}\cos{q_{y}}\cos{q_{z}}}\>.$ (99)

The integral here is one of the generalized Watson integrals first evaluated
by Maradudin [110, 111] and later studied by other authors [112, 113, 114].
Application of those results leads to the energy equation

$\frac{|U|}{|E_{0}|}\left(\frac{2}{\pi}\right)^{2}{\bf
K}^{2}\left[\sqrt{\frac{1}{2}-\frac{1}{2}\sqrt{1-\left(\frac{16t}{|E_{0}|}\right)^{2}}}\right]=1\>.$
(100)

The binding threshold is found by setting $E_{0}=-16t$, which results in

$|U^{\rm bcc}_{\rm
cr}|=\frac{4\pi^{2}\,t}{\left[K\\!\\!\left(\frac{1}{\sqrt{2}}\right)\right]^{2}}=\frac{64\pi^{3}\,t}{\left[\Gamma\\!\\!\left(\frac{1}{4}\right)\right]^{4}}=(11.484320\ldots)\,t\>.$
(101)

Expanding Eq. (100) near $E_{0}\approx-16t$ yields a quadratic dependence of
the binding energy near the threshold

$E_{0}(|U|\approx|U^{\rm bcc}_{\rm
cr}|)=-16\,t-\frac{\left[\Gamma\\!\\!\left(\frac{1}{4}\right)\right]^{16}}{2^{15}\pi^{10}}\frac{(|U|-|U^{\rm
bcc}_{\rm cr}|)^{2}}{t}\>.$ (102)

The numerical coefficient at the quadratic term is $=0.290501\ldots$. Formula
(102) is derived in E. Finally, expanding Eq. (100) at large $E_{0}$, one
obtains

$E_{0}(|U|\rightarrow\infty)=-|U|-\frac{32\,t^{2}}{|U|}+o\\!\left(\frac{t^{2}}{|U|}\right),$
(103)

which is consistent with strong-coupling perturbation theory.

Nonzero pair momenta are discussed next. Cases with one nonzero component, for
example ${\bf P}=(P_{x},0,0)$, can be reduced to the ground state case. The
energy denominator in $M^{\rm bcc}_{\bf 00}$ contains

$\varepsilon_{\bf q}+\varepsilon_{{\bf P}-{\bf
q}}=-16\,t\cos{\frac{P_{x}}{4}}\cos{\left(\frac{P_{x}}{4}-q_{x}\right)}\cos{q_{y}}\cos{q_{z}}\>.$
(104)

After shifting the integration variable $q_{x}$, $M^{\rm bcc}_{\bf 00}$ is
reduced to the ground state expression, Eq. (99), where $t$ is replaced with
$t\cos{(P_{x}/4)}$. Without rederiving all the results given above, let us
just mention a generalization of the pair binding condition:

$|U^{\rm bcc}_{\rm
cr}(P_{x},0,0)|=(11.484320\ldots)\,t\cos{\left(P_{x}/4\right)}\>.$ (105)

This expression should be compared with the simple cubic value, Eq. (87).

Lattice | $z$ | $|U_{\rm cr}|/t$ | $|U_{\rm cr}|/(zt)$
---|---|---|---
Simple cubic | 6 | $7.913552\ldots$ | $1.318925\ldots$
Body-centered cubic | 8 | $11.484320\ldots$ | $1.435540\ldots$
Face-centered cubic | 12 | $17.848362\ldots$ | $1.487363\ldots$
Table 1:  Pair binding thresholds for the attractive Hubbard model in the
three cubic lattices with isotropic nearest-neighbor hopping. $z$ is the
number of nearest neighbors. ${\bf P}=0$.

###  3.8  3D. Face-centered cubic (FCC) lattice

Similar to the FCC case, integration over the FCC Brillouin zone can be
replaced by integration over the cube with side length $(4\pi)$. In the ground
state, $\varepsilon_{\bf q}=\varepsilon_{{\bf P}-{\bf q}}$, and $M^{\rm
fcc}_{\bf 00}$ reduces to

$M^{\rm fcc}_{\bf
00}=\frac{1}{8t\pi^{3}}\int\limits^{\pi}_{0}\\!\\!\int\limits^{\pi}_{0}\\!\\!\int\limits^{\pi}_{0}\frac{{\rm
d}q_{x}\,{\rm d}q_{y}\,{\rm
d}q_{z}}{\frac{|E_{0}|}{8t}-\cos{q_{x}}\cos{q_{y}}-\cos{q_{y}}\cos{q_{z}}-\cos{q_{z}}\cos{q_{x}}}\>.$
(106)

This triple integral was first evaluated by Iwata [115] and later in a
different form by Joyce [92]. Joyce’s result reads

$M^{\rm fcc}_{\bf
00}(E_{0}<-24t)=\frac{1}{|E_{0}|}\frac{(1+3\zeta^{2})^{2}}{(1-\zeta)^{3}(1+3\zeta)}\left[\frac{2}{\pi}{\bf
K}(\kappa)\right]^{2},$ (107)
$\kappa^{2}(\zeta)=\frac{16\,\zeta^{3}}{(1-\zeta)^{3}(1+3\zeta)}\>,$ (108)
$\zeta=\zeta(w)=\frac{-1+\sqrt{1+\frac{w}{3}}}{1+\sqrt{1-w}}\>,\hskip
14.22636ptw=\frac{24\,t}{|E_{0}|}\>.$ (109)

Pair formation occurs at $E_{0}=-24\,t$ or $w=1$. The binding threshold reads

$|U^{\rm fcc}_{\rm cr}|=\frac{8\pi^{2}\,t}{\sqrt{3}\left[{\bf
K}\\!\\!\left(\frac{\sqrt{3}-1}{2\sqrt{2}}\right)\right]^{2}}=(17.848362\ldots)\,t\>.$
(110)

Table 1 summarizes threshold values for the three cubic lattices.

##  4  1D. $UV$ model on the one-dimensional chain

We now transition to $UV$ models with on-site repulsion $U$ and nearest-
neighbor attraction $V_{\bf b}\equiv-|V|$. A general feature of these models
is the existence of multiple bound states, whose number increases with lattice
dimensionality. Because of the complexity of general energy equation, Eq.
(17), it is advantageous to utilize the (anti)symmetrized formalism developed
in Section 2.3. We begin with the one-dimensional $UV$ model [87].

###  4.1  Singlet states

The symmetrized set of neighbor vectors consists of two elements: $\\{{\bf
b}_{+}\\}=\\{(0),(1)\\}$. Symmetrized Schrödinger equation, Eq. (23), reads:

$(E_{s}-\varepsilon_{k_{1}}-\varepsilon_{k_{2}})\,\phi^{+}_{k_{1}k_{2}}=U\frac{1}{N}\sum_{q}\phi^{+}_{q,\>k_{1}+k_{2}-q}-|V|\frac{1}{N}\sum_{q}\phi^{+}_{q,\>k_{1}+k_{2}-q}e^{iq}\left(e^{-ik_{1}}+e^{-ik_{2}}\right).$
(111)

Here $\varepsilon_{k}=-2t\cos{k}$, and subscript in $E_{s}$ indicates “spin-
singlet”. We introduce two auxiliary functions:

$\displaystyle\Phi^{+}_{0}(P)$
$\displaystyle=\frac{1}{N}\sum_{q}\phi^{+}_{q,\>k_{1}+k_{2}-q}\>,$ (112)
$\displaystyle\Phi^{+}_{1}(P)$
$\displaystyle=\frac{1}{N}\sum_{q}\phi^{+}_{q,\>k_{1}+k_{2}-q}\>e^{iq}\>,$
(113)

so that

$\phi^{+}_{k_{1}k_{2}}=\frac{U\Phi^{+}_{0}(P)-|V|\Phi^{+}_{1}(P)\left(e^{-ik_{1}}+e^{-ik_{2}}\right)}{E_{s}-\varepsilon_{k_{1}}-\varepsilon_{k_{2}}}\>.$
(114)

Substituting Eq. (114) back into the definitions, Eqs. (112) and (113), one
obtains:

$\displaystyle\Phi^{+}_{0}(P)$
$\displaystyle=U\Phi^{+}_{0}(P)\left(\frac{1}{N}\sum_{q}\frac{1}{E_{s}-\varepsilon_{q}-\varepsilon_{P-q}}\right)-|V|\Phi^{+}_{1}(P)\left(\frac{1}{N}\sum_{q}\frac{e^{-iq}+e^{-i(P-q)}}{E_{s}-\varepsilon_{q}-\varepsilon_{P-q}}\right),$
(115) $\displaystyle\Phi^{+}_{1}(P)$
$\displaystyle=U\Phi^{+}_{0}(P)\left(\frac{1}{N}\sum_{q}\frac{e^{iq}}{E_{s}-\varepsilon_{q}-\varepsilon_{P-q}}\right)-|V|\Phi^{+}_{1}(P)\left(\frac{1}{N}\sum_{q}\frac{e^{iq}(e^{-iq}+e^{-i(P-q)})}{E_{s}-\varepsilon_{q}-\varepsilon_{P-q}}\right).$
(116)

Next, a change of variables $q^{\prime}=q-\frac{P}{2}$ under the integrals
results in a $(2\times 2)$ matrix equation

$\displaystyle\Phi^{+}_{0}$
$\displaystyle=-UM_{0}\cdot\Phi^{+}_{0}+2|V|\,e^{-i\frac{P}{2}}M_{1}\cdot\Phi^{+}_{1}\>,$
(117) $\displaystyle\Phi^{+}_{1}$
$\displaystyle=-U\,e^{i\frac{P}{2}}M_{1}\cdot\Phi^{+}_{0}+|V|(M_{0}+M_{2})\cdot\Phi^{+}_{1}\>,$
(118)

where

$M_{n}=\frac{1}{N}\sum_{q}\frac{\cos{nq}}{|E_{s}|-4t\cos{(P/2)}\cos{q}}=\frac{1}{\sqrt{|E_{s}|^{2}-\alpha^{2}}}\left[\frac{\sqrt{|E_{s}|^{2}-\alpha^{2}}-|E_{s}|}{-\alpha}\right]^{n},$
(119)

and $\alpha\equiv 4t\cos{(P/2)}$. Note that

$\displaystyle M_{1}$
$\displaystyle=\frac{1}{\alpha}\left(|E_{s}|M_{0}-1\right)\>,$ (120)
$\displaystyle M_{2}$ $\displaystyle=\frac{2|E_{s}|}{\alpha}M_{1}-M_{0}\>.$
(121)

As a result, everything can be expressed via the basic integral $M_{0}$. The
bound state’s energy is determined by the consistency condition of Eqs. (117)
and (118). Expanding the determinant, one obtains

$(UM_{0}+1)+\frac{2}{\alpha^{2}}\,|V|(|E_{s}|+U)(1-|E_{s}|M_{0})=0\>.$ (122)

This form will be useful later in comparing with similar equations in higher
dimensions. Substitution of the explicit form of $M_{0}$ yields the final
expression:

$\left[U+\sqrt{|E_{s}|^{2}-16t^{2}\cos^{2}{\frac{P}{2}}}\right]\left[|E_{s}|+\sqrt{|E_{s}|^{2}-16t^{2}\cos^{2}{\frac{P}{2}}}\right]-2|V|\left(U+|E_{s}|\right)=0\>.$
(123)

This is a cubic equation for $E_{s}(P)$, which does not have a simple-form
analytical solution. Only at BZ boundary the equation simplifies to
$E_{s}(\pm\pi)=-|V|$. At strong coupling, $U,|V|,|E_{s}|\gg t$, Eq. (123)
yields for the ground state

$E_{s}(P=0)=-|V|-\frac{4t^{2}}{|V|}-\frac{8t^{2}}{U+|V|}+o(t^{2}/|V|)\>,$
(124)

which is consistent with second-order perturbation theory.

Figure 6: Singlet-pair dispersion $E_{s}$ in the 1D $UV$ model computed from
Eq. (123) (thick solid line). $U=15\,t$, $|V|=1.3\,t$. Note that $|V|$ is
below the $P=0$ binding threshold, $2Ut/(U+4t)$. The shaded area is the
continuum of two-particle scattering states. The circles mark threshold
momenta, $P=\pm 2\arccos{[UV/2/(U-2V)]}$, derived from Eq. (125). The triplet
pair energy $E_{t}$ given by Eq. (145) is very close to $E_{s}$ and therefore
not shown in the figure.

In order to obtain binding threshold, set $E_{s}$ equal to the lowest energy
of two free carriers, $E_{11}=-4t\cos{(P/2)}$, in Eq. (123). It results in

$|V^{s}_{\rm cr}|=\frac{2Ut\cos{(P/2)}}{U+4t\cos{(P/2)}}\>.$ (125)

This formula possesses several interesting properties. First of all,
$|V^{s}_{\rm cr}|$ is nonzero despite the model being one-dimensional. Here,
the attraction competes not only with kinetic energy but also with repulsion
$U$ which leads to a nonzero threshold. At weak repulsion, $U<t$, one has
$|V^{s}_{\rm cr}|\approx U/2$ which can be understood from the Born
approximation: there are two attractive sites for one repulsive site, hence
half as strong $V$ is needed to overcome $U$. In the opposite limit of strong
repulsion, $|V^{s}_{\rm cr}|$ approaches a finite limit. At large $U$, the on-
site wave function amplitude, $\Psi(m,m)\rightarrow 0$, which becomes a
boundary condition for the rest of $\Psi$. Once attraction is strong enough to
produce a bound state in the presence of this zero, further increase of $U$
has no effect. Finally, $|V^{s}_{\rm cr}|$ is a strong function of pair
momentum. Like in 3D Hubbard models, the pair becomes more stable at large
$P$. At BZ boundary, the pair is always stable for any $U$ however large, and
any $|V|$ however small.

A typical pair dispersion for $|V|<|V^{s}_{\rm cr}|$ is shown in Fig. 6.
Consider weakly bound pairs with $P=0$. Expanding the exact dispersion
relation, Eq. (123), for $E_{s}\approx-4t$ and utilizing Eq. (125), one
obtains

$E_{s}(P=0)\approx-4t-\frac{(|V|-|V^{s}_{{\rm cr},0}|)^{2}}{2t}\>,$ (126)

where $|V^{s}_{{\rm cr},0}|=2Ut/(U+4t)$. Thus the pair energy varies
quadratically near the threshold. Expanding Eq. (123) for small $P$, one
obtains, after transformations, the singlet’s effective mass

$\frac{m^{\ast}_{ps}}{m_{0}}\\!=\\!\frac{(U\\!-\\!2E_{s}\\!-\\!2|V|)\sqrt{E^{2}_{s}\\!-\\!(4t)^{2}}\\!+\\!2E^{2}_{s}\\!-\\!UE_{s}\\!-\\!16t^{2}}{(2t)\left[2\sqrt{E^{2}_{s}-(4t)^{2}}+U-E_{s}\right]}\>.$
(127)

At threshold, $E_{s}=-4t$, and the last formula yields
$m^{\ast}_{ps}/m_{0}=2$, as expected. At strong attraction, the mass generally
grows linearly with coupling, $m^{\ast}_{ps}/m_{0}\propto|V|/t$, with the
slope depending on $U$ and the $U/V$ ratio. As an example, Fig. 7 shows the
pair mass for $U=20\,t$.

Figure 7: Singlet pair mass in the 1D $UV$ model computed from Eqs. (123) and
(127). The thin solid line corresponds to a fixed $U=20\,t$. The mass grows
linearly with $|V|$. The circle marks threshold value, Eq. (125), for $P=0$.
The thick solid line corresponds to the light pairs regime with $U=-|V|$. The
mass stays of order $m_{0}$ for all $V$. At strong coupling,
$|V|\rightarrow\infty$, the mass approaches the theoretical limit,
$m^{\ast}_{ps}=\sqrt{8}\,m_{0}$ (the dashed line), see Eq. (135).

###  4.2  Light bound pairs

In this section, we introduce the topic of light bound pairs that will be
discussed later in several places of this paper. We use this term to describe
a situation when pairs are strongly bound, $E\rightarrow-\infty$, but at the
same time remain mobile with an effective mass of order $m_{0}$. It occurs
when a pair can move through the lattice without changing its energy, i.e.,
without breaking the most attractive bond. There are two primary reasons why
it can happen. First, because of geometry. In some lattices such as triangular
and FCC, one member of a pair can hop to another site while still remaining a
nearest neighbor to the second member, which keeps the configuration energy
unchanged at $-|V|$. This is followed by a similar hop by the second member.
The two particles hop in turns in a “crab-like” fashion, which results in
overall movement of the pair through the system. The second origin of light
pairs is a flat segment in the attractive part of the inter-particle
potential. In this case, the particles can also move in alternating order
without changing their energy.

In this section, we use the relative simplicity of the 1D $UV$ model to
illustrate the second mechanism. To this end, we set $U=-|V|$ in the formulas
of Section 4.1. Additionally, the strong-coupling limit,
$U=V\rightarrow-\infty$, can be treated analytically. Consider Fig. 8(b). It
is sufficient to include only two types of spin-singlet configurations:

$\displaystyle A_{m}$ $\displaystyle=|\uparrow\downarrow\rangle_{m}\>,$ (128)
$\displaystyle B_{m}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(|\uparrow\rangle_{m}|\downarrow\rangle_{m+1}+|\downarrow\rangle_{m}|\uparrow\rangle_{m+1}\right)\>.$
(129)

Hamiltonian action within this basis is

$\displaystyle\hat{H}A_{m}$
$\displaystyle=-\sqrt{2}tB_{m}-\sqrt{2}tB_{m-1}\>,$ (130)
$\displaystyle\hat{H}B_{m}$
$\displaystyle=-\sqrt{2}tA_{m}-\sqrt{2}tA_{m+1}\>.$ (131)

The Schrödinger equation in momentum space is

$\displaystyle\tilde{E}A_{P}$
$\displaystyle=-\sqrt{2}t\left(1+e^{-iPa}\right)B_{P}\>,$ (132)
$\displaystyle\tilde{E}B_{P}$
$\displaystyle=-\sqrt{2}t\left(1+e^{iPa}\right)A_{P}\>,$ (133)

where $\tilde{E}$ is pair energy counted from $-|V|$ and $a$ is the lattice
constant. Band dispersion is

$\tilde{E}_{1,2}(P)=\pm 2\sqrt{2}t\cos{\frac{Pa}{2}}\>,$ (134)

which corresponds to an effective mass

$m^{\ast}_{ps}=\sqrt{2}\,\frac{\hbar^{2}}{ta^{2}}=2\sqrt{2}\,m_{0}\>.$ (135)

Thus, pair mass remains of the order of free-particle mass $m_{0}$ even in the
limit of infinitely strong attraction. Figure 7 shows a numerical solution of
Eqs. (123) and (127) for the resonant potential $U=-|V|$. Indeed, the pair
mass never exceeds the strong-coupling limit $\sqrt{8}\,m_{0}$.

One might think of an attractive interaction with $U=-|V|$ as exotic, but it
can potentially be realized in cold gases where both $U$ and $V$ can be
independently controlled. In crystalline solids, one can envision more
realistic potentials comprising a strong repulsive core and a long-range
attractive tail. Such a potential will necessarily have a minimum at a finite
separation between particles. If the minimum is wide compared with the
interatomic distance, then with high probability there will be two separations
with equal attractive strengths. Such a situation is illustrated in Fig. 8(c),
where attraction on the third and fourth nearest neighbors are assumed equal,
$V_{3}=V_{4}$. (Importantly, other parts of the potential do not change the
argument given below because if the particles are allowed to access
configurations outside of the $-|V|$ basis, the pair mass will only decrease!)
A proper ground state basis in this example is

$A_{m}=|\bullet\rangle_{m}|\bullet\rangle_{m+3}\>;\hskip
14.22636ptB_{m}=|\bullet\rangle_{m}|\bullet\rangle_{m+4}\>.$ (136)

(Since the particles cannot really exchange, there is no need to consider spin
degrees of freedom. Singlet and triplet pairs will have the same mass.) The
Schrödinger equation reads

$\displaystyle\tilde{E}A_{P}$
$\displaystyle=-t\left(1+e^{-iPa}\right)B_{P}\>,$ (137)
$\displaystyle\tilde{E}B_{P}$ $\displaystyle=-t\left(1+e^{iPa}\right)A_{P}\>,$
(138)

which yields $\tilde{E}_{1,2}(P)=\pm 2t\cos{\frac{Pa}{2}}$ and

$m^{\ast}=\frac{2\hbar^{2}}{ta^{2}}=4\,m_{0}\>.$ (139)

Thus, the bound pair is no heavier than just four free particle masses. Note
that this conclusion does not depend on the separation distance at which the
flat section of the potential occurs.

Figure 8: Illustration of the light bound pair mechanism. (a) Conventional
pair movement. The intermediate configuration has a larger energy, which
results in a mass growing linearly with $|V|$. (b) In the resonant case,
$U=-|V|$, there are intermediate configurations with the same energy as the
initial and starting configurations. (c) An attractive potential with a flat
section at a nonzero separation between two particles. In this case,
$V_{3}=V_{4}$.

###  4.3  Triplet states

The antisymmetrized set of vectors consists of just one element $\\{{\bf
b}_{-}\\}=\\{(1)\\}$ and there is one basis function $\Phi^{-}_{1}(P)$. An
antisymmetrized Schrödinger equation, Eq. (40), reads:

$(E_{t}-\varepsilon_{k_{1}}-\varepsilon_{k_{2}})\,\phi^{-}_{k_{1}k_{2}}=-|V|\frac{1}{N}\sum_{q}\phi^{-}_{q,\>k_{1}+k_{2}-q}e^{iq}\left(e^{-ik_{1}}-e^{-ik_{2}}\right).$
(140)

In terms of the auxiliary function

$\Phi^{-}_{1}(P)=\frac{1}{N}\sum_{q}\phi^{-}_{q,\>k_{1}+k_{2}-q}\>e^{iq}\>,$
(141)

the pair wave function is expressed as

$\phi^{-}_{k_{1}k_{2}}=\frac{-|V|\Phi^{-}_{1}(P)\left(e^{-ik_{1}}-e^{-ik_{2}}\right)}{E-\varepsilon_{k_{1}}-\varepsilon_{k_{2}}}\>.$
(142)

Substituting Eq. (142) back in the definition, Eq. (141), one obtains:

$\Phi^{-}_{1}(P)=-|V|\Phi^{-}_{1}(P)\left(\frac{1}{N}\sum_{q}\frac{e^{iq}(e^{-iq}-e^{-i(P-q)})}{E-\varepsilon_{q}-\varepsilon_{P-q}}\right).$
(143)

Changing variables $q^{\prime}=q-\frac{P}{2}$ in the last equation yields

$\left\\{1-|V|\,(M_{0}-M_{2})\right\\}\Phi^{-}_{1}=0\>.$ (144)

Note that it is independent of $U$, as expected for a triplet pair. Direct
calculation results in

$E_{t}(P)=-|V|-\frac{4t^{2}\cos^{2}{\frac{P}{2}}}{|V|}\>,$ (145)

for the triplet energy, and in

$\frac{m^{\ast}_{pt}}{m_{0}}=\frac{|V|}{t}\>,$ (146)

for the triplet effective mass, where $m_{0}=\hbar^{2}/(2ta^{2})$. The triplet
pair is stable when

$|V|>|V^{t}_{\rm cr}|=2t\,\cos{\frac{P}{2}}\>.$ (147)

Notice that in the $U\rightarrow\infty$ limit, $|V^{s}_{\rm cr}|=|V^{t}_{\rm
cr}|$ and $E_{s}=E_{t}$. The phase diagram of the 1D $UV$ model at $P=0$ is
shown in Fig. 9(a).

Figure 9: Two-fermion phase diagrams in the 1D and 2D square $UV$ models at
${\bf P}=0$.

##  5  2D. $UV$ model on the square lattice

Two-dimensional lattice models at low carrier density have been popular in the
studies of HTSC because most high-$T_{c}$ superconductors including the copper
oxides are highly anisotropic. According to one point of view,
superconductivity in cuprates is essentially two-dimensional. The (repulsive)
2D Hubbard model [116] and its derivative, the 2D t-J model [59], were both
put forward as capturing the essential physics [117]. Although this simple
picture is being increasingly challenged [118, 119, 120, 121], pure 2D models
possess rich physics and remain popular in the fields of HTSC [122] and cold
gases [64]. For the purposes of this review, one should mention that the t-J
model “in the hole-rich regime” [33, 60, 61, 62] bears similarities with the
$UV$ model studied here and many of the results derived later in this section
apply equally to both models.

###  5.1  Singlet states. $\Gamma$-point

The symmetrized set of neighbor vectors consists of three elements: $\\{{\bf
b}_{+}\\}=\\{(0,0),(1,0),(0,1)\\}$. In writing down the $(3\times 3)$ system,
Eqs. (27) and (28), it is convenient to shift the inner variables in $M^{+}$:
$q^{\prime}_{x}=q_{x}-\frac{P_{x}}{2}$ and
$q^{\prime}_{y}=q_{y}-\frac{P_{y}}{2}$. That leads to a new set of functions:
$\tilde{\Phi}^{+}_{\bf 0}=\Phi^{+}_{00}$, $\tilde{\Phi}^{+}_{\bf
x}=e^{-i(P_{x}/2)}\Phi^{+}_{10}$, and $\tilde{\Phi}^{+}_{\bf
y}=e^{-i(P_{y}/2)}\Phi^{+}_{01}$. In terms of the new set, the consistency
condition reads

$\left|\begin{array}[]{ccc}1+UM_{00}&-2|V|M_{10}&-2|V|M_{01}\\\
UM_{10}&1-|V|(M_{20}+M_{00})&-2|V|M_{11}\\\
UM_{01}&-2|V|M_{11}&1-|V|(M_{02}+M_{00})\end{array}\right|=0\>,$ (148)

where

$M_{nm}=\int\limits^{\pi}_{-\pi}\\!\\!\int\limits^{\pi}_{-\pi}\frac{{\rm
d}q_{x}\,{\rm
d}q_{y}}{(2\pi)^{2}}\frac{\cos{nq_{x}}\cos{mq_{y}}}{|E|-\alpha\cos{q_{x}}-\beta\cos{q_{y}}}\>,$
(149)

$\alpha=4t\cos{(P_{x}/2)}$, and $\beta=4t\cos{(P_{y}/2)}$. $M_{00}$ was given
in Eq. (72). Other matrix elements in Eq. (148) can also be expressed via
complete elliptic integrals, see B.1 and B.5. The double integrals can also be
computed numerically.

At the $\Gamma$ point, $P_{x}=P_{y}=0$ and $\alpha=\beta=4t$. In this case,
$M_{10}=M_{01}$, $M_{20}=M_{02}$, and Eq. (148) acquires additional symmetry.
Introducing a new basis $\Phi^{+}_{\bf 0}$,
$\Phi^{+}_{s}=\frac{1}{2}(\Phi^{+}_{\bf x}+\Phi^{+}_{\bf y})$, and
$\Phi^{+}_{d}=\frac{1}{2}(\Phi^{+}_{\bf x}-\Phi^{+}_{\bf y})$, the equation
splits into $s$-symmetric and $d$-symmetric sectors. The $s$-sector involves
functions $\Phi^{+}_{\bf 0}$ and $\Phi^{+}_{s}$ and its consistency condition
reads

$\left|\begin{array}[]{cc}1+UM_{00}&-4|V|M_{10}\\\
UM_{10}&1-|V|(M_{20}+M_{00}+2M_{11})\end{array}\right|=0\>.$ (150)

The $d$-sector equation is obtained by subtracting the last two lines of Eq.
(148). It involves only one function $\Phi^{+}_{d}$ and does not include $U$:

$\left\\{1-|V|(M_{20}+M_{00}-2M_{11})\right\\}\Phi^{+}_{d}=0\>.$ (151)

We begin analysis with the $s$-symmetrical ground state described by Eq.
(150). First, we note that the combination
$M_{20}+M_{00}+2M_{11}=\frac{2|E|}{\alpha}M_{10}$ can be expressed via
$M_{10}$, and the latter can be expressed via $M_{00}$ as
$M_{10}=\frac{1}{2\alpha}(|E|M_{00}-1)$. Thus, all the matrix elements in Eq.
(150) are expressible via the base integral $M_{00}$. Expanding the
determinant, one obtains

$(UM_{00}+1)+\frac{1}{\alpha^{2}}|V|(|E_{s}|+U)(1-|E_{s}|M_{00})=0\>,$ (152)

where

$M_{00}=\frac{2}{\pi|E_{s}|}\,{\bf
K}\left(\frac{2\alpha}{|E_{s}|}\right);\hskip 14.22636pt\alpha=4t\>.$ (153)

Equation (152) determines the energy of $s$-states in the $\Gamma$ point.
Depending on the values of $U$ and $V$, there may be one, two, or no bound
states. Equation (152) should be compared with its 1D counterpart, Eq. (122).
The former has a factor 1 in the second term while the latter has a factor 2.
Otherwise, the two energy equations have similar structures.

Let us determine the pairing threshold for a positive $U$. To this end, set
$E_{s}=-8t-0$ in Eq. (152). Then $M_{00}$ logarithmically diverges. This
yields a critical coupling strength:

$|V^{s}_{\rm cr}|=\frac{2Ut}{U+8t}\>.$ (154)

This line separates the regions of “no pairs” and “s-pairs” in the 2D $UV$
phase diagram, see Fig. 9(b). Using the asymptotic behavior of the elliptic
integral

${\bf
K}\left(\frac{8t}{8t+\Delta}\right)\simeq\frac{1}{2}\log{\frac{64\,t}{\Delta}}\>;\hskip
14.22636pt\Delta\ll t\>,$ (155)

one obtains the binding energy near the threshold:

$\Delta=(64e^{-\pi})\,t\,\exp{\left(-\frac{2\pi t}{|V|-|V^{s}_{\rm
cr}|}\right)}\>.$ (156)

Note that the exponent is four times less than in the corresponding expression
in the attractive Hubbard model, Eq. (75). This is because the $UV$ model has
four attractive sites instead of one. In the $U\rightarrow\infty$ limit,
$|V^{s}_{\rm cr}|\rightarrow 2t$, and the general expression simplifies to

$\Delta(U=\infty)=64\,t\,\exp{\left(-\frac{\pi|V|}{|V|-2t}\right)}\>.$ (157)

In this form, the binding energy was given in Ref. [62].

Turning now to the $d$-symmetric state, Eq. (151), one observes that the
combination $M_{20}+M_{00}-2M_{11}$ converges in the limit $E\rightarrow-8t$.
Utilizing explicit expressions in B.3, one derives the $d$-pairing threshold

$|V^{d}_{\rm cr}|=\frac{2\pi t}{4-\pi}=(7.319584\ldots)\,t\>.$ (158)

It is independent of $U$, as expected for a $d$-symmetric wave function.
General expressions for $M_{00}$, $M_{11}$, and $M_{20}$ at arbitrary $E$ are
given in B.3. Using those, Eq. (151) defines the $d$-state energy as a
function of $V$.

###  5.2  Singlet states. Arbitrary momentum ${\bf P}$

Separation of the singlet dispersion relation, Eq. (148), into $s$ and $d$
sectors is also possible on the BZ diagonal. This is because at $P_{x}=P_{y}$,
the relations $M_{10}=M_{01}$ and $M_{20}=M_{02}$ continue to be valid.
Transformations described in Section 5.1 still apply, leading to the final
dispersion relations, Eqs. (152) and (151). The only difference is the
modified expression for $\alpha=4t\cos{(P_{x}/2)}$. This simple dependence on
$P$ along BZ diagonals can be used, for example, to extract the effective mass
separately for $s$\- and $d$-symmetrical pairs. Another consequence is a
simple modification of the binding thresholds:

$|V^{s}_{\rm
cr}(P_{x}=P_{y})|=\frac{2U\,t\cos{\frac{P_{x}}{2}}}{U+8\,t\cos{\frac{P_{x}}{2}}}\>.$
(159) $|V^{d}_{\rm
cr}(P_{x}=P_{y})|=\frac{2\pi\,t}{4-\pi}\,\cos{\frac{P_{x}}{2}}\>.$ (160)

Like in the 1D $UV$ model, see Eq. (125), these expressions indicate that the
pairing thresholds decrease at large lattice momenta. The energy of bound
states still grow with ${\bf P}$ but the lowest energy of two free particles
at the same ${\bf P}$ grows even faster. As a result, a bound pair may form at
a finite ${\bf P}$ even if it is unstable at ${\bf P}=0$. This physics is much
richer in 2D than in 1D. Below, we investigate it in some detail [87].

In order to determine the binding threshold at an arbitrary ${\bf P}$, energy
$E$ must be set equal to the minimal energy of two free particles
$E_{11}=-\alpha-\beta$. Upon substitution $E=E_{11}$, all $M_{nm}$ in Eq.
(148) diverge logarithmically. To regularize the determinant, express each
$M_{nm}$ as a sum of $M_{00}$ and remaining difference $L_{nm}$:

$M_{nm}=M_{00}+(M_{nm}-M_{00})\equiv M_{00}+L_{nm}\>.$ (161)

Note that all $L_{nm}$ converge to finite values in the $E=E_{11}$ limit.
Explicit analytical expressions are given in B.7. Insertion of Eq. (161) into
Eq. (148) and expansion of the determinant leads to a lengthy expression that
is a third-order polynomial in $M_{00}$. However, the $M^{3}_{00}$ and
$M^{2}_{00}$ terms cancel identically, and the determinant assumes the form

$A\cdot M_{00}+B=0\>,$ (162)

where

$\displaystyle A=$
$\displaystyle\;U|V|^{2}\left[(L_{20}-4L_{10})(L_{02}-4L_{01})-4(L_{10}+L_{01}-L_{11})^{2}\right]$
$\displaystyle+U|V|\left[(4L_{01}-L_{02})+(4L_{10}-L_{20})\right]+U+|V|^{2}\left[2(L_{02}+L_{20})-8L_{11}\right]-4|V|\>,$
(163)

and the specific form of $B$ is unimportant. In the limit $E\rightarrow
E_{11}$, both $A$ and $B$ remain finite whereas $M_{00}$ diverges. Thus, the
binding condition becomes

$A=0\>.$ (164)

Equations (163) and (164) determine the binding threshold sought. For a given
${\bf P}$, they define a function $|V_{\rm cr}|(U)$. Alternatively, for some
fixed $U$ and $|V|$, the threshold defines a line inside BZ that separates
unbound states at small ${\bf P}$ from bound pairs at large ${\bf P}$. An
example of such boundary lines is shown in Fig. 10. Properties of these lines
were studied in Ref. [87].

Figure 10: The pairing lines in the 2D square $UV$ model for $U=30\,t$ and
$|V|=1.1\,t$. Bound pairs form outside the respective lines. The $s$ and $d$
lines are solutions of Eqs. (163) and (164). The $p_{x}$ and $p_{y}$ lines are
solutions of Eqs. (167) and (168) at $E=-\alpha-\beta$. At these parameters,
both $p$ and $d$ lines terminate at the Brillouin zone boundaries. [87]

The existence of such lines poses an interesting question about pair formation
at a finite but low particle density. Figure 11 shows the same $s$ pairing
line as Fig. 10, together with the free-dispersion Fermi line corresponding to
a Fermi energy $E_{F}=-2.7\,t$ and filling factor 0.113. Consider states
within one of the four outlined circular segments, for example states with
$k_{x}=0$ and $k_{y}=k_{F}$ (the Fermi momentum). The combined momentum of two
such particles would be ${\bf P}=(0,2k_{F})$ which would land beyond the
$s$-pairing line. Therefore, the particles would tend to form a bound state.
Intriguingly, the states near the Fermi line close to BZ diagonals do not show
the same tendency because their combined momentum lies inside the pairing
surface, see Fig. 11. The entire Fermi line splits into eight disconnected
segments: four with pairing and four without. The segments without pairing
should produce sharp features in photoemission experiments (“Fermi arcs”),
while those with pairing will produce emission lines separated from the Fermi
energy by a gap. Such a behavior is consistent with photoemission signatures
of some cuprates superconductors [123].

Clearly, the central question is whether this empty-lattice picture retains
its qualitative features at small but finite densities. That requires a
reliable many-body method that can handle large systems such as Quantum Monte
Carlo [124, 125] or low-density $T$-matrix approximation [126, 127]. Such a
treatment is beyond the scope of the present work and is not analyzed here.

Figure 11: Comparison of the pairing and Fermi lines in the 2D square $UV$
model. The thick solid line is the same $s$ pairing line drawn in Fig. 10. The
hatched area in the center are free-particle states enclosed by a Fermi line
with energy $E_{F}=-2.7\,t$. It corresponds to a filling factor of 0.113 and
total particle density of $n=0.226$. The dashed line is the twice-Fermi
momentum line. The four circular segments outlined by thick lines enclose free
particle states that would tend to pair up if they were moving in an empty
lattice.

Next, we discuss full pair dispersion $E({\bf P})$ for all ${\bf P}$. Due to
complexity of the main dispersion relation, Eq. (148), very little can be done
analytically beyond the separation of $s$ and $d$ energies on BZ diagonals.
Some simplification occurs at BZ edges. Let us set for definitiveness,
$P_{x}=\pm\pi$. Then $M_{10}=M_{20}=M_{11}=0$ by symmetry. The remaining
matrix elements, $M_{00}$, $M_{01}$, and $M_{02}$, become one-dimensional
integrals given by Eq. (119). Upon expansion, the determinant splits in two
factors. One has an explicit solution

$E_{1}(\pm\pi,P_{y})=-\sqrt{|V|^{2}+\beta^{2}}\>,$ (165)

where $\beta=4t\cos{(P_{y}/2)}$. The second factor can be brought to the
following form

$2|V|(U+|E|)(|E|-\sqrt{|E|^{2}-\beta^{2}})-\beta^{2}(U+\sqrt{|E|^{2}-\beta^{2}})=0\>.$
(166)

This equation does not have a simple-form analytical solution. An example of
numerical solution of Eq. (148) is shown in Fig. 12. A bound pair behaves as a
single particle with a fairly complex dispersion. Notice degeneracies along
some high-symmetry lines.

Figure 12: Dispersion of bound pairs in the 2D square $UV$ model for $U=30\,t$
and $V=-14\,t$. The $s$ and $d$ energies are solutions of Eq. (148). The $s$
dispersion along the ${\rm X}-{\rm M}$ line is given by Eq. (165). $p_{x}$ and
$p_{y}$ energies are solutions of Eqs. (167) and (168), respectively. Notice
$p_{x}$, $p_{y}$ degeneracy along the ${\rm M}-\Gamma$ line and $s$, $p_{x}$
degeneracy along the ${\rm X}-{\rm M}$ line.

###  5.3  Triplet states

There are two antisymmetrized vectors $\\{{\bf b}_{-}\\}=\\{(1,0),$ $(0,1)\\}$
and two functions $\Phi^{-}$. Upon constructing the system, Eq. (43), the
integration variable in $M^{-}_{nm}$ is shifted as ${\bf q}^{\prime}={\bf
q}-\frac{\bf P}{2}$. Off-diagonal terms vanish by symmetry, and the $(2\times
2)$ system splits into two separate equations:

$\displaystyle\left\\{1-|V|\left(M_{00}-M_{20}\right)\right\\}\Phi^{-}_{p_{x}}$
$\displaystyle=0\>,$ (167)
$\displaystyle\left\\{1-|V|\left(M_{00}-M_{02}\right)\right\\}\Phi^{-}_{p_{y}}$
$\displaystyle=0\>.$ (168)

Note that such a decomposition takes place over the entire BZ. Since
$M_{20}\neq M_{02}$, $p_{x}$ and $p_{y}$ energies are different, as shown Fig.
12. Along the diagonals, $M_{20}=M_{02}$, and the dispersion becomes double-
degenerate.

Let us derive the binding condition along BZ diagonals. In the limit
$E\rightarrow E_{11}$, the difference $M_{00}-M_{20}$ converges. Making use of
Eq. (306), one obtains:

$|V^{p}_{\rm
cr}|=\frac{2\pi}{\pi-2}\,t\cos{\frac{P_{x}}{2}}=(5.503876\ldots)\,t\cos{\frac{P_{x}}{2}}\>.$
(169)

The two-particle phase diagram of the square $UV$ model for ${\bf P}=0$ is
shown in Fig. 9(b). The binding condition at arbitrary ${\bf P}$ can be
obtained from Eqs. (167) and (168) if analytical expressions for the $M$
differences, Eqs. (337) and (338), are utilized. An example of $p_{x}$ and
$p_{y}$ pairing lines is given in Fig. 10. Analytical properties of the $p$
pairing lines were studied in Ref. [87].

Consider now triplet energy at arbitrary ${\bf P}$. Again, certain
simplification takes place at BZ edges. At $P_{x}=\pm\pi$, $M_{20}=0$, and the
energy of $p_{x}$ pair is given by Eq. (165). Thus, the dispersion at BZ edges
is always double-degenerate. This can also be seen in Fig. 12. At the
$\Gamma$-point, the $p$ energy can be determined from the following equation

$\pi-2{\bf E}\\!\left(\frac{8t}{|E_{p}|}\right)=\frac{2\pi
t}{|V|}\frac{8t}{|E_{p}|}\>,$ (170)

where Eqs. (295) and (298) have been utilized. Equation (170) does not have a
simple analytical solution for $E_{p}$.

###  5.4  Longer-range attractions

The square $UV$ model is the rare case when the two-body problem was solved
for interactions beyond the nearest neighbors [99, 37, 128]. In Refs. [37,
128], next-nearest interaction and next-nearest hopping were considered. In
Ref. [99], attraction was extended up to the seventh nearest neighbors (with
the potential depth being constant within the radius of attraction $R$) but
the hopping was limited to the first neighbors only. Below, we provide the
reasoning behind the model analyzed in Ref. [99].

In order to form a bound pair, attraction must exceed a threshold to overcome
the strong repulsive core. For nearest-neighbor attraction, the threshold is
about $2t$, see Eq. (154). In physical units, it may be quite a large number.
Let us assume the effective hopping of holes in the cuprates to be $t_{h}\sim
0.1$ eV [129, 130]. Then, the attraction must be of order $V\simeq 0.1-0.2$
eV, which is arguably quite large. At the same time, the cuprates are
anisotropic polar solids with low carrier density and poorly screened
electron-phonon interactions. This favors formation of bipolarons [57, 58,
131, 132, 133]. It also results in longer-range shallow attractive potentials
within the copper–oxygen planes. Spreading attraction to many sites allows
lowering the threshold on individual sites. Since it is the “total power” of
the potential that matters for binding, one expects the threshold to scale
approximately inversely to the number of sites participating in attraction,
i.e., inversely with the potential area. (In continuous quantum mechanics the
scaling is also $V_{\rm cr}\propto 1/R^{2}$.) This is what was confirmed by
exact calculations [99]. The dependence of $|V_{cr}|(R)$ is shown in Fig. 13.
It approximately follows the $1/R^{2}$ scaling, as expected for a continuum
problem, with fluctuations around the line which reflects the discrete nature
of the lattice. One can see that by the 6th or 7th nearest neighbors, the
threshold falls by an order of magnitude to about $0.2\,t\simeq 0.02$ eV,
which would be easier to attain in real solids.

Figure 13: Binding threshold of the $UV$ model on the square lattice versus
the distance of $n$th nearest neighbors from the origin. $U=\infty$. Solid
symbols are exact values obtained from solving two-body Schrödinger equations.
(The values are 2.0, 0.90212, 0.57139, 0.31280, 0.25075, 0.20967, and
0.15584.) The dashed line is a phenomenological dependence $|V_{\rm
cr}|=2/R^{2}$ and is guide to the eye. Adapted from Ref. [99].

One can add to this argument the light pairs mechanism already discussed in
Section 4.2. In shallow long-range attractive potentials, two or more
neighbors will have equal or close attractive strengths with high probability.
That will enable resonant movement of pairs without breaking attractive bonds.
As a result, pair effective mass will remain of order $m_{0}$ even in the
strong coupling limit. We illustrate this mechanism for an attractive
potential extended to the second nearest neighbors with $V_{2}=V_{1}$. We
introduce the strong-coupling singlet dimer basis [128]

$D_{i,{\bf m}}=\frac{1}{\sqrt{2}}\left(|\uparrow\rangle_{\bf
m}|\downarrow\rangle_{{\bf m}+{\bf b}_{i}}+|\downarrow\rangle_{\bf
m}|\uparrow\rangle_{{\bf m}+{\bf b}_{i}}\right)\>,$ (171)

with ${\bf b}_{1}=(1,0)$, ${\bf b}_{2}=(1,1)$, ${\bf b}_{3}=(0,1)$, and ${\bf
b}_{4}=(-1,1)$. Hamiltonian action within this space is given by

$\displaystyle\hat{H}D_{1,{\bf m}}$ $\displaystyle=-t\left(D_{2,{\bf
m}}\\!+\\!D_{2,{\bf m}-{\bf b}_{3}}\\!+\\!D_{4,{\bf m}+{\bf
b}_{1}}\\!+\\!D_{4,{\bf m}-{\bf b}_{4}}\right),\makebox[5.69046pt]{}$
$\displaystyle\hat{H}D_{2,{\bf m}}$ $\displaystyle=-t\left(D_{1,{\bf
m}}\\!+\\!D_{1,{\bf m}+{\bf b}_{3}}\\!+\\!D_{3,{\bf m}}\\!+\\!D_{3,{\bf
m}+{\bf b}_{1}}\right),$ $\displaystyle\hat{H}D_{3,{\bf m}}$
$\displaystyle=-t\left(D_{4,{\bf m}}\\!+\\!D_{4,{\bf m}+{\bf
b}_{1}}\\!+\\!D_{2,{\bf m}}\\!+\\!D_{2,{\bf m}-{\bf b}_{1}}\right),$
$\displaystyle\hat{H}D_{4,{\bf m}}$ $\displaystyle=-t\left(D_{3,{\bf
m}}\\!+\\!D_{3,{\bf m}-{\bf b}_{1}}\\!+\\!D_{1,{\bf m}-{\bf
b}_{1}}\\!+\\!D_{1,{\bf m}+{\bf b}_{4}}\right).$ (172)

Converting the corresponding Schrödinger equation to momentum space and
expanding the determinant, one obtains dimer dispersion:

$\displaystyle\tilde{E}_{1,2}({\bf P})$
$\displaystyle=\pm(2t)\sqrt{2+\cos{(P_{x}a)}+\cos{(P_{y}a)}}\>,$ (173)
$\displaystyle\tilde{E}_{3,4}({\bf P})$ $\displaystyle=0\>,$ (174)

where $\tilde{E}$ is referenced from $-|V|$. At small ${\bf P}$,
$\tilde{E}_{1}\approx-4t+(ta^{2}P^{2})/4$, from where the pair mass is

$m^{\ast}_{p}=\frac{2\hbar^{2}}{ta^{2}}=4\,m_{0}\>,$ (175)

where $m_{0}=\hbar^{2}/(2ta^{2})$. Thus, in this model the mass enhancement is
limited by the same factor of 4 as in the 1D model with long-range attraction,
cf. Eq. (139). The dimer analysis can also be extended to nonzero second-
neighbor hopping [128].

Next, we rigorously solve the $UV$ model with attraction extended to second
nearest neighbors. In deviation from [99], the second neighbor attraction
$-|V_{2}|$ will be different from the first neighbor attraction $-|V_{1}|$. We
consider only singlet pairs and begin with ${\bf P}=0$. According to the
general theory, there are five symmetrized neighbor vectors. They can be
chosen, for example, as $\\{{\bf
b}_{+}\\}=\\{(0,0),(1,0),(0,1),(1,-1),(1,1)\\}$. Hence, there are five
equations, Eqs. (27) and (28), for five functions $\Phi^{+}$. The equations
mix $s$\- and $d$-symmetric pair states. To untangle them, add and subtract
the equations for $\Phi^{+}_{(1,0)}$ and $\Phi^{+}_{(0,1)}$, and then do the
same for the equations for $\Phi^{+}_{(1,-1)}$ and $\Phi^{+}_{(1,1)}$. The two
differences are

$\displaystyle\left\\{1-|V_{1}|(M_{20}+M_{00}-2M_{11})\right\\}\left(\Phi^{+}_{10}-\Phi^{+}_{01}\right)=0\,,$
(176)
$\displaystyle\left\\{1-|V_{2}|(M_{22}+M_{00}-2M_{20})\right\\}\left(\Phi^{+}_{1-\\!1}-\Phi^{+}_{11}\right)=0\,.$
(177)

The first equation is equivalent to Eq. (151) and describes a $d$-pair
previously discussed. It has lobes along $x$ and $y$ axes, and therefore may
be called a $d_{x^{2}-y^{2}}$ state. This state does not depend on $V_{2}$ and
has the same pairing threshold as before, Eq. (158). Equation (177) describes
another $d$ state with lobes along the square diagonals. It may be called a
$d_{xy}$ state. Note that the two $d$ states do not mix. The $d_{xy}$ state
does not involve potential $V_{1}$ because the wave function has nodes on the
first nearest neighbors. The linear combination of $M$’s in Eq. (177)
converges in the $E\rightarrow E_{11}$ limit. Applying Eqs. (306) and (312),
one derives the $d_{xy}$ pairing threshold

$|V^{d_{xy}}_{2,{\rm cr}}|=\frac{3\pi}{3\pi-8}\,t=(6.614909\ldots)\,t\>.$
(178)

It is smaller than the $d_{x^{2}-y^{2}}$ threshold, Eq. (158), by $\approx
10$%. Because of a node at ${\bf m}=(0,0)$, the pair wave function is larger
at the second nearest neighbors than at the first nearest neighbors. As a
result, potential $V_{2}$ has more “power” and produces a bound state at a
slightly lower value that $V_{1}$.

Returning to the general system, Eqs. (27) and (28), the two equation sums
combine with the equation for $\Phi^{+}_{(0,0)}$ and form the $s$-sector. The
full system reads

$\left[\begin{array}[]{ccc}1+UM_{00}&-4|V_{1}|M_{10}&-4|V_{2}|M_{11}\\\
UM_{10}&1-|V_{1}|(M_{20}+M_{00}+2M_{11})&-2|V_{2}|(M_{10}+M_{21})\\\
UM_{11}&-2|V_{1}|(M_{10}+M_{21})&1-|V_{2}|(M_{22}+M_{00}+2M_{20})\end{array}\right]\left[\begin{array}[]{c}\Phi^{+}_{00}\\\
\frac{1}{2}\Phi^{+}_{10}+\frac{1}{2}\Phi^{+}_{01}\\\
\frac{1}{2}\Phi^{+}_{1-1}+\frac{1}{2}\Phi^{+}_{11}\end{array}\right]=\left[\begin{array}[]{c}0\\\
0\\\ 0\end{array}\right].$ (179)

The binding threshold is derived by equating the determinant to zero and
setting $E=-8t$. Because all $M_{mn}$ logarithmically diverge, one needs to
apply the subtractive procedure, Eq. (161). After the substitution, the
determinant is expanded into the form $A\cdot M_{00}+B$, with complex
expressions for $A$ and $B$ that are both convergent. Again, terms with
$M^{2}_{00}$ and $M^{3}_{00}$ cancel to zero. Since only $M_{00}$ remains
divergent, the binding condition is simply $A=0$. Utilizing Eqs. (305)-(315),
the final result reads:

$(32-9\pi)U|V_{1}||V_{2}|-6\pi
U|V_{1}|t+(12\pi-64)U|V_{2}|t+(256-72\pi)|V_{1}||V_{2}|t+12\pi
Ut^{2}-48\pi|V_{1}|t^{2}-48\pi|V_{2}|t^{2}=0.$ (180)

For $V_{2}=0$, it reduces to Eq. (154). For $V_{1}=0$, Eq. (180) yields

$|V_{2,{\rm cr}}|=\frac{3\pi Ut}{(16-3\pi)U+12t}\>.$ (181)

For $V_{1}=V_{2}$ and $U\rightarrow\infty$, Eq. (180) reduces to

$(32-9\pi)|V|^{2}+(6\pi-64)|V|+12\pi=0\>,$ (182)

that was reported in Ref. [99]. The smallest root of this equation is $V_{\rm
cr}(2)=0.902120$, which is plotted in Fig. 13. $V_{\rm cr}(2)$ should be
compared with the $U\rightarrow\infty$ limit of Eq. (154), which is $V_{\rm
cr}(1)=2.0$.

Figure 14: Pair effective mass in the 2D square $UV_{1}V_{2}$ model with
attraction on the first and second nearest neighbors. The masses are obtained
by numerically solving Eq. (179) for $E$ and then applying Eq. (183).
$U=30\,t$. The dashed line marks the light-pair limit, Eq. (175).

Decomposition of the full dispersion relation into $s$ and $d$ sectors also
occurs on BZ diagonals, ${\bf P}=(P,P)$. Matrix elements $M_{nm}$ are still
given by their isotropic expressions listed in B.3 but with
$\alpha=4t\cos{(P/2)}$. This allows for numerical calculation of effective
masses using the formula:

$\frac{m^{\ast}_{p}}{m_{0}}=4t\left[\frac{\partial^{2}E(P)}{\partial
P^{2}}\right]^{-1}.$ (183)

Note that Eq. (179) contains more than one $s$ state [37], but we consider
only the lowest one. Figure 14 shows the effective mass for three different
cases: (i) Attraction on the first neighbors only; (ii) Attraction on the
second neighbors only; (iii) Attraction of equal strength on both the first
and the second neighbors, $V_{1}=V_{2}$. One can see that in the first two
cases the mass grows linearly with $V$ while in the third case the mass is
bounded by the light pair limit given by Eq. (175). We also note that for
$|V|<10\,t$, the mass enhancement is very modest, $m^{\ast}_{p}/m_{0}<4$, in
all three cases.

##  6  2D. $UV$ model on the rectangular lattice

In the rectangular $UV$ model, both hopping and interaction along $x$ and $y$
axes are different, as illustrated in Fig. 15. Such a model can potentially be
realized in cold gases. The model is rich in physical content and smoothly
interpolates between 1D and 2D square models investigated earlier. In this
section, we mostly focus on deriving binding conditions. In difference from
preceding sections, we consider the possibility of repulsive nearest-neighbor
interaction $V$. Therefore, the equations are written in terms of $V$ that can
be of any sign rather than $-|V|$.

###  6.1  Singlet states. $\Gamma$-point

The dispersion relation is given by Eq. (148) in which $-|V|$ is replaced with
$V_{x}$ and $V_{y}$ in the second and third columns, respectively:

$\left|\begin{array}[]{ccc}1+UM_{00}&2V_{x}M_{10}&2V_{y}M_{01}\\\
UM_{10}&1+V_{x}(M_{20}+M_{00})&2V_{y}M_{11}\\\
UM_{01}&2V_{x}M_{11}&1+V_{y}(M_{02}+M_{00})\end{array}\right|=0\>.$ (184)

Matrix elements $M_{nm}$ are still defined by Eq. (149) but with
$\alpha=4t_{x}\cos{(P_{x}/2)}$ and $\beta=4t_{y}\cos{(P_{y}/2)}$. In the
following, we only consider the ground state. Hence, $\alpha=4t_{x}$ and
$\beta=4t_{y}$. At the threshold, $E\rightarrow-\alpha-\beta$, and all
$M_{nm}$ logarithmically diverge. To obtain a finite result, we utilize the
subtractive procedure defined by Eq. (161). Substitution in Eq. (184) and
expansion of the determinant results in a third-order polynomial in $M_{00}$.
The $M^{3}_{00}$ and $M^{2}_{00}$ terms vanish identically, and the
determinant assumes the form of Eq. (162). From here, the binding condition is
$A=0$, or, in full form

$\displaystyle
UV_{x}V_{y}\left[(L_{20}-4L_{10})(L_{02}-4L_{01})-4(L_{10}+L_{01}-L_{11})^{2}\right]$
$\displaystyle+\>UV_{x}\left[L_{20}-4L_{10}\right]+UV_{y}\left[L_{02}-4L_{01}\right]+\>V_{x}V_{y}\left[2(L_{02}+L_{20})-8L_{11}\right]+\>U+2V_{x}+2V_{y}=0\>.$
(185)

Equation (185) is the general binding condition in the rectangular $UV$ model.
It reduces to Eq. (163) if $V_{x}=V_{y}=-|V|$.

Figure 15: The rectangular $UV$ model.

Let us analyze particular cases of Eq. (185). First, we investigate isotropic
hopping, $t_{y}=t_{x}=t$. Utilizing Eqs. (335)-(339), the binding condition
becomes

$\frac{4-\pi}{4\pi t^{2}}\,UV_{x}V_{y}+\frac{1}{\pi
t}\,U(V_{x}+V_{y})+\frac{2(4-\pi)}{\pi t}\,V_{x}V_{y}+U+2(V_{x}+V_{y})=0\>.$
(186)

For isotropic attraction, $V_{x}=V_{y}=-|V|$, Eq. (186) reduces to the product
of $s$\- and $d$\- thresholds, given by Eqs. (154) and (158), respectively. To
get a sense of the effect of strong nearest neighbor repulsion, we set
$V_{x}=U$ and then $U\rightarrow\infty$. Obviously, to form a bound state,
$V_{y}$ must be attractive. Equation (186) yields the threshold:

$|V_{y,{\rm cr}}|=\frac{4t}{4-\pi}=(4.569792\ldots)\,t\>.$ (187)

This should be compared with the $U\to\infty$ limit of the isotropic, four-
site attraction, $|V^{s}_{\rm cr}|=2t$, see Eq. (154). Making two of the
attractive sites strongly repulsive increases the necessary strength of the
two remaining attractions from $2t$ to $4.57\,t$.

Next, we consider the case of anisotropic attraction, $V_{x}=0$, $V_{y}<0$,
and arbitrary hopping anisotropy. Equation (186) reduces to

$|V_{y,{\rm cr}}|=\frac{U}{2+U(L_{02}-4L_{01})}\>.$ (188)

According to B.7,

$L_{02}-4L_{01}=\frac{1}{\pi
t_{y}}\left\\{\frac{t_{y}-t_{x}}{t_{y}}\arcsin{\sqrt{\frac{t_{y}}{t_{x}+t_{y}}}}+\sqrt{\frac{t_{x}}{t_{y}}}\right\\}.$
(189)

Several particular cases are of interest.

(i) $t_{x}=0$. In this case, the system splits into individual $y$ chains,
each of which hosts a 1D $UV$ model. $L_{02}-4L_{01}=(2t_{y})^{-1}$ and
$|V_{y,{\rm cr}}|=2Ut_{y}/(U+4t_{y})$, which coincides with the earlier
result, Eq. (125).

(ii) $t_{x}=t_{y}=t$. This is the case of isotropic hopping.
$L_{02}-4L_{01}=(\pi t)^{-1}$ and

$|V_{y,{\rm cr}}|=\frac{\pi Ut}{U+2\pi\,t}\>.$ (190)

This expression should be compared with its isotropic attraction counterpart,
Eq. (154). The attraction strength must be larger because there are only two
attractive sites rather than four. In the $U=\infty$ limit, the threshold
value is $\pi t$ rather than $2t$.

(iii) $t_{y}=0$. In this case, the system splits into individual $x$ chains.
The two particles reside on adjacent chains and feel an attraction $V_{y}$
when their $x$ coordinates coincide. The model is isomorphic to the 1D
attractive Hubbard model. In the $t_{y}\rightarrow 0$ limit, $L_{02}-4L_{01}$
diverges as $\simeq 4(3\pi)^{-1}(t_{x}t_{y})^{-1/2}$. Hence, $|V_{y,{\rm
cr}}|\rightarrow 0$, as expected for the 1D attractive Hubbard model.

###  6.2  Triplet states. $\Gamma$-point

Dispersions of $p_{x}$ and $p_{y}$ pairs are given by Eqs. (167) and (168)
with $V$ replaced by $V_{x}$ and $V_{y}$, respectively. The binding thresholds
at ${\bf P}=(0,0)$ are

$\displaystyle|V_{x,{\rm cr}}|$ $\displaystyle=\frac{\pi
t_{x}}{\frac{t_{x}+t_{y}}{t_{x}}\arcsin{\sqrt{\frac{t_{x}}{t_{x}+t_{y}}}}-\sqrt{\frac{t_{y}}{t_{x}}}}\>,$
(191) $\displaystyle|V_{y,{\rm cr}}|$ $\displaystyle=\frac{\pi
t_{y}}{\frac{t_{x}+t_{y}}{t_{y}}\arcsin{\sqrt{\frac{t_{y}}{t_{x}+t_{y}}}}-\sqrt{\frac{t_{x}}{t_{y}}}}\>.$
(192)

Consider the $p_{y}$ pair. In the $t_{x}=0$ limit, $|V_{y,{\rm cr}}|=2t_{y}$,
which is the correct threshold for the 1D $UV$ model, see Eq. (147). In the
isotropic hopping case, $t_{x}=t_{y}=t$, $|V_{y,{\rm cr}}|=2\pi t/(\pi-2)$,
which coincides with Eq. (169). Finally, in the $t_{y}\rightarrow 0$ limit,
$|V_{y,{\rm cr}}|=\frac{3\pi}{2}\sqrt{t_{x}t_{y}}\rightarrow 0$. The $p_{y}$
pair forms with a zero threshold.

##  7  2D. $UV$ model on the triangular lattice

The triangular $UV$ model possesses another qualitative feature: light pairs
are formed due to lattice topology rather than degeneracy of the attractive
potential. Consider Fig. 16. In a tightly bound pair, the members can move
through the lattice while remaining nearest neighbors to each other and
without breaking the main attractive bond. As a result, the pair remains light
even when attraction is limited to first nearest neighbors. Analysis of dimer
motion [134, 135, 136] predicts that the mass of the lowest singlet pair is
$m^{\ast}_{p}<6\,m_{0}$.

This situation is more general than it might seem. There exist other lattices
that can support light pairs with first neighbor attraction. One example is
the staggered ladder [57, 134], where the two ladder chains are shifted
relative to each other by half of lattice constant. Similarly, two square
lattices stuck one on top of each other and shifted by half diagonal of
elementary plaquette will support “crab motion” as long as the two constituent
particles reside on different layers. This system will be considered in
Section 10.1. In 3D, the face-centered cubic lattice supports light pairs
[90]. Furthermore, if the range of primary hopping is extended, then more
lattices are added to the list. For example, the square lattice with next
nearest neighbor hopping (across the elementary plaquette) produces crab
motion and light pairs as a result [57, 58, 133].

Figure 16: In the triangular $UV$ model, the pair moves in first order in
hopping even when the attraction is limited to first nearest neighbors.

The light pair effect has implications for the formation of phonon bipolarons.
It is by now well understood [137] that long-range electron–ion interactions
exponentially reduce the effective mass of polarons, and consequently of
polaron pairs, bipolarons. On triangular-like lattices, the bipolarons acquire
additional lightness due to crab motion and become “superlight” [133, 134,
138, 101]. Based on that and on the BEC formula, Eq. (1), it was suggested
that triangular-like lattices could provide even higher-$T_{c}$
superconductivity than square-like lattices [138].

The $UV$ model on triangular lattice was solved by Bak [88] using a method
similar to our unsymmetrized solution, and then in Ref [101]. Here we analyze
the model with emphasis on the pair mass, dispersion, and binding conditions.

###  7.1  General dispersion relations

In the singlet sector, there are four symmetrized vectors: $\\{{\bf
b}_{+}\\}=\\{(0,0),(1,0),(\frac{1}{2},\frac{\sqrt{3}}{2}),(-\frac{1}{2},\frac{\sqrt{3}}{2})\\}\equiv\\{{\bf
0},{\bf 1},{\bf 2},{\bf 3}\\}$. Changing momentum variable, ${\bf
q}^{\prime}={\bf q}-\frac{\bf P}{2}$, and redefining amplitudes as
$\tilde{\Phi}^{+}_{\bf b}=e^{-i({\bf P}/2){\bf b}}\Phi^{+}_{\bf b}$, the
master equations, Eqs. (27) and (28), take the form

$\left[\begin{array}[]{cccc}1+UM^{+}_{\bf 00}&-|V|M^{+}_{\bf
01}&-|V|M^{+}_{\bf 02}&-|V|M^{+}_{\bf 03}\\\ UM^{+}_{\bf 10}&1-|V|M^{+}_{\bf
11}&-|V|M^{+}_{\bf 12}&-|V|M^{+}_{\bf 13}\\\ UM^{+}_{\bf 20}&-|V|M^{+}_{\bf
21}&1-|V|M^{+}_{\bf 22}&-|V|M^{+}_{\bf 23}\\\ UM^{+}_{\bf 30}&-|V|M^{+}_{\bf
31}&-|V|M^{+}_{\bf 32}&1-|V|M^{+}_{\bf
33}\end{array}\right]\left[\begin{array}[]{c}\tilde{\Phi}^{+}_{0}\\\
\tilde{\Phi}^{+}_{1}\\\ \tilde{\Phi}^{+}_{2}\\\
\tilde{\Phi}^{+}_{3}\end{array}\right]=0\>.$ (193)

Here

$\displaystyle M^{+}_{\bf 00}$ $\displaystyle=\frac{1}{N}\sum_{\bf
q}\frac{1}{W}\>,$ (194) $\displaystyle M^{+}_{{\bf 0}{\bf b}^{\prime}_{+}}$
$\displaystyle=\frac{1}{N}\sum_{\bf q}\frac{2\cos{({\bf q}{\bf
b}^{\prime}_{+})}}{W}\>,$ (195) $\displaystyle M^{+}_{{\bf b}_{+}{\bf 0}}$
$\displaystyle=\frac{1}{N}\sum_{\bf q}\frac{e^{i{\bf q}{\bf b}_{+}}}{W}\>,$
(196) $\displaystyle M^{+}_{{\bf b}_{+}{\bf b}^{\prime}_{+}}$
$\displaystyle=\frac{1}{N}\sum_{\bf q}\frac{2\cos{({\bf q}{\bf
b}^{\prime}_{+})}\,e^{i{\bf q}{\bf b}_{+}}}{W}\>,$ (197)
$W=|E|-\alpha\cos{q_{x}}-\beta\cos{\frac{q_{x}}{2}}\cos{\frac{\sqrt{3}q_{y}}{2}}-\gamma\sin{\frac{q_{x}}{2}}\sin{\frac{\sqrt{3}q_{y}}{2}}\>,$
(198)

and $\alpha$, $\beta$, and $\gamma$ are defined in Eqs. (345)-(347).
Integration over BZ can be replaced with integration over the rectangle $0\leq
q_{x}\leq 2\pi$, $0\leq q_{y}\leq(4\pi)/\sqrt{3}$. The expressions for $M^{+}$
are given in C.2. The following relations hold for any ${\bf P}$:

$\displaystyle M^{+}_{\bf 01}$ $\displaystyle=2M^{+}_{\bf 10}\>,\hskip
2.84544ptM^{+}_{\bf 02}=2M^{+}_{\bf 20}\>,\hskip 2.84544ptM^{+}_{\bf
03}=2M^{+}_{\bf 30}\>,$ (199) $\displaystyle M^{+}_{\bf 12}$
$\displaystyle=M^{+}_{\bf 21}\>,\hskip 2.84544ptM^{+}_{\bf 13}=M^{+}_{\bf
31}\>,\hskip 2.84544ptM^{+}_{\bf 23}=M^{+}_{\bf 32}\>.$ (200)

Additional relations hold along BZ symmetry lines.

Figure 17: Dispersion of bound pairs in the triangular $UV$ model for
$U=30\,t$ and $V=-16\,t$. $E=-12\,t$ is the lowest energy of two free
particles at the $\Gamma$-point.

In the triplet sector, there are three anti-symmetrized vectors: $\\{{\bf
b}_{-}\\}=\\{{\bf 1},{\bf 2},{\bf 3}\\}$. Performing the same variable change,
${\bf q}^{\prime}={\bf q}-\frac{\bf P}{2}$, and amplitude redefinition
$\tilde{\Phi}^{-}_{{\bf b}_{-}}=e^{-i({\bf P}/2){\bf b}_{-}}\Phi^{-}_{{\bf
b}_{-}}$, the master equations yield

$\left[\begin{array}[]{ccc}1-|V|M^{-}_{\bf 11}&-|V|M^{-}_{\bf
12}&-|V|M^{-}_{\bf 13}\\\ -|V|M^{-}_{\bf 21}&1-|V|M^{-}_{\bf
22}&-|V|M^{-}_{\bf 23}\\\ -|V|M^{-}_{\bf 31}&-|V|M^{-}_{\bf
32}&1-|V|M^{-}_{\bf
33}\end{array}\right]\left[\begin{array}[]{c}\tilde{\Phi}^{-}_{1}\\\
\tilde{\Phi}^{-}_{2}\\\ \tilde{\Phi}^{-}_{3}\end{array}\right]=0\>,$ (201)

where

$M^{-}_{{\bf b}_{-}{\bf b}^{\prime}_{-}}=\frac{1}{N}\sum_{\bf
q}\frac{(-2i)\sin{({\bf q}{\bf b}^{\prime}_{-})}\,e^{i{\bf q}{\bf
b}_{-}}}{W}\>.$ (202)

The expressions for $M^{-}$ are also given in C.2, Eqs. (380)-(388). One has:

$M^{-}_{\bf 12}=M^{-}_{\bf 21}\>,\hskip 8.5359ptM^{-}_{\bf 13}=M^{-}_{\bf
31}\>,\hskip 8.5359ptM^{-}_{\bf 23}=M^{-}_{\bf 32}\>.$ (203)

Thus, the triplet dispersion is described by a symmetric matrix for all ${\bf
P}$. A typical pair dispersion is shown in Fig. 17. Of note is the fact that
attraction $V=-16\,t$ is barely above the pairing threshold for the highest,
$f$, state given by Eq. (217). Thus, all six pair states are well defined in
the entire BZ. At weaker attractions, the bound states begin to disappear into
the two-particle continuum one-by-one.

The effective mass of the lowest, $s$-symmetric pair is shown in Fig. 18. One
can observe that the mass indeed approaches the strong-coupling limit [134,
135, 136], $m^{\ast}_{p}/m_{0}=6$, for all $U$. However, the approach is slow.
For most realistic attractions, $|V|<20\,t$, the pair mass is no heavier than
just 4 free-particle masses.

###  7.2  Pairing thresholds at the $\Gamma$-point

Determination of pairing thresholds from the general dispersion relations,
Eqs. (193) and (201), at arbitrary ${\bf P}$ can be done only numerically. At
the $\Gamma$-point, the systems acquire additional symmetries and the
thresholds can be derived analytically. Analysis is greatly aided by the group
theory. All the necessary information is given in C.4.

Figure 18: Effective mass of the lowest, $s$-symmetric bound pair in the
triangular $UV$ model. The strong coupling limit is $m^{\ast}_{p}/m_{0}=6$.

For the singlet states, we combine $\Phi_{0}$, which is unchanged, with
symmetrized combinations of basis functions given by Eqs. (402) and (405) into
a new basis:

$\left[\\!\begin{array}[]{c}\Phi_{0}\\\ \Phi_{s}\\\ \Phi_{d_{1}}\\\
\Phi_{d_{2}}\end{array}\\!\right]=\left[\begin{array}[]{rrrr}1&0&0&0\\\
0&1&1&1\\\ 0&2&-1&-1\\\
0&0&1&-1\end{array}\right]\\!\left[\\!\begin{array}[]{c}\Phi_{0}\\\
\Phi^{+}_{1}\\\ \Phi^{+}_{2}\\\
\Phi^{+}_{3}\end{array}\\!\right]\equiv\hat{A}_{S}\left[\\!\begin{array}[]{c}\Phi_{0}\\\
\Phi^{+}_{1}\\\ \Phi^{+}_{2}\\\ \Phi^{+}_{3}\end{array}\\!\right].$ (204)

In terms of the new basis, dispersion equation, Eq. (193), transforms into a
block-diagonal form:

$\hat{A}_{S}\cdot[\ldots]\cdot\hat{A}^{-1}_{S}=\left[\\!\begin{array}[]{cccc}1+UM^{+}_{\bf
00}&-2|V|M^{+}_{\bf 10}&0&0\\\ 3UM^{+}_{\bf 10}&1-|V|(M^{+}_{\bf
11}+2M^{+}_{\bf 12})&0&0\\\ 0&0&1-|V|(M^{+}_{\bf 11}-M^{+}_{\bf 12})&0\\\
0&0&0&1-|V|(M^{+}_{\bf 11}-M^{+}_{\bf
12})\end{array}\right]\left[\begin{array}[]{c}\Phi_{0}\\\ \Phi_{s}\\\
\Phi_{d_{1}}\\\ \Phi_{d_{2}}\end{array}\right]=0\>.$ (205)