# Beyond mean field: condensate coupled with pair excitations

Stephen Sorokanich  Department of Mathematics, University of Maryland, College
Park<EMAIL_ADDRESS>
###### Abstract

We prove existence results for a system of partial differential equations
describing the approximate condensate wavefunction and pair-excitation kernel
of a dilute (T=0) Bose gas in the stationary setting, in the presence of a
trapping potential and repulsive pairwise atomic interactions. Notably, the
Hartree-type equation for the condensate in this system contains contributions
from non-condensate particles, and the pair excitation kernel satisfies a
nonlinear operator equation. The equations studied here are inspired by the
work of Griffin, who derived this system in the study of finite temperature
condensates. The techniques employed include a variational principle, which
exploits the connection between unitary Bogoliubov rotations and a nonlinear
operator equation for the pair excitation kernel. An iterative procedure for
constructing solutions is also included.

Quantum many-body dynamics, nonlinear integro-partial differential equation,
non-self-adjoint operators, variational principle, Fock space, quasiparticle,
operator Riccati equation

††preprint: AIP/123-QED

##  I Introduction

In this paper, we study a nonlinear system of integro-partial-differential
equations which arise in the study of a trapped Bose-Einstein condensate (BEC)
in equilibrium with particles outside the condensate. The equations are
derived (in a physically-inspired although heuristic fashion) from quantum
mechanical first principles via a model for a finite-but-large number of
interacting bosons, inspired by several preeminent works in physics
Griffin1996 ; fetter72 ; wu61 . The quantities described by this system are
the macroscopic condensate wavefunction,
$\phi(x):\mathbb{R}^{3}\to\mathbb{C}$, and the pair-excitation kernel,
$k(x,y):\mathbb{R}^{3\times 3}\to\mathbb{C}$. Together, these quantities
provide a description of the quantum ground state of a dilute Bose gas which
generalizes the mean-field theory MGM ; GrillakisMargetisSorokanich ;
Margetis2012 ; GMM2010 .

The mean-field theory for BEC specifies a nonlinear Schrödinger equation for
the single-particle condensate wavefunction, $\phi(x)$. Assuming that the
effective ground state of the system is a tensor product involving only this
mean-field condensate function yields an approximation to the exact many-body
Hamiltonian which is quadratic and diagonal in field operators for states
orthogonal to the condensate. This is not true for models that account for the
effects of pair-excitation. Such models can be described by two functions,
referred to as $\phi(x)$ and $k(x,y)$ in this work, and are associated with an
approximate Hamiltonian that is quadratic, but not necessarily diagonal, in
the field operators. The heuristic approximation of the exact many-body
Hamiltonian considered in this work takes includes contributions from all two-
point correlation functions associated with non-condensate particles.

The models and techniques presented here, as well as the use of fundamental
quantities $\phi,\,k$ to describe the pair-excited ground state, are inspired
by several well-known works in the physical modeling of Bose-Einstein
condensation. The quantum field theory of BEC originated in the work of
Bogoliubov, Gross, and Pitaevskii Bogoliubov1947 ; Gross61 ; Pitaevskii61 ,
who formulated the mean-field approximation and the concept of macroscopic
wavefunction for the condensate (proportional to $\phi$ in our terminology).
Lee, Huang, Yang, and Luttinger began a comprehensive study of the Bose gas in
the periodic box in the 1950s leehuangyang ; leeyang ; HuangYang1957 ;
huangyangluttinger  and derived some of the first corrections to the ground
state energy predicted by the mean field theory; these corrections are
produced by correlations between the condensate and pairs of non-condensate
particles with equal-and-opposite momenta. This early work for the periodic
box was adapted to the non-periodic setting in the 1960s by T. T. Wu wu58 ;
wu61 , who was the first to formulate the pair-excitation kernel (the quantity
$k$), and derived the fundamental correction to the mean-field ground state
represented schematically by:

$\Psi_{0}=\exp(P[k])\Psi_{\mathrm{mf}}[\phi].$ (I.1)

Here $\Psi_{\mathrm{mf}}[\phi]$ is the tensor-product ground state of the
mean-field theory, and $P[k]$ is the (non-Hermitian) pair-excitation operator.

The primary motivations for this work are the quadratic Hamiltonian models of
Fetter (1971) fetter72 , and Griffin (1996) Griffin1996 ; GriffinBook . These
authors formulate the unitary (Bogoliubov) rotation of quadratic, non-
translation-invariant Hamiltonians describing systems in the presence of a
trapping potential. By proper choice of this unitary transformation (which
amounts to the choice of a non-orthogonal basis of single-particle states), a
diagonal transformed Hamiltonian is obtained, and the excitation spectrum can
be directly inferred. The operator basis in which the Hamiltonian is
diagonalized is commonly referred to as the basis of quasiparticles
huang2008statistical . We instead approach the mentioned quadratic models from
the formalism of pair-excitation. In previous work GrillakisMargetisSorokanich
; Sorokanich-Thesis , we showed that there is an intimate connection between
the quasiparticle basis which appears in the work of Fetter fetter72 , and the
pair-excitation kernel that appears in the second-order correction to the
ground state described by Wu wu61 , and represented by the expression (I.1).
The connection between this kernel and the single-particle quasiparticle basis
parallels results in operator theory (see, e.g., the works of Tretter Tretter-
book ; Tretter2016 , Albeverio, et. al. AlbeverioMotovilov2010 ;
AlbeverioMotovilov2019 ; Kostrykin03-chap ). An existence theory for this
kernel in the specific context of the Bose gas was developed in
GrillakisMargetisSorokanich , based on a variational principle.

The purpose of this paper is twofold. First, we apply the method of pair-
excitation to a family of quadratic models for the Bose gas which are not
described by our previous work. This extension is nontrivial, since the model
in this work contains nonlinear couplings between the condensate and excited
states (in either the quasiparticle representation or in the representation
using the pair-excitation kernel). Second, we utilize the new representation
of this system in terms of the excitation kernel to give an existence theory
for the resulting time-independent nonlinear equations. The operator equation
for $k$ in this model will contain nonlinear terms which represent non-
condensate particle densities. We will show how it is possible to solve the
nonlinear system for $k$ and $\phi$ using a variational approach. Two methods
for solving the variational problem are given; the first method is general and
is a direct extension of our previous method. The second method is an
iteration scheme, which gives more detailed information about the solution
found.

Our primary source in the physics literature for the model studied here is
Griffin Griffin1996 ; GriffinBook , although it is sufficiently general to be
called the Bogoliubov de Gennes system Bach2022 . The excitations of a Bose
gas in non-translation invariant settings has been studied by many authors;
see, e.g., Dalfovoetal1999 ; Gardiner1997 ; Leggett2001 ; Ozeri2005 ;
Rovenchak2016 . While time-dependent systems related to the one studied here
have been the subject of past work Bach2022 ; Nam2017 ; GM2013-a ; GM2013-b ,
we believe that our focus on the stationary case can potentially shed new
light onto the problem. As far as the author is aware, this is the first
rigorous treatment of the explicitly stationary problem, and the first time
these equations have been linked to a pair excitation formalism. In addition
to these works, we refer the reader to recent efforts to place the quadratic
approximation of bosonic Hamiltonians on a rigorous basis Seiringer ;
NamNapiorkowski2017 ; NamNapiorkowski2017-II ; NamSolovej2016 ;
NamSeiringer2015 ; BrenneckeSchlein2019 ; the model presented here has
potential to be rigorously treated as an approximation to the exact
Schrödinger dynamics over quasifree states Zagrebnov2001 , although we do not
follow this path in this work.

###  I.1 Mathematical prerequisites and Hamiltonian model

We now introduce the exact Hamiltonian for interacting bosons, and the
heuristic reduction of this Hamiltonian to a quadratic expression in field
operators inspired by the work of Griffin Griffin1996 ; GriffinBook . The
Bogoliubov rotation of the resulting reduced Hamiltonian is introduced, and
the mean field equations for the quantities
$\big{(}\phi(x),u_{j}(x),p_{j}(x)\big{)}_{j=1}^{\infty}$ are derived. We then
describe the equivalent description of these quantities in terms of the
condensate wavefunction and the pair-excitation kernel,
$\big{(}\phi(x),k(x,y)\big{)}$. This subsection will also serve to introduce
the notation and terminology used throughout the rest of the work.

Our starting point is the bosonic Fock space, denoted $\mathbb{F}$, which is
defined as the direct sum of all nn-particle symmetric Hilbert spaces,

$\mathbb{F}:=\mathbb{C}\oplus\bigoplus_{n=1}^{\infty}{L^{2}_{\mathrm{sym}}(\mathbb{R}^{3n})}.$

The operators $a_{x},a_{x}^{\ast}$ are the usual bosonic field operators on
$\mathbb{F}$, which satisfy the canonical commutation relations:

$[a_{x},a_{y}]=0=[a_{x}^{\ast},a_{y}^{\ast}],\quad[a_{x},a_{y}^{\ast}]=\delta(x-y).\quad
x,y\in\mathbb{R}^{3}.$ (I.2)

Here, as in the rest of the work, $\delta(x-y)$ denotes the Dirac mass of the
delta function. The exact many-body Hamiltonian, $\mathcal{H}$, for $N$ bosons
with a repulsive, mean-field, pairwise interaction $\upsilon_{N}$, in the
presence of trapping potential $V_{\mathrm{trap}}(x)$, is given in second-
quantized form by the Fock space operator:

$\mathcal{H}=\int{{\rm
d}x\,\\{a_{x}^{\ast}\big{(}-\Delta+V_{\mathrm{trap}}(x)\big{)}a_{x}\\}}+\frac{1}{2}\iint{{\rm
d}x{\rm
d}y\,\big{\\{}a_{x}^{\ast}a_{y}^{\ast}\,\frac{1}{N}\upsilon_{N}(x-y)a_{x}a_{y}\big{\\}}}.$
(I.3)

The Hamiltonian $\mathcal{H}$ includes the mean-field type interaction,
$\frac{1}{N}\upsilon_{N}(x)$, which is required for the analysis that follows.
Specifically, we take $\upsilon_{N}(x)=N^{3\beta}\upsilon(N^{\beta}x)$, where
$\upsilon$ is a smooth, positive, radially symmetric function, $0<\beta\leq
1$, with $\int{{\rm d}x\\{\upsilon(x)\\}}=g$, for $g>0$. As $N\to\infty$,
$\upsilon_{N}$ is assumed to approximate the delta function. The trapping
potential $V_{\mathrm{trap}}(x)$ is assumed to be positive, and moreover, we
assume that the single-particle Schrödinger operator
$(-\Delta_{x}+V_{\mathrm{trap}}(x))$ has a positive spectrum with a gap
between it’s lowest and first excited states. It is useful to take the
concrete example $V_{\mathrm{trap}}(x)=|x|^{2}$ for $x\in\mathbb{R}^{3}$ as
the prototypical trapping potential. We also introduce the number operator,
$\mathcal{N}:\mathbb{F}\to\mathbb{F}$, defined by

$\mathcal{N}:=\int{{\rm d}x\\{a_{x}^{\ast}a_{x}\\}}.$ (I.4)

Since the following approximation scheme will not conserve the total number of
particles, we consider instead the operator
$\mathcal{K}:=\mathcal{H}-\mu\mathcal{N}$, where $\mu>0$ is the chemical
potential, a Lagrange multiplier constraining the average total number of
particles.

The assumptions on the interaction potential yield the following $L^{\infty}$
estimate

$\frac{1}{N}\upsilon_{N}(x-y)\leq\|\upsilon\|_{L^{\infty}},\quad
0\leq\beta\leq 1/3,$ (I.5)

and by Young’s convolution inequality (where the operation $\ast$ denotes
convolution), the estimate

$\|\frac{1}{N}(\upsilon_{N}\ast f)\|_{L^{2}}\leq g\|f\|_{L^{2}},\quad f\in
L^{2}(\mathbb{R}^{3}).$ (I.6)

The exact operator $\mathcal{K}$, with $\mathcal{H}$ coming from equation
(I.3), is a quartic expression in field operators $a,\ a^{\ast}$ due to the
presence of the two-particle interaction. The main question that motivates the
quadratic reduction of this exact Hamiltonian is: what is the low-lying
spectrum of $\mathcal{K}$, and can this spectrum be approximated by the
spectrum of an effective quadratic operator? We do not take up this question
in rigorous detail here; our goal is rather to study a simplified effective
model which incorporates features of the many-body system for states close to
the ground state, and we take the lead from the significant work in physics
already mentioned. We emphasize the utility of an approximation which is
quadratic in field operators orthogonal to the condensate, since quadratic
Hamiltonians can be diagonalized by unitary rotation. Two systematic, although
heuristic, approximation schemes are now summarized, motivated by the works
Griffin1996  and fetter72  respectively, and the unitary rotation of the
resulting quadratic operators. For the sake of completeness, elements of the
reduction of operator $\mathcal{K}$ are given. Details of specific
calculations are omitted, and we refer to Griffin1996 ; GriffinBook ; fetter72
for the ultimate justifications that these authors give. The approximation
schemes are heuristic, since they do not address, e.g., how accurately the
ground state of the reduced Hamiltonian approximates the ground state of
Hamiltonian (I.3).

Formal reduction (Griffin, 1996 Griffin1996 ): First, the heuristic Bogoliubov
approximation is employed by which the field operators $a_{x}^{\ast},\,a_{x}$
appearing in the operator $\mathcal{K}=\mathcal{H}-\mu\mathcal{N}$ are
decomposed into a macroscopic scalar field proportional to the condensate
wavefunction $\phi(x)$, in addition to field operators,
$b^{\ast}_{x},\,b_{x}$, which create/annihilate particles in states orthogonal
to the condensate, viz.,

$a_{x}\approx\sqrt{N}\,\,\overline{\phi(x)}+b_{x},\quad
a_{x}^{\ast}\approx\sqrt{N}\phi(x)+b^{\ast}_{x}.$ (I.7)

The symbol “$\approx$" is used here, and throughout this section, to denote
the fact that (I.7) is a physically-motivated, but ultimately formal
substitution in the current treatment (referring, e.g., to the references
Seiringer ; Zagrebnov2001 ; NamNapiorkowski2017 ; NamNapiorkowski2017-II  for
a nuanced treatment of this canonical approximation). The operators
$b_{x},\,b_{x}^{\ast}:\mathbb{F}\to\mathbb{F}$ are understood to satisfy the
following commutation relations:

$\begin{split}[b_{x}^{\ast},b_{y}^{\ast}]&=0=[b_{x},b_{y}],\\\
[b_{x},\,b^{\ast}_{y}]&=\delta(x-y)-\phi(x)\overline{\phi(y)}:=\widehat{\delta}(x-y),\quad
x,y\in\mathbb{R}^{3}.\end{split}$ (I.8)

After making the Bogoliubov approximation, operator $\mathcal{K}$ is further
reduced by replacing cubic and quartic terms in $b_{x},b_{x}^{\ast}$ with
effective quadratic operators. This replacement depends on the three (as yet
unspecified) functions, which we denote
$m^{\mathrm{pair}}(x,y),\,n^{\mathrm{pair}}(x,y),\,\rho^{\mathrm{pair}}(x)$,
via,

$\begin{split}b_{x}^{\ast}b_{y}^{\ast}b_{x}b_{y}&\approx\,\,\,\,2n^{\mathrm{pair}}(x,y)b^{\ast}_{y}b_{x}+2\rho^{\mathrm{pair}}(x)b^{\ast}_{y}b_{y}+\overline{m^{\mathrm{pair}}(x,y)}b_{x}b_{y}+m^{\mathrm{pair}}(x,y)b^{\ast}_{x}b^{\ast}_{y}\,\,,\\\
b^{\ast}_{x}b_{x}b_{y}&\approx\,\,\rho^{\mathrm{pair}}(x)b_{y}+n^{\mathrm{pair}}(x,y)b_{x}+m^{\mathrm{pair}}(x,y)b^{\ast}_{x}\,\,.\end{split}$
(I.9)

The specific combinations which appear in (I.9) are consistent with their role
in Griffin1996 , and will soon be clarified. With the substitutions (I.7) and
(I.9), operator $\mathcal{K}$ becomes a quadratic expression involving the
undetermined quantities
$\phi(x),\,m^{\mathrm{pair}}(x,y),\,n^{\mathrm{pair}}(x,y),\,\rho^{\mathrm{pair}}(x)$,
which we rewrite as the following ‘approximate’ Hamiltonian:

$\begin{split}&\mathcal{H}_{\mathrm{approx}}:=\mathrm{const}+\int{dx\big{\\{}\overline{B(x)}b_{x}+B(x)b_{x}^{\ast}\big{\\}}}\\\
&+\int{dxdy\Big{\\{}(b_{x},b_{x}^{\ast})\begin{pmatrix}-H^{T}(x,y)&(\upsilon_{N}\overline{m})(x,y)\\\
-(\upsilon_{N}m)(x,y)&H(x,y)\end{pmatrix}\begin{pmatrix}-b_{y}^{\ast}\\\
b_{y}\end{pmatrix}\Big{\\}}},\end{split}$ (I.10)

with the operators $H(x,y)$, $B(x)$ in this expression defined by

$\begin{split}H(x,y)&:=\\{-\Delta_{x}+V_{\mathrm{trap}}(x)+(\upsilon_{N}\ast\rho^{\mathrm{pair}})(x)-\mu\\}\delta(x-y)+(\upsilon_{N}n)(x,y),\\\
B(x)&:=\big{(}-\Delta_{x}+V_{\mathrm{trap}}(x)+\frac{1}{N}(\upsilon_{N}\ast\rho)(x)-\mu\big{)}\phi(x)\\\
&+\int{dy\\{(\upsilon_{N}n)(x,y)\phi(y)+(\upsilon_{N}m)(x,y)\overline{\phi(y)}\\}},\end{split}$
(I.11)

for,

$\begin{split}n(x,y)&:=\phi(x)\overline{\phi(y)}+\frac{1}{N}n^{\mathrm{pair}}(x,y),\quad(\upsilon_{N}n)(x,y):=\upsilon_{N}(x-y)n(x,y)\\\
\quad
m(x,y)&:=\phi(x)\phi(y)+\frac{1}{N}m^{\mathrm{pair}}(x,y),\quad(\upsilon_{N}m)(x,y):=\upsilon_{N}(x-y)m(x,y)\\\
\rho(x)&:=n(x,x),\quad\rho^{\mathrm{pair}}(x):=n^{\mathrm{pair}}(x,x).\end{split}$
(I.12)

Observe that both the operator kernels $B(x)$ and $H(x,y)$ contain the single-
particle Schrödinger operator
$\big{(}-\Delta_{x}+V_{\mathrm{trap}}(x)-\mu\big{)}$, along with additional
terms related to the density of particles in the condensate, as well as
$n^{\mathrm{pair}},\,m^{\mathrm{pair}},\,\rho^{\mathrm{pair}}$, which,
following Griffin, are identified with marginal densities of particles outside
the condensate (equation (I.14) below). The Hartree-type equation for the
mean-field condensate $\phi(x)$ results by enforcing the condition $B(x)\equiv
0$,

$\begin{split}\mu\phi(x)&=(-\Delta_{x}+V_{\mathrm{trap}}(x))\phi(x)+(\upsilon_{N}\ast|\phi|^{2})(x)\phi(x)+\frac{1}{N}(\upsilon_{N}\ast\rho^{\mathrm{pair}})(x)\phi(x)\\\
&+\frac{1}{N}\int{dy\\{\upsilon_{N}(x-y)n^{\mathrm{pair}}(x,y)\phi(y)\\}}+\frac{1}{N}\int{dy\\{\upsilon_{N}(x-y)m^{\mathrm{pair}}(x,y)\overline{\phi(y)}\\}}.\end{split}$
(I.13)

The mean field condensate $\phi$ thus depends on the quantities
$n^{\mathrm{pair}},\,m^{\mathrm{pair}},\,\rho^{\mathrm{pair}}$. Griffin
proposes a closure to this system, i.e., a determination of all functions
$\\{\phi,\,m^{\mathrm{pair}},n^{\mathrm{pair}},\rho^{\mathrm{pair}}\\}$, via
the hypothesis that these functions are ground-state correlations in the
second-order approximation to $\mathcal{K}$. We write these correlations using
bra-ket notation:

$\langle b_{x}^{\ast}b_{y}\rangle:=n^{\mathrm{pair}}(x,y),\quad\langle
b_{x}b_{y}\rangle:=m^{\mathrm{pair}}(x,y),\quad\langle
b_{x}^{\ast}b_{x}\rangle:=\rho^{\mathrm{pair}}(x).$ (I.14)

Under this assumption, diagonalization of the effective operator provides the
algebraic closure of the system. We now briefly explain how this is done by
introducing the Fock space quasiparticle operator basis
$\alpha_{j},\,\alpha_{j}^{\ast}:\mathbb{F}\to\mathbb{F}$, for $j\in\mathbb{N}$
and a basis of single-particle states
$\\{u_{j}(x),\,p_{j}(x)\\}_{j=1}^{\infty}$ for the Hilbert space orthogonal to
the condensate, denoted $(\phi)_{\perp}$. These quantities are defined so that
the quasiparticle creation and annihilation operators
$\alpha^{\ast}_{j},\,\alpha_{j}$ satisfy the canonical commutation relations

$[\alpha_{j},\,\alpha_{k}^{\ast}]=\delta_{jk},\quad[\alpha_{j},\alpha_{k}]=[\alpha_{j}^{\ast},\alpha_{k}^{\ast}]=0,$

and such that the field operators $b_{x},b_{x}^{\ast}$ have the following
decomposition

$\begin{split}b_{x}=\sum_{j}{\\{u_{j}(x)\alpha_{j}-\overline{p_{j}(x)}\alpha_{j}^{\ast}\\}},\quad
b^{\ast}_{x}=\sum_{j}{\\{\overline{u_{j}(x)}\alpha_{j}^{\ast}-p_{j}(x)\alpha_{j}\\}}.\end{split}$
(I.15)

The commutation relations (I.8) for $b_{x},b_{x}^{\ast}$ imply the following
completeness relations for the basis of single-particle states (for all
$j,j^{\prime}\in\mathbb{N}$ and $\delta_{jj^{\prime}}$ denoting Kronecker’s
delta):

$\begin{split}\int{\rm d}x\
\\{u_{j}(x)p_{j^{\prime}}(x)-p_{j}(x)u_{j^{\prime}}(x)\\}&=0~{},\\\ \int{\rm
d}x\
\\{u_{j}(x)\overline{u_{j^{\prime}}(x)}-p_{j}(x)\overline{p_{j^{\prime}}(x)}\\}&=\delta_{jj^{\prime}}~{}.\end{split}$
(I.16)

The reader may verify that the reduced Hamiltonian becomes diagonal in the
$\alpha$-operators, provided the operator matrix equation for the basis holds:

$\int{{\rm d}y\Big{\\{}\begin{pmatrix}H(x,y)&-(\upsilon_{N}m)(x,y)\\\
\overline{(\upsilon_{N}m)(x,y)}&-H^{T}(x,y)\end{pmatrix}\begin{pmatrix}u_{j}(y)\\\
p_{j}(y)\end{pmatrix}\Big{\\}}}=E_{j}\begin{pmatrix}u_{j}(x)\\\
p_{j}(x)\end{pmatrix},\quad\|u_{j}\|^{2}-\|p_{j}\|^{2}=1.$ (I.17)

In this setting, the effective Hamiltonian $\mathcal{H}_{\mathrm{approx}}$
reads,

$\mathcal{H}_{\mathrm{approx}}=\sum_{j}{E_{j}\alpha_{j}^{\ast}\alpha_{j}},\quad
E_{j}<E_{j+1},\quad j\in\mathbb{N}.$

The normalization condition for $\\{u_{j}(x),\,p_{j}(x)\\}_{j=1}^{\infty}$ in
(I.17) is natural due to the constraint (I.16), which is a consequence of the
Fock space commutation relations. Using (I.15), Griffin finds that the
functions $n^{\mathrm{pair}}(x,y),\,m^{\mathrm{pair}}(x,y)$ have the following
expansion in the single particle basis
$\\{u_{j}(x),\,p_{j}(x)\\}_{j=1}^{\infty}$ :

$\begin{split}n^{\mathrm{pair}}(x,y)=\sum_{j}{p_{j}(x)\overline{p_{j}(y)}},\quad
m^{\mathrm{pair}}(x,y)=-\sum_{j}{u_{j}(x)\overline{p_{j}(y)}}.\end{split}$
(I.18)

Equations (I.13), (I.17), and (I.18) constitute the fundamental system of
integro-partial-differential equations which are the subject of this work. We
emphasize that (I.13) and (I.17) contain nonlinear terms which come from the
pair correlation functions $m^{\mathrm{pair}}$ and
$n^{\mathrm{pair}},\,\rho^{\mathrm{pair}}$ given in (I.18).

######  Remark I.1.

(Reduction of Fetter, 1971 fetter72 ) The above reduction can be simplified to
yield a regularized version of the approximate Hamiltonian used by Fetter in
fetter72 . Specifically, setting
$m^{\mathrm{pair}}(x,y)=n^{\mathrm{pair}}(x,y)=0$ in the above derivation
results in a reduced Hamiltonian which is quadratic in the field operators
$b_{x},\,b_{x}^{\ast}$. In fetter72 , this approximation is accomplished
without reference to two-particle correlation functions by explicitly dropping
cubic and quartic terms in $b_{x},\,b_{x}^{\ast}$ that result from applying
the Bogoliubov approximation (I.7) to $\mathcal{K}$. The Hartree-type equation
for the mean-field condensate $\phi(x)$ which results is given by

$(-\Delta_{x}+V_{\mathrm{trap}}(x))\phi(x)+(\upsilon_{N}\ast|\phi|^{2})(x)\phi(x)=\mu\phi(x).$
(I.19)

Diagonalization of the corresponding quadratic operator results in a linear
matrix system analogous to (I.17). In GrillakisMargetisSorokanich , we
identify the spectrum of this linear eigenvalue problem with quasiparticle
excitations of the Bose gas.

###  I.2 Connection to the method of pair excitation

We now describe how the infinite nonlinear coupled system (I.17) for basis
$\\{u_{j}(x),\,p_{j}(x)\\}_{j=1}^{\infty}$ is related to a single nonlinear
equation for a kernel $k(x,y)$. This is summarized in the following
Proposition:

######  Proposition 1.

  1. 1.

Let $\phi(x),\,m^{\mathrm{pair}}(x,y),\,n^{\mathrm{pair}}(x,y)$ be given as
elements of $L^{2}(\mathbb{R}^{3})$ and
$L^{2}(\mathbb{R}^{3}\times\mathbb{R}^{3})$ respectively, and let
$m(x,y),\,n(x,y)$ be defined as in (I.11). If
$k:(\phi_{\perp})\to(\phi_{\perp})$ is a symmetric Hilbert-Schmidt solution to
the following equation

$H\circ k+k\circ H^{T}+(\upsilon_{N}m)+k\circ(\upsilon_{N}\overline{m})\circ
k=0,\quad\|k\|_{\mathrm{op}}<1,$ (I.20)

then there exists a basis $\\{u_{j}(x),p_{j}(x)\\}_{j=1}^{\infty}$ for
$(\phi_{\perp})$, which solves (I.17) with normalization
$\|u_{j}\|^{2}-\|p_{j}\|^{2}=1$. The spectrum $\\{E_{j}\\}_{j=1}^{\infty}$ in
this case is positive and the following integral relation holds:

$\int{{\rm d}y\\{k(x,y)u_{j}(y)\\}}=p_{j}(x),\quad j\in\mathbb{N}.$ (I.21)

  2. 2.

Let $m^{\mathrm{pair}}(x,y),\,n^{\mathrm{pair}}(x,y)$ in
$L^{2}(\mathbb{R}^{3}\times\mathbb{R}^{3})$ be given and suppose the single-
particle states $\phi(x),\,\\{u_{j}(x),p_{j}(x)\\}_{j=1}^{\infty}$ solve the
system (I.13), and (I.17), where $\\{u_{j}\\}_{j=1}^{\infty}$is a basis of
$(\phi_{\perp})$, and furthermore suppose that the integral relation (I.21)
holds for some $k:(\phi_{\perp})\to(\phi_{\perp})$ with
$\|k\|_{\mathrm{op}}<1$ and all $j$. Then $k(x,y)$ solves the nonlinear
operator equation (I.20).

The proof of this theorem is an adaptation of the spectral theory of J-self-
adjoint operator matrices Tretter-book ; the proof is the same here, after
substituting the appropriate operators, as it was in the case
$\rho^{\mathrm{pair}},n^{\mathrm{pair}},m^{\mathrm{pair}}=0$, which was
handled in GrillakisMargetisSorokanich . We include it here for sake of
completeness.

###### Proof.

The matrix in equation (I.17) is diagonalized by an operator matrix on
$\phi_{\perp}\oplus\phi_{\perp}$ constructed using $k$ which solves (I.20). In
particular, define the operator matrices

$W:=\begin{pmatrix}\widehat{\delta}&\quad k\\\
\overline{k}&\quad\widehat{\delta}\end{pmatrix},\quad
W^{-1}=\begin{pmatrix}(\widehat{\delta}-k\circ\overline{k})^{-1}&\quad-k(\widehat{\delta}-\overline{k}\circ
k)^{-1}\\\
-\overline{k}(\widehat{\delta}-k\circ\overline{k})^{-1}&\quad(\widehat{\delta}-\overline{k}\circ
k)^{-1}\end{pmatrix}$

(recalling the definition
$\widehat{\delta}(x,y):=\delta(x-y)-\phi(x)\overline{\phi(y)}$), and

$D:=\begin{pmatrix}H^{T}+k\circ\overline{(\upsilon_{N}m)}&0\\\
0&-H-\overline{k}\circ{(\upsilon_{N}m)}\end{pmatrix}~{}.$

A direct calculation, which uses the fact that $k$ satisfies (I.20), shows
that the following holds:

$\begin{pmatrix}H^{T}+k\circ\overline{(\upsilon_{N}m)}&0\\\
0&-H-\overline{k}\circ{(\upsilon_{N}m)}\end{pmatrix}~{}=W\circ\begin{pmatrix}H&-(\upsilon_{N}m)\\\
\overline{(\upsilon_{N}m)}&-H^{T}\end{pmatrix}\circ W^{-1}.$

In order for both $W$ and $W^{-1}$ to be bounded, we require that
$(\delta-k\circ\overline{k})^{-1}$ be a bounded operator, which is true if
$\|k\|_{\mathrm{op}}<1$. The operator
$H^{T}+k\circ\overline{(\upsilon_{N}m)}$, while not self-adjoint, is similar
to the self-adjoint operator:

$(\delta-k\circ\overline{k})^{-1/2}\circ\big{(}H^{T}+k\circ\overline{(\upsilon_{N}m)}\,\big{)}\circ(\delta-k\circ\overline{k})^{1/2},$
(I.22)

(this is a consequence of the equation (I.20) for $k$). The spectrum
$\sigma\big{(}H^{T}+k\circ\overline{(\upsilon_{N}m)}\big{)}$ is therefore
real, discrete, and bounded below (as a compact perturbation of
$-\Delta+V_{\mathrm{trap}}(x)-\mu$), and the system of eigenvectors of (I.22),
denoted $\\{\eta_{j}(x)\\}_{j=1}^{\infty},$ is orthogonal and complete in the
space $(\phi)_{\perp}$. Given the solutions $\eta_{j}(x)$, we can construct
$u_{j}(x),\,p_{j}(x)$ which solve the matrix system via

$u_{j}(x)=(\delta-k\circ\overline{k})^{-1/2}(x,\eta_{j}),\quad
p_{j}(x)=k(x,u_{j}).$

The spectrum of $D$ consists of two discrete and semi-bounded parts:
$\sigma(H^{T}+k\circ\overline{(\upsilon_{N}m)})$ is bounded below (since
$\sigma(H^{T})$ is bounded below), and
$-\sigma(H^{T}+k\circ\overline{(\upsilon_{N}m)})$ is bounded above. The
normalization condition for $\\{u_{j},p_{j}\\}$, which comes from the fact
that these functions represent Bogoliubov rotations, implies that the
Bogoliubov spectrum of the quadratic Hamiltonian consists of one part only,
namely
$\sigma(\mathcal{K})=\sigma\big{(}H^{T}+k\circ(\upsilon_{N}\overline{m})\big{)}$.
∎

######  Corollary 1.

Using (I.18) and (I.21), the following formulas hold for the density
operators:

$\begin{split}m^{\mathrm{pair}}(x,y)&=-\sum_{j}{u_{j}(x)\overline{p_{j}(y)}}=\Big{(}k\circ(\delta-\overline{k}\circ
k)^{-1}\Big{)}(x,y)\\\
n^{\mathrm{pair}}(x,y)&=\sum_{k}{p_{j}(x)\overline{p_{j}(y)}}=\Big{(}(k\circ\overline{k})\circ(\delta-\overline{k}\circ
k)^{-1}\Big{)}(x,y).\end{split}$ (I.23)

This Corollary shows that the terms of (I.20) which involve the density
operators $m^{\mathrm{pair}},\,n^{\mathrm{pair}}$ can be rewritten exclusively
using the kernel $k$ when system (I.17) is satisfied. In particular, the
operator kernels $H(x,y)$ and $(\upsilon_{N}m)(x,y)$ depend on the quantities
$\phi(x),\,k(x,y)$. Motivated by Corollary 1, we define:

$\begin{split}&H[\phi,k;\mu](x,y):=\big{(}-\Delta_{x}+V_{\mathrm{trap}}(x)-\mu\big{)}\delta(x-y)\\\
&\quad\quad\quad\quad\quad\quad+\overline{\phi(y)}\upsilon_{N}(x-y)\phi(x)+(\upsilon_{N}*|\phi|^{2})(y)\delta(x-y)\\\
&+\frac{1}{N}\upsilon_{N}(x-y)\Big{\\{}\frac{k\circ\overline{k}}{\delta-k\circ\overline{k}}\Big{\\}}(x,y)+\frac{1}{N}\int{{\rm
d}x^{\prime}\upsilon_{N}(x-x^{\prime})\Big{\\{}\frac{k\circ\overline{k}}{\delta-k\circ\overline{k}}\Big{\\}}(x^{\prime},x^{\prime})}\delta(x-y).\end{split}$
(I.24)

We also define the operator kernel $\Theta[\phi,k](x,y)$,

$\begin{split}\Theta[\phi,k](x,y)&:=\upsilon_{N}(x-y)\Big{\\{}\phi(x)\phi(y)+\frac{1}{N}\big{(}k\circ(\delta-\overline{k}\circ
k)^{-1}\big{)}(x,y)\Big{\\}},\end{split}$ (I.25)

so that the operator equation (I.20) for $k$ under this condition reads:

$H[\phi,k;\mu]\circ k+k\circ
H^{T}[\phi,k;\mu]+\Theta[\phi,k]+k\circ\overline{\Theta[\phi,k]}\circ
k=0,\quad k:(\phi)_{\perp}\to(\phi)_{\perp}.$ (I.26)

######  Remark I.2.

  1. 1.

For simplicity, we have written the equation for $k$ as it holds on the space
$\phi_{\perp}$. This allows us to avoid including the Lagrange multiplier term
$\lambda\otimes_{\mathrm{s}}\phi$ in the equation (I.20) for $k$. Also, the
fraction notation $k\circ\overline{k}/(\delta-k\circ\overline{k})$ is
justified here because of the commutation

$k\circ\overline{k}\circ(\delta-k\circ\overline{k})^{-1}=(\delta-k\circ\overline{k})^{-1}\circ
k\circ\overline{k}.$

  2. 2.

An alternative approach to the approximation of $\mathcal{K}$, which also
results in the equation (I.26) for the kernel $k(x,y)$, involves the
transformation of $\mathcal{K}$ by the non-Hermitian pair-excitation operator,
which takes the form

$\mathcal{K}\mapsto\exp\big{(}P[k]\big{)}\mathcal{K}\exp\big{(}-P[k]\big{)},\quad
P[k]:=-\frac{1}{2}\iint{{\rm d}x{\rm d}yk(x,y)b_{x}^{\ast}b_{y}^{\ast}}.$

The transformation, and therefore $k(x,y)$, is chosen so that there are no
terms in the resulting non-Hermitian operator which are proportional to the
product $b_{x}b_{y}$. This criterion for determining $k$ is attributed to Wu
wu58 ; wu61 ; wu98 . Non-Hermitian quantum mechanical Hamiltonians with real
spectra have more recently been studied by Bender Bender2007 .

##  II Main results

Proposition 1 and its Corollary show that the system (I.13), (I.17), (I.18)
can be recast as a coupled system for the quantities $\phi(x)$ and $k(x,y)$,
which is the subject of the current analysis. We adopt a variational approach
to find a unique solution to this system.

The following function space notation will be used throughout the rest of the
work. For trapping potential $V_{\mathrm{trap}}(x)$, we define the space

$L^{2}_{V}(\mathbb{R}^{d})=\big{\\{}f(x):\,\int{{\rm
d}x\\{V_{\mathrm{trap}}(x)|f(x)|^{2}\\}}<\infty\big{\\}}.$

Function spaces on $\mathbb{R}^{d}$ (e.g., $d=3$) are denoted by lowercase
gothic letters, viz.,

$\mathfrak{h}(\mathbb{R}^{d}):=L^{2}(\mathbb{R}^{d})~{},\quad\mathfrak{h}^{1}(\mathbb{R}^{d}):=H^{1}(\mathbb{R}^{d})~{},\quad\mathfrak{h}_{V}^{1}(\mathbb{R}^{d}):=H^{1}(\mathbb{R}^{d})\cap
L_{V}^{2}(\mathbb{R}^{d})~{}.$

We write $\mathfrak{h}$, $\mathfrak{h}^{1}$, $\mathfrak{h}_{V}^{1}$ for these
spaces if $d=3$. As an exception to this notation, as before, we define
$(\phi_{\perp}):=\big{\\{}e\in\mathfrak{h}^{1}_{V}\ \big{|}\
e\perp\phi\big{\\}}$ where $\phi\in\mathfrak{h}^{1}$ is the condensate wave
function. The space of bounded linear operators on $\mathfrak{h}$ is denoted
$\mathfrak{B}(\mathfrak{h})$, with norm $\|\cdot\|_{\mathrm{op}}$. Also, the
space of trace-class operators on $\mathfrak{h}$ is denoted
$\mathfrak{B}_{1}(\mathfrak{h})$ with norm

$\|A\|_{\mathfrak{B}_{1}(\mathfrak{h})}=\|A\|_{1}=\mathrm{tr}|A|~{},\quad\forall\,\,A\in\mathfrak{B}_{1}(\mathfrak{h})~{}.$

Similarly, the space of Hilbert-Schmidt operators on $\mathfrak{h}$ is
$\mathfrak{B}_{2}(\mathfrak{h})$ with norm

$\|A\|_{\mathfrak{B}_{2}(\mathfrak{h})}=\|A\|_{2}=(\mathrm{tr}|A^{\ast}A|)^{1/2}~{},\quad\forall\,\,A\in\mathfrak{B}_{2}(\mathfrak{h})~{}.$

The space of compact operators on $\mathfrak{h}$ is
$\mathfrak{B}_{0}(\mathfrak{h})$. Note the inequalities

$\|A\|_{\mathrm{op}}\leq\|A\|_{2}\leq\|A\|_{1}~{},$

and the inclusions
$\mathfrak{B}_{1}(\mathfrak{h})\subseteq\mathfrak{B}_{2}(\mathfrak{h})\subseteq\mathfrak{B}_{0}(\mathfrak{h})\subseteq\mathfrak{B}(\mathfrak{h})$.

###  II.1 Main theorem and variational formulation

######  Theorem II.1.

Suppose the spectrum of the single-particle operator
$\big{(}-\Delta+V_{\mathrm{trap}}(x)\big{)}$ contains a gap between its ground
state and first excited state. Also, suppose the interaction potential is
given by $\upsilon_{N}=N^{3\beta}\upsilon_{N}(N^{\beta}x)$ for $0\leq\beta\leq
1$ where $\upsilon(x)$ is smooth and bounded, and $\widehat{\upsilon}\geq 0$,
where $\widehat{\upsilon}$ is the Fourier transform of $\upsilon$. Then for
$g>0$ small enough and $0\leq\beta\leq 1$, there exist
$\phi\in\mathfrak{h}^{1}_{V}$, with $\|\phi\|_{L^{2}}=1$, and
$k\in\mathfrak{B}_{2}(\phi_{\perp},\phi_{\perp})$, with
$\|k\|_{\mathrm{op}}<1$ which solve the coupled system consisting of the
equation for $\phi$:

$\begin{split}(&-\Delta_{x}+V_{\mathrm{trap}}(x))\phi(x)+(\upsilon_{N}\ast|\phi|^{2})(x)\phi(x)\\\
&+\frac{1}{N}\big{(}\upsilon_{N}\ast\Big{\\{}\frac{k\circ\overline{k}}{\delta-k\circ\overline{k}}\Big{\\}}(x^{\prime},x^{\prime})\big{)}\phi(x)+\frac{1}{N}\int{dy\\{\upsilon_{N}(x-y)\Big{\\{}\frac{k\circ\overline{k}}{\delta-k\circ\overline{k}}\Big{\\}}(x,y)\phi(y)\\}}\\\
&+\frac{1}{N}\int{dy\\{\upsilon_{N}(x-y)\Big{\\{}\frac{k}{\delta-k\circ\overline{k}}\Big{\\}}(x,y)\overline{\phi(y)}\\}}=\mu\phi(x),\end{split}$
(II.1)

and the nonlinear operator equation for
$k(x,y):(\phi)_{\perp}\to(\phi)_{\perp},$

$H[\phi,k;\mu]\circ k+k\circ
H^{T}[\phi,k;\mu]+\Theta[\phi,k]+k\circ\overline{\Theta[\phi,k;\mu]}\circ
k=0,$ (II.2)

where the operators $H[\phi,k;\mu](x,y)$ and $\Theta[\phi,k](x,y)$ are given
by formulas (LABEL:gN) and (I.25) respectively.

The functions $\phi$ and $k$ are sought as minimizers of the total energy,
$\mathcal{E}_{\mathrm{tot}}[\phi,k]$, which is now introduced.

######  Definition II.2.

The Hartree functional, denoted $E_{\mathrm{H}}$, is defined for all
$\phi\in\mathfrak{h}^{1}_{V}$, via

$E_{\mathrm{H}}[\phi]:=\iint{{\rm d}x{\rm
d}y\big{\\{}\overline{\phi(x)}\big{(}-\Delta+V_{\mathrm{trap}}(x)\big{)}\delta(x-y)\phi(y)+\frac{1}{2}|\phi(x)|^{2}\upsilon_{N}(x-y)|\phi(y)|^{2}\big{\\}}}.$
(II.3)

Given the operators
$H[\phi,k],\,\Theta[\phi,k]:\mathfrak{h}_{V}^{1}\to\mathfrak{h}^{1}_{V}$
(defined by equations (LABEL:gN), (I.25)), let

$H_{0}[\phi,k]:=H[\phi,k;\mu]+\mu\delta(x-y),$ (II.4)

(so that $H_{0}[\phi,k]$ contains no contribution proportional to $\mu$).
Given $\phi\in\mathfrak{h}_{V}^{1}$, the functional

$\mathcal{E}\big{[}k;H_{0}[\phi,k],\Theta[\phi,k]\big{]},\quad\mathrm{dom}(\mathcal{E}):=\big{\\{}k\in\mathfrak{B}_{2}(\mathfrak{h}^{1}_{V}),\quad\|k\|_{\mathrm{op}}<1\big{\\}},$

is defined by the operator trace formula,

$\mathcal{E}[k;\,H_{0},\Theta]:=\mathrm{tr}\Big{\\{}(\delta-\overline{k}\circ
k)^{-1}\circ\Big{(}\overline{k}\circ H_{0}[\phi,k]\circ
k+\frac{1}{2}\overline{k}\circ\Theta[\phi,k]+\frac{1}{2}\overline{\Theta[\phi,k]}\circ
k\Big{)}\Big{\\}}.$ (II.5)

The total energy, $\mathcal{E}_{\mathrm{tot}}$ is defined for all
$\phi\in\mathfrak{h}_{V}^{1}$ and
$\big{\\{}k\in\mathfrak{B}_{2}(\mathfrak{h}_{V}^{1}),\quad\|k\|_{\mathrm{op}}<1\big{\\}}$
by the sum

$\mathcal{E}_{\mathrm{tot}}[\phi,k;N]:=\Big{\\{}E_{\mathrm{H}}[\phi]+\frac{1}{N}\mathcal{E}\big{[}k;\,H_{0}[\phi,k],\Theta[\phi,k]\big{]}\Big{\\}},$
(II.6)
$\mathrm{dom}(\mathcal{E}_{\mathrm{tot}})=\\{\phi\in\mathfrak{h}^{1}_{V},\quad
k\in\mathrm{dom}(\mathcal{E})\\}.$

It is worth commenting on the manner in which the energy
$\mathcal{E}[k;H_{0},\Theta]$ is written in equation (II.6). Specifically, it
might be claimed that our use of operators $H_{0}[\phi,k]$ and
$\Theta[\phi,k]$ obscure important nonlinear dependences of the energy on $k$.
We choose to write $\mathcal{E}$ in this way in order to draw the connection
to the variational method of our previous work GrillakisMargetisSorokanich ,
in which the analogues of operators $H_{0},\,\Theta$ did not depend on the
argument $k$.

The following lemma describes how the coupled system (II.1) and (II.2) can be
derived from the total energy functional. The proof is direct, but rather
involved. It is therefore omitted.

######  Lemma II.3.

Let $H[\phi,k;\mu],\Theta[\phi,k]$ be given by (LABEL:gN), (I.25), and let the
total energy $\mathcal{E}_{\mathrm{tot}}$ be given by formula (II.6), i.e.,

$\mathcal{E}_{\mathrm{tot}}[\phi,k;N]:=\Big{\\{}E_{\mathrm{H}}[\phi]+\frac{1}{N}\mathcal{E}\big{[}k;\,H_{0}[\phi,k],\Theta[\phi,k]\big{]}\Big{\\}}.$

Then (i) For fixed $k,\,\overline{k}$, critical points of the total energy
with respect to constrained variations $\overline{\phi}\in
L^{2}(\mathbb{R}^{3})$, $\|\phi\|^{2}=1$, i.e., points satisfying
$\,\delta\mathcal{E}_{\mathrm{tot}}/\delta\overline{\phi}=\mu\phi\,$, where
$\mu$ is the Lagrange multiplier enforcing the constraint, are solutions to
the nonlinear Hartree-type equation (II.1).

(ii) For fixed $\phi$ and $\mu$ satisfying (II.1), critical points of

$\mathcal{E}_{\mathrm{tot}}+\mu\,\mathrm{tr}\Big{\\{}(\delta-\overline{k}\circ
k)^{-1}\circ\Big{(}\overline{k}\circ
k\Big{)}\Big{\\}}:=E_{\mathrm{H}}[\phi]+\frac{1}{N}\mathcal{E}\big{[}k;\,H[\phi,k;\mu],\Theta[\phi,k]\big{]},$

with respect to symmetric variations of
$\overline{k}\in\mathfrak{B}_{2}(\mathfrak{h}_{V})$, i.e., points satisfying

$\delta/\delta\overline{k}\Bigg{[}\mathcal{E}_{\mathrm{tot}}+\mu\,\mathrm{tr}\Big{\\{}(\delta-\overline{k}\circ
k)^{-1}\circ\Big{(}\overline{k}\circ k\Big{)}\Big{\\}}\Bigg{]}=0,$

are solutions to the nonlinear equation (II.2).

######  Remark II.4.

Integrating (II.1) against $\overline{\phi(x)}$ gives the formula for $\mu$:

$\begin{split}\mu&=\iint{{\rm d}x{\rm
d}y\,\overline{\phi(x)}\Big{(}-\Delta_{x}+V_{\mathrm{trap}}(x)\Big{)}\delta(x-y)\phi(y)+|\phi(x)|^{2}\upsilon_{N}(x-y)|\phi(y)|^{2}}\\\
&+\frac{1}{N}\iint{{\rm d}x{\rm
d}y\,\Big{\\{}\frac{k\circ\overline{k}}{\delta-\overline{k}\circ
k}\Big{\\}}(y,y)\upsilon_{N}(x-y)|\phi(x)|^{2}}\\\ &+\frac{1}{N}\iint{{\rm
d}x{\rm d}y\Big{\\{}\frac{k\circ\overline{k}}{\delta-\overline{k}\circ
k}\Big{\\}}(x,y)\upsilon_{N}(x-y)\overline{\phi(x)}\phi(y)}\\\
&+\frac{1}{N}\iint{{\rm d}x{\rm d}y\Big{\\{}\frac{k}{\delta-\overline{k}\circ
k}\Big{\\}}(x,y)\upsilon_{N}(x-y)\overline{\phi(x)}\,\,\overline{\phi(y)}}.\end{split}$
(II.7)

###  II.2 Proof of Theorem II.1 by direct minimization

The existence of a weak minimizer of $\mathcal{E}_{\mathrm{tot}}$, denoted by
the triple $\\{\phi(x),k(x,y);\mu\\}$, can be demonstrated by showing that the
functional $\mathcal{E}_{\mathrm{tot}}$ is bounded below, and extracting the
weak limit of a minimizing sequence. The proof which follows applies to the
full range $0\leq\beta\leq 1$ of the mean-field interaction potential.
However, the method of proof does not give much information about the nature
of the minimizer. It also does not deal directly with the quantity $k(x,y)$,
but rather a kernel denoted $\psi(x,y)$, which is an analytic function of
$k(x,y)$.

###### Proof.

Given $k\in\mathfrak{B}_{2}(\mathfrak{h}^{1}_{V}(\phi_{\perp}))$, with
$\|k\|_{\mathrm{op}}<1$, we define the operator kernel $\psi(x,y)$ via a
change-of-variables

$\psi(x,y):=\Big{(}k\circ{(\delta-k\circ\overline{k})^{-1/2}}\Big{)}(x,y).$
(II.8)

Inverting this relationship gives

$k=\psi\circ(\delta+\overline{\psi}\circ\psi)^{-1/2}=(\delta+\psi\circ\overline{\psi})^{-1/2}\circ\psi,$
(II.9)

and observe that $\delta-\overline{k}\circ
k=(\delta+\overline{\psi}\circ\psi)^{-1}.$ This constitutes a change of
variables for minimization problem; the above relationship implies that we can
solve the problem for $\psi$ and transform the solution to a solution $k$. The
functional $\mathcal{E}_{\mathrm{tot}}\big{[}\phi,k\big{]}$ defined in (II.6),
whose minimization produces the nonlinear system for $\phi$ and $k$ is
rewritten using the operator $\psi$. It can be verified that
$\mathcal{E}_{\mathrm{tot}}[\phi,k]=\widetilde{\mathcal{E}}[\phi,\psi]$, where

$\begin{split}\widetilde{\mathcal{E}}\big{[}\phi,\psi\big{]}&:=N\int{{\rm
d}x\\{|\nabla\phi(x)|^{2}+V_{\mathrm{trap}}(x)|\phi(x)|^{2}\\}}\\\
&+N\iint{{\rm d}x{\rm
d}y\\{\upsilon_{N}(x-y)\big{|}\phi(x)\phi(y)+\frac{1}{N}\psi\circ(\delta+\overline{\psi}\circ\psi)^{1/2}(x,y)\big{|}^{2}\\}}\\\
&+\frac{1}{2}\iint{{\rm d}x{\rm
d}y\\{|\nabla_{x}\psi(x,y)|^{2}+|\nabla_{y}\psi(x,y)|^{2}+(V_{\mathrm{trap}}(x)+V_{\mathrm{trap}}(y))|\psi(x,y)|^{2}\\}}\\\
&+\frac{1}{2}\iint{{\rm d}x{\rm
d}y\\{\upsilon_{N}(x-y)[\phi(x)\overline{\phi(y)}(\psi\circ\overline{\psi})(x,y)+\overline{\phi(x)}\phi(y)(\overline{\psi}\circ\psi)(x,y)]\\}}\\\
&+\iint{{\rm d}x{\rm
d}y\\{|\phi(x)|^{2}(\overline{\psi}\circ\psi)(y,y)+|\phi(y)|^{2}(\overline{\psi}\circ\psi)(x,x)\\}}\\\
&+\frac{1}{2N}\iint{{\rm d}x{\rm
d}y\\{\upsilon_{N}(x-y)\big{(}|(\overline{\psi}\circ\psi)(x,y)|^{2}+(\overline{\psi}\circ\psi)(x,x)(\overline{\psi}\circ\psi)(y,y)\big{)}\\}}.\end{split}$
(II.10)

Recall that we assume $\widehat{\upsilon}\geq 0$ as in the statement of the
Theorem. This is not necessary but it simplifies the argument and makes more
transparent the main idea.

The constraint for the minimization is modified in a way that accommodates
both the constraints on $\phi$ and $\psi$. Namely, we solve

$\min\widetilde{\mathcal{E}}[\phi,\psi],\quad\mathrm{over}\quad\big{\\{}\phi\in\mathfrak{h}^{1}_{V},\,\psi\in\mathfrak{B}_{2}(\mathfrak{h}^{1}_{V}),\quad\|\phi\|_{L^{2}(\mathbb{R}^{3})}^{2}+\frac{1}{N}\|\psi\|_{L^{2}(\mathbb{R}^{3}\times\mathbb{R}^{3})}^{2}=1\big{\\}}.$

In the new coordinate system the system can be written

$\begin{split}\frac{\delta\widetilde{\mathcal{E}}}{\delta\overline{\phi}}&=\mu\phi\\\
\frac{\delta\widetilde{\mathcal{E}}}{\delta\overline{\psi}}&=\mu\psi.\end{split}$

We will prove that the minimum exists and is non-trivial. Let us define the
norm for $\psi(x,y)$,

$\|\psi\|^{2}_{H^{1}_{\mathrm{trap}}}:=\iint{{\rm d}x{\rm
d}y\\{|\nabla_{x}\psi(x,y)|^{2}+|\nabla_{y}\psi(x,y)|^{2}+(V_{\mathrm{trap}}(x)+V_{\mathrm{trap}}(y))|\psi(x,y)|^{2}\\}}.$

We know by Sobolev inequalities in $\mathbb{R}^{6}$,

$\begin{split}\|\psi\|_{L^{3}(\mathbb{R}^{3}\times\mathbb{R}^{3})}&\leq
C\|\psi\|_{H^{1}(\mathbb{R}^{3}\times\mathbb{R}^{3})},\\\
\|\psi\|_{L^{2}(\mathbb{R}^{3}\times\mathbb{R}^{3})}&\leq
C\|\psi\|_{H^{1}_{\mathrm{trap}}(\mathbb{R}^{3}\times\mathbb{R}^{3})}.\end{split}$
(II.11)

Moreover, the embedding $L^{p}(\mathbb{R}^{3}\times\mathbb{R}^{3})\subset
H^{1}_{\mathrm{trap}}(\mathbb{R}^{3}\times\mathbb{R}^{3})$ is compact for
$2\leq p\leq 3$. Defining the density
$\rho^{\mathrm{pair}}(x):=(\overline{\psi}\circ\psi)(x,x)$, we have

$\|\nabla\rho^{\mathrm{pair}}\|_{L^{2}(\mathbb{R}^{3})}^{3/2}\leq
C\|\psi\|^{3/2}_{H^{1}(\mathbb{R}^{3}\times\mathbb{R}^{3})}\|\rho^{\mathrm{pair}}\|_{L^{3}(\mathbb{R}^{3})}^{3/4}$

(see Lemma III.5 below for a detailed proof of this). Sobolev again implies

$\|\rho^{\mathrm{pair}}\|_{L^{3}(\mathbb{R}^{3})}\leq
C\|\nabla\rho^{\mathrm{pair}}\|_{W^{1,3/2}(\mathbb{R}^{3})}.$

We also know

$\|\rho^{\mathrm{pair}}\|_{L^{1}(\mathbb{R}^{3})}=\|\psi\|_{L^{2}(\mathbb{R}^{3}\times\mathbb{R}^{3})}^{2}.$

We have the inequalities:

$|(\overline{\psi}\circ\psi)(x,y)|\leq\sqrt{\rho^{\mathrm{pair}}(x)\rho^{\mathrm{pair}}(y)}$

and

$|\nabla_{x}(\overline{\psi}\circ\psi)(x,y)|\leq\Big{(}\int{{\rm
d}y|\nabla_{x}\psi(x,y)|^{2}}\Big{)}^{1/2}(\rho^{\mathrm{pair}}(y))^{1/2}$

hence:

$\|\overline{\psi}\circ\psi\|_{H^{1}(\mathbb{R}^{3}\times\mathbb{R}^{3})}\leq\|\psi\|^{2}_{H^{1}_{\mathrm{trap}}(\mathbb{R}^{3}\times\mathbb{R}^{3})}$
$\|\overline{\psi}\circ\psi\|_{L^{3}(\mathbb{R}^{3}\times\mathbb{R}^{3})}\leq\|\psi\|^{2}_{H^{1}_{\mathrm{trap}}(\mathbb{R}^{3}\times\mathbb{R}^{3})}.$

Recall the Hartree energy $E_{\mathrm{H}}[\phi]$, and notice that
$\widetilde{\mathcal{E}}(\phi,0)=NE_{\mathrm{H}}[\phi]$. The trapping norm for
$\phi$ is

$\|\phi\|_{H^{1}_{\mathrm{trap}}(\mathbb{R}^{3})}^{2}=\int{{\rm
d}x\\{|\nabla_{x}\phi(x)|^{2}+V_{\mathrm{trap}}(x)|\phi(x)|^{2}\\}}$

Notice that if we set $\rho^{\mathrm{cond}}(x):=|\phi(x)|^{2}$ for the
condensate density then

$\|\rho^{\mathrm{cond}}\|_{L^{3}(\mathbb{R}^{3})}\leq
C\|\phi\|_{H^{1}_{\mathrm{trap}}(\mathbb{R}^{3})}.$

Now we have all the ingredients in order to prove the existence of a minimum.
We start with a minimizing sequence $(\phi^{(j)},\psi^{(j)})$ which without
loss of generality we may assume converges weakly to
$(\phi_{\mathrm{min}},\psi_{\mathrm{min}})$ in
$H^{1}_{\mathrm{trap}}(\mathbb{R}^{3})\times
H^{1}_{\mathrm{trap}}(\mathbb{R}^{3}\times\mathbb{R}^{3})$. We know that

$N\|\phi_{\mathrm{min}}\|_{H^{1}_{\mathrm{trap}}(\mathbb{R}^{3})}+\frac{1}{2}\|\psi_{\mathrm{min}}\|_{H^{1}_{\mathrm{trap}}(\mathbb{R}^{3}\times\mathbb{R}^{3})}\leq
NE_{\mathrm{H}},$

where $E_{\mathrm{H}}:=\min_{\|\phi\|=1}E_{\mathrm{H}}[\phi]$. We also know
that $E_{\mathrm{H}}$ is finite and this fact implies that the
$H^{1}_{\mathrm{trap}}$ norms of $(\phi^{(j)},\psi^{(j)})$ are bounded (for
fixed $N$). Here we used the fact that the sum of lines 4 and 5 in the energy
functional $\widetilde{\mathcal{E}}$ is positive. The compact embedding
implies that

$\phi^{(j)}\to\phi_{\mathrm{min}},\quad\mathrm{strongly}\,\,\mathrm{in}\quad
L^{2}(\mathbb{R}^{3}).$

This guarantees that the solution is non-trivial. The weak convergence is in
fact strong, i.e.,

$\phi^{(j)}\to\phi_{\mathrm{min}},\quad\mathrm{strongly}\,\,\mathrm{in}\quad
H^{1}_{\mathrm{trap}}(\mathbb{R}^{3}),$

and

$\psi^{(j)}\to\psi_{\mathrm{min}},\quad\mathrm{strongly}\,\,\mathrm{in}\quad
H^{1}_{\mathrm{trap}}(\mathbb{R}^{3}\times\mathbb{R}^{3}).$

We explain this conclusion for $\psi^{(j)}$ (the conclusion for $\phi^{(j)}$
follows by a similar argument). The role of $V_{\mathrm{trap}}$ is to
guarantee that $\|\psi_{\mathrm{min}}\|_{L^{2}}$ is dominated by
$\|\psi_{\mathrm{min}}\|_{H^{1}_{\mathrm{trap}}}$. In fact, all the terms in
$\widetilde{\mathcal{E}}[\phi^{(j)},\psi^{(j)}]$ converge strongly except
those involving $\|\phi^{(j)}\|_{H^{1}_{\mathrm{trap}}}$ and
$\|\psi^{(j)}\|_{H^{1}_{\mathrm{trap}}}$. For these terms we have, e.g.,

$\|\psi_{\mathrm{min}}\|^{2}_{H^{1}_{\mathrm{trap}}(\mathbb{R}^{3}\times\mathbb{R}^{3})}\leq\liminf_{j}\|\psi^{(j)}\|^{2}_{H^{1}_{\mathrm{trap}}(\mathbb{R}^{3}\times\mathbb{R}^{3})}.$

If the inequality were strict it would contradict the minimization because,
e.g.,

$\widetilde{\mathcal{E}}[\phi,\psi_{\mathrm{min}}]<\min_{\psi}\widetilde{\mathcal{E}}[\phi,\psi]$

Therefore,

$\|\psi_{\mathrm{min}}\|_{H^{1}_{\mathrm{trap}}}^{2}=\liminf_{j}\|\psi^{(j)}\|_{H^{1}_{\mathrm{trap}}}^{2}$

which implies up to a subsequence $\psi^{(j)}\to\psi_{\mathrm{min}}$ strongly
in $H^{1}_{\mathrm{trap}}(\mathbb{R}^{3}\times\mathbb{R}^{3})$. Thus the
minimum is attained at $\psi_{\mathrm{min}}$ and therefore it is a weak
solution. ∎

##  III Iterative Minimization

The starting point of the iteration scheme is the system
$\\{\phi_{0},\,k_{0}=0;\,\mu_{0}\\}$ which is the solution to (II.1), (II.2)
by forcing pair, $\mathrm{pair},\,n^{\mathrm{pair}}\equiv 0$. In particular,
the equation for $\phi_{0}$ is the usual Hartree equation. The equation for
$k_{1}$ is the Riccati equation which we studied extensively in our previous
work GrillakisMargetisSorokanich . The existence result for this method is
stated in Theorem III.1. The solution discovered by iteration corresponds to
the restricted range of values $0<\beta\leq 1/6$. The result nonetheless seems
interesting in its own right, since the kernel found by this method will be a
small perturbation of the kernel $k_{1}$, which has been studied before.

######  Theorem III.1.

Suppose the single-particle operator
$\big{(}-\Delta+V_{\mathrm{trap}}(x)\big{)}$ contains a gap between its ground
state and first excited state, and the interaction potential is given by
$\upsilon_{N}=N^{3\beta}\upsilon_{N}(N^{\beta}x)$ for $0\leq\beta\leq 1/6$
where $\upsilon(x)$ is smooth and bounded. Also assume that
$\widehat{\upsilon}\geq 0$, where $\widehat{\upsilon}$ is the Fourier
transform of $\upsilon$. Then for $g>0$ small enough, there exist functions
$\phi\in\mathfrak{h}^{1}_{V}(\phi_{\perp})$, with $\|\phi\|_{L^{2}}=1$, and
$k\in\mathfrak{B}_{2}\big{(}\mathfrak{h}^{1}_{V}(\phi_{\perp})\big{)}$, with
$\|k\|_{\mathrm{op}}<1$ which solve the coupled system (II.1), (II.2), for
$\mu>0$.

The iteration provides solutions $\\{\phi_{n},k_{n};\mu_{n}\\}_{n=0}^{\infty}$
to a sequence of simpler minimization problems. It will be shown that this
sequence converges strongly (in $\mathfrak{h}$ and
$\mathfrak{B}(\mathfrak{h})$) to weak solutions $\\{\phi,\,k;\mu\\}$ of the
nonlinear system. To this end, define the following shorthand for operators in
this sequence:

$H[\phi_{n},k_{n};\mu_{n}]:=H[n],$ (III.1)

as well as

$\Theta[\phi_{n},k_{n}]:=\Theta[n],$ (III.2)

where $H[\phi,k;\mu],\Theta[\phi,k]$ are given by the formulas (LABEL:gN),
(I.25). The iteration scheme is now summarized:  Iterative Minimization for
$\\{\phi,\,k;\,\mu\\}$:

  * Zero iterate: Set $k_{0}=0$, and solve the constrained minimization problem

$\min_{\|\phi\|=1}E_{\mathrm{H}}[\phi],\quad\mathrm{for}\quad\phi_{0}\in\mathfrak{h}.$

This entails that $\phi_{0}$ satisfies $\delta
E_{\mathrm{H}}/\delta\overline{\phi}=\mu_{0}\phi$, which gives the Hartree
equation

$\Big{(}-\Delta_{x}+V_{\mathrm{trap}}(x)\Big{)}\phi_{0}(x)+(\upsilon_{N}\ast|\phi_{0}|^{2})(x)\phi_{0}(x)=\mu_{0}\phi_{0}(x).$

For
$\epsilon_{0}(x,y):=\big{(}-\Delta_{x}+V_{\mathrm{trap}}(x)\big{)}\delta(x-y)$,
the constant $\mu_{0}$ is given by (II.7),

$\mu_{0}=\iint{{\rm d}x{\rm
d}y\\{\overline{\phi_{0}(x)}\epsilon_{0}(x,y)\phi_{0}(y)+|\phi_{0}(x)|^{2}\upsilon_{N}(x-y)|\phi_{0}(y)|^{2}\\}}.$
(III.3)

  * ($n+1$) iterate: Assume that the functions $\\{\phi_{n},\,k_{n}\\}$ as well as the constant $\mu_{n}$ are given. The operators $H[n](x,y),\,\Theta[n](x,y)$ are then defined by the formulas (III.1) and (III.2). We solve the variational problem for $k=k_{n+1}$:

$\min_{\|k\|_{\mathrm{op}}<1}\Big{\\{}E_{\mathrm{H}}[\phi_{n}]+\frac{1}{N}\mathcal{E}\big{[}k;\,\,H[n],\,\Theta[n]\big{]},\quad
k\in\mathfrak{B}_{2}(\mathfrak{h}^{1}_{V})\Big{\\}},$ (III.4)

and next the variational problem for $\phi=\phi_{n+1}$:

$\min_{\|\phi\|=1}\Big{\\{}E_{\mathrm{H}}[\phi]+\frac{1}{N}\mathcal{E}\big{[}k_{n+1};\,H_{0}[\phi,k_{n+1}],\,\Theta[\phi,k_{n+1}]\big{]},\quad\phi\in\mathfrak{h}^{1}_{V}\Big{\\}}.$
(III.5)

The minimization (III.4) implies that
$k_{n+1}:(\phi_{n})_{\perp}\to(\phi_{n})_{\perp}$ solves the Riccati equation

$H[n]\circ k_{n+1}+k_{n+1}\circ
H[n]^{T}+\Theta[n]+k\circ\overline{\Theta[n]}\circ k=0,$ (III.6)

and the minimization (III.5) yields the function $\phi_{n+1}$ satisfying a
Hartree-type equation of the form (II.1) under the substitutions
$\phi\equiv\phi_{n+1},\,k\equiv k_{n+1}$, and the constant
$\mu\equiv\mu_{n+1}$ given by (II.7).

######  Remark III.2.

(Solutions to the Riccati equation for $k_{n}$). Suppose the functions
$\phi_{n},\,k_{n}$ which solve (III.5), (III.4) are given, in addition to the
constant $\mu_{n}$ and the operators $H[n],\,\Theta[n]$. We review some
conclusions about minimizers of the energy functional
$\mathcal{E}\big{[}k;\,H[n],\Theta[n]\big{]}$ by means of which $k_{n+1}$ is
generated. Namely, if the following operator estimate holds for all
$e\in(\phi_{n})_{\perp}$:

$\Big{\\{}H[n](e,\overline{e})-\big{|}\Theta[n](\overline{e},\overline{e})\big{|}\Big{\\}}\geq
c\|e\|^{2},\quad c>0,\quad e\in(\phi_{n})_{\perp},$ (III.7)

then a Hilbert-Schmidt minimizer of $\mathcal{E}$, denoted
$k_{n+1}\in\mathfrak{B}_{2}(\mathfrak{h}^{1}_{V})$, with
$\|k_{n+1}\|_{\mathrm{op}}<1$, exists. The proof of this statement consists in
constructing an orthonormal basis
$\\{e_{j}(n+1)\\}_{j=1}^{\infty}\subset\mathfrak{h}_{V}^{1}$ of
$(\phi_{n})_{\perp}$, and a sequence
$\\{z_{j}(n+1)\\}_{j=1}^{\infty}\subset\ell^{2}$, such that

$k_{n+1}(x,y)=\sum_{j}{z_{j}(n+1)\\{e_{j}(x)e_{j}(y)\\}}.$ (III.8)

(We will refer to the basis as $\\{e_{j}\\}_{j=1}^{\infty}$ since its
dependence on $n+1$ is clear from the context and does not play a role in the
current argument). Plugging this expansion into the trace formula which
defines $\mathcal{E}$, and setting
$\partial\mathcal{E}/\partial{\overline{z_{j}}}=0,$ yields two possible
choices of coefficient, $z_{j}^{(\pm)}(n+1)$ for every $j$, where
$|z_{j}^{+}(n+1)|<1$ and $|z_{j}^{-}(n+1)|>1$. The root $z_{j}^{+}(n+1)$ is
the one of interest since this is the solution which guarantees
$\|k\|_{\mathrm{op}}<1$:

$z_{j}^{+}(n+1)=\frac{\Theta[n](e_{j},e_{j})}{H[n](e_{j},\overline{e_{j}})+\sqrt{H[n](e_{j},\overline{e_{j}})^{2}-|\Theta[n](e_{j},e_{j})|^{2}}}.$
(III.9)

This allows us to write the energy
$\mathcal{E}\big{[}k_{n+1};H[n],\Theta[n]\big{]}$ as

$\begin{split}\mathcal{E}\big{[}k_{n+1};H[n],\Theta[n]\big{]}&=\mathcal{E}\Big{(}\\{z_{j}^{+}\\},\\{e_{j}\\}\Big{)}\\\
&=-\frac{1}{2}\sum_{j=1}^{\infty}{\Big{\\{}H[n](\overline{e_{j}},e_{j})-\sqrt{H^{2}[n](\overline{e_{j}},e_{j})-|\Theta[n](\overline{e_{j}},\overline{e_{j}})|^{2}}\Big{\\}}}\\\
&:=-\frac{1}{2}\sum_{j}{\mathcal{F}(e_{j})}.\end{split}$ (III.10)

The existence of $k_{n+1}$ proceeds by maximizing
$\sum_{j}{\mathcal{F}(e_{j})}$ over all orthogonal frames $\\{e_{j}\\}$. The
expression for $\mathcal{F}$ given above implies that $H[n](\overline{e},e)$
is bounded for any maximizing sequence in
$(\mathrm{span}\\{e_{i}\\}_{i=1}^{j})_{\perp}$, which allows us to extract a
strongly convergent subsequence in $L^{2}(\mathbb{R}^{3})$ to produce
$e_{j}(x)$. Having constructed the orthonormal basis and coefficients, the
proof of existence for $k_{n+1}$ is then complete.

######  Remark III.3.

(Solutions to the nonlinear Schrödinger-type equation for $\phi_{n+1}$) The
existence of constrained minimizer $\phi_{n+1}$ of

$\min_{\|\phi\|=1}\Big{\\{}E_{\mathrm{H}}[\phi]+\frac{1}{N}\mathcal{E}\big{[}k_{n+1};\,H[\phi,k_{n+1}],\,\Theta[\phi,k_{n+1}]\big{]},\quad\phi\in\mathfrak{h}^{1}_{V}\Big{\\}},$
(III.11)

will follow the usual proof for minimizing the Hartree functional
$E_{\mathrm{H}}[\phi]$ (see, for example, LiebSeiringer-book ), noting that
the additional terms introduced by the functional $\mathcal{E}$ are bounded
when $k_{n+1}$ is fixed with $\|k_{n+1}\|_{\mathrm{op}}<1$.

###  III.1 Preliminary lemmas

Several preliminary results are required to carry out the proof of Theorem
III.1 via iteration, and are presented here. Lemma III.4 is a general result
on the stability of spectral gaps for relatively-bounded perturbations of
self-adjoint operators. It is an extension of the work of Kato
kato2012perturbation ; a concise proof can be found in the work of Cuenin and
Tretter Cuenin . Lemma III.5 and its Corollary prove that the kernel

$\rho^{\mathrm{pair}}[n](x):=\Big{\\{}\frac{k_{n}\circ\overline{k_{n}}}{\delta-
k_{n}\circ\overline{k_{n}}}\Big{\\}}(x,x),$

defines $(\upsilon_{N}\ast\rho^{\mathrm{pair}}[n])(x)$ as a bounded operator
on $L^{2}(\mathbb{R}^{3})$ for all $n$. The discussion which follows Lemma
III.5 provides the existence of the sequence of functions
$\\{\phi_{n},\,k_{n}\\}_{n=0}^{\infty}$ which solve (III.5), (III.4) in the
iteration scheme. To show that the induction can be carried out, it must be
shown that the sequence $\Theta[n]$ remains uniformly bounded in Hilbert-
Schmidt and operator norm, and that the gap condition for the operators
$H[n],\,\Theta[n]$ (i.e., inequality (III.7)) holds uniformly for all $n$.
These estimates allow us to take strong limits of the operators
$H[n],\,\Theta[n]$ as $n\to\infty$ in (III.5), (III.4) to construct solutions
to the coupled system of Theorem III.1.

######  Lemma III.4.

(Cuenin and Tretter, 2016 Cuenin ) Let $T$ be a self-adjoint operator in a
Hilbert space $\mathfrak{h}$ and let $A$ be bounded, i.e., $\|Ax\|\leq a\|x\|$
for all $x\in\mathfrak{h}$ with $a\geq 0$. Then if $T$ has a spectral gap
$(\alpha_{T},\beta_{T})\subset\mathbb{R}$, i.e.,
$\sigma(T)\cap(\alpha_{T},\beta_{T})=\emptyset$ with
$\alpha_{T},\beta_{T}\in\sigma(T),$ and if

$2a<\beta_{T}-\alpha_{T},$

then $T+A$ has a stable spectral free strip
$(\alpha_{T+A},\beta_{T+A})+i\mathbb{R}\subset\mathbb{C}$, i.e.,

$\sigma(T+sA)\cap\\{z\in\mathbb{C}:\alpha_{T+A}<\mathrm{Re}(z)<\beta_{T+A}\\}=\emptyset,\quad
s\in[0,1],$

with

$\alpha_{T+A}:=\alpha_{T}+a,\quad\beta_{T+A}:=\beta_{T}-a.$

######  Lemma III.5.

Suppose that $k_{n}\in\mathfrak{B}_{2}(\mathfrak{h}^{1}_{V})$ solves (III.4)
with $\|k_{n}\|_{\mathrm{op}}<1$, where the operators $H[n-1](x,y)$ and
$\Theta[n-1](x,y)$ are given by formula (III.20). Also, suppose that the
operator $\Theta[n-1](x,y)$ is Hilbert-Schmidt. Then the density

$\rho^{\mathrm{pair}}[n](x):=\Big{(}\frac{k_{n}\circ\overline{k_{n}}}{\delta-
k_{n}\circ\overline{k_{n}}}\Big{)}(x,x),$ (III.12)

satisfies $\rho^{\mathrm{pair}}[n](x)\in L^{3}(\mathbb{R}^{3})\cap
L^{1}(\mathbb{R}^{3})$. In particular, the following estimate holds in
Hilbert-Schmidt norm:

$\|\rho^{\mathrm{pair}}[n]\|_{L^{3}}\leq C\Big{\|}\frac{k_{n}}{\delta-
k_{n}\circ\overline{k_{n}}}\Big{\|}_{\mathrm{HS}}\cdot\Big{(}\Big{\|}\Theta[n-1]\Big{\|}_{\mathrm{HS}}+\mu[n-1]\Big{)},$
(III.13)

where $C$ is some constant independent of $n$.

###### Proof.

Define the operator kernel $\psi(x,y)$ as in Section II.2 by

$\psi(x,y):=\Big{(}\frac{k_{n}}{(\delta-
k_{n}\circ\overline{k_{n}})^{1/2}}\Big{)}(x,y),$ (III.14)

so that $\rho^{\mathrm{pair}}[n](x)=(\psi\circ\overline{\psi})(x,x)$. Then
using Hölder’s inequality,

$\begin{split}|\nabla\rho^{\mathrm{pair}}[n](x)|^{(3/2)}&=\Big{|}\nabla_{x}\int{{\rm
d}z\\{\psi(x,z)\overline{\psi(z,x)}\\}\Big{|}}^{3/2}\\\ &\leq\Big{(}\int{{\rm
d}z|\nabla_{x}\psi(x,z)|^{2}}\Big{)}^{3/4}\Big{(}\int{{\rm
d}z|\psi(x,z)|^{2}}\Big{)}^{3/4}\\\ &=\Big{(}\int{{\rm
d}z|\nabla_{x}\psi(x,z)|^{2}}\Big{)}^{3/4}\Big{(}\rho^{\mathrm{pair}}[n](x)\Big{)}^{3/4}.\end{split}$
(III.15)

Using Holder’s inequality again, ($p=4/3,\,q=4$) gives:

$\begin{split}\int{{\rm
d}x|\nabla\rho^{\mathrm{pair}}[n](x)|^{3/2}}&\leq\int{{\rm d}x\Big{(}\int{{\rm
d}z|\nabla\psi(x,z)|^{2}}\Big{)}^{3/4}\Big{(}\rho^{\mathrm{pair}}[n](x)\Big{)}^{3/4}}\\\
&\leq\Big{(}\iint{{\rm d}x{\rm
d}z\big{|}\nabla_{x}\psi(x,z)\big{|}^{2}}\Big{)}^{3/4}\Big{(}\int{{\rm
d}x\big{(}\rho^{\mathrm{pair}}[n](x)\big{)}^{3}}\Big{)}^{1/4},\end{split}$
(III.16)

The Sobolev embedding $W^{1,\,3/2}(\mathbb{R}^{3})\subseteq
L^{3}(\mathbb{R}^{3})$ gives (up to a constant)

$\Big{(}\int{{\rm
d}x\rho^{\mathrm{pair}}[n](x)^{3}}\Big{)}^{1/3}\leq\|\nabla\rho^{\mathrm{pair}}[n]\|_{L^{3/2}}.$

Combining this with the previous equation, we have

$\Big{(}\int{{\rm
d}x\rho^{\mathrm{pair}}[n](x)^{3}}\Big{)}^{1/3}\leq\Big{(}\iint{{\rm d}x{\rm
d}z\big{|}\nabla_{x}\psi(x,z)\big{|}^{2}}\Big{)}.$

Now, we know that

$\begin{split}&\mathcal{E}(k_{n})=\\\ &\mathrm{tr}\Big{\\{}\big{(}\delta-
k_{n}\circ\overline{k_{n}})^{-1}\circ\Big{(}k_{n}\circ
H[n-1]\circ\overline{k_{n}}+\frac{1}{2}k_{n}\circ\overline{\Theta[n-1]}+\frac{1}{2}\Theta[n-1]\circ\overline{k_{n}}\Big{)}\Big{\\}}\leq
0.\end{split}$ (III.17)

This completes the proof. ∎

######  Corollary 2.

Under the assumptions of Lemma III.5, $\rho^{\mathrm{pair}}[n](x)\in L^{2}$
with

$\|\rho^{\mathrm{pair}}[n]\|_{L^{2}}\leq\Big{\|}\frac{k_{n}}{\delta-
k_{n}\circ\overline{k_{n}}}\Big{\|}_{\mathrm{HS}}\cdot\Big{(}\Big{\|}\Theta[n-1]\Big{\|}_{\mathrm{HS}}+\mu[n-1]\Big{)}^{3/4}.$
(III.18)

###### Proof.

Interpolation gives

$\|\rho^{\mathrm{pair}}[n]\|_{L^{2}}\leq\|\rho^{\mathrm{pair}}[n]\|^{3/4}_{L^{3}}\|\rho^{\mathrm{pair}}[n]\|_{L^{1}}^{1/4},$

and

$\|\rho^{\mathrm{pair}}[n]\|_{L^{1}}=\mathrm{tr}\Big{(}\frac{k_{n}\circ\overline{k_{n}}}{\delta-
k_{n}\circ\overline{k_{n}}}\Big{)}\leq\|k_{n}\|_{\mathrm{op}}\Big{\|}\frac{k_{n}}{\delta-
k_{n}\circ\overline{k_{n}}}\Big{\|}_{\mathrm{HS}},$

which proves the statement. ∎

###  III.2 Proof of Theorem III.1

###### Proof.

We first demonstrate the existence of solutions
$\\{\phi_{n},\,k_{n};\,\mu_{n}\\}_{n=0}^{\infty}$ to the iteration scheme. We
also verify the gap condition for the operators $H[n],\,\Theta[n]$ at every
step, showing that there exists $c>0$ independent of $n$ such that

$H[n](e,\overline{e})-|\Theta[n](\overline{e},\overline{e})|\geq c\|e\|,\quad
e\in(\phi_{n})_{\perp},\quad n=0,1,2,\dots.$ (III.19)

We show this inductively, beginning with $\phi_{0}(x)$, which is a solution to
the Hartree equation:

$\big{(}-\Delta_{x}+V_{\mathrm{trap}}(x)\big{)}\phi_{0}(x)+(\upsilon_{N}\ast|\phi_{0}|^{2})(x)\phi_{0}(x)=\mu_{0}\phi_{0}(x).$

The constant $\mu_{0}$ is given by

$\mu_{0}=\iint{{\rm d}x{\rm
d}y\\{\overline{\phi_{0}(x)}\epsilon_{0}(x,y)\phi_{0}(y)+|\phi_{0}(x)|^{2}\upsilon_{N}(x-y)|\phi_{0}(y)|^{2}\\}},$

and the operators $H[0],\,\Theta[0]$ are defined by

$\begin{split}H[0]&:=\big{(}-\Delta_{x}+V_{\mathrm{trap}}(x)\big{)}\delta(x-y)+\overline{\phi_{0}(y)}\upsilon_{N}(x-y)\phi_{0}(x)\\\
&\quad\quad+(\upsilon_{N}*|\phi_{0}|^{2})(y)\delta(x-y)-\mu_{0},\\\
\Theta[0]&:=\upsilon_{N}(x,y)\Big{\\{}\phi_{0}(x)\phi_{0}(y)\Big{\\}}.\end{split}$
(III.20)

In particular, the operator $\Theta[0]$ is Hilbert-Schmidt, with norm

$\|\Theta[0]\|_{\mathrm{HS}}\leq
gN^{3\beta}\|\upsilon\|_{L^{\infty}}\|\phi_{0}\|_{L^{2}}^{2}.$

It is not required that $\Theta[0]$ must be small in this scheme, only that it
is bounded (see equation (III.37)). Using the assumption
$\widehat{\upsilon}\geq 0$ for the Fourier transform of $\upsilon$,

$\begin{split}H[0](e,\overline{e})-|\Theta[0](\overline{e},\overline{e})|&\geq\int{{\rm
d}x\\{e(x)\Big{(}-\Delta_{x}+V_{\mathrm{trap}}(x)\Big{)}\overline{e(x)}\\}}\\\
&+\iint{{\rm d}x{\rm
d}y\\{\upsilon_{N}(x-y)|\phi_{0}(x)|^{2}|e(y)|^{2}\\}}-\mu_{0}\|e\|^{2}.\end{split}$
(III.21)

The gap condition for $H[0],\,\Theta[0]$ follows from this inequality, because
$\phi_{0}(x)$, as the minimizer of $E_{\mathrm{H}}(\phi)$, is also the ground
state of the linear Hartree operator, $H_{\mathrm{Hartree}}$, defined by

$H_{\mathrm{Hartree}}[\phi_{0}]:=\Big{(}-\Delta_{x}+V_{\mathrm{trap}}(x)\Big{)}\delta(x-y)+(\upsilon_{N}\ast|\phi_{0}|^{2})(x)\delta(x-y).$

$H_{\mathrm{Hartree}}$ has a discrete spectrum, with a gap between its lowest
and first eigenvalues, since it is a compact perturbation of
$\Big{(}-\Delta_{x}+V_{\mathrm{trap}}(x)\Big{)}$, which follows by Young’s
convolution inequality, and the compact Sobolev embedding
$L^{4}(\mathbb{R}^{3})\subset\mathfrak{h}^{1}_{V}$,

$\|(\upsilon_{N}\ast|\phi_{0}|^{2})\|_{L^{2}}\leq\|\upsilon_{N}\|_{L^{1}}\|~{}|\phi_{0}|^{2}~{}\|_{L^{2}}=\|\upsilon_{N}\|_{L^{1}}\|\phi_{0}\|_{L^{4}}\leq\|\upsilon_{N}\|_{L^{1}}\|\phi_{0}\|_{\mathfrak{h}^{1}_{V}},$
(III.22)

i.e.,

$\|(\upsilon_{N}\ast|\phi_{0}|^{2})\|_{L^{2}}\leq\|\upsilon_{N}\|_{L^{1}}\cdot
E_{\mathrm{H}}(\phi_{0})\leq g\cdot E_{\mathrm{H}}(\phi_{0}).$ (III.23)

The perturbation to $(\upsilon_{N}\ast|\phi_{0}|^{2})(x)$ can be made
arbitrarily small in $L^{2}$ by assuming that $g$ is small enough at the
beginning of the induction.

For the induction step, assume that there exists $c>0$ and $C>0$ such that for
$s=0,1,2,\dots,n$,

$\begin{cases}\big{\|}\Theta[s]\big{\|}_{\mathrm{HS}}\leq CN^{3\beta},\\\
\mu[s]\leq C,\\\
H[s](e,\overline{e})-|\Theta[s](\overline{e},\overline{e})|\geq
c\|e\|^{2},\quad e\in(\phi_{s})_{\perp},\end{cases}$ (III.24)

We prove the inequalities of (III.24) for $s=n+1$.

The equation for $k_{n+1}$ on the space $\mathfrak{h}(\phi_{\perp})$ is

$H[n]\circ k_{n+1}+k_{n+1}\circ
H[n]^{T}+\Theta[n]+k_{n+1}\circ\overline{\Theta[n]}\circ k_{n+1}=0,$ (III.25)

and the equation for $\phi_{n+1}$ is given by (II.1) with $\phi=\phi_{n+1}$,
$k=k_{n+1}$ and $\mu=\mu_{n+1}$.

To show the first inequality of (III.24) for $s=n+1$, the induction assumption
allows us to construct $k_{n+1}$ as a minimizer of
$\mathcal{E}\big{[}k;H[n],\Theta[n]\big{]}$, with the formula:

$\Big{(}\frac{k_{n+1}}{\delta-
k_{n+1}\circ\overline{k_{n+1}}}\Big{)}(x,y)=\sum_{j}{\Big{(}\frac{z_{j}}{1-|z_{j}|^{2}}\Big{)}e_{j}(x)e_{j}(y)}.$
(III.26)

Here $\\{e_{j}(n+1)\\}_{j=1}^{\infty}$ is an orthonormal basis of
$(\phi_{n})_{\perp}$, which is determined as part of the construction of
$k_{n+1}$. The formula for the coefficient $z_{j}(n+1)$ gives

$\Big{|}\frac{z_{j}(n+1)}{1-|z_{j}(n+1)|^{2}}\Big{|}^{2}=\frac{1}{4}\Big{(}\frac{|\Theta[n](e_{j},e_{j})|^{2}}{H[n]^{2}(e_{j},\overline{e_{j}})-|\Theta[n](e_{j},e_{j})|^{2}}\Big{)}\leq\frac{1}{4c^{2}}|\Theta[n](e_{j},e_{j})|^{2},$
(III.27)

where the inequality follows since we assume the gap condition for the
operators $H[n],\,\Theta[n]$. Thus,

$\|\Big{(}\frac{k_{n+1}}{\delta-
k_{n+1}\circ\overline{k_{n+1}}}\Big{)}\|_{\mathrm{HS}}\leq\frac{1}{2c}\|\Theta[n]\|_{\mathrm{HS}}\leq\frac{C}{2c}N^{3\beta}.$

From this (recalling that $\upsilon_{N}(x)=gN^{3\beta}\upsilon(N^{\beta}x)$
for $\beta\leq 1/6$), we have

$\|\frac{1}{N}\upsilon_{N}\Big{(}\frac{k_{n+1}}{\delta-
k_{n+1}\circ\overline{k_{n+1}}}\Big{)}\|_{\mathrm{HS}}\leq\frac{gC^{2}}{2c}\|\upsilon\|_{L^{\infty}}\Big{(}\frac{N^{6\beta}}{N}\Big{)},$
(III.28)

as well as

$\|\frac{1}{N}(\upsilon_{N}\Big{(}\frac{k_{n+1}\circ\overline{k_{n+1}}}{\delta-
k_{n+1}\circ\overline{k_{n+1}}}\Big{)}\|_{\mathrm{HS}}\leq\frac{gC^{2}}{2c}\|\upsilon\|_{L^{\infty}}\Big{(}\frac{N^{6\beta}}{N}\Big{)},$
(III.29)

and by Corollary (2) with Young’s convolution inequality:

$\begin{split}\|\frac{1}{N}(\upsilon_{N}\ast\rho^{\mathrm{pair}}[n+1])\|_{L^{2}}&\leq\frac{gC}{2c}\frac{N^{3\beta}}{N}\cdot\Big{(}\Big{\|}\Theta[n-1]\Big{\|}_{\mathrm{HS}}+\mu[n-1]\Big{)}^{3/4}\\\
&\leq\frac{gC^{2}}{c}\Big{(}\frac{N^{6\beta}}{N}\Big{)}.\end{split}$ (III.30)

It is now clear that for $\beta\leq 1/6$ the constant $g>0$ can be chosen at
the start of the induction so that

$\|\Theta[n+1]\|_{\mathrm{HS}}^{2}<CN^{3\beta}.$

This is because, e.g.,

$\begin{split}\|\Theta[n+1]\|_{\mathrm{HS}}&\leq
gN^{3\beta}\|\upsilon\|_{L^{\infty}}\|\phi_{n+1}\|_{L^{2}}^{2}+\frac{1}{N}\|(\Theta^{\mathrm{pair}}[n+1])\|_{\mathrm{HS}}\\\
&\leq
g\|\upsilon\|_{L^{\infty}}\Big{\\{}N^{3\beta}+\frac{C^{2}}{2c}\Big{\\}},\end{split}$
(III.31)

where

$\Theta^{\mathrm{pair}}[n+1](x,y):=\frac{1}{N}\upsilon_{N}(x-y)\Big{(}\frac{k_{n+1}}{\delta-
k_{n+1}\circ\overline{k_{n+1}}}\Big{)}(x,y).$ (III.32)

We remind the reader that the inductive argument only requires $\Theta[n+1]$
to be bounded (not small), and that the particle number $N$ is large but fixed
throughout the entire work.

We use the above conclusions to show that $\mu[n+1]<C$. Indeed, $\phi_{n+1}$
minimizes the functional $\tilde{\mathcal{E}}[n+1](\phi)$, which consists of
only the terms of $\mathcal{E}_{\mathrm{tot}}[\phi,k_{n+1}]$ which have a
dependence on $\phi$, and satisfies the relation (for $n>0$),

$\begin{split}\tilde{\mathcal{E}}[n+1](\phi_{n+1})&=\mu_{n+1}-\frac{1}{2}\iint{{\rm
d}x{\rm d}y|\phi_{n+1}(x)|^{2}\upsilon_{N}(x-y)|\phi_{n+1}(y)|^{2}}\\\
&-\frac{1}{2N}(\upsilon_{N}\Big{(}\frac{k_{n+1}}{\delta-
k_{n+1}\circ\overline{k_{n+1}}}\Big{)})(\overline{\phi_{n+1}},\overline{\phi_{n+1}}).\end{split}$
(III.33)

Since $\phi_{n+1}$ is the minimizer of $\tilde{\mathcal{E}}[n+1]$, we have
(recall the definition of $\rho^{\mathrm{pair}}[n+1](x)$ in (III.12))

$\begin{split}\tilde{\mathcal{E}}[n+1](\phi_{n+1})&\leq\tilde{\mathcal{E}}[n+1](\phi_{0})=\mu_{0}-\frac{1}{2}\iint{{\rm
d}x{\rm d}y|\phi_{0}(x)|^{2}\upsilon_{N}(x-y)|\phi_{0}(y)|^{2}}\\\
&+\frac{1}{N}\iint{{\rm d}x{\rm
d}y|\phi_{0}(x)|^{2}\upsilon_{N}(x-y)\rho^{\mathrm{pair}}[n+1](y)}\\\
&+\frac{1}{N}(\upsilon_{N}\Big{(}\frac{k_{n+1}\circ\overline{k_{n+1}}}{\delta-
k_{n+1}\circ\overline{k_{n+1}}}\Big{)})(\phi_{0},\overline{\phi_{0}})\\\
&+\frac{1}{2}\frac{1}{N}(\upsilon_{N}\Big{(}\frac{k_{n+1}}{\delta-
k_{n+1}\circ\overline{k_{n+1}}}\Big{)})(\overline{\phi_{0}},\overline{\phi_{0}}).\end{split}$
(III.34)

Combining (III.33) and (III.34) with inequalities (III.28)–(III.30), we
conclude that

$\begin{split}\mu_{n+1}&\leq\mu_{0}+\frac{g}{2}\big{(}E_{\mathrm{H}}(\phi_{0})\big{)}^{2}+\frac{g}{2}\big{(}\tilde{\mathcal{E}}[n+1](\phi_{n+1})\big{)}^{2}\\\
&+\frac{2gC}{c}\|\upsilon\|_{L^{\infty}}+\frac{gC^{2}}{Nc}\mathcal{E}_{\mathrm{Hartree}}(\phi_{0}),\end{split}$
(III.35)

where we have made use of Sobolev’s and Young’s inequalities to write

$\begin{split}\iint{{\rm d}x{\rm
d}y|\phi_{0}(x)|^{2}\upsilon_{N}(x-y)|\phi_{0}(y)|^{2}}&\leq
E_{\mathrm{H}}(\phi_{0})\cdot\|(\upsilon_{N}\ast|\phi_{0}|^{2})\|_{L^{2}}\\\
\iint{{\rm d}x{\rm
d}y|\phi_{n+1}(x)|^{2}\upsilon_{N}(x-y)|\phi_{n+1}(y)|^{2}}&\leq\tilde{\mathcal{E}}[n+1](\phi_{n+1})\cdot\|(\upsilon_{N}\ast|\phi_{n+1}|^{2})\|_{L^{2}},\end{split}$
(III.36)

along with the following $L^{2}$ bound which is analogous to (III.22)

$\|(\upsilon_{N}\ast|\phi_{n+1}|^{2})\|_{L^{2}}\leq\|\upsilon_{N}\|_{L^{1}}\cdot\|\phi_{n+1}\|_{\mathfrak{h}^{1}_{V}}\leq
g\cdot\tilde{\mathcal{E}}[n+1](\phi_{n+1}),$

and finally,

$\frac{1}{N}\iint{{\rm d}x{\rm
d}y|\phi_{0}(x)|^{2}\upsilon_{N}(x-y)\rho^{\mathrm{pair}}[n+1])(y)}\leq\frac{g}{N}\cdot
E_{\mathrm{H}}(\phi_{0})\|\rho^{\mathrm{pair}}[n+1]\|_{L^{2}}.$

This shows that $\mu_{n+1}$ modifies $\mu_{0}$ by a correction which is
controlled by the constant $g$, which allows us to conclude the uniform
boundedness of the constants $\mu[s]$ for $s\leq n+1$.

Finally, we show the gap condition (III.19) for $s=n+1$. Observe that for any
$e(x)\in(\phi_{n})_{\perp}$,

$\begin{split}H[n+1](e,\overline{e})-|\Theta[n+1](\overline{e},\overline{e})|&\geq
h[n+1](e,\overline{e})-\mu_{n+1}\\\
&-\frac{1}{N}|(\upsilon_{N}\Big{(}\frac{k_{n+1}}{\delta-
k_{n+1}\circ\overline{k_{n+1}}}\Big{)})(\overline{e},\overline{e})|,\end{split}$
(III.37)

where the operator $h[n+1]$ is defined, using $\rho^{\mathrm{pair}}[n+1]$ of
(III.12), by

$\begin{split}h[n+1](x,y)&:=\big{(}-\Delta_{x}+V_{\mathrm{trap}}(x)\big{)}\delta(x-y)+(\upsilon_{N}\ast|\phi_{n+1}|^{2})(x)\delta(x-y)\\\
&+\frac{1}{N}(\upsilon_{N}\Big{(}\frac{k_{n+1}\circ\overline{k_{n+1}}}{\delta-
k_{n+1}\circ\overline{k_{n+1}}}\Big{)})(x,y)+\frac{1}{N}(\upsilon_{N}\ast\rho^{\mathrm{pair}}[n+1])(x)\delta(x-y).\end{split}$
(III.38)

Also observe that $\phi_{n+1}$ satisfies the self-adjoint operator matrix
equation

$\begin{pmatrix}h[n+1]&\Theta^{\mathrm{pair}}[n+1]\\\
\overline{\Theta^{\mathrm{pair}}[n+1]}&h[n+1]\end{pmatrix}\circ\begin{pmatrix}\phi_{n+1}\\\
\overline{\phi_{n+1}}\end{pmatrix}=\mu_{n+1}\begin{pmatrix}\phi_{n+1}\\\
\overline{\phi_{n+1}}\end{pmatrix},$ (III.39)

for

$\Theta^{\mathrm{pair}}[n+1](x,y):=\frac{1}{N}\upsilon_{N}(x-y)\Big{(}\frac{k_{n+1}}{\delta-
k_{n+1}\circ\overline{k_{n+1}}}\Big{)}(x,y).$ (III.40)

The rest of the argument follows by showing that: (i) $\phi_{n+1}$ which
minimizes $\tilde{\mathcal{E}}[n+1](\phi)$ is the ground state of the operator
in (III.39); and (ii) the higher eigenvalues of the problem (III.39), which
solve

$\begin{pmatrix}h[n+1]&\Theta^{\mathrm{pair}}[n+1]\\\
\overline{\Theta^{\mathrm{pair}}[n+1]}&h[n+1]\end{pmatrix}\circ\begin{pmatrix}e\\\
\overline{e}\end{pmatrix}=\lambda(e)\begin{pmatrix}e\\\
\overline{e}\end{pmatrix},\quad\lambda(e)\geq\mu_{n+1},\quad
e(x)\in(\phi_{n+1})_{\perp},$ (III.41)

are distinct from $\mu_{n+1}$, which has no degeneracy, i.e.,

$\lambda(e)-\mu_{n+1}\geq\delta,\quad
e\in(\phi_{n+1})_{\perp},\quad\mathrm{for}\quad\delta>0.$

The $\lambda(e)$ are given by the formula

$\lambda(e)=h[n+1](e,\overline{e})+\Theta^{\mathrm{pair}}[n+1](\overline{e},\overline{e}).$
(III.42)

Both of these arguments follow from the fact that the matrix in equations
(III.39), (III.41) is a perturbation of the diagonal matrix
$\mathrm{diag}(\epsilon_{0},\,\epsilon_{0})$, where $\epsilon_{0}$ is the
single-particle Schrödinger operator,
$\epsilon_{0}:=-\Delta+V_{\mathrm{trap}}$. The perturbation can be written as
the matrix $M^{\prime}[n+1]=M_{1}[n+1]+M_{2}[n+1]$ with

$\begin{split}M_{1}[n+1]&:=\mathrm{diag}\Big{(}\\{(\upsilon_{N}\ast|\phi_{n+1}|^{2})(x)+\frac{1}{N}(\upsilon_{N}\ast\rho^{\mathrm{pair}})(x)\\}\delta(x-y)\\\
&+\frac{1}{N}(\upsilon_{N}\Big{(}\frac{k_{n+1}\circ\overline{k_{n+1}}}{\delta-
k_{n+1}\circ\overline{k_{n+1}}}\Big{)})(x,y)\Big{)},\end{split}$ (III.43)

and

$M_{2}[n+1]:=\begin{pmatrix}0&\Theta^{\mathrm{pair}}[n+1]\\\
\overline{\Theta^{\mathrm{pair}}[n+1]}&0\end{pmatrix}.$

The perturbation $M^{\prime}[n+1]$ is bounded in operator norm by

$\|M^{\prime}[n+1]\|_{\mathrm{op}}\leq
g\cdot\tilde{\mathcal{E}}[n+1](\phi_{n+1})+g\|\upsilon\|_{L^{\infty}}C,\quad\mathrm{for}\quad
C<\infty.$ (III.44)

The first term of this operator bound comes from Young’s inequality (i.e.,
equation (III.35)). The second term in (III.44) meanwhile comes from (III.28).
The important feature is that the operator bound for $M[n+1]$ can be made as
small as desired by choosing $g$ small enough.

We invoke the perturbation lemma, Lemma III.4, in order to claim that the
spectrum of the problem (III.41) is a perturbation of the spectrum
$\sigma(\epsilon_{0})$. In particular, for $g$ small enough we can guarantee
that $\mu_{n+1}$ is the ground state energy of this system, since it is the
only element which is in a neighborhood of size $gC$ of the ground state of
$\epsilon_{0}$ (where $C$ is some constant). It follows that

$\lambda(e)-\mu_{n+1}\geq\delta-2\|M^{\prime}[n+1]\|_{\mathrm{op}},$

for all $e\in(\phi_{n+1})_{\perp}$, where $\delta$ is the size of the gap for
the operator $\epsilon_{0}.$ Since $\|M^{\prime}[n+1]\|_{\mathrm{op}}$ can be
made arbitrarily small by making $g$ small enough, we can choose the constants
$g,\,c$ to be small enough so that $\lambda(e)-\mu_{n+1}>c$. This completes
the proof of existence for the sequence
$\\{\phi_{n},\,k_{n};\,\mu_{n}\\}_{n=0}^{\infty}$, and the uniform gap
condition for the operators $H[n],\,\Theta[n]$.

We finally show that the sequence
$\\{\phi_{n},\,k_{n};\,\mu_{n}\\}_{n=0}^{\infty}$ converges to the appropriate
solution as described in the statement of Theorem III.1. Since the gap
condition holds uniformly for every step of the iteration, the sequence
$\\{\phi_{n},k_{n}\\}$ determined by the minimization problems (III.5) and
(III.4) exists, with $\|\phi_{n}\|=1$ and $\|k_{n}\|_{\mathrm{HS}}<C$,
$\|(\delta-k_{n}\circ\overline{k_{n}})^{-1}\|_{\mathrm{op}}<C$ for all $n$ and
some $C$. The sequence $\\{\mu_{n}\\}$ is also uniformly bounded. Since we
also have

$\begin{split}\mathrm{tr}\Big{\\{}\overline{k_{n}}\circ\epsilon_{0}\circ
k_{n}\Big{\\}}&\leq\Big{|}\mathrm{tr}\Big{\\{}(\delta-\overline{k_{n}}\circ
k_{n})^{-1}\circ(\overline{k_{n}}\circ\epsilon_{0}\circ
k_{n})\Big{\\}}\Big{|}\\\
&\leq\frac{1}{2}\Big{|}\mathrm{tr}\Big{\\{}(\delta-\overline{k_{n}}\circ
k_{n})^{-1}\circ\\{\overline{k_{n}}\circ\Theta[n-1]+\overline{\Theta[n-1]}\circ
k_{n}\\}\Big{\\}}\Big{|}\\\
&+\mu_{n}\mathrm{tr}\Big{\\{}(\delta-\overline{k_{n}}\circ
k_{n})^{-1}\circ(\overline{k_{n}}\circ k_{n})\Big{\\}}\Big{|},\end{split}$
(III.45)

the sequence $\\{k_{n}\\}$ is uniformly bounded in
$\mathfrak{B}_{2}(\mathfrak{h}^{1}_{V},\mathfrak{h}^{1}_{V})$. Because
$\\{\phi_{n}\\}$ and $\\{k_{n}\\}$ are uniformly bounded in
$\mathfrak{h}^{1}_{V}$, $\mathfrak{B}_{2}(\mathfrak{h}^{1}_{V})$, we can find
some subsequence (again denoted $\\{\phi_{n}\\}$ and $\\{k_{n}\\}$) which
converges weakly in $\mathfrak{h}^{1}_{V}$,
$(\mathfrak{B}_{2}(\mathfrak{h}^{1}_{V})$ respectively). Since
$\mathfrak{h}^{1}_{V}$ is compactly embedded in $L^{2}_{V}\cap L^{4}_{V}$, we
conclude that $\phi_{n}$ converges strongly in $L^{2}_{V}\cap L^{4}_{V}$ to
some $\phi$, and $k_{n}$ converges strongly in $\mathfrak{B}_{2}(L^{2}_{V})$
to some $k$.

The following strong limits in $\mathfrak{B}_{2}(\mathfrak{h})$ therefore
hold:

$\begin{split}\upsilon_{N}(x-y)\phi_{n}(x)\phi_{n}(y)&\to\upsilon_{N}(x-y)\phi(x)\phi(y),\\\
\upsilon_{N}(x-y)\Big{\\{}\frac{k_{n}}{\delta-\overline{k_{n}}\circ
k_{n}}\Big{\\}}&\to\upsilon_{N}(x-y)\Big{\\{}\frac{k}{\delta-\overline{k}\circ
k}\Big{\\}},\\\
\upsilon_{N}(x-y)\Big{\\{}\frac{k_{n}\circ\overline{k_{n}}}{\delta-
k_{n}\circ\overline{k_{n}}}\Big{\\}}&\to\upsilon_{N}(x-y)\Big{\\{}\frac{k\circ\overline{k}}{\delta-k\circ\overline{k}}\Big{\\}}.\end{split}$
(III.46)

We also have the strong limits in $L^{2}(\mathbb{R}^{3})$:

$(\upsilon_{N}\ast|\phi_{n}|^{2})(y)\delta(x-y)\to(\upsilon_{N}\ast|\phi|^{2})(y)\delta(x-y),$
(III.47)

and

$\begin{split}&\int{{\rm
d}x^{\prime}\upsilon_{N}(x-x^{\prime})\Big{\\{}\frac{k_{n}\circ\overline{k_{n}}}{\delta-
k_{n}\circ\overline{k_{n}}}\Big{\\}}(x^{\prime},x^{\prime})}\delta(x-y)\\\
&\to\int{{\rm
d}x^{\prime}\upsilon_{N}(x-x^{\prime})\Big{\\{}\frac{k\circ\overline{k}}{\delta-k\circ\overline{k}}\Big{\\}}(x^{\prime},x^{\prime})}\delta(x-y).\end{split}$
(III.48)

The first limit in (III.47) follows by noticing (via Young’s convolution
inequality),

$\|(\upsilon_{N}\ast|\phi_{n}|^{2})\|_{L^{2}}\leq\|\upsilon_{N}\|_{L^{1}}\cdot\||\phi_{n}|^{2}\|_{L^{2}}=\|\upsilon_{N}\|_{L^{1}}\cdot\|\phi_{n}\|_{L^{4}},$

and so $(\upsilon_{N}\ast|\phi_{n}|^{2})(x)\to(\upsilon_{N}\ast|\phi|^{2})(x)$
strongly in $L^{2}$, since $\phi_{n}\to\phi$ strongly in
$L^{4}(\mathbb{R}^{3})$. The second of these limits follows by the fact that
$\rho^{\mathrm{pair}}[n](x)\in L^{2}(\mathbb{R}^{3})$ for all $n$.

Taking $n\to\infty$ in the equations satisfied by $\phi_{n},k_{n}$ therefore
give equations (II.1), (II.2) for $\phi,k$. Finally, we must check that
$k\perp\phi$. Since $k_{n}(\phi_{n},x)=0$ by construction, we have

$\|k(\phi,x)\|_{L^{2}}\leq\|k\|_{\mathrm{HS}}\cdot\|\phi-\phi_{n}\|_{L^{2}}+\|k-k_{n}\|_{\mathrm{HS}}\cdot\|\phi_{n}\|_{L^{2}}.$

By picking $n$ sufficiently large, we can bound $\|k(\phi,x)\|$ by an
arbitrarily small quantity. This completes the proof. ∎

##  IV Conclusion

The result here represents an extension of the pair-excitation formalism to
quadratic models for the dilute Bose gas which include the coupling of the
condensate to excited states. The pair-excitation kernel introduced here is
related to the diagonalization of the quadratic model in a fundamental way;
the solution $k(x,y)$ can be used to construct a basis
$(u_{j}(x),\,p_{j}(x))_{j}$ which describes the unitary rotation of the
quadratic Hamiltonian.

The fact that the mathematical formalism of GrillakisMargetisSorokanich  can
be adapted in the presence of more complicated mean-field couplings perhaps
lends credit to understanding pair-excitation as a fundamental physical
mechanism. The model presented here is moreover amenable to further
generalization. In Griffin1996 , Griffin describes the contribution of thermal
mean-field effects through finite-temperature averages, i.e., replacing the
ground state expectations with thermal expectations: $\langle
b_{x}b_{y}\rangle_{T>0}$. Such expectations involve the Bose-Einstein
statistics for indistinguishable bosons. A rigorous treatment of this system
would necessitate, along with the basis $\\{u_{j},p_{j}\\}$, a procedure to
determine the thermal mean-field averages given by

$\langle\alpha_{j}^{\ast}\alpha_{k}\rangle=N_{0}(E_{j})\delta_{jk}\,,\quad
N_{0}(E_{j}):=\Big{(}\frac{1}{e^{\beta(T)E_{j}}-1}\Big{)},\quad\beta(T):=1/k_{B}T,\quad
T\geq 0.$ (IV.1)

The results presented in this work thus constitute the T=0 case.

## Acknowledgements

The author would like to thank Professor Manoussos Grillakis for his guidance,
useful discussions, and for his help on the proof of Theorem II.1. The author
would also like to acknowledge Professor Dionisios Margetis for inspiring the
study of this system.

## References

  * (1) S. Albeverio, A. K. Motovilov, and C. Tretter, _Bounds on the spectrum and reducing subspaces of a $J$-self-adjoint operator_, Indiana Univ. Math. J. 59 (2010), no. 5, 1737–1776.
  * (2) Sergio Albeverio and Alexander K. Motovilov, _Solvability of the operator Riccati equation in the Feshbach case_ , Math. Notes 105 (2019), no. 4, 485–502.
  * (3) Volker Bach, Sebastien Breteaux, Thomas Chen, Jürg Fröhlich, and Israel Michael Sigal, _The time-dependent hartree–fock–bogoliubov equations for bosons_ , Journal of Evolution Equations 22 (2022), no. 2, 46.
  * (4) Carl M. Bender, _Making sense of non-Hermitian Hamiltonians_ , Rep. Prog. Phys. 70 (2007), no. 6, 947–1018.
  * (5) N. Bogoliubov, _On the theory of superfluidity_ , J. Phys. 11 (1947), 23–32.
  * (6) C. Brennecke, P. T. Nam, M. Napiórkowski, and B. Schlein, _Fluctuations of N-particle quantum dynamics around the nonlinear Schrödinger equation_ , Ann. Inst. Henri Poincaré 36 (2019), no. 5, 1201–1235.
  * (7) Jean-Claude Cuenin and Christiane Tretter, _Non-symmetric perturbations of self-adjoint operators_ , J. Math. Anal. Appl. 441 (2016), 235–258.
  * (8) , _Non-symmetric perturbations of self-adjoint operators_ , Journal of Mathematical Analysis and Applications 441 (2016), no. 1, 235–258.
  * (9) F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, _Theory of Bose-Einstein condensation in trapped gases_ , Rev. Mod. Phys. 71 (1999), no. 3, 463–512.
  * (10) A. L. Fetter, _Nonuniform states of an imperfect bose gas_ , Annals Phys. 70 (1972), no. 1, 67–101.
  * (11) C. W. Gardiner, _Particle-number-conserving Bogoliubov method which demonstrates the validity of the time-dependent Gross-Pitaevskii equation for a highly condensed Bose gas_ , Phys. Rev. A 56 (1997), no. 2, 1414–1423.
  * (12) A. Griffin, _Conserving and gapless approximations for an inhomogeneous Bose gas at finite temperatures_ , Phys. Rev. B 53 (1996), no. 14, 9341–9347.
  * (13) A. Griffin, T. Nikuni, and E. Zaremba, _Bose-condensed gases at finite temperatures._ , Cambridge University Press, 2009.
  * (14) M. Grillakis and M. Machedon, _Beyond mean field: On the role of pair excitations in the evolution of condensates_ , J. Fixed Point Theory Appl. 14 (2013), 91–111.
  * (15) , _Pair excitations and the mean field approximation of interacting Bosons.I._ , Comm. Math. Phys. 324 (2013), 601–636.
  * (16) Manoussos Grillakis, Matei Machedon, and Dionisios Margetis, _Second-order corrections to mean field evolution of weakly interacting Bosons. I._ , Comm. Math. Phys. 294 (2010), 273–301.
  * (17) , _Evolution of the Boson gas at zero temperature: Mean-field limit and second-order correction_ , Quart. Appl. Math. 75 (2017), no. 1, 69–104.
  * (18) E. P. Gross, _Structure of a quantized vortex in boson systems_ , Nuovo Cim. 20 (1961), 454–477.
  * (19) K. Huang, _Statistical mechanics, 2nd ed_ , Wiley India Pvt. Limited, 2008.
  * (20) K. Huang and C. N. Yang, _Quantum-mechanical many-body problem with hard-sphere interaction_ , Phys. Rev. 105 (1957), no. 3, 767–775.
  * (21) K. Huang, C. N. Yang, and J. M. Luttinger, _Imperfect Bose gas with hard-sphere interaction_ , Phys. Rev. 105 (1957), no. 3, 776–784.
  * (22) T. Kato, _Perturbation Theory for Linear Operators_ , 2nd ed., Springer-Verlag, Berlin, 1995.
  * (23) V. Kostrykin, K.. A. Makarov, and A. K Motovilov, _Existence and uniqueness of solutions to the operator Riccati equation. A geometric approach_ , Advances in Differential Equations and Mathematical Physics (E. Carlen, E. M. Harrell, and M. Loss, eds.), Contemporary Mathematics, vol. 327, Amer. Math. Soc., Providence, RI, 2003, pp. 181–198.
  * (24) T. D. Lee, K. Huang, and C. N. Yang, _Eigenvalues and eigenfunctions of a Bose System of hard spheres and its low-temperature properties_ , Phys. Rev. 106 (1957), no. 6, 1135–1145.
  * (25) T. D. Lee and C. N. Yang, _Low-temperature behavior of a dilute Bose system of hard spheres. I. Equilibrium properties_ , Phys. Rev. 112 (1958), no. 5, 1419–1429.
  * (26) Anthony J. Leggett, _Bose-Einstein condensation in the alkali gases: Some fundamental concepts_ , Rev. Mod. Phys. 73 (2001), no. 2, 307–358.
  * (27) E. H. Lieb, R. Seiringer, J. P. Solovej, and J. Yngvanson, _The Mathematics of the Bose gas and its Condensation_ , Birkhäuser, Basel, Switzerland, 2005.
  * (28) S. Sorokanich M. G. Grillakis, D. Margetis, _Many-body excitations in trapped bose gas: A non-hermitian view_ , preprint (2021).
  * (29) D. Margetis, _Bose-Einstein condensation beyond mean field: Many-body state of periodic microstructure_ , Multiscale Model. Simul. 10 (2012), no. 2, 383–417.
  * (30) P. T. Nam and M. Napiórkowski, _Bogoliubov correction to the mean-field dynamics of interacting Bosons_ , Adv. Theor. Math. Phys. 108 (2017), no. 3, 683–738.
  * (31) , _A note on the validity of Bogoliubov correction to mean-field dynamics_ , J. Math. Pure Appl. 108 (2017), no. 5, 662–688.
  * (32) P. T. Nam, M. Napiórkowski, and J. P. Solovej, _Diagonalization of bosonic quadratic Hamiltonians by Bogoliubov transformations_ , J. Funct. Anal. 270 (2016), no. 11, 4340–4368.
  * (33) P. T. Nam and R. Seiringer, _Collective excitations of Bose gases in the mean-field regime_ , Arch. Rat. Mech. Anal. 215 (2015), no. 2, 381–417.
  * (34) Phan Thành Nam and Marcin Napiórkowski, _Norm approximation for many-body quantum dynamics and bogoliubov theory_ , pp. 223–238, Springer International Publishing, Cham, 2017.
  * (35) R. Ozeri, N. Katz, J. Steinhauer, and N. Davidson, _Colloquium: Bulk Bogoliubov excitations in a Bose-Einstein condensate_ , Rev. Mod. Phys. 77 (2005), no. 1, 187–205.
  * (36) L. P. Pitaevskii, _Vortex lines in an imperfect Bose gas_ , Soviet Phys. JETP 13 (1961), no. 2, 451–454.
  * (37) A. Rovenchak, _Effective Hamiltonian and excitation spectrum of harmonically trapped bosons_ , Low Temp. Phys. 42 (2016), no. 1, 36–41.
  * (38) R. Seiringer, _Bose gases, Bose-Einstein condensation, and the Bogoliubov approximation_ , J. Math. Phys. 55 (2014), no. 7, 075209.
  * (39) Stephen Sorokanich, _Non-Hermitian approaches for pair-excitation in quantum Boson dynamics_ , 2022, PhD Thesis.
  * (40) Christiane Tretter, _Spectral Theory of Block Operator Matrices_ , 1st ed., Imperial College Press, London, UK, 2008.
  * (41) T. T. Wu, _Ground state of a Bose system of hard spheres_ , Phys. Rev. 115 (1959), no. 6, 1390–1404.
  * (42) , _Some nonequilibrium properties of a Bose system of hard spheres at extremely low temperatures_ , J. Math. Phys. 2 (1961), no. 1, 105–123.
  * (43) , _Bose-Einstein condensation in an external potential at zero temperature: General theory_ , Phys. Rev. A 58 (1998), no. 2, 1465–1474.
  * (44) V. A. Zagrebnov and J.-B. Bru, _The Bogoliubov model of weakly imperfect Bose gas_ , Phys. Rep. 350 (2001), 291–434.