# Bose-Hubbard models with on-site and nearest-neighbor interactions: Exactly
solvable case

Saidakhmat Lakaev, Shokhrukh Kholmatov, Shakhobiddin Khamidov  Samarkand State
University, University boulevard, 15, 140104 Samarkand, Uzbekistan
<EMAIL_ADDRESS>University of Vienna, Oskar-Morgenstern-Platz 1, 1090 Vienna,
Austria<EMAIL_ADDRESS>Samarkand State University,
University boulevard, 15, 140104 Samarkand, Uzbekistan
<EMAIL_ADDRESS>
###### Abstract.

We study the discrete spectrum of the two-particle Schrödinger operator
$\widehat{H}_{\mu\lambda}(K),$ $K\in\mathbb{T}^{2},$ associated to the Bose-
Hubbard Hamiltonian $\widehat{\mathbb{H}}_{\mu\lambda}$ of a system of two
identical bosons interacting on site and nearest-neighbor sites in the two
dimensional lattice $\mathbb{Z}^{2}$ with interaction magnitudes
$\mu\in\mathbb{R}$ and $\lambda\in\mathbb{R},$ respectively. We completely
describe the spectrum of $\widehat{H}_{\mu\lambda}(0)$ and establish the
optimal lower bound for the number of eigenvalues of
$\widehat{H}_{\mu\lambda}(K)$ outside its essential spectrum for all values of
$K\in\mathbb{T}^{2}.$ Namely, we partition the $(\mu,\lambda)$-plane such that
in each connected component of the partition the number of bound states of
$\widehat{H}_{\mu\lambda}(K)$ below or above its essential spectrum cannot be
less than the corresponding number of bound states of
$\widehat{H}_{\mu\lambda}(0)$ below or above its essential spectrum.

###### Key words and phrases:

Two-particle system, discrete Schrödinger operator, essential spectrum, bound
states, Fredholm determinant

##  1\. Introduction

In this paper we consider the family $\widehat{H}_{\mu\lambda}(K)$ Schrödinger
operators, associated to the Bose-Hubbard Hamiltonian of a system of two
identical bosons on the two dimensional lattice $\mathbb{Z}^{2}$ with on-site
interaction $\mu\in\mathbb{R}$ and nearest-neighbor interaction
$\lambda\in\mathbb{R}.$

Lattice Bose-Hubbard models have become popular in recent years since they
represent a minimal, natural Hamiltonian in ultracold atoms in optical
lattices, systems with highly controllable parameters such as lattice geometry
and dimensionality, particle masses, tunneling, two-body potentials,
temperature etc. (see e.g., [2, 6, 7, 17] and references therein). Unlike the
traditional condensed matter systems, where stable composite objects are
usually formed by attractive forces and repulsive forces separate particles in
free space, the controllability of collision properties of ultracold atoms has
enabled to experimentally observe a stable _repulsive_ bound pair of ultracold
atoms in the optical lattice $\mathbb{Z}^{3}$, see e.g., [4, 20, 23, 24]. In
all these observations Bose-Hubbard Hamiltonians became a link between the
theoretical basis and experimental results.

The main difficulty in solving Bose-Hubbard Hamiltonian even with the on-site
interaction is due to the tunneling, i.e., the kinetic energy necessary for a
boson to hop from site to site, since it is highly nonlocal. Moreover, unlike
its continuous counterpart, the lattice Hamiltonian, corresponding to short-
range interacting sytems of particle pairs, is not rotationally invariant,
hence, the separation of lattice Hamiltonian related to the center of mass-
motion is not possible. However, the translation-invariance of the Hamiltonian
(in $d$-dimensional lattice $\mathbb{Z}^{d}$) allows to use the Floquet-Bloch
decomposition: the underlying Hilbert space $\ell^{2}(\mathbb{Z}^{d})^{2}$ and
the total two-particle Hamiltonian $\widehat{\mathbb{H}}$ are represented as
the von Neumann direct integral

$\ell^{2}[(\mathbb{Z}^{d})^{2}]\simeq\int\limits_{K\in\mathbb{T}^{d}}\oplus\ell^{2}(\mathbb{Z}^{d})\mathrm{d}K,\qquad\widehat{\mathbb{H}}\simeq\int\limits_{K\in\mathbb{T}^{d}}\oplus\widehat{H}(K)\mathrm{d}K,$

associated to the representation of the discrete group $\mathbb{Z}^{d}$ by the
shift operators, here $\mathbb{T}^{d}$ is the $d$-dimensional torus, a
Pontryagin dual group to $\mathbb{Z}^{d}.$ The fiber Hamiltonians
$\widehat{H}_{\mu}(K)$ nontrivially depend on the so-called quasi-momentum
$K\in\mathbb{T}^{d}$ (see e.g., [1, 3, 18, 19]).

In what follows any nonzero solution $\widehat{\psi}_{K}$ to the Schrödinger
equation

$\widehat{H}(K)\widehat{\psi}_{K}=e_{K}\widehat{\psi}_{K},\qquad\widehat{\psi}_{K}\in\ell^{2}(\mathbb{Z}^{d}),$

will be called bound state with energy $e_{K}$.

In the current paper, we study the family

$\widehat{H}_{\mu\lambda}(K):=\widehat{H}_{0}(K)+\widehat{V}_{\mu\lambda},\qquad
K\in\mathbb{T}^{2},$

discrete Schrödinger operators associated to the Hamiltonian
$\widehat{\mathbb{H}}_{\mu\lambda}$ of a system of two identical bosons on the
two dimensional lattice $\mathbb{Z}^{2}$ with zero-range on-site interaction
$\mu\neq 0$ and nearest-neighbor interaction $\lambda\neq 0$ (see (2.5) in
Subsection 2.3). To the best of our knowledge, this is a new, exactly solvable
model, for which the exact number of eigenvalues and their locations as well
as exact lower and upper bounds for all values of the pair interactions
$\mu,\lambda\in\mathbb{R}$ can be found.

First we observe that the essential spectrum of $\widehat{H}_{\mu\lambda}(K)$
consists of a segment
$[\mathcal{E}_{\min}(K),\allowbreak\mathcal{E}_{\max}(K)],$ where

$\mathcal{E}_{\min}(K):=2\sum\limits_{i=1}^{2}\Big{(}1-\cos\tfrac{K_{i}}{2}\Big{)},\qquad\mathcal{E}_{\max}(K):=2\sum\limits_{i=1}^{2}\Big{(}1+\cos\tfrac{K_{i}}{2}\Big{)}$

(see Subsection 2.4).

Now we study the discrete spectrum of $\widehat{H}_{\mu\lambda}(K).$ The
eigenvalue problem for $\widehat{H}_{\mu\lambda}(K)$ with general $K$ is not
an easy problem; note that $\widehat{H}_{\mu\lambda}(K)$ becomes a small
perturbation of $\widehat{V}_{\mu\lambda}$ if $K$ is close to
$\vec{\pi}:=(-\pi,\pi),$ and $\widehat{V}_{\mu\lambda}$ is a small pertubation
of $\widehat{H}_{0}$ provided that $\mu,\lambda$ is small. However, it turns
out that for any $\mu,\lambda\in\mathbb{R},$ if $e_{n}(K;\mu,\lambda)$ and
$E_{n}(K;\mu,\lambda)$ are the $n$-th eigenvalues of
$\widehat{H}_{\mu\lambda}(K)$ below and above the essential spectrum,
respectively, then the functions $\mathcal{E}_{\min}(K)-e_{n}(K;\mu,\lambda)$
and $E_{n}(K;\mu,\lambda)-\mathcal{E}_{\max}(K)$ as a function of $K$ will
have a minimum at $K=0$ (see Lemma 4.9). Subsequently, the number of
eigenvalues $\mathrm{n}_{+}(\widehat{H}_{\mu\lambda}(K))$ (resp.
$\mathrm{n}_{-}(\widehat{H}_{\mu\lambda}(K))$) of
$\widehat{H}_{\mu\lambda}(K)$ above (resp. below) the essential spectrum is
not less than that of $\widehat{H}_{\mu\lambda}(0)$ (Theorem 3.1). This result
is a generalization of [1, Theorems 1 and 2] which contain the assertion
related only to the ground state of $\widehat{H}_{\mu\lambda}(K).$

To study the discrete spectrum of $\widehat{H}_{\mu\lambda}(0)$ we introduce
the Fredholm determinants $\Delta_{\mu\lambda}^{(s)}(z)$ and
$\Delta_{\lambda}^{(a)}(z),$ associated to the restriction of
$\widehat{H}_{\mu\lambda}(0)$ onto symmetric and antisymmetric functions in
$\ell^{2,e}(\mathbb{Z}^{2})=\ell^{2,e}(\mathbb{Z}\times\mathbb{Z})$ w.r.t.
coordinate permutations111Note that the terms ’symmetric’ and ’antisymmetric’
are not related to the particle symmetry, but rather to the invariant
subspaces of $\widehat{H}_{\mu\lambda}(0)$ related to the permutation of
coordinates in $\mathbb{Z}^{2}:=\mathbb{Z}\times\mathbb{Z}$ of one-particle.,
respectively. It is well-known that eigenvalues of $H_{\mu\lambda}(0)$ are
zeros of these determinants, and vice versa (see also Lemma 4.5). Moreover,
the number of zeros of $\Delta_{\mu\lambda}^{(s)}(z)$ and
$\Delta_{\lambda}^{(a)}(z)$ can change if and only if their asymptotics as $z$
approaches to the threholds of the essential spectrum vanish.

Therefore, we partition $(\mu,\lambda)$-plane of interactions into connected
components by means of hyperbolas $\lambda\mu+4\lambda+2\mu=0$ and
$\lambda\mu-4\lambda-2\mu=0$ and straight lines $\lambda=\pm\frac{\pi}{4-\pi}$
(see Figure 1), where $\lambda\mu\pm 4\lambda\pm 2\mu$ and
$\lambda=\frac{\pi}{4-\pi}$ appear as constant in front of the main term in
the asymptotics of the Fredholm determinants $\Delta_{\mu\lambda}^{(s)}(z)$
and $\Delta_{\lambda}^{(a)}(z),$ respectively, as $z$ converges to the edges
of the essential spectrum (Proposition 4.4). As we noticed above, the number
of eigenvalues of $H_{\mu\lambda}(0)$ changes if and only if one of the
’constants’ $\lambda\mu\pm 4\lambda\pm 2\mu$ and $\lambda\pm\frac{\pi}{4-\pi}$
changes its sign (see also Lemmas 4.6-4.8).

Figure 1. A partition of the $(\mu,\lambda)$-plane using the hyperbolas
$\lambda\mu\pm 4\lambda\pm 2\mu=0$ and straight lines
$\lambda=\pm\frac{\pi}{4-\pi}.$ The figure shows also the dynamics of change
in the number of eigenvalues: in the region $[i,j]$ the operator
$H_{\mu\lambda}(0)$ has exactly $i$ eigenvalues below the essential spectrum
and exactly $j$ eigenvalues above the essential spectrum. If $(\mu,\lambda)$
moves from the region $[i,j]$ to the region $[i-1,j]$ resp. $[i,j-1]$ one
eigenvalue of $H_{\mu\lambda}(0)$ is absorbed into the essential spetrum at
the lower resp. upper edge of the essential spectrum. Conversely, if
$(\mu,\lambda)$ moves from the region $[i,j]$ to the region $[i+1,j]$ resp.
$[i,j+1]$ one eigenvalue of $H_{\mu\lambda}(0)$ is released the essential
spetrum at the lower resp. upper edge of the essential spectrum. This
absoption and release of eigenvalues can occur simultaneously: for instance,
if $(\mu,\lambda)$ moves from the region $[1,0]$ to $[0,1]$ then one
eigenvalue is absorbed at the bottom and another releases at the top.

Hence, while $(\mu,\lambda)$ runs in $\mathbb{R}^{2}$ and does not cross any
of these curves, no qualitative or quantitative changes occur in the discrete
spectrum of $\widehat{H}_{\mu\lambda}(0);$ however, as soon as $(\mu,\lambda)$
crosses any of those hyperbolas resp. straight lines, the essential spectrum
of $\widehat{H}_{\mu\lambda}(0)$ either ’gives birth’ or ’absorb’ a bound
state of $\widehat{H}_{\mu\lambda}(0)$ which is symmetric resp. antisymmetric
w.r.t. permutation of coordinates (Theorem 3.4).

This fact is tightly connected to _coupling constant threshold phenomenon_
[10, 14, 15]: if $H(t),$ $t\geq 0,$ is a one-parameter family of Schrödinger
operators (in continuum or lattice) and $0$ is a lower edge of the essential
spectrum, then $t_{0}\geq 0$ is a coupling constant threshold if and only if
there exists a negative eigenvalue $E(t)$ of $H(t)$ for $t>t_{0}$ such that
$E(t)\nearrow 0$ as $t\searrow t_{0},$ i.e., as $t\searrow t_{0}$ an
eigenvalue is absorbed at the threshold of continuum, and conversely, as
$t\nearrow t_{0}+\varepsilon$ an eigenvalue is released from the continuum.
Indeed, considering $\widehat{H}(t):=\widehat{H}_{t\mu,t\lambda}$ for $t\geq
0,$ from Figure 1 below we immediately obtain that the only coupling constant
thresholds are $t_{0}=0,$
$t_{0}=\big{|}\frac{4\lambda+2\mu}{\lambda\mu}\big{|}$ (when $\lambda\mu\neq
0$) and $t_{0}=\frac{\pi}{4-\pi}.$

Surprisingly, the maximum number of isolated eigenvalues is achieved only in
four connected components in which both $\mu$ and $\lambda$ run on infinite
intervals.

In general, some results such as existence of an eigenvalue and also the
finiteness of the number of eigenvalues can be obtained for wide classes of
operators (see e.g., [8, 9, 10, 22]). However, Figure 1 shows that the study
of a qualitative change in the number of eigenvalues of
$\widehat{H}_{\mu\lambda}(K),$ even for $K=0,$ is very delicate: there are
discs in the $(\mu,\lambda)$-plane with arbitrarily small radius in which the
number of eigenvalues has jump (see Theorem 3.4).

Recall that in [8] authors extend the results of [9, 22] to the lattice case.
Namely, the spectral properties of one-particle discrete Schrödinger operator

$\widehat{\bf h}_{t}=\widehat{\bf h}_{0}+t\widehat{\bf v},\qquad t\geq 0,$

in $\mathbb{Z}^{1}$ and $\mathbb{Z}^{2}$ have been studied, here $\widehat{\bf
h}_{0}$ is a self-adjoint Laurent-Toeplitz-type operator generated by a
dispersion relation $\widehat{\mathcal{E}}:\mathbb{Z}^{d}\to\mathbb{C}$ of the
particle and the potential $\widehat{\bf v}$ is the multiplication operator by
$\widehat{v}:\mathbb{Z}^{d}\to\mathbb{R}.$ Under certain regularity assumption
on $\widehat{\mathcal{E}}$ and a decay assumption on $\widehat{v}$, authors
show that if $\sum\widehat{v}(x)$ is nonnegative resp. nonpositive, then the
discrete spectrum of $\widehat{\bf h}_{t}$ above resp. below its essential
spectrum is non-empty for any $t>0.$ Moreover, authors prove the existence of
$t_{0}>0$ (depending only on $\widehat{v}$ and $\widehat{\mathcal{E}}$) such
that

  * (a)

if $\sum\widehat{v}(x)>0,$ then for any $t\in(0,t_{0})$ the operator
$\widehat{\bf h}_{t}$ has no eigenvalues below the essential spectrum and has
a unique eigenvalue above the essential spectrum;

  * (b)

if $\sum\widehat{v}(x)<0,$ then for any $t\in(0,t_{0})$ the operator
$\widehat{\bf h}_{t}$ has no eigenvalues above the essential spectrum and has
a unique eigenvalue below the essential spectrum;

  * (c)

if $\sum\widehat{v}(x)=0,$ then for any $t\in(0,t_{0})$ the operator
$\widehat{\bf h}_{t}$ has unique eigenvalues in both below and above the
essential spectrum.

Figure 2. The definition of $t_{0}$ given by [8].

Note that for our model $\sum\widehat{v}(x)=\mu+2\lambda.$ Hence, our results
for the discrete spectrum of $\widehat{H}_{\mu\lambda}(0)$ show the (exact!)
dependence of $t_{0}$ on $\mu$ and $\lambda:$ namely, $t_{0}$ will the
intersection of the half-line $\\{(t\mu,t\lambda):\,\,t\geq 0\\}$ with the
branches of hyperbolas $\lambda\mu\pm 4\lambda\pm 2\mu=0$ passing through the
origin (see Figure 2).

Moreover, from Figure 2 we observe that the unique eigenvalue of [8] appears
from the symmetric space and eigenvalues in the antisymmetric space come out
later (see Theorem 4.2).

In general, it is well-known that in the case of $d\geq 3$ or in the case of
fermions with $d\geq 1$ the bound states appear from the essential spectrum
either as a threshold bound state or as a threshold resonance [11, 12].
However, our results show that in the case $d=2,$ even though antisymmetric
bound states come out from the threshold eigenvalues, all symmetric bound
states of $\widehat{H}_{\mu\lambda}(0)$ appear (only!) from the singularity of
the associated Fredholm determinant at the thresholds, namely, associated
Fredholm determinant cannot be analytically (even continuously) extended to
the edges of the essential spectrum (see Proposition 4.4). Such a result holds
also in $d=1$ (see [16]). This is a strict mathematical explanation in the
difference of the appearance of eigenvalues for $d=1,2$ and $d\geq 3,$ or for
bosons and fermions.

As far as we know for continuous two-body Schrödinger operators
$\mathbb{R}^{2}$ there are no analogous examples and results. In the discrete
case similar results for the number of eigenvalues of one-particle Schrödinger
operators in $\mathbb{Z}^{d}$ with zero-range on-site and nearest-neighbor
interactions have been obtained, for instance, in [12] for $d=3$ with
attractive potential field, in [16] for $d=1,$ and in [5] for all $d\geq 1$
considering only negative eigenvalues.

The paper is organized as follows: In Section 2 we introduce the Hamiltonian
$\widehat{\mathbb{H}}_{\mu\lambda}$ of a system of two bosons in the position
and momentum space and also the Schrödinger operator
$\widehat{H}_{\mu\lambda}(K)$ associated to
$\widehat{\mathbb{H}}_{\mu\lambda}$. Main results of the paper are stated in
Section 3 and their proofs are contained in Section 4. We conclude the paper
with an appendix containing the proof of Proposition 4.4.

##  2\. Discrete Schrödinger operators on lattices

###  2.1. The two-particle Hamiltonian: the position-space representation

Let $\mathbb{Z}^{2}$ be the two dimensional lattice and let
$\ell^{2,s}(\mathbb{Z}^{2}\times\mathbb{Z}^{2})$ be the Hilbert space of
square-summable symmetric functions on $\mathbb{Z}^{2}\times\mathbb{Z}^{2}.$

In the position-space representation the two particle Hamiltonian
$\widehat{\mathbb{H}}_{\mu\lambda},$ associated to a system of two bosons
interacting via zero-range and nearest-neighbor potential
$\widehat{v}_{\mu\lambda}$ is a bounded self-adjoint operator acting in
$\ell^{2,s}(\mathbb{Z}^{2}\times\mathbb{Z}^{2})$ as

$\widehat{\mathbb{H}}_{\mu\lambda}=\widehat{\mathbb{H}}_{0}+\widehat{\mathbb{V}}_{\mu\lambda},\qquad\mu,\lambda\in\mathbb{R}.$

Here the free Hamiltonian $\widehat{\mathbb{H}}_{0}$ of a system of two
identical particles (bosons) is a bounded self–adjoint operator acting in
$\ell^{2,s}(\mathbb{Z}^{2}\times\mathbb{Z}^{2})$ as

$\widehat{\mathbb{H}}_{0}\widehat{f}(x,y)=\sum_{n\in\mathbb{Z}^{2}}\widehat{\varepsilon}(x-n)\widehat{f}(n,y)+\sum\limits_{n\in\mathbb{Z}^{2}}\widehat{\varepsilon}(y-n)\widehat{f}(x,n),$

where

$\widehat{\varepsilon}(s)=\begin{cases}2&\text{if $|s|=0,$}\\\
-\frac{1}{2}&\text{if $|s|=1,$}\\\ 0&\text{if $|s|>1,$}\end{cases}$ (2.1)

and $|s|=|s_{1}|+|s_{2}|$ for $s=(s_{1},s_{2})\in\mathbb{Z}^{2}.$

The interaction $\widehat{\mathbb{V}}_{\mu\lambda}$ is the multiplication
operator

$\widehat{\mathbb{V}}_{\mu\lambda}\widehat{f}(x,y)=\widehat{v}_{\mu\lambda}(x-y)\widehat{f}(x,y),$

given by the function

$\widehat{v}_{\mu\lambda}(s)=\begin{cases}\mu&\text{if $|s|=0,$}\\\
\frac{\lambda}{2}&\text{if $|s|=1,$}\\\ 0&\text{if $|s|>1.$}\end{cases}$ (2.2)

###  2.2. The two-particle Hamiltonian: the momentum-space representation

Let
$\mathbb{T}^{2}=\allowbreak(\mathbb{R}/2\pi\mathbb{Z})^{2}\equiv[-\pi,\pi)^{2}$
be the two dimensional torus, the Pontryagin dual group of $\mathbb{Z}^{2}$,
equipped with the Haar measure $\mathrm{d}p,$and let
$L^{2,s}(\mathbb{T}^{2}\times\mathbb{T}^{2})$ be the Hilbert space of square-
integrable symmetric functions on $\mathbb{T}^{2}\times\mathbb{T}^{2}.$ Let
$\mathcal{F}:\ell^{2}(\mathbb{Z}^{2})\rightarrow L^{2}(\mathbb{T}^{2})$ be the
standard Fourier transform

$\mathcal{F}\widehat{f}(p)=\frac{1}{2\pi}\sum_{x\in\mathbb{Z}^{2}}\widehat{f}(x)e^{ip\cdot
x},$

where $p\cdot x:=p_{1}x_{1}+p_{2}x_{2}$ for $p=(p_{1},p_{2})\in\mathbb{T}^{2}$
and $x=(x_{1},x_{2})\in\mathbb{Z}^{2}.$

In the momentum-space via $\mathcal{F}\otimes\mathcal{F}$ the two-particle
Hamiltonian is represented in $L^{2,s}(\mathbb{T}^{2}\times\mathbb{T}^{2})$ as

$\mathbb{H}_{\mu\lambda}:=(\mathcal{F}\otimes\mathcal{F})\widehat{\mathbb{H}}_{\mu\lambda}(\mathcal{F}\otimes\mathcal{F})^{*}:=\mathbb{H}_{0}+\mathbb{V}_{\mu\lambda}.$

Here the free Hamiltonian
$\mathbb{H}_{0}=(\mathcal{F}\otimes\mathcal{F})\widehat{\mathbb{H}}_{0}(\mathcal{F}\otimes\mathcal{F})^{*}$
is the multiplication operator:

$\mathbb{H}_{0}f(p,q)=[\varepsilon(p)+\varepsilon(q)]f(p,q),$

where

$\varepsilon(p):=\sum\limits_{i=1}^{2}\big{(}1-\cos p_{i}),\quad
p=(p_{1},p_{2})\in\mathbb{T}^{2},$

is the _dispersion relation_ of a single boson. The interaction
$\mathbb{V}_{\mu\lambda}=(\mathcal{F}\otimes\mathcal{F})\widehat{\mathbb{V}}_{\mu\lambda}(\mathcal{F}\otimes\mathcal{F})^{*}$
is the (partial) integral operator

$\mathbb{V}_{\mu\lambda}f(p,q)=\frac{1}{(2\pi)^{2}}\int_{\mathbb{T}^{2}}v_{\mu\lambda}(p-u)f(u,p+q-u)\mathrm{d}u,$

where

$v_{\mu\lambda}(p)=\mu+\lambda\sum_{i=1}^{2}\cos p_{i},\quad
p=(p_{1},p_{2})\in\mathbb{T}^{2}.$

###  2.3. The Floquet-Bloch decomposition of $\mathbb{H}_{\mu\lambda}$ and
discrete Schrödinger operators

Since $\widehat{H}_{\mu\lambda}$ commutes with the representation of the
discrete group $\mathbb{Z}^{2}$ by shift operators on the lattice, we can
decompose the space $L^{2,s}(\mathbb{T}^{2}\times\mathbb{T}^{2})$ and
$\mathbb{H}_{\mu\lambda}$ into the von Neumann direct integral as

$L^{2,s}(\mathbb{T}^{2}\times\mathbb{T}^{2})\simeq\int\limits_{K\in\mathbb{T}^{2}}\oplus
L^{2,e}(\mathbb{T}^{2})\,\mathrm{d}K$ (2.3)

and

$\mathbb{H}_{\mu\lambda}\simeq\int\limits_{K\in\mathbb{T}^{2}}\oplus
H_{\mu\lambda}(K)\,\mathrm{d}K,$ (2.4)

where $L^{2,e}(\mathbb{T}^{2})$ is the Hilbert space of square-integrable even
functions on $\mathbb{T}^{2}$ (see, e.g., [1]).

The fiber operator $H_{\mu\lambda}(K),$ $K\in\mathbb{T}^{2},$ is a self-
adjoint operator defined in $L^{2,e}(\mathbb{T}^{2})$ as

$H_{\mu\lambda}(K):=H_{0}(K)+V_{\mu\lambda},$

where the unperturbed operator $H_{0}(K)$ is the multiplication operator by
the function

$\mathcal{E}_{K}(p):=2\sum_{i=1}^{2}\Big{(}1-\cos\tfrac{K_{i}}{2}\,\cos
p_{i}\Big{)},$

and the perturbation $V_{\mu\lambda}$ is defined as

$V_{\mu\lambda}f(p)=\frac{1}{(2\pi)^{2}}\int_{\mathbb{T}^{2}}\Big{(}\mu+\lambda\sum\limits_{i=1}^{2}\cos
p_{i}\cos q_{i}\Big{)}f(q)\mathrm{d}q.$

In the literature the parameter $K\in\mathbb{T}^{2}$ is called the _two-
particle quasi-momentum_ and the fiber $H_{\mu\lambda}(K)$ is called the
_discrete Schrödinger operator_ associated to the two-particle Hamiltonian
$\widehat{\mathbb{H}}_{\mu\lambda}.$

Using the Fourier transform (2.3) and (2.4) can be represented as

$\ell^{2,s}(\mathbb{Z}^{2}\times\mathbb{Z}^{2})\simeq\int\limits_{K\in\mathbb{T}^{2}}\oplus\ell^{2,e}(\mathbb{Z}^{2})\,\mathrm{d}K$

and

$\widehat{\mathbb{H}}_{\mu\lambda}\simeq\int\limits_{K\in\mathbb{T}^{2}}\oplus\widehat{H}_{\mu\lambda}(K)\,\mathrm{d}K,$

where $\ell^{2,e}(\mathbb{Z}^{2})$ is the Hilbert space of square-summable
even functions on $\mathbb{Z}^{2}$ and

$\widehat{H}_{\mu\lambda}(K):=\mathcal{F}^{*}H_{\mu\lambda}(K)\mathcal{F},$

where

$\widehat{H}_{\mu\lambda}(K)=\widehat{H}_{0}(K)+\widehat{V}_{\mu\lambda}(K),$
(2.5)

and

$\widehat{H}_{0}(K)f(x)\ =\
\sum_{s\in\mathbb{Z}^{2}}\widehat{\mathcal{E}}_{K}(x-s)\widehat{f}(s),\qquad
f\in\ell^{2,e}\mathbb{Z}^{2}$

with

$\widehat{\mathcal{E}}_{K}(x)=2\widehat{\varepsilon}(x)\,\cos\frac{K\cdot
x}{2}$

and the operator $\widehat{V}_{\mu\lambda}$ acts in
$\ell^{2,e}(\mathbb{Z}^{2})$ as

$\widehat{V}_{\mu\lambda}f(x)=\widehat{v}_{\mu\lambda}(x)f(x),$

and the functions $\widehat{\varepsilon}$ and $\widehat{v}_{\mu\lambda}$ are
given by (2.1) and (2.2), respectively.

###  2.4. The essential spectrum of discrete Schrödinger operators

Since $V_{\mu\lambda}$ has rank at most three, by Weyl’s Theorem for any
$K\in\mathbb{T}^{2}$ the essential spectrum
$\sigma_{\mathrm{ess}}(H_{\mu\lambda}(K))$ coincides with the spectrum of
$H_{0}(K),$ i.e.,

$\sigma_{\mathrm{ess}}(H_{\mu\lambda}(K))=\sigma(H_{0}(K))=[\mathcal{E}_{\min}(K),\mathcal{E}_{\max}(K)],$
(2.6)

where

$\displaystyle\mathcal{E}_{\min}(K):=$
$\displaystyle\min_{p\in\mathbb{T}^{2}}\,\mathcal{E}_{K}(p)=2\sum\limits_{i=1}^{2}\Big{(}1-\cos\tfrac{K_{i}}{2}\Big{)}\geq
0=\mathcal{E}_{\min}(0),$ $\displaystyle\mathcal{E}_{\max}(K):=$
$\displaystyle\max_{p\in\mathbb{T}^{2}}\,\mathcal{E}_{K}(p)=2\sum\limits_{i=1}^{2}\Big{(}1+\cos\tfrac{K_{i}}{2}\Big{)}\leq
8=\mathcal{E}_{\max}(0).$

##  3\. Main results

Our first main result is the following generalization of [1, Theorems 1 and
2].

######  Theorem 3.1.

Suppose that $H_{\mu\lambda}(0)$ has $n$ eigenvalues below resp. above the
essential spectrum for some $\mu,\lambda\in\mathbb{R}.$ Then for every
$K\in\mathbb{T}^{2}$ the operator $H_{\mu\lambda}(K)$ has at least $n$
eigenvalues below resp. above its essential spectrum.

Next we find exact lower bound for the number of eigenvalues of
$H_{\mu\lambda}(K)$ depending only on $\mu$ and $\lambda.$

In the $(\mu,\lambda)$-plane let us define the following nine sets:

$\displaystyle\textit{S}_{01}:=$
$\displaystyle\Big{\\{}(\mu,\lambda)\in\mathbb{R}^{2}:\,\,\lambda>\frac{\pi}{4-\pi}\Big{\\}},$
(3.1) $\displaystyle\textit{S}_{00}:=$
$\displaystyle\Big{\\{}(\mu,\lambda)\in\mathbb{R}^{2}:\,\,|\lambda|<\frac{\pi}{4-\pi}\Big{\\}},$
$\displaystyle\textit{S}_{10}:=$
$\displaystyle\Big{\\{}(\mu,\lambda)\in\mathbb{R}^{2}:\,\,\lambda<\frac{\pi}{\pi-4}\Big{\\}},$
$\displaystyle\allowbreak\textit{C}_{0}^{+}:=$
$\displaystyle\Big{\\{}(\mu,\lambda)\in\mathbb{R}^{2}:\,\,\lambda\mu-4\lambda-2\mu>0,\,\,\lambda<2\Big{\\}},$
$\displaystyle\textit{C}_{1}^{+}:=$
$\displaystyle\Big{\\{}(\mu,\lambda)\in\mathbb{R}^{2}:\,\,\lambda\mu-4\lambda-2\mu<0\Big{\\}},$
$\displaystyle\textit{C}_{2}^{+}:=$
$\displaystyle\Big{\\{}(\mu,\lambda)\in\mathbb{R}^{2}:\,\,\lambda\mu-4\lambda-2\mu>0,\,\,\lambda>2\Big{\\}},$
$\displaystyle\textit{C}_{0}^{-}:=$
$\displaystyle\Big{\\{}(\mu,\lambda)\in\mathbb{R}^{2}:\,\,\lambda\mu+4\lambda+2\mu>0,\,\,\lambda>-2\Big{\\}},$
$\displaystyle\textit{C}_{1}^{-}:=$
$\displaystyle\Big{\\{}(\mu,\lambda)\in\mathbb{R}^{2}:\,\,\lambda\mu+4\lambda+2\mu<0\Big{\\}},$
$\displaystyle\textit{C}_{2}^{-}:=$
$\displaystyle\Big{\\{}(\mu,\lambda)\in\mathbb{R}^{2}:\,\,\lambda\mu+4\lambda+2\mu>0,\,\,\lambda<-2\Big{\\}}$

(see Figures 3-3 below).

(a) (a)

(b) (b)

(c) (c)

(d) (d)

Figure 3. Schematic locations of sets $\textit{S}_{ij}$ and
$\textit{C}_{i}^{\pm}$ in (3.1).

Let $\mathrm{n}_{+}(H_{\mu\lambda}(K))$ resp.
$\mathrm{n}_{-}(H_{\mu\lambda}(K))$ be the number of $H_{\mu\lambda}(K)$ above
resp. below its essential spectrum.

######  Theorem 3.2.

Let $K\in\mathbb{T}^{2}$ and $(\mu,\lambda)\in\mathbb{R}^{2}.$ Then

$\displaystyle(\mu,\lambda)\in\textit{C}_{2}^{+}\cap\textit{S}_{01}$
$\displaystyle\Longrightarrow\qquad\mathrm{n}_{+}(H_{\mu\lambda}(K))=3,$ (3.2)
$\displaystyle(\mu,\lambda)\in\textit{C}_{2}^{+}\Delta\textit{S}_{01}$
$\displaystyle\Longrightarrow\qquad\mathrm{n}_{+}(H_{\mu\lambda}(K))\geq 2,$
$\displaystyle(\mu,\lambda)\in\textit{C}_{1}^{+}\setminus\textit{S}_{01}$
$\displaystyle\Longrightarrow\qquad\mathrm{n}_{+}(H_{\mu\lambda}(K))\geq 1,$
$\displaystyle(\mu,\lambda)\in\overline{\textit{C}_{0}^{+}}$
$\displaystyle\Longrightarrow\qquad\mathrm{n}_{+}(H_{\mu\lambda}(K))=0,$

and

$\displaystyle(\mu,\lambda)\in\textit{C}_{2}^{-}\cap\textit{S}_{10}$
$\displaystyle\Longrightarrow\qquad\mathrm{n}_{-}(H_{\mu\lambda}(K))=3,$ (3.3)
$\displaystyle(\mu,\lambda)\in\textit{C}_{2}^{-}\Delta\textit{S}_{10}$
$\displaystyle\Longrightarrow\qquad\mathrm{n}_{-}(H_{\mu\lambda}(K))\geq 2,$
$\displaystyle(\mu,\lambda)\in\textit{C}_{1}^{-}\setminus\textit{S}_{10}$
$\displaystyle\Longrightarrow\qquad\mathrm{n}_{-}(H_{\mu\lambda}(K))\geq 1,$
$\displaystyle(\mu,\lambda)\in\overline{\textit{C}_{0}^{-}}$
$\displaystyle\Longrightarrow\qquad\mathrm{n}_{-}(H_{\mu\lambda}(K))=0,$

where $A\Delta B:=(A\setminus B)\cup(B\setminus A)$ is the symmetric
difference of sets and $\overline{A}$ is the closure of a set $A.$

Theorem 3.2 provides a lower bound for the number of eigenvalues for both
sides of the essential spectrum of $H_{\mu\lambda}(K)$. Recall that by the
min-max principle $H_{\mu\lambda}(K)$ can have at most three bound states
outside its essential spectrum.

######  Remark 3.3.

Theorem 3.2 implies that in some subsets of the $(\mu,\lambda)$-plane (for
example in $\textit{C}_{1}^{+}\cap\textit{C}_{1}^{-}$ ) the eigenvalues of
$H_{\mu\lambda}(K)$ can appear simultaneously on both sides of the essential
spectrum.

The following result shows that the estimates for
$\mathrm{n}_{\pm}(H_{\mu\lambda}(K))$ given by Theorem 3.2 are sharp.

######  Theorem 3.4.

Let $K=0.$ Then all inequalities in (3.2) and (3.3) are in fact equalities.
Moreover, the spaces

$L^{2,e,s}(\mathbb{T}^{2}):=\Big{\\{}f\in
L^{2,e}(\mathbb{T}^{2}):\,\,f(p_{1},p_{2})=f(p_{2},p_{1}),\,\,p_{1},p_{2}\in\mathbb{T}\Big{\\}}$

and

$L^{2,e,a}(\mathbb{T}^{2}):=\Big{\\{}f\in
L^{2,e}(\mathbb{T}^{2}):\,\,f(p_{1},p_{2})=-f(p_{2},p_{1}),\,\,p_{1},p_{2}\in\mathbb{T}\Big{\\}}$

of symmetric and antisymmetric even functions are invariant with respect to
$H_{\mu\lambda}(0)$ and:

  * (1)

if $(\mu,\lambda)\in\overline{\textit{S}_{00}},$ then $H_{\mu\lambda}(0)$ has
no antisymmetric bound states outside the essential spectrum;

  * (2)

if $(\mu,\lambda)\in\textit{S}_{01}$ resp. $(\mu,\lambda)\in\textit{S}_{10},$
then $H_{\mu\lambda}(0)$ has a unique antisymmetric bound state above resp.
below the essential spectrum;

  * (3)

if $(\mu,\lambda)\in\textit{C}_{2}^{+}$ resp.
$(\mu,\lambda)\in\textit{C}_{2}^{-},$ then $H_{\mu\lambda}(0)$ has a exactly
two symmetric bound states above resp. below the essential spectrum;

  * (4)

if $(\mu,\lambda)\in\textit{C}_{1}^{+}\cup\partial\textit{C}_{2}^{+}$ resp.
$(\mu,\lambda)\in\textit{C}_{1}^{-}\cup\partial\textit{C}_{2}^{-},$ then
$H_{\mu\lambda}(0)$ has a unique symmetric bound state above resp. below the
essential spectrum, where $\partial A$ is the topological boundary of a set
$A;$

  * (5)

if $(\mu,\lambda)\in\overline{\textit{C}_{0}^{+}}$ resp.
$(\mu,\lambda)\in\overline{\textit{C}_{0}^{-}},$ then $H_{\mu\lambda}(0)$ has
no symmetric bound states above resp. below the essential spectrum.

##  4\. Proofs of the main results

###  4.1. Discrete spectrum of $H_{\mu\lambda}(0)$

Unlike the case $K\neq 0$ in the case $K=0$ the Fredholm determinant
$\Delta_{\mu\lambda}(0,z)$ is easier to study. Notice that the operator
$H_{0}(0)$ is the multiplication operator by the symmetric function
$\mathcal{E}_{0}(p)=2\varepsilon(p)$ in $L^{2,e}(\mathbb{T}^{2}).$ Hence, both
$L^{2,e,s}(\mathbb{T}^{2})$ and $L^{2,e,a}(\mathbb{T}^{2})$ are invariant with
respect to $H_{0}(0).$ Moreover, since

$\displaystyle 2\cos p_{1}\cos q_{1}+2\cos p_{2}\cos q_{2}$ $\displaystyle=$
$\displaystyle(\cos p_{1}+\cos p_{2})(\cos q_{1}+\cos q_{2})+(\cos p_{1}-\cos
p_{2})(\cos q_{1}-\cos q_{2}),$

the spaces of symmetric and antisymmetric even functions in $\mathbb{T}^{2}$
are invariant also with respect $V_{\mu\lambda}.$ Thus,

######  Lemma 4.1.

$\sigma(H_{\mu\lambda}(0))=\sigma(H_{\mu\lambda}^{s})\cup\sigma(H_{\lambda}^{a}),$

where

$H_{\mu\lambda}^{s}:=H_{0}(0)+V_{\mu\lambda}^{s}\qquad\text{and}\qquad
H_{\lambda}^{a}:=H_{0}(0)+V_{\lambda}^{a},$

with

$V_{\mu\lambda}^{s}f(p)=\frac{\mu}{4\pi^{2}}\int_{\mathbb{T}^{2}}f(q)\,\mathrm{d}q+\frac{\lambda}{8\pi^{2}}\,(\cos
p_{1}+\cos p_{2})\int_{\mathbb{T}^{2}}(\cos q_{1}+\cos
q_{2})f(q)\,\mathrm{d}q$

and

$V_{\lambda}^{a}f(p)=\frac{\lambda}{8\pi^{2}}\,(\cos p_{1}-\cos
p_{2})\int_{\mathbb{T}^{2}}(\cos q_{1}-\cos q_{2})f(q)\,\mathrm{d}q,$

are the restrictions of $H_{\mu\lambda}(0)$ onto $L^{2,e,s}(\mathbb{T}^{2})$
and $L^{2,e,a}(\mathbb{T}^{2}).$

######  Theorem 4.2.

Fix $(\mu,\lambda)\in\mathbb{R}^{2}.$

  * (a)

If $2\lambda+\mu\geq 0,$ then $H_{\mu\lambda}^{s}$ has at least one eigenvalue
greater than $8.$

  * (b)

If $2\lambda+\mu\leq 0,$ then $H_{\mu\lambda}^{s}$ has at least one negative
eigenvalue.

###### Proof.

Note that the trace of $V_{\mu\lambda}^{s}$ is exactly $2\lambda+\mu.$ Hence,
the proof can be done along the essentially same lines of [8, Theorem 1.3]
introducing the non-symmetric Birman-Schwinger operator in
$L^{2,e,s}(\mathbb{T}^{2})$ in place of $L^{2}(\mathbb{T}^{2})$ and rewriting
it as a small perturbation of rank-one projection (which is not identically
zero in $L^{2,e,s}(\mathbb{T}^{2})$). ∎

######  Remark 4.3.

A statement in $L^{2,e,a}(\mathbb{T}^{2})$ analogous to Theorem 4.2 cannot be
proven using the arguments of [8, Theorem 1.3] since in this case the rank-one
projection (which is norm-close to the Birman-Schwinger operator) is
identically zero in $L^{2,e,a}(\mathbb{T}^{2})$.

Note that the essential spectrum of all Hamiltonians $H_{\mu\lambda}(0),$
$H_{\mu\lambda}^{s}$ and $H_{\mu\lambda}^{a}$ coincide with the segment
$[0,8].$ Let us find an (implicit) equation for the discrete eigenvalues of
$H_{\mu\lambda}(0).$ By Lemma 4.1 it is enough to solve

$H_{\mu\lambda}^{s}f=zf\qquad\text{resp.}\qquad H_{\lambda}^{a}f=zf$

in $z\in\mathbb{C}\setminus[0,8]$ and nonzero $f.$ Here we apply Fredholm’s
determinants theory (see, e.g., [21]). For shortness, writing

$\displaystyle
a(z):=\frac{1}{4\pi^{2}}\int_{\mathbb{T}^{2}}\frac{\mathrm{d}p}{\mathcal{E}_{0}(p)-z},\qquad$
$\displaystyle b(z):=\frac{1}{8\pi^{2}}\int_{\mathbb{T}^{2}}\frac{(\cos
p_{1}+\cos p_{2})^{2}\,\mathrm{d}p}{\mathcal{E}_{0}(p)-z},$ $\displaystyle
c(z):=\frac{1}{8\pi^{2}}\int_{\mathbb{T}^{2}}\frac{(\cos p_{1}+\cos
p_{2})\,\mathrm{d}p}{\mathcal{E}_{0}(p)-z},\qquad$ $\displaystyle
d(z):=\frac{1}{8\pi^{2}}\int_{\mathbb{T}^{2}}\frac{(\cos p_{1}-\cos
p_{2})^{2}\,\mathrm{d}p}{\mathcal{E}_{0}(p)-z},$

define the Fredholm determinants associated to $H_{\mu\lambda}^{s}$ and
$H_{\lambda}^{a},$ respectively, as

$\Delta_{\mu\lambda}^{(s)}(z):=\,\Delta_{\mu
0}(z)\Delta_{0\lambda}(z)-2\mu\lambda c(z)^{2},$

and

$\Delta_{\lambda}^{(a)}(z):=\,1+\lambda d(z),$

where

$\displaystyle\Delta_{\mu 0}(z):=1+\mu
a(z),\qquad\Delta_{0\lambda}(z):=1+\lambda b(z).$

Note that the functions $a(\cdot),$ $b(\cdot)$ and $d(\cdot)$ are analytic in
$\mathbb{C}\backslash[0,8]$, strictly increasing in
$\mathbb{R}\backslash[0,\,8],$ negative in $(8,+\infty)$ and positive in
$(-\infty,0).$ Their behaviour near $z=0$ and $z=8$ are established in the
following proposition.

######  Proposition 4.4.

There exists $\delta>0$ such that for every $\omega\in\\{a,b,c,d\\}$ there
exist functions $P_{\omega}^{0},$ and $Q_{\omega}^{0},$ analytic in the disc
$\\{z\in\mathbb{C}:\,\,|z|<\delta\\}$, and $P_{\omega}^{1},$ and
$Q_{\omega}^{1},$ analytic in the disc
$\\{z\in\mathbb{C}:\,\,|z-8|<\delta\\}$, such that

  * •

for every $z\in(-\delta,0)$

$\omega(z)=P_{\omega}^{0}(z)\ln(-z)+Q_{\omega}^{0}(z),$ (4.1)

where

$P_{a}^{0}(0)=-\frac{1}{4\pi},\quad P_{b}^{0}(0)=-\frac{1}{2\pi},\quad
P_{c}^{0}(0)=-\frac{1}{4\pi},\quad P_{d}^{0}(0)=0$

and

$Q_{a}^{0}(0)=\frac{5\ln 2}{4\pi},\quad Q_{b}^{0}(0)=\frac{5\ln
2-\pi}{2\pi},\quad Q_{c}^{0}(0)=\frac{5\ln 2-\pi}{4\pi},\quad
Q_{d}^{0}(0)=\frac{4-\pi}{\pi};$

  * •

for every $z\in(8,8+\delta)$

$\omega(z)=P_{\omega}^{1}(z)\ln(z-8)+Q_{\omega}^{1}(z),$ (4.2)

where

$P_{a}^{1}(8)=\frac{1}{4\pi},\quad P_{b}^{1}(8)=\frac{1}{2\pi},\quad
P_{c}^{1}(8)=\frac{1}{4\pi},\quad P_{d}^{1}(8)=0$

and

$Q_{a}^{1}(8)=-\frac{5\ln 2}{4\pi},\quad Q_{b}^{1}(8)=-\frac{5\ln
2-\pi}{2\pi},\quad Q_{c}^{1}(8)=-\frac{5\ln 2-\pi}{4\pi},\quad
Q_{d}^{1}(8)=\frac{\pi-4}{\pi}.$

The proof of Proposition 4.4 is postponed to the Appendix A.

######  Lemma 4.5.

A number $z\in\mathbb{C}\setminus[0,8]$ is an eigenvalue of
$H_{\mu\lambda}^{s}$ resp. $H_{\lambda}^{a}$ of multiplicity $m\geq 1$ if and
only if it is a zero of $\Delta_{\mu\lambda}^{(s)}(\cdot)$
resp.$\Delta_{\lambda}^{(a)}(\cdot)$ of multiplicity $m.$ Moreover, in
$\mathbb{R}\setminus[0,8]$ the function $\Delta_{\mu\lambda}^{(s)}(\cdot)$ has
at least one and at most two zeros and the function
$\Delta_{\lambda}^{(a)}(\cdot)$ has at most one zero.

###### Proof.

The first assertion follows from the Fredholm determinants theory. By Theorem
4.2 $H_{\mu\lambda}^{s}$ has at least one eigenvalue outside the essential
spectrum. Moreover, since $H_{\mu\lambda}^{s}$ is of rank two resp. rank one,
by the min-max principle, it has at most two eigenvalues outside the essential
spectrum. Hence, by the first part of the proposition
$\Delta_{\mu\lambda}^{(s)}(\cdot)$ has at least one and at most two zeros in
$\mathbb{R}\setminus[0,8].$

The last assertion follows from the rank-one property of $H_{\lambda}^{(a)}$
and the first part of the proposition. ∎

Next we study zeros of $\Delta_{\lambda}^{(a)}.$

######  Lemma 4.6.

Let $\lambda\in\mathbb{R}.$

  * (a)

If $\lambda\in[\frac{\pi}{\pi-4},\frac{\pi}{4-\pi}],$ then
$\Delta_{\lambda}^{(a)}(\cdot)>0$ in $\mathbb{R}\setminus[0,8].$

  * (b)

If $\lambda<\frac{\pi}{\pi-4},$ then $\Delta_{\lambda}^{(a)}(\cdot)$ has a
unique zero in $(-\infty,0)$ and $\Delta_{\lambda}^{(a)}(\cdot)>1$ in
$(8,+\infty).$

  * (c)

If $\lambda>\frac{\pi}{4-\pi},$ then $\Delta_{\lambda}^{(a)}(\cdot)$ has a
unique zero in $(8,\infty)$ and $\Delta_{\lambda}^{(a)}(\cdot)>1$ in
$(-\infty,0).$

###### Proof.

Note that the map $z\mapsto d(z)$ strictly increases in both connected
components of $\mathbb{R}\setminus[0,8].$ Hence, the equation
$\Delta_{\lambda}^{(a)}(z)=0$ has at most one zero in
$\mathbb{R}\setminus[0,8].$ Moreover,

$\lim\limits_{z\to\pm\infty}\Delta_{\lambda}^{(a)}(z)=1$ (4.3)

and by Proposition 4.4

$\lim\limits_{z\nearrow
0}\Delta_{\lambda}^{(a)}(z)=1+\frac{(4-\pi)\lambda}{\pi},\qquad\lim\limits_{z\searrow
8}\Delta_{\lambda}^{(a)}(z)=1+\frac{(\pi-4)\lambda}{\pi}.$ (4.4)

Now the assertions of the lemma follow from the strict monotonicity and
continuity of $z\mapsto\Delta_{\lambda}^{(a)}(z),$ and (4.3)-(4.4). ∎

The following lemma provides the dependence of the number of zeros of
$\Delta_{\mu\lambda}^{(s)}$ in $(8,+\infty)$ on $\mu$ and $\lambda.$

######  Lemma 4.7.

Let $(\mu,\lambda)\in\mathbb{R}^{2}.$

  * (a)

If $\lambda\mu-4\lambda-2\mu\geq 0$ and $\lambda<2,$ then
$\Delta_{\mu\lambda}^{(s)}$ has no zeros in $(8,+\infty).$

  * (b)

If $\lambda\mu-4\lambda-2\mu<0$ or $\lambda\mu-4\lambda-2\mu=0$ with
$\lambda>2,$ then $\Delta_{\mu\lambda}^{(s)}$ has a unique zero in
$(8,+\infty).$

  * (c)

If $\lambda\mu-4\lambda-2\mu>0$ and $\lambda>2,$ then
$\Delta_{\mu\lambda}^{(s)}$ has two zeros in $(8,+\infty).$

###### Proof.

Note that

$\lim\limits_{z\to+\infty}\Delta_{\mu\lambda}^{(s)}(z)=1$ (4.5)

and by Proposition 4.4

$\displaystyle\Delta_{\mu\lambda}^{(s)}(z)=$
$\displaystyle-\frac{\lambda\mu-4\lambda-2\mu}{8\pi}\,\ln(z-8)$
$\displaystyle+\Big{(}1-\frac{5\ln 2}{4\pi}\mu-\frac{5\ln
2-\pi}{2\pi}\lambda+\frac{5\ln 2-\pi}{8\pi}\lambda\mu\Big{)}+o(1)$

as $z\searrow 8$ so that

$\lim\limits_{z\searrow
8}\Delta_{\mu\lambda}^{(s)}(z)=\begin{cases}+\infty&\text{if
$\lambda\mu-4\lambda-2\mu>0,$}\\\ -\infty&\text{if
$\lambda\mu-4\lambda-2\mu<0,$}\\\ 1-\frac{\mu}{4}&\text{if
$\lambda\mu-4\lambda-2\mu=0.$}\end{cases}$ (4.6)

Notice that $\lambda\mu-4\lambda-2\mu=0$ and $\lambda>2,$ then $1-\mu/4<0.$

(a) We observe that if $\lambda\mu-4\lambda-2\mu\geq 0$ and $\lambda<2,$ then
either $\lambda\leq 0$ or $\mu\leq 0.$ Hence, by the minmax principle,
$H_{\mu\lambda}^{s}$ can have at most one eigenvalue above the essential
spectrum. By Lemma 4.5 $\Delta_{\mu\lambda}^{(s)}$ has at most one zero in
$(8,+\infty).$ Then (4.5) and (4.6) imply that
$\Delta_{\mu\lambda}^{(s)}(z)>0$ in $(8,+\infty),$ otherwise
$\Delta_{\mu\lambda}^{(s)}(\cdot)$ would cross $(8,+\infty)$ at least two
times.

(b) If $\lambda\mu-4\lambda-2\mu<0$ or $\lambda\mu-4\lambda-2\mu=0$ with
$\lambda>2,$ then by (4.5) and (4.6) the continuous function
$\Delta_{\mu\lambda}^{(s)}(\cdot)$ changes sign in $(8,+\infty)$ so that the
equation $\Delta_{\mu\lambda}^{(s)}(z)=0$ has at least one zero in
$(8,+\infty).$ If this equation had at least two zeros, then
$\Delta_{\mu\lambda}^{(s)}(\cdot)$ should have at least three zeros since
$\Delta_{\mu\lambda}^{(s)}(\cdot)$ has different signs at the endpoints of
$(8,+\infty).$ This contradicts to Lemma 4.5.

(c) If $\lambda\mu-4\lambda-2\mu>0$ and $\lambda>2,$ then $\mu+2\lambda>0,$
and hence, by Theorem 4.2 $H_{\mu\lambda}^{s}$ has at least one eigenvalue
above the essential spectrum. Then by Lemma 4.5
$\Delta_{\mu\lambda}^{(s)}(\cdot)$ has at least one zero in $(8,+\infty).$ On
the other hand, by (4.5) and (4.6) $\Delta_{\mu\lambda}^{(s)}(\cdot)$ has the
same signs at the endpoints of $(8,+\infty),$ hence, by continuity the
equation $\Delta_{\mu\lambda}^{(s)}(z)=0$ has at least two solutions. Now
Lemma 4.5 implies that $\Delta_{\mu\lambda}^{(s)}(\cdot)$ has two zeros in
$(8,+\infty).$ ∎

The zeros of $\Delta_{\mu\lambda}^{(s)}$ in $(-\infty,0)$ are studied in the
following lemma whose proof can be done along the lines of Lemma 4.7.

######  Lemma 4.8.

Let $(\mu,\lambda)\in\mathbb{R}^{2}.$

  * (a)

If $\lambda\mu+4\lambda+2\mu\geq 0$ and $\lambda>-2,$ then
$\Delta_{\mu\lambda}^{(s)}$ has no zeros in $(-\infty,0).$

  * (b)

If $\lambda\mu+4\lambda+2\mu<0$ or $\lambda\mu+4\lambda+2\mu=0$ with
$\lambda<-2,$ then $\Delta_{\mu\lambda}^{(s)}$ has a unique zero in
$(-\infty,0).$

  * (c)

If $\lambda\mu+4\lambda+2\mu>0$ and $\lambda<-2,$ then
$\Delta_{\mu\lambda}^{(s)}$ has two zeros in $(-\infty,0).$

###### Proof of Theorem 3.4.

The assertions of the theorem follow from Lemmas 4.1, 4.5, 4.6, 4.7 and 4.8. ∎

###  4.2. The discrete spectrum of $H_{\mu\lambda}(K)$

For every $n\geq 1$ define

$e_{n}(K;\mu,\lambda):=\sup\limits_{\phi_{1},\ldots,\phi_{n-1}\in
L^{2,e}(\mathbb{T}^{2})}\,\,\inf\limits_{\psi\in[\phi_{1},\ldots,\phi_{n-1}]^{\perp},\,\|\psi\|=1}(H_{\mu\lambda}(K)\psi,\psi)$

and

$E_{n}(K;\mu,\lambda):=\inf\limits_{\phi_{1},\ldots,\phi_{n-1}\in
L^{2,e}(\mathbb{T}^{2})}\,\,\sup\limits_{\psi\in[\phi_{1},\ldots,\phi_{n-1}]^{\perp},\,\|\psi\|=1}(H_{\mu\lambda}(K)\psi,\psi).$

By the minmax principle, $e_{n}(K;\mu,\lambda)\leq\mathcal{E}_{\min}(K)$ and
$E_{n}(K;\mu,\lambda)\geq\mathcal{E}_{\max}(K).$ Moreover, choosing
$\phi_{1}\equiv 1,$ $\phi_{2}(p)=\cos p_{1}$ and $\phi_{3}(p)=\cos p_{2}$ we
immediately see that $e_{n}(K;\mu,\lambda)=\mathcal{E}_{\min}(K)$ and
$E_{n}(K;\mu,\lambda)=\mathcal{E}_{\max}(K)$ for all $n\geq 4.$

######  Lemma 4.9.

Let $n\geq 1$ and $i\in\\{1,2\\}.$ Then the map

$K_{i}\in\mathbb{T}\mapsto\mathcal{E}_{\min}(K)-e_{n}(K;\mu,\lambda)$

is non-increasing in $(-\pi,0]$ and non-decreasing in $[0,\pi]$. Similarly,
the map

$K_{i}\in\mathbb{T}\mapsto E_{n}(K;\mu,\lambda)-\mathcal{E}_{\max}(K)$

is non-increasing in $(-\pi,0]$ and non-decreasing in $[0,\pi]$.

###### Proof.

Without loss of generality we assume that $i=1.$ Given $\psi\in
L^{2}(\mathbb{T}^{2})$ consider

$((H_{0}(K)-\mathcal{E}_{\min}(K))\psi,\psi)=\int_{\mathbb{T}^{2}}\sum\limits_{i=1}^{2}\cos\tfrac{K_{i}}{2}\,\big{(}1-\cos
q_{i}\big{)}|\psi(q)|^{2}\,\mathrm{d}q.$

Thus, the map
$K_{1}\in\mathbb{T}\mapsto((H_{0}(K)-\mathcal{E}_{\min}(K))\psi,\psi)$ is non-
decreasing in $(-\pi,0]$ and is non-increasing in $[0,\pi].$ Since
$V_{\mu\lambda}$ is independent of $K,$ from the definition of
$e_{n}(K;\mu,\lambda)$ the map $K_{1}\in\mathbb{T}\mapsto
e_{n}(K;\mu,\lambda)-\mathcal{E}_{\min}(K)$ is non-decreasing in $(-\pi,0]$
and is non-increasing in $[0,\pi].$

The case of $K_{i}\mapsto E_{n}(K;\mu,\lambda)-\mathcal{E}_{\max}(K)$ is
similar. ∎

###### Proof of Theorem 3.1.

From Lemma 4.9 for any $K\in\mathbb{T}^{2}$ and $m\geq 1$

$0\leq\mathcal{E}_{\min}(0)-e_{m}(0;\mu,\lambda)\leq\mathcal{E}_{\min}(K)-e_{m}(K;\mu,\lambda),$
(4.7)

and

$E_{m}(K;\mu,\lambda)-\color[rgb]{0,0,1}\mathcal{E}_{\max}(K)\color[rgb]{0,0,0}\geq
E_{m}(0;\mu,\lambda)-\color[rgb]{0,0,1}\mathcal{E}_{\max}(0)\color[rgb]{0,0,0}\geq
0.$

By assumption, $e_{n}(0;\mu,\lambda)$ is a discrete eigenvalue of
$H_{\mu\lambda}(0)$ for some $\mu,\lambda\in\mathbb{R}.$ Thus,
$\mathcal{E}_{\min}(0)-e_{n}(0;\mu,\lambda)>0,$ and hence, by (4.7) and (2.6)
$e_{n}(K;\mu,\lambda)$ is a discrete eigenvalue of $H_{\mu\lambda}(K)$ for any
$K\in\mathbb{T}^{2}.$ Since $e_{1}(K;\mu,\lambda)\leq\ldots\leq
e_{n}(K;\mu,\lambda)<\mathcal{E}_{\min}(K),$ it follows that
$H_{\mu\lambda}(K)$ has at least $n$ eigenvalue below the essential spectrum.

The case of $E_{n}(K;\mu,\lambda)$ is similar. ∎

###### Proof of Theorem 3.2.

Just combine Theorem 3.1 with Theorem 3.4. ∎

##  Appendix A Proof of Proposition 4.4

The general form of expansions (4.1) and (4.2) of $a,b,c,d$ can be established
following, for example, to [13]. The first two terms of these expansions,
i.e., the pairs $(P_{\omega}^{0}(0),Q_{\omega}^{0}(0))$ and
$(P_{\omega}^{1}(8),Q_{\omega}^{1}(8))$ are given by the following lemma.

######  Lemma A.1.

$\displaystyle\int_{\mathbb{T}^{2}}\frac{\mathrm{d}q_{1}\mathrm{d}q_{2}}{2-\cos
q_{1}-\cos q_{2}-z}$ $\displaystyle\hskip
34.1433pt=\begin{cases}-2\pi\ln(-z)+8\pi\ln 2+o(1)&\text{as $z\nearrow 0,$}\\\
2\pi\ln(z-4)-8\pi\ln 2+o(1)&\text{as $z\searrow 4,$}\end{cases}$ (A.1)
$\displaystyle\int_{\mathbb{T}^{2}}\frac{\cos q_{1}\cos
q_{2}\,\mathrm{d}q_{1}\mathrm{d}q_{2}}{2-\cos q_{1}-\cos q_{2}-z}$
$\displaystyle\hskip 34.1433pt=\begin{cases}-2\pi\ln(-z)+8\pi(\ln
2-1)+o(1)&\text{as $z\nearrow 0,$}\\\ 2\pi\ln(z-4)-8\pi(\ln 2-1)+o(1)&\text{as
$z\searrow 4.$}\end{cases}$ (A.2)
$\displaystyle\int_{\mathbb{T}^{2}}\frac{(\cos q_{1}+\cos
q_{2})\mathrm{d}q_{1}\mathrm{d}q_{2}}{2-\cos q_{1}-\cos q_{2}-z}$
$\displaystyle\hskip 34.1433pt=\begin{cases}-4\pi\ln(-z)+16\pi\ln
2-4\pi^{2}+o(1)&\text{as $z\nearrow 0,$}\\\ 4\pi\ln(z-4)-16\pi\ln
2+4\pi^{2}+o(1)&\text{as $z\searrow 4,$}\end{cases}$ (A.3)
$\displaystyle\int_{\mathbb{T}^{2}}\frac{(\cos q_{1}+\cos
q_{2})^{2}\mathrm{d}q_{1}\mathrm{d}q_{2}}{2-\cos q_{1}-\cos q_{2}-z}$
$\displaystyle\hskip 34.1433pt=\begin{cases}-8\pi\ln(-z)+32\pi\ln
2-8\pi^{2}+o(1)&\text{as $z\nearrow 0,$}\\\ 8\pi\ln(z-4)-32\pi\ln
2+8\pi^{2}+o(1)&\text{as $z\searrow 4,$}\end{cases}$ (A.4)
$\displaystyle\int_{\mathbb{T}^{2}}\frac{(\cos q_{1}-\cos
q_{2})^{2}\mathrm{d}q_{1}\mathrm{d}q_{2}}{2-\cos q_{1}-\cos
q_{2}-z}=\begin{cases}32\pi-8\pi^{2}+o(1)&\text{as $z\nearrow 0,$}\\\
-32\pi+8\pi^{2}+o(1)&\text{as $z\searrow 4.$}\end{cases}$ (A.5)

###### Proof.

Since the idea is the same, we prove (A.1)-(A.5) only for $z<0.$ By
simmetricity, it suffices to prove only (A.1) and (A.2). First we observe that
if $|u|>1,$ then

$\int_{\mathbb{T}}\frac{\mathrm{d}t}{u-\cos
t}=\begin{cases}\frac{-2\pi}{\sqrt{u^{2}-1}}&\text{if $u<-1,$}\\\
\frac{2\pi}{\sqrt{u^{2}-1}}&\text{if $u>1.$}\end{cases}$ (A.6)

Let us prove (A.1). Since $2-z-\cos q_{2}>1$ for any
$z\in\mathbb{R}\setminus[0,4]$ and $q_{2}\in\mathbb{T},$ by (A.6)

$\int_{\mathbb{T}^{2}}\frac{\mathrm{d}q_{1}\mathrm{d}q_{2}}{2-\cos q_{1}-\cos
q_{2}-z}=\int_{\mathbb{T}}\frac{\mathrm{d}q_{2}}{\sqrt{(2-z-\cos
q_{2})^{2}-1}}.$

Then using the change of variables
$q_{2}:=2\arctan\big{(}\frac{z-2}{z-4}\big{)}^{1/2}v$ we get

$\displaystyle\int_{\mathbb{T}^{2}}\frac{\mathrm{d}q_{1}\mathrm{d}q_{2}}{2-\cos
q_{1}-\cos q_{2}-z}=$
$\displaystyle\frac{8\pi}{|2-z|}\int_{0}^{\infty}\frac{\mathrm{d}v}{\sqrt{\big{(}\tfrac{z^{2}-4z}{(z-2)^{2}}+v^{2}\big{)}\big{(}1+v^{2}\big{)}}}$
$\displaystyle=$ $\displaystyle J_{1}(z)+J_{2}(z)+J_{3}(z),$

where

$\displaystyle J_{1}(z):=$
$\displaystyle\frac{8\pi}{|2-z|}\int_{0}^{1}\frac{\mathrm{d}v}{\sqrt{\tfrac{z^{2}-4z}{(z-2)^{2}}+v^{2}}},$
$\displaystyle J_{2}(z):=$
$\displaystyle\frac{8\pi}{|2-z|}\int_{0}^{1}\frac{\big{(}(1+v^{2})^{-1/2}-1\big{)}\,\mathrm{d}v}{\sqrt{\tfrac{z^{2}-4z}{(z-2)^{2}}+v^{2}}},$
$\displaystyle J_{3}(z):=$
$\displaystyle\frac{8\pi}{|2-z|}\int_{1}^{\infty}\frac{\mathrm{d}v}{\sqrt{\big{(}\tfrac{z^{2}-4z}{(z-2)^{2}}+v^{2}\big{)}\big{(}1+v^{2}\big{)}}}.$

Note that

$J_{1}(z)=\frac{8\pi}{|2-z|}\Big{(}\ln\Big{[}1+\sqrt{1+\tfrac{z^{2}-4z}{(z-2)^{2}}}\Big{]}-\ln\tfrac{\sqrt{z^{2}-4z}}{|z-2|}\Big{)}=-2\pi\ln(-z)+4\pi\ln
2+o(1)$

as $z\nearrow 0.$ Moreover, $J_{2}$ and $J_{3}$ are continuous at $z=0$ and
$J_{2}(0)=4\pi\ln\frac{2}{1+\sqrt{2}}$ and $J_{3}(0)=4\pi\ln(1+\sqrt{2}).$
Thus, (A.1) follows.

Now we prove (A.2). Using $\int_{\mathbb{T}}\cos q_{2}\mathrm{d}q_{2}=0$ and
(A.6)

$\displaystyle\int_{\mathbb{T}^{2}}\frac{\cos q_{1}\cos
q_{2}\mathrm{d}q_{1}\mathrm{d}q_{2}}{2-\cos q_{1}-\cos q_{2}-z}=$
$\displaystyle\int_{\mathbb{T}}\frac{(2-z-\cos q_{2})\cos
q_{2}\,\mathrm{d}q_{2}}{\sqrt{(2-z-\cos q_{2})^{2}-1}}.$

Thus changing variables as $q_{2}=2\arctan u$ we get

$\int_{\mathbb{T}^{2}}\frac{\cos q_{1}\cos
q_{2}\mathrm{d}q_{1}\mathrm{d}q_{2}}{2-\cos q_{1}-\cos
q_{2}-z}=I_{1}(z)+I_{2}(z),$

where

$\displaystyle I_{1}(z):=$ $\displaystyle
8\pi(1-z)\int_{0}^{\infty}\frac{\mathrm{d}u}{(1+u^{2})^{2}\sqrt{-z+(2-z)u^{2}}\sqrt{2-z+(4-z)u^{2}}}$
$\displaystyle I_{2}(z):=$ $\displaystyle
8\pi\int_{0}^{\infty}\frac{\big{(}2u^{2}-(3-z)u^{4}\big{)}\,\mathrm{d}u}{(1+u^{2})^{2}\sqrt{-z+(2-z)u^{2}}\sqrt{2-z+(4-z)u^{2}}}.$

Obviously, $I_{2}$ is continuous at $z=0$ and $I_{2}(0)=2\pi(\pi-5).$ We
represent $I_{1}(z)$ as

$I_{1}(z)=I_{11}(z)+I_{12}(z)+I_{13}(z),$

where

$\displaystyle I_{11}(z):=$
$\displaystyle\frac{8\pi(1-z)}{\sqrt{2-z}}\int_{0}^{1}\frac{\mathrm{d}u}{\sqrt{-z+(2-z)u^{2}}}$
$\displaystyle I_{12}(z):=$ $\displaystyle
8\pi(1-z)\int_{0}^{1}\frac{\big{[}(1+u^{2})^{-2}(2-z+(4-z)u^{2})^{-1/2}-(2-z)^{-1/2}\big{]}\,\mathrm{d}u}{\sqrt{-z+(2-z)u^{2}}}$
$\displaystyle I_{13}(z):=$ $\displaystyle
8\pi(1-z)\int_{1}^{\infty}\frac{\mathrm{d}u}{(1+u^{2})^{2}\sqrt{-z+(2-z)u^{2}}\sqrt{2-z+(4-z)u^{2}}}.$

Then

$\displaystyle I_{11}(z)=$
$\displaystyle\frac{8\pi(1-z)}{2-z}\Big{[}\ln\Big{(}1+\sqrt{1+\frac{z}{z-2}}\Big{)}-\ln\sqrt{\frac{z}{z-2}}\Big{]}$
$\displaystyle=$ $\displaystyle-2\pi\ln(-z)+6\pi\ln 2+o(1)$

as $z\nearrow 0.$ Note that the functions $I_{12}(z)$ and $I_{13}(z)$ are
continuous at $z=0$ and

$I_{12}(0)=4\pi\Big{[}\ln\frac{2}{1+\sqrt{3}}-\frac{\pi}{6}+\frac{2-\sqrt{3}}{4}\Big{]}$

and

$I_{13}(0)=4\pi\Big{[}\ln\frac{\sqrt{3}+1}{\sqrt{2}}-\frac{\pi}{3}+\frac{\sqrt{3}}{4}\Big{]}.$

Hence

$I_{1}(z)=-2\pi\ln(-z)+8\pi\ln 2-2\pi^{2}+2\pi+o(1)$

and (A.2) follows. ∎

### Acknowledgment

The first and third author acknowledge support from the Foundation for Basic
Research of the Republic of Uzbekistan (Grant No. OT-F4-66). The second author
acknowledges support from the Austrian Science Fund (FWF) project M 2571-N32.

## References

  * [1] S. Albeverio, S.N. Lakaev, K.A. Makarov, Z.I. Muminov: The threshold effects for the two-particle Hamiltonians on lattices. Commun. Math. Phys. 262 (2006), 91–115.
  * [2] I. Bloch: Ultracold quantum gases in optical lattices. Nat. Phys. 1 (2005), 23–30.
  * [3] P.A. Faria Da Veiga, L. Ioriatti, M. O’Carroll: Energy-momentum spectrum of some two-particle lattice Schrödinger Hamiltonians. Phys. Rev. E 66 (2002), 016130.
  * [4] M.-S. Heo, T.T. Wang, C.A. Christensen, T.M. Rvachov, D.A. Cotta, J.-H. Choi, Y.-R. Lee, W. Ketterle: Formation of ultracold fermionic NaLi Feshbach molecules. Phys. Rev. A 86 (2012), 021602.
  * [5] F. Hiroshima, Z. Muminov, U. Kuljanov: Threshold of discrete Schrödinger operators with delta-potentials on $N$-dimensional lattice. Linear and Multilinear Algebra (2020).
  * [6] D. Jaksch, C. Bruder, J. Cirac, C.W. Gardiner, P. Zoller: Cold bosonic atoms in optical lattices. Phys. Rev. Lett. 81 (1998), 3108–3111.
  * [7] D. Jaksch, P. Zoller: The cold atom Hubbard toolbox. Ann. Phys. 315 (2005), 52–79.
  * [8] Sh.Yu. Kholmatov, S.N. Lakaev, F. Almuratov: Bound states of discrete Schrödinger operators on one and two dimensional lattices. Submitted, arXiv:2007.04035 (2020).
  * [9] M. Klaus: On the bound state of Schrödinger operators in one dimension. Ann. Phys. 108 (1977), 288-300.
  * [10] M. Klaus, B. Simon: Coupling constant thresholds in non-relativistic quantum mechanics. I. Short range two body case. Ann. Phys. 130 (1980), 251–281.
  * [11] S. N. Lakaev, S. Kh. Abdukhakimov: Threshold effects in a two-fermion system on an optical lattice. Theor. Math. Phys. 203 (2020), 251–268.
  * [12] S.N. Lakaev, I.N. Bozorov: The number of bound states of a one-particle Hamiltonian on a three-dimensional lattice. Theoret. and Math. Phys. 158 (2009), 360–376.
  * [13] S.N. Lakaev, A.M. Khalkhuzhaev, Sh.S. Lakaev: Asymptotic behavior of an eigenvalue of the two-particle discrete Schrödinger operator. Theoret. Math. Phys. 171 (2012), 800–811.
  * [14] S.N. Lakaev, Sh.Yu. Kholmatov: Asymptotics of the eigenvalues of a discrete Schrödinger operator with zero-range potential. Izvestiya Math. 76 (2012), 946–966.
  * [15] S.N. Lakaev, Sh.Yu. Kholmatov: Asymptotics of eigenvalues of two-particle Schrödinger operators on lattices with zero-range interaction. J. Phys. A: Math. Theor. 44 (2011).
  * [16] S.N. Lakaev, E. Özdemir: The existence and location of eigenvalues of the one particle Hamiltonians on lattices. Hacettepe J. Math. Stat. 45 (2016), 1693–1703.
  * [17] M. Lewenstein, A. Sanpera, V. Ahufinger: Ultracold Atoms in Optical Lattices: Simulating Quantum Many-body Systems. Oxford University Press, Oxford, 2012.
  * [18] D. Mattis: The few-body problem on a lattice. Rev. Mod. Phys. 58 (1986), 361–379.
  * [19] A. Mogilner: Hamiltonians in solid-state physics as multiparticle discrete Schrödinger operators: Problems and results. Advances in Societ Math. 5 (1991), 139–194.
  * [20] C. Ospelkaus, S. Ospelkaus, L. Humbert, P. Ernst, K. Sengstock, K. Bongs: Ultracold heteronuclear molecules in a 3d optical lattice. Phys.Rev. Lett. 97 (2006).
  * [21] B. Simon: Notes on infinite determinants of Hilbert space operators. Advances in Math. 24 (1977), 244-273.
  * [22] B. Simon: The bound state of weakly coupled Schrödinger operators in one and two dimensions. Ann. Phys. 97 (1976), 279–288.
  * [23] K. Winkler, G. Thalhammer, F. Lang, R. Grimm, J. Hecker Denschlag, A.J. Daley, A. Kantian, H.P. Büchler, P. Zoller: Repulsively bound atom pairs in an optical lattice. Nature 441 (2006), 853–856.
  * [24] J.J. Zirbel, K.-K. Ni, S. Ospelkaus, T.L. Nicholson, M.L. Olsen, P.S. Julienne, C.E. Wieman, J. Ye, D.S. Jin: Heteronuclear molecules in an optical dipole trap. Phys. Rev. A 78 (2008).


# Stabilizer Approximation

Xinying Li  Jianan Wang  Chuixiong Wu  Fen Zuo111Email<EMAIL_ADDRESS>Shanghai MiQro Era Digital Technology Ltd, Shanghai, China

###### Abstract

We propose a heuristic method to obtain the approximate groundstate for a
Hamiltonian in the qubit form, based on the stabilizer formalism. These states
may serve as proper initial states for further refined computation. It would
be interesting to assess the efficiency and scalability of the method.

###### Contents

  1. I INTRODUCTION
  2. II THE METHOD
    1. II.1 H2
    2. II.2 LiH
    3. II.3 BeH2
  3. III DISCUSSION

##  I INTRODUCTION

Recently the stabilizer formalism Gottesman96 ; Gottesman97 ; Gottesman98  is
utilized to obtain approximate groundstates of chemical and manybody physical
systems Shang2021 ; CAFQA ; CAFQA2 . However, these studies focus on the
normalizers, namely the Clifford circuits, rather than the stabilizers
themselves, and thus are not as efficient as expected Gottesman98 . To fully
capture the essence of the stabilizer formalism, it would be better to work
with the stabilizers directly. For the chemical and manybody problems, they
could be engineered naturally from the qubit form of the Hamiltonians. In the
following we will illustrate such a procedure using some examples. In some
sense, our method could be considered as a follow-up to the qubit-tapering
algorithm Tapering .

##  II THE METHOD

We will take the H2, LiH and BeH2 molecules as examples, and work in the
STO-3G basis. The distance unit is chosen to be Å, and energy unit hartree.
The chemical framework PySCF PySCF  is utilized to generate the second-
quantized fermionic Hamitonians, and do the restricted/unrestricted Hartree-
Fock calculations for comparison. The open source package Qiskit Qiskit  is
then employed to transform the fermonic Hamiltonians to the qubit form. We
choose the parity transformation Parity , but the other transformations work
as well. The nuclear-nuclear potentials are not included, and will be added to
the final results by hand.

###  II.1 H2

Using two-qubit reduction technique Tapering , the H2 molecule can be
described with 2 qubits. Let us start from the equilibrium point $d=0.74$,
where we know the Hartree-Fock solution approximates the exact result well.
The Hamiltonian reads

$\displaystyle H_{\text{H}_{2}}$ $\displaystyle=$
$\displaystyle-1.0534210769165204*II$ $\displaystyle+0.39484436335590356*IZ$
$\displaystyle-0.39484436335590367*ZI$ $\displaystyle+0.1812104620151969*XX$
$\displaystyle-0.011246157150821112*ZZ.$

To approach the groundstate as close as possible, the stabilizers should be
chosen to make Hamiltonian as small as possible. It is straightforward to
figure out the best choice is

$Z_{1},\quad-Z_{0}.$

Here the qubits are labeled from right to left as $0,1,2...$. And for $M/2$
spatial orbitals, the first $M/2$ qubits correspond to the spin-up spin
orbitals, and the latter $M/2$ qubits correspond to the spin-down ones. So the
stabilizer state is fixed to be $|0_{1}1_{0}\rangle$, which is simply the
parity transformation of the corresponding Hartree-Fock state. The
corresponding total energy is about $-1.12$, while the exact value is $-1.14$.

Now we proceed to study the region away from the equilibrium point. Taking
$d=2.8$, the Hamiltonian reads:

$\displaystyle H^{\prime}_{\text{H}_{2}}$ $\displaystyle=$
$\displaystyle-0.8284676561247681*II$ $\displaystyle+0.2930431286727852*XX$
$\displaystyle+0.016170000066607376*IZ$
$\displaystyle-0.016170000066607328*ZI$
$\displaystyle-0.0001469354633982234*ZZ.$

It is straightforward to choose the first stabilizer as $-X_{1}X_{0}$. Then
single $Z$ operators are not allowed to be stabilizers, and we have to choose
$Z_{1}Z_{0}$. The stabilizer state is then

$\frac{|0_{1}0_{0}\rangle-|1_{1}1_{0}\rangle}{\sqrt{2}}.$ (1)

The corresponding total energy is about $-0.93$, nearly the same as the exact
result. However, the coefficient of the $Z_{1}Z_{0}$ term is very close to
zero, which indicates a possible two-fold degeneracy. If so, $-Z_{1}Z_{0}$
could be an alternative choice for the second stabilizer, and the
corresponding state is

$\frac{|0_{1}1_{0}\rangle-|1_{1}0_{0}\rangle}{\sqrt{2}}.$ (2)

The splitting of the two states we see here could be an artificial effect of
the qubit reduction procedure. Also we find that the particle number
explanation is a little tricky in such a situation.

So let us go back to the original $4$-qubit form. The Hamiltonian at $d=2.8$
is given by

$\displaystyle H^{\prime\prime}_{\text{H}_{2}}$ $\displaystyle=$
$\displaystyle-0.7340910665455848*IIII$
$\displaystyle+0.12194654795642478*ZZZZ$
$\displaystyle+0.12044907695778825*ZZIZ$
$\displaystyle+0.12044907695778825*IZIZ$
$\displaystyle+0.11909854142254998*IZZZ$
$\displaystyle+0.07326078216819633*ZXIX$
$\displaystyle-0.07326078216819633*IXZX$
$\displaystyle-0.07326078216819633*ZXZX$
$\displaystyle+0.07326078216819633*IXIX$
$\displaystyle+0.04718829478959189*IIZI$
$\displaystyle+0.04718829478959189*ZIZI$
$\displaystyle+0.047527808451266154*IIIZ$
$\displaystyle+0.04752780845126615*IZZI$
$\displaystyle+0.03135780838465877*IIZZ$
$\displaystyle+0.03135780838465882*ZZII.$

After some work one will find out the first three stabilizers could be chosen
as

$-Z_{2}Z_{0},Z_{1},Z_{3}.$

But then the Hamiltonian gives no more information for us to fix the last
stabilizer. So we conclude that the stabilizer groundstate is indeed 2-fold
degenerate, and could be taken to be

$\frac{|0_{3}0_{2}0_{1}1_{0}\rangle+|0_{3}1_{2}0_{1}0_{0}\rangle}{\sqrt{2}},\quad\frac{|0_{3}0_{2}0_{1}1_{0}\rangle-|0_{3}1_{2}0_{1}0_{0}\rangle}{\sqrt{2}}.$
(3)

This means the last stabilizer has been chosen as $X_{2}X_{0}$ and
$-X_{2}X_{0}$, respectively. However, the degeneracy actually destroys the
entanglement structure in the above states, making them equivalent to product
states, or Hartree-Fock states. A unrestricted Hartree-Fock calculation
confirms such a conclusion.

One could sample more points on the potential energy curve. It is not
difficult to show that the potential curve at small distance could be well
approximated by the Hartree-Fock/product states, while at large distance is
well approximated by the entangling states (1,2,3). This is in accordance with
the numerical results in CAFQA .

###  II.2 LiH

To avoid the degeneracy of H2, we go ahead to study the asymmetric LiH
molecule. One can estimate the required qubits by counting the atomic
orbitals. We take the $1s,2s,2p$ orbitals of Li and the $1s$ orbital of H into
account, so we have in total $6$ spatial orbitals, or $12$ spin orbitals.
After the two-qubit reduction, we would get a $10$-qubit Hamiltonian, which
contains as many as $631$ Pauli terms.

We start from the point $d=1.5$ near the equilibrium. The Hamiltonian reads:

$\displaystyle H_{\text{LiH}}$ $\displaystyle=$
$\displaystyle-5.161946304557396*IIIIIIIIII$
$\displaystyle+1.010986985827016*IIIIIIIIIZ$
$\displaystyle+1.0109869858270164*IIIIZIIIII$
$\displaystyle+0.4145416937880457*IIIIZIIIIZ$
$\displaystyle-0.3982629120981239*IIIIIZIIII$
$\displaystyle-0.39826291209812403*ZIIIIIIIII$
$\displaystyle-0.22878247783636188*ZZIIIIIIII$
$\displaystyle-0.2287824778363618*IZZIIIIIII$
$\displaystyle-0.22878247783636196*IIIIIZZIII$
$\displaystyle-0.2287824778363618*IIIIIIZZII$
$\displaystyle-0.19731447155414064*IIIIIIIZZI$
$\displaystyle-0.19731447155414092*IIZZIIIIII$
$\displaystyle-0.11495682604662383*IIIIIIIIZZ$
$\displaystyle-0.11495682604662373*IIIZZIIIII$ $\displaystyle+...$

From these dominant terms one can easily fix the stabilizers to be

$Z_{9},Z_{8},Z_{7},Z_{6},-Z_{5},Z_{4},Z_{3},Z_{2},Z_{1},-Z_{0}.$

So the stabilizer state is the product state

$|0_{9}0_{8}0_{7}0_{6}1_{5}0_{4}0_{3}0_{2}0_{1}1_{0}\rangle,$ (4)

or Hartree-Fock state, as expected. The total energy can be estimated to be
$-7.86$, close to the exact value $-7.88$.

To see the entangling structure in the faraway region, we generate the
Hamiltonian at $d=5.0$, which is given by :

$\displaystyle H^{\prime}_{\text{LiH}}$ $\displaystyle=$
$\displaystyle-4.712365947275493*IIIIIIIIII$
$\displaystyle+1.0043414021528052*IIIIIIIIIZ$
$\displaystyle+1.0043414021528052*IIIIZIIIII$
$\displaystyle+0.415119697523212*IIIIZIIIIZ$
$\displaystyle-0.23366481986596774*IIIIIZIIII$
$\displaystyle-0.2336648198659678*ZIIIIIIIII$
$\displaystyle-0.23360151362602502*IIIIIIZZII$
$\displaystyle-0.23360151362602524*IIIIIZZIII$
$\displaystyle-0.23360151362602505*IZZIIIIIII$
$\displaystyle-0.23360151362602527*ZZIIIIIIII$
$\displaystyle-0.12132590572285508*IIZZIIIIII$
$\displaystyle-0.12132590572285518*IIIIIIIZZI$
$\displaystyle-0.11613981631234098*IIIIIIIIZZ$
$\displaystyle-0.11613981631234092*IIIZZIIIII$ $\displaystyle+...$

From these dominant terms one can already get $8$ stabilizers:

$Z_{9},Z_{8},Z_{7},-Z_{5},Z_{4},Z_{3},Z_{2},-Z_{0}.$

A detailed analysis of some subdominant terms suggests the remaining two
stabilizers be chosen as $-Z_{6}Z_{1}$ and $-X_{6}X_{1}$. So the resulting
state is the entangling state:

$\frac{|0_{9}0_{8}0_{7}1_{6}1_{5}0_{4}0_{3}0_{2}0_{1}1_{0}\rangle-|0_{9}0_{8}0_{7}0_{6}1_{5}0_{4}0_{3}0_{2}1_{1}1_{0}\rangle}{\sqrt{2}}.$
(5)

The corresponding total energy is roughly $-7.76$, again close to the exact
value $-7.78$. Now we do not have any degeneracy any more, and the entangling
structure do not collapse. So in principle, the above entangling state can not
be directly obtained through a Hartree-Fock procedure, which only samples the
product states. However, by treating the spin-up and spin-down sectors
separately as in the unrestricted Hartree-Fock approach, one could somehow
imitate the entangling effect. Numerical calculation shows that unrestricted
Hartree-Fock approach gives similar results for the total energy in this
region.

Same as the H2 molecule, one can check that the product state (4) and the
entangling state (5) together give the best stabilizer approximation for the
potential energy curve in the whole region. Again this is in accordance with
the numerical results in CAFQA .

###  II.3 BeH2

For molecules with more than two atoms, one would expect more dominant
stabilizer configurations, and configurations with complicated entangling
structure. So we apply the method to the BeH2 molecule. Near the equilibrium
point the atoms are on a line, with Be at the middle. We take the Be-H
distance $d$ as the varying parameter, consider the equilibrium point
$d=1.32$, and the faraway point $d=5.0$. The procedure is nearly the same as
LiH, so we simply list the results.

At $d=1.32$, the stabilizers are:

$Z_{11},Z_{10},Z_{9},Z_{8},-Z_{7},Z_{6},-Z_{5},-Z_{4},-Z_{3},-Z_{2},Z_{1},-Z_{0}.$

The corresponding state is

$|0_{11}0_{10}0_{9}0_{8}1_{7}0_{6}1_{5}1_{4}1_{3}1_{2}0_{1}1_{0}\rangle.$ (6)

The total energy is about $-15.56$, while the exact energy is $-15.60$.

At $d=5.0$, the first $11$ stabilizers can be easily chosen as:

$Z_{11},Z_{10},Z_{9},-Z_{7},Z_{6},-Z_{5},-Z_{4},-Z_{3},Z_{1},-Z_{0},-X_{8}X_{2}.$

Again here the Hamiltonian gives no more information for the last stabilizer.
This is because, just as H2, BeH2 is also symmetric, and the groundstate is
doubly degenerate. For simplicity, let us choose the last stabilizer to be
$Z_{8}Z_{2}$, or $-Z_{8}Z_{2}$. The corresponding stabilizer states are

$\displaystyle\frac{|0_{11}0_{10}0_{9}0_{8}1_{7}0_{6}1_{5}1_{4}1_{3}0_{2}0_{1}1_{0}\rangle-|0_{11}0_{10}0_{9}1_{8}1_{7}0_{6}1_{5}1_{4}1_{3}1_{2}0_{1}1_{0}\rangle}{\sqrt{2}},$
(7)
$\displaystyle\frac{|0_{11}0_{10}0_{9}0_{8}1_{7}0_{6}1_{5}1_{4}1_{3}1_{2}0_{1}1_{0}\rangle-|0_{11}0_{10}0_{9}1_{8}1_{7}0_{6}1_{5}1_{4}1_{3}0_{2}0_{1}1_{0}\rangle}{\sqrt{2}}.$
(8)

Notice the similarity with (1) and (2). Again the particle number explanation
is ambiguous. The total energy is about $-15.28$, to be compared with the
exact value $-15.34$. The unrestricted Hartree-Fock procedure gives an energy
value -$15.22$, which is close to but not exactly the same as ours.

Unfortunately, BeH2 does not give rise to any new structure, probably due to
the symmetry of the system.

##  III DISCUSSION

As a generation of the traditional Hartree-Fock states, stabilizer states
could approximate the true groundstates reasonably well. We find that the
corresponding stabilizers could be naturally engineered from the qubit-form
Hamiltonians. Importantly, a specific stabilizer configuration seems to
dominate over a large region, making the calculation of the potential energy
surface much easier. For the simple molecules analyzed here, the stabilizer
configuration exhibits an interesting pattern, a transition from the product
state/Hartree-Fock state near the equilibrium to the entangling state in the
faraway region.

However, since the method is heuristic, it is not clear whether the method
will always be efficient. Also it is not verified if it performs well in much
larger systems.

We will continue to explore the properties of the method to see if it
qualifies as a satisfactory initial-state preparation procedure, which would
be important for further refined computation GKC2208 .

## CODE AVAILABILITY

A preliminary version of the algorithm is developed at
https://github.com/MiqroEra/Stabilizer.

## References

  * (1) Daniel Gottesman, 1996,  Class of quantum error-correcting codes saturating the quantum Hamming bound, Phys. Rev. A54 (1996) 1862, [https://arxiv.org/abs/quant-ph/9604038].
  * (2) Daniel Gottesman, 1997,  Stabilizer Codes and Quantum Error Correction, Ph.D. thesis, California Institute of Technology, Pasadena, CA, 1997, [https://arxiv.org/abs/quant-ph/9705052].
  * (3) Daniel Gottesman, 1998,  The Heisenberg Representation of Quantum Computers, [https://arxiv.org/abs/quant-ph/9807006].
  * (4) Zhong-Xia Shang, Ming-Cheng Chen, Xiao Yuan, Chao-Yang Lu and Jian-Wei Pan,  Schrödinger-Heisenberg Variational Quantum Algorithms, [arXiv: 2112.07881].
  * (5) Gokul Subramanian Ravi, Pranav Gokhale, Yi Ding, William M. Kirby, Kaitlin N. Smith, Jonathan M. Baker, Peter J. Love, Henry Hoffmann, Kenneth R. Brown, and Frederic T. Chong, 2022,  CAFQA: A classical simulation bootstrap for variational quantum algorithms, In Proceedings of 28th Conference on Architectural Support for Programming Languages and Operating Systems (ASPLOS). ACM, New York, NY, USA, 15 pages. [https://doi.org/10.1145/nnnnnnn.nnnnnnn], [arXiv: 2202.12924].
  * (6) M. H. Cheng, K. E. Khosla, C. N. Self, M. Lin, B. X. Li, A. C. Medina, M. S. Kim, 2022,  Clifford Circuit Initialisation for Variational Quantum Algorithms, [arXiv: 2207.01539].
  * (7) Sergey Bravyi, Jay M. Gambetta, Antonio Mezzacapo, Kristan Temme, 2017,  Tapering off qubits to simulate fermionic Hamiltonians, [https://arxiv.org/abs/1701.08213].
  * (8) Qiming Sun, Timothy C. Berkelbach, Nick S. Blunt, George H. Booth, Sheng Guo, Zhendong Li, Junzi Liu, James McClain, Elvira R. Sayfutyarova, Sandeep Sharma, Sebastian Wouters, and Garnet Kin-Lic Chan. 2017,  The Python-based Simulations of Chemistry Framework (PySCF). [arXiv: 1701.08223].
  * (9) Jacob T. Seeley, Martin J. Richard, Peter J. Love, 2012,  The Bravyi-Kitaev transformation for quantum computation of electronic structure, The Journal of Chemical Physics 137 (22), 224109.
  * (10) Héctor Abraham et al., 2019,  Qiskit: An Open-source Framework for Quantum Computing, [https://doi.org/10.5281/zenodo.2562110].
  * (11) Seunghoon Lee, Joonho Lee, Huanchen Zhai, Yu Tong, Alexander M. Dalzell, Ashutosh Kumar, Phillip Helms, Johnnie Gray, Zhi-Hao Cui, Wenyuan Liu, Michael Kastoryano, Ryan Babbush, John Preskill, David R. Reichman, Earl T. Campbell, Edward F. Valeev, Lin Lin, Garnet Kin-Lic Chan, 2022, Is there evidence for exponential quantum advantage in quantum chemistry? [arXiv: 2208.02199].


0pt0pt

# The Landau Collision Integral in the Particle Basis in the PETSc Library

Joseph V. Pusztay, Filippo Zonta, Matthew G. Knepley, Mark F. Adams

###### Abstract

The Landau collision integral is often considered the gold standard in the
context of kinetic plasma simulation due to its conservative properties,
despite challenges involved in its discretization. The primary challenge when
implementing an efficient computation of this operator is conserving physical
properties of the continuum equation when the system is discretized. Recent
work has achieved continuum discretizations using the method of Finite
Elements which maintain conservation of mass, momentum, and energy, but which
lacks monotonic entropy production. More recently, a particle discretization
has been introduced which conserves mass, momentum, and energy, but maintains
the benefit of monotonic entropy production necessary for the metriplecticity
of the system. We present here an implementation of the particle basis Landau
collision integral in the Portable Extensible Toolkit for Scientific Computing
in 2 and 3V for the construction of a full geometry solver with a novel
approach to computation of the entropy functional gradients. Verification of
the operator is achieved with thermal equilibration and isotropization tests.
All examples are available, open source, in the PETSc repository for
reproduction.

##  1 Introduction

The metriplectic representation of a kinetic plasma [1, 2] poses the system in
terms of two brackets; the Poisson bracket where Hamiltonian dynamics may be
evolved symplecticly, and the metric bracket which evolves the system
according to collisional (dissipitive) effects. This system, by construction,
has favorable mathematical properties such as conservation of mass, momentum,
and energy of the modified Hamiltonian in the symplectic structure, as well as
the monotonic production of entropy in the dissipitive bracket. Discretization
of the continuum in the Poisson bracket is relatively straightforward and has
been well studied [1] with conservative time stepping schemes well understood
[3, 4], where conservative discretizations of the metric bracket require
successful discretization of some collision operator. The Landau collision
integral [5] is often considered the gold standard in this regard.

Extensive efforts have been expended on the Landau collision operator, casting
it in terms of the metric bracket [6, 7] for continuum and particle
representations. Recent work has discretized the operator in its finite
element representation whilst maintaining conservation [8], which was
implemented in PETSc [9, 10, 11], and demonstrates decent performance at
exascale [12, 13]. For particle codes, this demands a projection operation
between particle and grid representations, and then a projection of the
evolved distribution function weights back to the particles. Suitable,
scalable operators exist within PETSc [14] for this purpose, however there are
caveats to this method. Firstly, the necessity of projection to a grid in
codes that are purely particle based. The second being the lack of monotonic
entropy production in the finite element discretization. A natural remedy to
the first is to construct a collision operator which exists purely in the
particle basis, which recent theoretical advancements have made possible in a
manner that alleviates the second point of concern.

Carillo et. al [15] recently proved the Maxwellian steady state of a
regularized particle basis representation of the Landau collision integral, as
well as the conservation of mass, momentum, and energy. This discretization
however suffered from lack of conservation in energy without a suitable time
discretization. A solution to this problem was presented by Hirvijoki in [16],
with a simple conservative time discretization alongside a discrete gradient
time stepping scheme which are conservative in mass, momentum, and energy, as
well as monotonic in entropy production. Tests of these conservative
properties were additionally performed in [17] with a GPU implementation of
this operator in 2 dimensions. These efforts form the necessary foundation to
develop a rigorous community implementation of the particle basis Landau
collision integral for the construction of, or implementation into, general
Particle-in-Cell (PIC) codes with support for full geometry solvers, with
potential for an implementation of the Discrete-Gradient operator.

We present here an implementation of the particle basis Landau collision
integral with the time discretization presented in [16] with 2 particle
species of realistic mass ratios in 2 and 3 dimensions. To begin, we will
reiterate the metriplectic formulation of the Vlasov-Poisson-Landau system
alongside the new particle basis approach to its evaluation. We present a
novel approach to evaluation of the entropy functionals via high order
quadrature on implicitly defined domains[18] and compare the quality of such
methods against the standard Gauss-Hermite quadratures used in [17].
Verification of the methods are demonstrated via basic tests such as thermal
equilibration and isotropization of multiple species.

##  2 Landau Collisions in the Particle Basis

Consider the evolution of a distribution function, $f$, with respect to time.
In a system of charged particles with self consistent and external forces,
this reads as

$\displaystyle\frac{\partial f}{\partial t}$ $\displaystyle=\frac{\partial
f}{\partial t}+\frac{\partial\vec{x}}{\partial
t}\cdot\nabla_{x}f+\frac{\partial\vec{v}}{\partial t}\cdot\nabla_{v}f$
$\displaystyle=\frac{\partial f}{\partial
t}+\vec{v}\cdot\nabla_{x}f+\frac{e}{m}(\vec{E}+\vec{v}\crossproduct\vec{B})\cdot\nabla_{v}f$
$\displaystyle=C$ (1)

which is the Vlasov-Maxwell-Boltzmann model with velocity $\vec{v}$, position
$\vec{x}$, electric field $\vec{E}$, magnetic field $\vec{B}$, mass $m$,
charge $e$, and diffusive collisional term C. This system contains two primary
constituents; namely the Poisson bracket corresponding to the electromagnetic
terms with $\frac{\partial f}{\partial t}=0$, and the diffusive metric bracket
incorporating collisional terms as

$\displaystyle\frac{\partial f}{\partial t}$
$\displaystyle=\frac{e}{m}\vec{E}\cdot\nabla f$ $\displaystyle=C$ (2)

where C is expressed in the continuum as the Landau collision integral

$\displaystyle
C=\nu_{\alpha\beta}\frac{m_{0}}{m_{\alpha}}\int_{\Omega}U(\vec{v},\bar{v})\cdot(\frac{m_{0}}{m_{\alpha}}\bar{f}_{\beta}\nabla
f_{\alpha}-\frac{m_{0}}{m_{\beta}}f_{\alpha}\bar{\nabla}\bar{f}_{\beta})$ (3)

for species $\alpha$ and $\beta$.

Assume now we have a PIC discretization of the distribution function, $f$,
with the distribution function sampled via Dirac delta functions. This
generates a phase space distribution of the form

$\displaystyle f=\sum_{p}w_{p}\delta(z-z_{p})$ (4)

where p represents the particle indices, $w$ the particle weights, and the
configuration space variable $z_{p}=(x_{p},v_{p})$. It is clear that this
basis doesn’t allot evaluation of the collision operator of 3 due to 4 being
neither smooth nor differentiable. One option is to deposit particles to the
grid and evaluate the collision integral with finite elements and synethesize
the new weights from the grid at the particles. An in depth conservative
discretization of such a method is given in [8] with accuracy and performance
benchmarks in [12] and the conservative projection procedure discussed in
[14]. A more natural approach for codes which do not inherently employ a grid
based method would be to formulate 3 such that the basis of particles are
represented sufficiently smooth. Recent theoretical advances have proven the
necessary qualities of such a formualation [15].

Recall the evolution in terms of the brackets

$\displaystyle\frac{dU}{dt}=\\{U,F\\}+(U,F)$ (5)

with $F=H-S$ representing the free energy, and $S$ the entropy. The first
term, $\\{U,F\\}$ represents Hamiltonian dynamics in the Poisson bracket, with
the right term, $(U,F)$ representing collisional dynamics in the metric
bracket. The mass, momentum, and energy are Casimirs of both brackets with the
Hamiltonian being a Casimir of the Poisson and metric bracket. The entropy
functional is a Casimir of the Poisson bracket, but the metric bracket
enforces $\frac{dF}{dt}\geq 0$ corresponding to monotonic entropy production.
For an in depth discussion on these Casimir invariants, see [6, 7, 1].

We can think of the collision integral as advective by introducing the entropy
functional

$\displaystyle S[f]=-\int flnf$ (6)

and using the identity

$\displaystyle\partial_{v}f=f\partial_{v}lnf=-f\partial_{v}\frac{\delta
S}{\delta f}$ (7)

with

$\displaystyle\frac{\partial f}{\partial t}+\frac{\partial f}{\partial
v}\cdot(Ff)=0$ (8)

and

$\displaystyle F[f](v_{p})=\gamma\int
f(v^{\prime})Q(\xi)\cdot(\frac{\partial}{\partial v_{p}}\frac{\delta S}{\delta
f}-\frac{\partial}{\partial v_{\bar{p}}}\frac{\delta S}{\delta f})$ (9)

Hirvijoki demonstrates for arbitrary functionals, we have the identity

$\displaystyle\frac{\partial}{\partial v_{p}}\frac{\delta A}{\delta
f}=\frac{1}{w_{p}}\frac{\partial A}{\partial v_{p}}$ (10)

for distributions, $f$, of the form in 4, and

$\displaystyle Q(\xi)=\frac{1}{|\xi|^{2}}(I-\frac{\xi^{T}\xi}{|\xi|})$ (11)

being the scaled projection matrix with

$\displaystyle\xi=v_{p}-v_{\bar{p}}$ (12)

This formulation demonstrates the evolution of the particle velocities in
terms of the entropy functional, 6, but the caveat remains the lack of
continuity in $f$.

The remaining component for the evaluation of this operator comes in the form
of a mollification procedure [15] with the mollification process in [19].
Recall again the entropy functional 6 over the particle distribution. We can
generate a smooth, continuous distribution function in terms of particles by
convolution with some smooth radial basis function. For instance, for a
Gaussian given by

$\displaystyle\psi_{\epsilon}(v-v_{p})=\frac{1}{2\epsilon}e^{-(v-v_{p})^{2}/\epsilon}$
(13)

such that

$\displaystyle S_{\epsilon}=\int\psi_{\epsilon}*fln\psi_{\epsilon}*f$ (14)

In practice, this results in replacing the delta functions in the evaluation
of the entropy functional with Gaussians regularized by the parameter
$\epsilon$, controlling how "wide" particle interactions are in terms of its
entropy. The resultant particle basis entropy functional reads

$\displaystyle
S_{\epsilon}=\sum_{v_{p}}\int\psi_{\epsilon}(v-v_{p})ln\sum_{v_{\bar{p}}}\psi_{\epsilon}(v-v_{\bar{p}})$
(15)

Now we can express the evolution of the velocity of a particle in terms of the
entropy functional gradients around that particles as

$\displaystyle\frac{dv_{p}}{dt}=\sum_{v_{\bar{p}}}1(p,\bar{p})Q(\xi)\Gamma(S_{\epsilon},p,\bar{p})$
(16)

where we have maintained Hirvijoki’s convention of $1(p,\bar{p})$ evaluating
to 1 if particles are within the same collision cell, and zero otherwise, and
the term $\Gamma(A,p,\bar{p})=\frac{1}{w_{p}}\frac{\partial A}{\partial
v_{p}}-\frac{1}{w_{\bar{p}}}\frac{\partial A}{\partial v_{\bar{p}}}$ from the
identity in 10.

In the first paper, [15], the mathematical rigor was focused primarily on
proving the Maxwellian steady state and conservative properties of this
integral and forward Euler was used for the numerical tests. This time
discretization was not conservative in energy, and the solution was put
forward by Hirvijoki in [16] where he evaluates the scaled projection matrix
in terms of half time steps to maintain the conservative properties of the
time continuous form given in 16. This time discretization reads

$\displaystyle\frac{v_{p}^{n+1}-v_{p}^{n}}{\Delta
t}=\sum_{v_{\bar{p}}}Q(\xi^{(n+\frac{1}{2})})\Gamma(S_{\epsilon}^{n},p,\bar{p})$
(17)

and grants the conservation of mass, momentum, and energy. The final component
remains the evaluation of the gradient with respect to 15 as the
$\Gamma(S,p,\bar{p})$ term requires it in relation to particle velocity at
index p. The derivation is straightforward so long as we are careful of
particle indices.

$\displaystyle\frac{\partial S_{\epsilon}}{\partial v_{p}}$
$\displaystyle=-\int\nabla_{v_{p}}\sum_{p}w_{p}\psi_{\epsilon}(v-v_{p})ln\sum_{k}w_{k}\psi_{\epsilon}(v-v_{k})$
$\displaystyle=-\int\sum_{p}w_{p}\psi^{{}^{\prime}}(v-v_{p})\delta_{pq}ln\sum_{k}w_{k}\psi(v-v_{k})$
$\displaystyle-\int\sum_{q}w_{q}\psi(v-v_{q})\frac{1}{\sum_{k}w_{k}\psi(v-v_{k})}\sum_{k}w_{k}\psi(v-v_{k})\delta_{kp}$
$\displaystyle=-\int
w_{p}\psi^{{}^{\prime}}(v-v_{p})ln\sum_{k}w_{k}\psi(v-v_{k})$
$\displaystyle-\int\sum_{q}w_{q}\psi(v-v_{q})\frac{w_{p}\psi^{\prime}(v-v_{p})}{\sum_{k}w_{k}\psi(v-v_{k})}$
$\displaystyle=-\int
w_{p}\psi^{\prime}(v-v_{p})ln\sum_{k}w_{k}\psi(v-v_{k})-\int
w_{p}\psi^{\prime}(v-v_{p})$ $\displaystyle=-\int
w_{p}\psi^{\prime}(v-v_{p})(1+ln\sum_{k}w_{k}\psi(v-v_{k}))$ (18)

An important point to note before moving on to numerical implementation and
testing is what the equation in 17 is evolving. This formulation results in an
evolution of the particle state in phase space (acts directly on velocity),
whereas classic interpretations of the integral capture the transfer of
momentum through redistribution of the particle weights according to the
evolution of the distribution function on the grid.

Another nuance here is the sampling regime these methods converge for, and the
importance of the parameter $\epsilon$, as this has stringent implications on
the convergence of the particle method to physically meaningful solutions. For
a detailed discussion on the convergence of the methods, see [19], however it
is important enough to briefly discuss again here.

With the particle spacing

$\displaystyle h=\frac{2L}{N}$ (19)

with N being the number of particles along any given slice of an axis and
total number of particles given by $N^{d}$. This gives and expression for
$\epsilon$ which has been experimentally demonstrated to converge well as
$\epsilon=2h^{1.98}$. These grid methods typically converge with
$\epsilon=o(h)$ [15]. This, for the moment, precludes an expression for
$\epsilon$ in terms of the discrepency of the set used to sample $f$. In other
words, a Klimontovich representation with Monto-Carlo or quasi Monte-Carlo
initializations of $f$ with evenly distributed weights are infeasible with
this expression. An $\epsilon$ can be found to produce physical results in
this regime, but not reliably in experimentation. An expression cast in terms
of the discrepency is an active area of research by the authors.

##  3 Numerical Implementation and Experiment

We present in this section the details on the implementation of the collision
operator and some numerical experiments. To start, we will discuss the details
of time discretization as well as how the quadrature is handled to evaluate
the velocity space gradients of 6 in 17. All code for these examples is
publicly available in the main branch of PETSc (at the time of writing,
version 3.19).

###  3.1 Numerical Tools

It is obvious by the form of 17 in relation to the entropy functional that the
particle loop constitutes an $O(N^{2})$ computation per time step,
necessitating careful handling. This is further compounded by the implicit
solve necessary to compute the particle evolution at each time step, or fixed
point method used in [17]. For this task, PETSc’s time stepping library is
used to construct a sufficiently accurate and fast solver such that we can
evaluate collisions in a reasonable amount of time, which are additionally
straightforward to parallelize by construction.

To that purpose, the main particle loop is handled via Implicit Midpoint to
evaluate the $Q(\xi^{n+\frac{1}{2}})$, whereas $\Gamma(S,p,\bar{p})$ terms are
evaluated at time $n$ and are thus precomputed prior to each time step to
reduce overhead in each solver iteration. The Jacobian is approximated using
either a standard finite difference, or matrix free approach. The matrix free
approach drastically improves runtime, particularly due preallocations of the
Jacobian matrix being nonoptimal at the time of writing. Particle management
is performed by the PETSc object DMSwarm, which initializes the distribution
functions, and provides data structure for particle management as well as
algorithms for particle migration between MPI ranks in parallel runs.

The initial particle distribution is the regular sampling scheme. Here,
particles are laid out in a grid with the spacing determined by
$h=\frac{2L}{N}$, where $[-L,L]^{d}$ is the domain represented by the grid,
and N the number of particles along a line in one of the dimensions. Particle
weights are assigned by performing the integral

$\displaystyle
w_{p}=\int_{v_{p}-h/2}^{v_{p}+h/2}\frac{1}{\sqrt{2\pi\theta}}e^{(-\vec{v_{p}}^{2}/\theta)}$
(20)

with $\theta=k_{B}T_{s}/m_{s}$ in the isotropic case. Note that the spatial
component is omitted as this exists within one spatial cell. Extension to a
global maxwellian normalized to all space within spatial cells would merely
require regularization of the particle number to the global sum

$\displaystyle f_{p}=\sum_{p}w_{p}\delta(v-v_{p})=1$ (21)

extension to which is trivial in this scheme for the collision operator.

####  3.1.1 Quadrature

Given the form of 2, a sufficiently accurate quadrature scheme is necessary to
perform the integration. This can naturally be a Gauss-Hermite quadrature
given the form of the entropy functional, however we use a different
technique. Saye presents an approach to generating high order qudratures on
implicit domains defined by polynomials in [18], which is publicly available
in the package Algoim. The benefits of using such a package are the ability to
tailor a high order quadrature for the computational domain of varying degree
and fidelity.

To compute integrals over $S_{\epsilon}$, we generate a high order quadrature
over a disc in 2V, and a ball in 3V. The domain is selected to sufficiently
encapsulate the particles within the quadrature. Figure 1 presents a
visualization of these quadrature schemes in Paraview. Domains are updated and
new quadratures formed as $\epsilon$ evolves. This ensures, without making a
change of variables at each step, that the quadrature accurately encapsulates
evaluations of the entropy functional gradients.

(a)

(b)

Figure 1: Quadratures in 2 (left) and 3 (right) dimensions generated by
algoim. Quadratures are third order with 6 points per level set.

###  3.2 Equilibration

The most basic test for the collision operator is the appropriate thermal
convergence of two thermal populations. The first population is given a
temperature of 0.35KeV, whereas the second population is given a temperature
of 0.2KeV. The distribution then relaxes to an equilibrium state. In this
case, an initial electron-positron distribution is selected to demonstrate
conservative properties of the operator in a trivial case of $m_{1}/m_{2}=1$.
Larger mass ratios are investigated in later tests such as isotropization.

(a)

(b)

(c)

(d)

Figure 2: Convergence to equilibrium for electron-positron collisions with
both Gauss-Hermite and Algoim evaluations of the entropy functional gradient
quadrature. An initial jump in energy error occurs on the order of $10^{-}12$
and remains consistent with numerical noise.

The appropriate behavior is observed with the two populations effectively
reaching an equilibrium state, with error on the order of $10^{-12}$ in
kinetic energy of the system. Effective error reduction can be achieved by
utilizing a smaller time step, as well as using the dense finite difference
approximation to the Jacobian. The matrix free approach is leveraged here with
relatively large time steps, where the solver is easily capable of converging.
This is not the case with species presenting relatively large mass ratios on
the order of realistic electron-ion collisions, where the step size must be
significantly reduced. This test case is also light on the quadrature used,
with a 3rd order quadrature generated over the disc using 10 points per level
set.

###  3.3 Two Species Isotropization

In isotropization tests, each velocity space dimension is given a different
value for its temperature for each species in 2V. In the 3V case, we have
$T_{x}=T_{y}$ to generate a uniform parallel and transverse temperature. The
relaxation rate is given in the Naval Research Laboratory Plasma Formulary
(NRLPF) as [20]

$\displaystyle\frac{dT_{\perp}}{dt}=\frac{1}{2}\frac{dT_{||}}{dt}=-v_{T}^{\alpha}(T_{\perp}-T_{||})$
(22)

for $A=T_{\perp}/T_{||}-1>0$ we have

$\displaystyle
v_{T}^{\alpha}=\frac{2\sqrt{\pi}e_{\alpha}^{2}e_{\beta}^{2}n_{\alpha}\lambda_{\alpha\beta}}{m_{\alpha}^{1/2}(kT_{||})^{1/2}}A^{-2}[-3+(A+3)\frac{tan^{-1}(A^{1/2})}{A^{1/2}}]$
(23)

When $A<0$, the arctan is replaced with hyperbolic arctan and its arguments
negated. Additional parameters include the mass, in which case we use the
electron mass for species 0, and two times the proton mass for species 1 (an
approximate Deuteron). Isotropization results are shown to agree well with
empirical results from the NRL.

(a)

(b)

Figure 3: Two species isotropization with quadrature over the disc (right) and
gauss hermit quadrature (left) in 2V. Variation in rate from the NRL is
attributed to the NRL assuming each successive step is Maxwellian where there
is no guarantee in the intermediate steps towards convergence.

Decent agreement is seen with a minor deviation in the rate. This is
attributed to the assumptions in the empirical equations that the intermediate
states of the distribution are Maxwellian, where the guarantee of closness to
the Maxwellian of intermediate states is not enforced. Adams et al. [13] and
Hager [21] present qualitatively similar results.

##  4 Concluding Remarks and Future Developments

We have reiterated the metriplectic formulation, and the recent theoretical
advancements that have made evaluation of the Landau collision integral
possible. Demonstration of the appropriate behavior of the collision integral
was performed on common tests for collision operators and compared to data
from the NRLPF, providing cross verification of the collision integral done by
Zonta et al. in a completely separate implementation, for a different machine
architecture, with different methods for computing its integrals. All methods
demonstrate considerable agreement. In addition to cross verification, we
presented and provided a community software implementation of this collision
integral for use in the construction of, or supplementation of, large scale
plasma codes. Naturally, some additional studies and optimizations are a
subject of ongoing development.

The stringent structure of the sampling method to obtain stability of this
operator is of first consideration. In the case of simulations which use free
streaming particles, there is no guarantee of the structure at each collision
step being sufficiently regular, as such, two solutions are posed to be
investigated, firstly the method of coarse graining. In this case, particles
are projected from the input basis to a sufficiently regular particle basis.
Weights would then be determined according to the projected moments of the
initial basis using the particle mass matrix presented in [14]. The second
method, which is more mathematically rigorous, is the cast the epsilon
parameter in terms of the discrepency of sampling function. This is a more
natural approach but would incur the additional overhead of measuring the
discrepency of the sample prior to the collision step. As stated, both of
these avenues are left to future work.

The final optimization to be examined is the extension of the PETSc
implementation to GPUs, similar to the work performed in [17] by Zonta et al,
but in the PETSc Vec GPU back end. This work is underway and similarly to the
operator provided in the continuum will leverage a Kokkos [22] backend for
portability requirements. This additional component sets the stage for real
world physics studies and comparisons between collision operators on exascale
problems with a particle basis full geometry solver in conjunction with
additional verification test of the Discrete Gradients implementation underway
in PETSc as well.

##  5 Acknowledgements

This work was supported in part by the U.S. Department of Energy, Office of
Science, Office of Fusion Energy Sciences and Office of Advanced Scientific
Computing Research, Scientific Discovery through Advanced Computing (SciDAC)
program through the FASTMath Institute under Contract No. DE-AC02-05CH11231 at
Lawrence Berkeley National Laboratory. Additionally, this work was funded in
part by the Exascale Computing Project (17-SC-20-SC), a collaborative effort
of the U.S. Department of Energy Office of Science and the National Nuclear
Security Administration. We gratefully acknowledge assistance on the Geosolver
cluster at the UB Center for Computational Research by CCR staff, and the
author of Algoim, Robert Saye, for insight into effective usage of the package
and its algorithms.

## References

  * [1] P.. Morrison  “Structure and structure-preserving algorithms for plasma physics”  In _Physics of Plasmas_ 24.5, 2017, pp. 055502  DOI: 10.1063/1.4982054
  * [2] Philip J. Morrison  “A paradigm for joined Hamiltonian and dissipative systems”  In _Physica D: Nonlinear Phenomena_ 18.1, 1986, pp. 410–419  DOI: https://doi.org/10.1016/0167-2789(86)90209-5
  * [3] Ernst Hairer, Christian Lubich and Gerhard Wanner  “Geometric numerical integration” Structure-preserving algorithms for ordinary differential equations 31, Springer Series in Computational Mathematics  Springer-Verlag, Berlin, 2006, pp. xviii+644
  * [4] Ernst Hairer, Christian Lubich and Gerhard Wanner  “Geometric numerical integration illustrated by the Störmer–Verlet method”  In _Acta Numerica_ 12 Cambridge University Press, 2003, pp. 399–450  DOI: 10.1017/S0962492902000144
  * [5] L D Landau  “Kinetic equation for the case of Coulomb interaction”  In _Phys. Z. Sowjet_ , 1936
  * [6] Michael Kraus and Eero Hirvijoki  “Metriplectic integrators for the Landau collision operator”  In _Physics of Plasmas_ 24.10, 2017, pp. 102311  DOI: 10.1063/1.4998610
  * [7] Eero Hirvijoki, Michael Kraus and Joshua W Burby  “Metriplectic particle-in-cell integrators for the Landau collision operator”  In _arXiv_ , 2018  eprint: 1802.05263
  * [8] Eero Hirvijoki and Mark Adams  “Conservative discretization of the Landau collision integral”  In _Physics of Plasmas_ 24, 2016  DOI: 10.1063/1.4979122
  * [9] Satish Balay et al.  “PETSc Users Manual”, 2020  URL: https://www.mcs.anl.gov/petsc
  * [10] Satish Balay et al.  “PETSc Web page”, https://www.mcs.anl.gov/petsc, 2019  URL: https://www.mcs.anl.gov/petsc
  * [11] Shrirang Abhyankar, Jed Brown, Emil M Constantinescu, Debojyoti Ghosh, Barry F Smith and Hong Zhang  “PETSc/TS: A Modern Scalable ODE/DAE Solver Library”  In _arXiv preprint arXiv:1806.01437_ , 2018
  * [12] Mark F Adams, Dylan P Brennan, Matthew G Knepley and Peng Wang  “Exascale Landau collision operator in the CUDA programming model applied to thermal quench plasmas”  In _arXiv preprint arXiv:2104.10000_ , 2021
  * [13] Mark F. Adams, Peng Wang, Jacob Merson and Matthew G. Knepley  “A performance portable, fully implicit Landau collision operator with batched linear solvers”, 2023  arXiv:2209.03228 [physics.plasm-ph]
  * [14] Joseph V. Pusztay, Matthew G. Knepley and Mark F. Adams  “Conservative Projection Between Finite Element and Particle Bases”  In _SIAM Journal on Scientific Computing_ 44.4, 2022, pp. C310–C319  DOI: 10.1137/21M1454079
  * [15] J. Carrillo, Jingwei Hu, Li Wang and Jeremy Wu  “A particle method for the homogeneous Landau equation”  In _Journal of Computational Physics: X_ 7, 2020, pp. 100066  DOI: 10.1016/j.jcpx.2020.100066
  * [16] Eero Hirvijoki  “Structure-preserving marker-particle discretizations of Coulomb collisions for particle-in-cell codes”  In _Plasma Physics and Controlled Fusion_ 63.4 IOP Publishing, 2021, pp. 044003  DOI: 10.1088/1361-6587/abe884
  * [17] Filippo Zonta, Joseph V. Pusztay and Eero Hirvijoki  “Multispecies structure-preserving particle discretization of the Landau collision operator”  In _Physics of Plasmas_ 29.12, 2022, pp. 123906  DOI: 10.1063/5.0105182
  * [18] Robert I. Saye  “High-order quadrature on multi-component domains implicitly defined by multivariate polynomials”  In _Journal of Computational Physics_ 448, 2022, pp. 110720  DOI: https://doi.org/10.1016/j.jcp.2021.110720
  * [19] J. Carrillo, Katy Craig and Francesco Saverio Patacchini  “A blob method for diffusion”  In _Calculus of Variations and Partial Differential Equations_ 58, 2019  DOI: 10.1007/s00526-019-1486-3
  * [20] J.. Huba  “NRL plasma formulary”  Washington, DC: Naval Research Laboratory], 2009
  * [21] Robert Hager, Eisung Yoon, Seung-Hoe Ku, Ed F. D’Azevedo, Patrick H. Worley and Choong-Seock Chang  “A fully non-linear multi-species Fokker–Planck–Landau collision operator for simulation of fusion plasma”  In _Journal of Computational Physics_ 315 Elsevier BV, 2016, pp. 644–660  DOI: 10.1016/j.jcp.2016.03.064
  * [22] H. Edwards, Christian R. Trott and Daniel Sunderland  “Kokkos: Enabling manycore performance portability through polymorphic memory access patterns” Domain-Specific Languages and High-Level Frameworks for High-Performance Computing  In _Journal of Parallel and Distributed Computing_ 74.12, 2014, pp. 3202–3216  DOI: https://doi.org/10.1016/j.jpdc.2014.07.003