# Permanent variational wave functions for bosons

J. M. Zhang<EMAIL_ADDRESS>Fujian Provincial Key Laboratory of Quantum
Manipulation and New Energy Materials, College of Physics and Energy, Fujian
Normal University, Fuzhou 350007, China  Fujian Provincial Collaborative
Innovation Center for Optoelectronic Semiconductors and Efficient Devices,
Xiamen, 361005, China  H. F. Song<EMAIL_ADDRESS>Laboratory of
Computational Physics, Institute of Applied Physics and Computational
Mathematics, Beijing 100088, China  Y. Liu<EMAIL_ADDRESS>Laboratory of
Computational Physics, Institute of Applied Physics and Computational
Mathematics, Beijing 100088, China

###### Abstract

We study the performance of permanent states (the bosonic counterpart of the
Slater determinant state) as approximating functions for bosons, with the
intention to develop variational methods based upon them. For a system of $N$
identical bosons, a permanent state is constructed by taking a set of $N$
arbitrary (not necessarily orthonormal) single-particle orbitals, forming
their product and then symmetrizing it. It is found that for the one-
dimensional Bose-Hubbard model with the periodic boundary condition and at
unit filling, the exact ground state can be very well approximated by a
permanent state, in that the permanent state has high overlap (at least 0.96
even for 12 particles and 12 sites) with the exact ground state and can
reproduce both the ground state energy and the single-particle correlators to
high precision. For a generic model, we have devised a greedy algorithm to
find the optimal set of single-particle orbitals to minimize the variational
energy or maximize the overlap with a target state. It turns out that quite
often the ground state of a bosonic system can be well approximated by a
permanent state by all the criterions of energy, overlap, and correlation
functions. And even if the error is apparent, it can often be remedied by
including more configurations, i.e., by allowing the variational wave function
to be a combination of multiple permanent states. The algorithm is used to
study the stability of a two-particle system, with great success. All these
suggest that permanent states are very effective as variational wave functions
for bosonic systems, and hence deserve further studies.

##  I Introduction

The Hartree-Fock approximation for fermions is a paradigm in quantum mechanics
hartree ; fock ; gaunt ; slater1 ; slater2 ; ostlund . Conceptually, it is
very simple. It is a variational method. For an $N$-fermion system, one just
takes $N$ orthonormal single-particle orbitals $\\{\phi_{i}|1\leq i\leq N\\}$,
constructs the product state
$\phi_{1}(x_{1})\phi_{2}(x_{2})\ldots\phi_{N}(x_{N})$, and then anti-
symmetrizes it to obtain the Slater determinant wave function

$\displaystyle\Xi(x_{1},x_{2},\ldots,x_{N})$ (1) $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N!}}\sum_{P\in
S_{N}}(-1)^{P}\phi_{P_{1}}(x_{1})\phi_{P_{2}}(x_{2})\ldots\phi_{P_{N}}(x_{N})$
$\displaystyle=$ $\displaystyle\frac{1}{\sqrt{N!}}\det(\phi_{i}(x_{j})).$

Here $S_{N}$ denotes the symmetric group of degree $N$. By construction, the
determinant state satisfies the anti-symmetry condition and constitutes a
legitimate wave function for a collection of identical fermions. With the
variational wave function built in this way, the rest work is an optimization
problem. One has to choose the $N$ orthonormal orbitals optimally so as to
minimize the energy expectation value of the $N$-body variational state.

It is a natural idea to generalize this approach to bosons. One can take $N$
single-particle orbitals $\\{\phi_{i}|1\leq i\leq N\\}$, form their product,
but then symmetrize it to obtain the following state,

$\displaystyle\Phi(x_{1},x_{2},\ldots,x_{N})$ (2) $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N!}}\sum_{P\in
S_{N}}\phi_{P_{1}}(x_{1})\phi_{P_{2}}(x_{2})\ldots\phi_{P_{N}}(x_{N})$
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N!}}\operatorname{per}(\phi_{i}(x_{j})),$

which we shall refer to as a permanent state. Unlike the fermionic case, here
because of the symmetry instead of anti-symmetry condition, the single-
particle orbitals are not necessarily orthogonal to each other fermi , but
could even be identical. In the extremal case in which all the orbitals are
constrained to be the same, we have the Gross-Pitaevskii approximation gross ;
pitaevskii , which has been proven to be very successful for weakly
interacting bose gases rmp1 ; rmp2 . However, for more general systems, such
as the Bose-Hubbard model which we shall study below, the Gross-Pitaevskii
approximation is too restrictive and we had better allow more freedom for the
$N$ orbitals.

The idea seems very simple. However, probably because mathematically the
permanent of a matrix lacks many of the nice properties of the determinant,
such an approach has rarely been put into practice. As far as we know, the
very limited literature starts with two papers of Romanovsky _et al._ in 2004
and 2006 igor1 ; igor2 . They employed the permanent state as variational wave
functions for some few-boson systems in two-dimensional harmonic traps. They
went beyond the Gross-Pitaevskii approximation by allowing each particle to
occupy a different orbital, for which they coined the term unrestricted Bose-
Hartree-Fock approximation. However, because of the perceived high complexity
of the self-consistency equations, they did not seek self-consistent orbitals,
but prescribed them as displaced Gaussians. Subsequently, the self-consistency
equations for the orbitals were derived by Heimsoth martin1 ; martin2 .
Unfortunately, the formalism was still unnecessarily complicated, and he did
not even implemented the Ryser algorithm for permanent computation.
Consequently, he could handle at most six particles.

In this paper, we resume research in this vein but with different
perspectives. For some reason, people almost always look at orthonormal
single-particle orbitals, and henceforth construct orthonormal Fock states as
basis functions for a multi-particle system. But working with orthonormal
orbitals is _artificial_ —The bosons themselves have no notion of
orthogonality and are totally happy to reside in non-orthogonal orbitals.
Therefore, for studying the structure of a bosonic wave function, one should
give up the obsession of orthogonal orbitals. It is absolutely possible that a
bosonic wave function looking complex with respect to a Fock basis is actually
of a simple and compact structure, namely, is equal or close to a permanent
state built of non-orthogonal orbitals. As we shall see below, this is indeed
the case for many models, in particular, for the one-dimensional Bose-Hubbard
model at unit filling. An immediate implication of such fortunate facts is
that permanent states can be used as building blocks for constructing
variational wave functions and it is worthwhile to develop variational methods
based upon them. Of course, besides this pragmatic purpose, permanent states
are of interest in their own right. They are the simplest bosonic wave
functions, yet of rich structures, with information (energy, correlation
functions, etc.) not easy to extract sometimes.

This paper is organized as follows. First in Sec. II, we review the connection
between the first and second quantization formalisms, and we shall establish
some analytic facts about the permanent state. Then in Sec. III, we show that
for the one-dimensional Bose-Hubbard model with periodic boundary condition
and at unit filling, which is the standard setting for studying the
superfluid-Mott insulator transition, the permanent state can be a very good
approximation of the exact ground state. It is good not only by the usual
energy criterion, but also by the more stringent criterions of overlap and
correlation functions. For a Bose-Hubbard model with 12 particles on 12 sites,
the energy-minimizing permanent state with prescribed orbitals has an overlap
with the exact ground state as large as $0.96$ in the worst case. Of course
one should not be satisfied with prescribed orbitals. It is desirable to have
more flexibility and presumably the numbers could be further improved if the
orbitals are really unrestricted. We thus propose an iterative algorithm in
Sec. IV for searching for the optimal set of orbitals minimizing the energy.
The equations are equivalent to what Heimsoth derived martin1 ; martin2 .
However, because of the different point of view, our derivation is more
elementary and straightforward, and the formulation is more amenable for
numerical implementation. Moreover, our formalism can handle easily the multi-
configuration case, i.e., the case when the variational wave function is a
combination of multiple permanent states. Note that this was not considered
previously, but is necessary and effective for improving accuracy. The
algorithm can actually be employed to solve another optimization problem,
namely, for a given wave function, finding the single- or multi-
configurational variational wave function most close to it, i.e., having the
largest possible overlap with it. Although this problem is rarely studied in
the literature and is not our focus in this paper, it should be a meaningful
question for studying the structure of a bosonic wave function. With the
optimization algorithm, we can tackle more general models. This is what we do
in Sec. V. We shall see that in many cases, a single- or multi-configurational
variational wave function is a very good approximation of the exact ground
state of the system. This enables us to use the algorithm to study the
stability of a two-boson system, which is analogous to the negative ion of
hydrogen. Finally, we conclude in Sec. VI with some open problems. We would
like to mention that the whole paper is actually a by-product of studying
these open problems.

##  II Permanent wave functions

For the sake of simplicity, let us assume a finite-dimensional single-particle
Hilbert space

$\displaystyle\mathcal{H}=\operatorname{span}\\{|x\rangle,1\leq x\leq L\\},$
(3)

where $|x\rangle$ are orthonormal basis vectors. The associated creation
(annihilation) operators will be denoted as $a_{x}^{\dagger}$ ($a_{x}$). They
satisfy the usual commutation relations. A generic (not necessarily
normalized) single-particle state or a single-particle orbital in this space
is $|\phi\rangle=\sum_{x=1}^{L}|x\rangle\langle
x|\phi\rangle=\sum_{x=1}^{L}\phi(x)|x\rangle$. The associated creation
operator is $a_{\phi}^{\dagger}=\sum_{x=1}^{L}\phi(x)a_{x}^{\dagger}$.

For an $N$-boson system, the many-body Hilbert space is spanned by the
orthonormal Fock states

$\displaystyle|\textbf{n}\rangle$ $\displaystyle=$
$\displaystyle\frac{(a_{1}^{\dagger})^{n_{1}}(a_{2}^{\dagger})^{n_{2}}\ldots(a_{L}^{\dagger})^{n_{L}}}{\sqrt{n_{1}!n_{2}!\ldots
n_{L}!}}|vac\rangle.$ (4)

where $\textbf{n}\equiv(n_{1},n_{2},\ldots,n_{L})$ is an $L$-tuple with
$n_{x}\geq 0$ and $\sum_{x=1}^{L}n_{x}=N$. The number of such Fock states or
the dimension of the many-body Hilbert space is

$\displaystyle\mathcal{D}$ $\displaystyle=$
$\displaystyle\binom{N+L-1}{N}=\frac{(N+L-1)!}{N!(L-1)!}.$ (5)

A generic $N$-boson state $|\Psi\rangle$ expands as
$|\Psi\rangle=\sum_{\textbf{n}}C(\textbf{n})|\textbf{n}\rangle$. In first
quantization, the same state is expressed by the wave function
$\Psi(x_{1},x_{2},\ldots,x_{N})$, with $1\leq x_{i}\leq L$. Under the action
of particle permutations, the coordinate tuples
$\textbf{x}\equiv(x_{1},x_{2},\ldots,x_{N})$ break into different equivalent
classes labelled by the occupation tuple n. The wave function should be
constant on each class. We say a coordinate tuple x belongs to n and denote it
as $\textbf{x}\vdash\textbf{n}$ if in x, the value $x$ appears $n_{x}$ times.
The cardinality of the class n is $N!/\prod_{x=1}^{L}n_{x}!$, and thus by
considering the norm of $\Psi$ in both the first and second quantization form,
we have bogoliubov

$\displaystyle\Psi(\textbf{x})$ $\displaystyle=$
$\displaystyle\sqrt{\frac{n_{1}!n_{2}!\ldots
n_{L}!}{N!}}C(\textbf{n}),\quad\textbf{x}\in\textbf{n}.$ (6)

With this formula, one can convert a wave function in the first quantization
form to the second quantization form, and vice versa.

Now suppose we have a set of $N$ arbitrary orbitals $\\{\phi_{i}(x),1\leq
i\leq N\\}$. The simplest symmetric $N$-particle wave function one can
construct out of them is the permanent state in (2). We shall use the notation

$\displaystyle\Phi$ $\displaystyle=$
$\displaystyle\hat{\mathcal{S}}(\phi_{1},\phi_{2},\ldots,\phi_{N})$ (7)

to indicate that $\Phi$ is built out of the orbitals $\\{\phi_{i}(x),1\leq
i\leq N\\}$ according to the product and symmetrization procedure in (2). By
the correspondence (6), it is easy to show that in second quantization, this
state has the expression

$\displaystyle|\Phi\rangle$ $\displaystyle=$
$\displaystyle\prod_{i=1}^{N}a_{\phi_{i}}^{\dagger}|vac\rangle.$ (8)

This form should reminds us of the standard Fock states. They are also
permanent states, but with orthonormal orbitals.

Here we emphasize that in this paper _we abandon the orthogonality and
normalization of the orbitals_. This is a fundamental difference between the
current approach and the multiconfigurational Hartree theory for bosons (MCHB)
cederbaum , which although also treats the single-particle orbitals as
variational parameters, insists on their orthogonality and normalization. One
apparent advantage of using non-orthogonal orbitals is that the expression of
the wave function is more _compact_ farid —A single permanent state in the
form of (8) would expand into a multitude of Fock states if the orbitals
$\phi_{i}$ are expanded in terms of an orthonormal basis. The downside is that
we have to compute the permanents of overlap matrices, which is expensive in
CPU time.

###  II.1 Five simple propositions

At least five simple facts about a permanent state can be easily established.
Although some of them are not much used in this paper, we collect them all
here for completeness.

In the fermionic case, because of a well-known property of the determinant,
the determinant state (1) vanishes identically if the $N$ orbitals are not
linearly independent. However, such a perfect cancellation cannot occur for
the permanent state in (2) or (8) in the case of bosons. Although this might
not be surprising as the terms in (2) are all of the same sign, we have
formulated it as a proposition, which will be referred to later.

######  Proposition 1.

The $N$-particle permanent state $\Phi$ constructed with $N$ nonzero orbitals
$\\{\phi_{1\leq i\leq N}\\}$ according to (2) is necessarily non-vanishing.

###### Proof.

Consider the inner product of $\Phi$ with the condensate-type symmetric state

$\displaystyle\Omega$ $\displaystyle=$ $\displaystyle v(x_{1})v(x_{2})\ldots
v(x_{N}),$ (9)

where $v\in\mathcal{H}$ is an arbitrary single-particle state. We have

$\displaystyle\frac{1}{\sqrt{N!}}\langle\Phi|\Omega\rangle=\prod_{i=1}^{N}\langle\phi_{i}|v\rangle.$
(10)

The product on the right hand vanishes if and only if $v$ is orthogonal to
some $\phi_{i}$, or $v\in\bigcup_{i=1}^{N}\ker(\phi_{i})$. Here by an abuse of
notation, we have also used $\phi_{i}$ to denote the linear functional
$\langle\phi_{i}|\cdot\rangle$. The kernel of this functional $\ker(\phi_{i})$
is simply the hyperplane orthogonal to $\phi_{i}$. By the well-known folklore
in mathematics that a vector space over $\mathbb{C}$ cannot be the union of a
finite number of proper subspaces rotman , we know the union
$\bigcup_{i=1}^{N}\ker(\phi_{i})$ cannot cover the whole space $\mathcal{H}$,
and for some $v$ the product is non-vanishing, which in turn means that $\Phi$
must be non-vanishing. ∎

The second one is about uniqueness. Given a set of orbitals $\\{\phi_{i},1\leq
i\leq N\\}$, one can construct a permanent state according to (2) or (8). One
might ask whether the same permanent state can be built with a different set
of orbitals. Here of course, two sets of orbitals should be deemed equivalent
if they differ just by a permutation or some linear scaling. The answer is no
as we have

######  Proposition 2.

Suppose a non-vanishing $N$-particle permanent state $\Phi$ can be constructed
with two sets of orbitals $\\{f_{1\leq i\leq N}\\}$ and $\\{g_{1\leq i\leq
N}\\}$, i.e.,

$\displaystyle\Phi$ $\displaystyle=$
$\displaystyle\hat{\mathcal{S}}(f_{1},f_{2},\ldots,f_{N})=\hat{\mathcal{S}}(g_{1},g_{2},\ldots,g_{N}),$
(11)

then for some permutation $\sigma\in S_{N}$, $f_{i}\propto g_{\sigma(i)}$ for
all $1\leq i\leq N$.

###### Proof.

We provide two proofs. The first one is elementary but lengthy. Again, let us
consider the inner product of $\Phi$ with the condensate-type state in (9). We
have

$\displaystyle\frac{1}{\sqrt{N!}}\langle\Phi|\Omega\rangle=\prod_{i=1}^{N}\langle
f_{i}|v\rangle$ $\displaystyle=$ $\displaystyle\prod_{i=1}^{N}\langle
g_{i}|v\rangle.$ (12)

Now suppose $v$ is orthogonal to $f_{1}$, i.e., $\langle f_{1}|v\rangle=0$, or
$v\in\ker(f_{1})$. The equality above implies
$v\in\bigcup_{i=1}^{N}\ker(g_{i})$. As this holds for any $v\in\ker(f_{1})$,
we obtain $\ker(f_{1})\subseteq\bigcup_{i=1}^{N}\ker(g_{i})$, which in turn
means

$\displaystyle\ker(f_{1})$ $\displaystyle=$
$\displaystyle\ker(f_{1})\bigcap\left(\bigcup_{i=1}^{N}\ker(g_{i})\right)$
(13) $\displaystyle=$
$\displaystyle\bigcup_{i=1}^{N}\left(\ker(f_{1})\bigcap\ker(g_{i})\right).$

Again by the well-known folklore that a vector space cannot be the union of a
finite number of proper subspaces rotman , we know there must be some $i$ such
that $\ker(f_{1})=\ker(f_{1})\bigcap\ker(g_{i})$, or
$\ker(f_{1})=\ker(g_{i})$, which means $f_{1}\propto g_{i}$. By permutation or
relabeling, we can assume $i=1$ and simply $f_{1}=g_{1}$.

We then wish to factor out $f_{1}$ and $g_{1}$ in (12) to get

$\displaystyle\prod_{i=2}^{N}\langle f_{i}|v\rangle$ $\displaystyle=$
$\displaystyle\prod_{i=2}^{N}\langle g_{i}|v\rangle$ (14)

for all $v\in\mathcal{H}$. This already holds for $v\notin\ker(f_{1})$. To
show that it actually holds also for $v\in\ker(f_{1})$, consider
$v\in\ker(f_{1})$ and $w\notin\ker(f_{1})$. For any nonzero $t\in\mathbb{C}$,
$v+tw\notin\ker(f_{1})$, otherwise $w=[(v+tw)-v]/t$ would be in $\ker(f_{1})$.
By (14), we have then

$\displaystyle\prod_{i=2}^{N}\langle f_{i}|v+tw\rangle$ $\displaystyle=$
$\displaystyle\prod_{i=2}^{N}\langle g_{i}|v+tw\rangle$ (15)

for all $t\neq 0$. That is, two polynomials of $t$ evaluate to the same value
for all $t\neq 0$. This means that the two polynomials are actually the same.
In particular, their constant terms are identical. We have thus proven that
(14) holds for all $v\in\ker(f_{1})$ too. The proposition is then proven by
induction.

The second proof is more direct. Let
$|f_{i}\rangle=\sum_{j=1}^{L}f_{ij}|j\rangle$ and similarly
$|g_{i}\rangle=\sum_{j=1}^{L}g_{ij}|j\rangle$. By (8) and (11),

$\displaystyle\Phi$ $\displaystyle=$
$\displaystyle\prod_{i=1}^{N}\left(\sum_{j=1}^{L}f_{ij}a_{j}^{\dagger}\right)|vac\rangle=\prod_{i=1}^{N}\left(\sum_{j=1}^{L}g_{ij}a_{j}^{\dagger}\right)|vac\rangle.\quad$
(16)

The fact that the $a_{j}^{\dagger}$ operators commute and the Fock basis
states are linearly independent implies the polynomial equality

$\displaystyle\prod_{i=1}^{N}\left(\sum_{j=1}^{L}f_{ij}z_{j}\right)=\prod_{i=1}^{N}\left(\sum_{j=1}^{L}g_{ij}z_{j}\right),$
(17)

where $z_{j}$ are indeterminates. The proposition is proven by the unique
factorization theorem of multivariate polynomials over $\mathbb{C}$ fraleigh .
∎

To appreciate Proposition 2, one should note that for fermions, a Slater
determinant state constructed out of $N$ orbitals is determined by the
subspace (a point on the so-called Grassmannian manifold aoto ) spanned by the
orbitals, not by the orbitals themselves. The orbitals are just a basis of the
subspace. Another basis would yield the same Slater determinant state up to a
global constant.

By Proposition 2, the permanent state $\Phi$ is invariant under the
permutation $f_{i}\rightarrow f_{\sigma_{i}}$ with $\sigma$ being an arbitrary
permutation, and the scaling $f_{i}\rightarrow\lambda_{i}f_{i}$ with the
scaling factors satisfying the condition $\prod_{i=1}^{N}\lambda_{i}=1$. This
allows us to count the degrees of freedom manifold  of a permanent state as

$\displaystyle d=NL-(N-1)=N(L-1)+1.$ (18)

Apparently, like not every $N$-fermion state is a Slater determinant state,
not every $N$-boson state is a permanent state. However, in the special case
of $\dim\mathcal{H}=L=2$, this is indeed the case. We have

######  Proposition 3.

Suppose the single-particle Hilbert space is of dimension 2, i.e.,
$\dim\mathcal{H}=L=2$, then every $N$-boson state is a permanent state.

###### Proof.

In this case, the $N$-boson Hilbert space is of dimension $N+1$ and a basis is
the Fock states
$\\{(a_{1}^{\dagger})^{i}(a_{2}^{\dagger})^{N-i}|vac\rangle,0\leq i\leq N\\}$.
But by (18), the number of degrees of freedom of a permanent state is also
$N+1$, so the proposition is anticipated by dimension counting. To prove it
rigorously, we expand an arbitrary state $\Psi$ as

$\displaystyle\Psi$ $\displaystyle=$
$\displaystyle(a_{2}^{\dagger})^{m}\sum_{k=0}^{N-m}c_{k}(a_{1}^{\dagger})^{N-m-k}(a_{2}^{\dagger})^{k}|vac\rangle.$
(19)

Here $m$ is an integer between $0$ and $N$, and $c_{0}\neq 0$. The fact that
$a_{1}^{\dagger}$ and $a_{2}^{\dagger}$ commute motivates us to consider the
polynomial $P(z)=\sum_{k=0}^{N-m}c_{k}z^{N-m-k}$. Let it factorize as

$\displaystyle P(z)$ $\displaystyle=$ $\displaystyle
c_{0}\prod_{j=1}^{N-m}(z-z_{j}).$ (20)

Then it is easy to see that $\Psi$ factorize as

$\displaystyle\Psi$ $\displaystyle=$ $\displaystyle
c_{0}(a_{2}^{\dagger})^{m}\prod_{j=1}^{N-m}(a_{1}^{\dagger}-z_{j}a_{2}^{\dagger})|vac\rangle.$
(21)

By (8), we see it is a permanent state and can read off the single-particle
orbitals. ∎

Proposition 3 means that for a two-site Bose-Hubbard model foerster , which is
a canonical model for studying the Bose-Josephson effect, any state is a
permanent state. The proposition also reminds us of a similar proposition for
fermions zhang1 ; zhang2 . Although an arbitrary fermionic wave function is
not generally a Slater determinant state, for the case of $N$ fermions in
$L=N+1$ orbitals, the wave function is necessarily a Slater determinant.

Proposition 3 is about the case when the dimension of the single-particle
Hilbert space is minimal (but still nontrivial). Similarly, when the particle
number is minimal (but still nontrivial), i.e., when $N=2$, we have

######  Proposition 4.

Suppose $N=2$ and $\dim\mathcal{H}=L$, then every $N$-boson state can be
written as the sum of at most $\lfloor(L+1)/2\rfloor$ permanent states. Here
$\lfloor x\rfloor$ is the floor function denoting the greatest integer less
than or equal to $x$.

###### Proof.

For $N=2$, the bosonic wave function $\Psi(x_{1},x_{2})$ with $1\leq x_{i}\leq
L$ can be considered as a complex symmetric matrix. By the Autonne-Takagi
theorem horn , it can be factorized as

$\displaystyle\Psi(x_{1},x_{2})$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{L}\sqrt{D_{j}}f_{j}(x_{1})f_{j}(x_{2}),$ (22)

where $f_{j}$ are a set of orthonormal functions and $D_{j}$ are a set of non-
negative numbers in decreasing order. The functions $f_{j}$ are actually the
so-called natural orbitals, i.e., eigenvectors of the one-body reduced density
matrix $\rho=2\Psi\Psi^{\dagger}$ associated with $\Psi$, and $2D_{j}$ are the
occupation numbers. In the following, we shall assume that $\Psi$ is
normalized, i.e., $\langle\Psi|\Psi\rangle=1$, so that
$\sum_{j=1}^{L}D_{j}=1$.

In the form of (22), the wave function $\Psi$ is already a sum of $L$
permanent states. Now for any pair $j\neq k$, we can combine
$\sqrt{D_{j}}f_{j}(x_{1})f_{j}(x_{2})$ and
$\sqrt{D_{k}}f_{k}(x_{1})f_{k}(x_{2})$ into a single permanent state, i.e.,

$\displaystyle\hat{\mathcal{S}}(u,v)=\frac{1}{\sqrt{2}}\left[u(x_{1})v(x_{2})+v(x_{1})u(x_{2})\right]$

with (note that by Proposition 2, this is essentially the only solution)

$\displaystyle u$ $\displaystyle\equiv$
$\displaystyle\frac{1}{\sqrt[4]{2}}(\sqrt[4]{D_{j}}f_{j}+i\sqrt[4]{D_{k}}f_{k}),$
$\displaystyle v$ $\displaystyle\equiv$
$\displaystyle\frac{1}{\sqrt[4]{2}}(\sqrt[4]{D_{j}}f_{j}-i\sqrt[4]{D_{k}}f_{k}).$

The proposition is then proven by noting that in (22), when $L$ is even, we
have $L/2$ pairs and when $L$ is odd, we have $(L-1)/2$ pairs and an extra
unpaired term, which is already a permanent state. ∎

It is easy to see that the number $\lfloor(L+1)/2\rfloor$ in Proposition 4 is
not only sufficient but also necessary for a generic $2$-boson state. For a
generic state, the matrix $\Psi$ is nonsingular or full-ranked. On the other
hand, a permanent state is at most 2-ranked. Hence, by the subadditivity
property of the rank of a matrix
[$\operatorname{rank}(A+B)\leq\operatorname{rank}(A)+\operatorname{rank}(B)$],
we need at least $\lfloor(L+1)/2\rfloor$ permanent states to fully recover the
original state. However, in practice, one might just need a sufficiently good
approximation of the original state. The question is then, by taking the sum
of $M<\lfloor(L+1)/2\rfloor$ permanent states, i.e., by constructing a state
in the form of

$\displaystyle\Omega$ $\displaystyle=$
$\displaystyle\sum_{\alpha=1}^{M}\hat{\mathcal{S}}(u^{(\alpha)},v^{(\alpha)}),$
(23)

where $u^{(\alpha)}$ and $v^{(\alpha)}$ are arbitrary orbitals, to what extent
can we approximate a target function $\Psi$? Quantitatively, what is the
maximal value of the overlap

$\displaystyle O$ $\displaystyle=$
$\displaystyle\frac{|\langle\Omega|\Psi\rangle|^{2}}{\langle\Omega|\Omega\rangle\langle\Psi|\Psi\rangle}$
(24)

achievable with such an $M$-configuration state? For this problem, we have

######  Proposition 5.

Suppose $N=2$, $\dim\mathcal{H}=L$, and $M\leq\lfloor(L+1)/2\rfloor$. For the
target state $\Psi$ in (22), the largest possible value of the overlap of an
$M$-configuration state in the form of (23) with it is

$\displaystyle O_{max}$ $\displaystyle=$
$\displaystyle\max_{\Omega}\frac{|\langle\Omega|\Psi\rangle|^{2}}{\langle\Omega|\Omega\rangle\langle\Psi|\Psi\rangle}=\sum_{j=1}^{2M}D_{j}.$
(25)

###### Proof.

Like $\Psi$ in (22), the $M$-configuration state $\Omega$ in (23) is also
symmetric and can also be factorized as

$\displaystyle\Omega$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{2M}\sqrt{C_{j}}\varphi_{j}(x_{1})\varphi_{j}(x_{2}),$
(26)

where $\varphi_{1\leq j\leq 2M}$ are a set of orthonormal vectors which can be
extended into a complete orthonormal basis $\varphi_{1\leq j\leq L}$, and
$C_{j}\geq 0$ are ordered in decreasing order. Without loss of generality, let
us assume $\Omega$ is normalized so that $\sum_{j=1}^{2M}C_{j}=1$. Note that
by construction $\Omega$ is at most of rank $2M$, and thus here in (26) we
have at most $2M$ nonzero terms. For the overlap between $\Omega$ and $\Psi$,
we have

$\displaystyle O$ $\displaystyle=$
$\displaystyle\left|\sum_{j=1}^{2M}\sqrt{C_{j}}\langle\varphi_{j}\varphi_{j}|\Psi\rangle\right|^{2}$
$\displaystyle\leq$
$\displaystyle\left(\sum_{j=1}^{2M}C_{j}\right)\left(\sum_{j=1}^{2M}|\langle\varphi_{j}\varphi_{j}|\Psi\rangle|^{2}\right)$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{2M}|\langle\varphi_{j}\varphi_{j}|\Psi\rangle|^{2}\leq\sum_{j=1}^{2M}\sum_{k=1}^{L}|\langle\varphi_{j}\varphi_{k}|\Psi\rangle|^{2}=\operatorname{Tr}(AA^{\dagger}),$

where $A\equiv V^{\dagger}\Psi W^{*}$, with $V$ the $L\times 2M$ matrix whose
columns are $\varphi_{1\leq j\leq 2M}$ and $W$ the $L\times L$ matrix whose
columns are $\varphi_{1\leq k\leq L}$. Note that $W$ is unitary. We have

$\displaystyle\operatorname{Tr}(AA^{\dagger})$ $\displaystyle=$
$\displaystyle\operatorname{Tr}(V^{\dagger}\Psi W^{*}W^{T}\Psi^{\dagger}V)$
(27) $\displaystyle=$
$\displaystyle\operatorname{Tr}(V^{\dagger}\Psi\Psi^{\dagger}V)\leq\sum_{j=1}^{2M}D_{j}.$

Here we used the Ky-Fan inequality bhatia  for the matrix
$\Psi\Psi^{\dagger}$, which states the for an $L\times L$ hermitian matrix,
the sum of its $m$ diagonal elements is less than or equal to the sum of its
$m$ largest eigenvalues, for all $1\leq m\leq L$. We have thus proven that the
overlap is upper bounded by $\sum_{j=1}^{2M}D_{j}$. That this upper bound can
be achieved is obvious, as we can just take the first $2M$ terms in (22) and
by Proposition 4 it is an $M$-configuration state. ∎

In general, the occupation numbers $D_{j}$ decrease fast and a truncation of
(22) with a very limited number of terms can yield a good approximation of the
original state.

###  II.2 Basic formulae

A generic many-body Hamiltonian is of the form

$\displaystyle H=H_{1}+H_{2}=\sum_{i=1}^{N}K(i)+\sum_{1\leq i<j\leq N}U(i,j).$
(28)

Here the first sum is over each particle, with $K(i)$ denoting the sum of the
kinetic energy and the external potential of the $i$th particle, while the
second sum is over each pair, with $U(i,j)$ denoting the interaction between
the $i$th and $j$th particle. The quantity of primary interest is the
expectation value of $H$ with respect to a permanent state (2), i.e.,

$\displaystyle E_{var}$ $\displaystyle=$
$\displaystyle\frac{\langle\Phi|H|\Phi\rangle}{\langle\Phi|\Phi\rangle}=\frac{\langle\Phi|H_{1}|\Phi\rangle+\langle\Phi|H_{2}|\Phi\rangle}{\langle\Phi|\Phi\rangle}.$
(29)

It is straightforward to calculate the denominator and the numerator here. But
for the sake of generality and in view of the extension to the
multiconfigurational case below, let us consider the off-diagonal matrix
elements instead of the diagonal matrix elements. Suppose $\Phi^{(1)}$ and
$\Phi^{(2)}$ are two permanent states constructed with two sets of orbitals
$\\{\phi^{(\alpha)}_{1\leq i\leq N},\alpha=1,2\\}$, i.e.,

$\displaystyle\Phi^{(\alpha)}$ $\displaystyle=$
$\displaystyle\hat{\mathcal{S}}(\phi_{1}^{(\alpha)},\phi_{2}^{(\alpha)},\ldots,\phi_{N}^{(\alpha)}),\quad\alpha=1,2.$
(30)

Let us first define the $N\times N$ overlap matrix of the orbitals,

$\displaystyle A_{ij}$ $\displaystyle=$
$\displaystyle\langle\phi_{i}^{(1)}|\phi_{j}^{(2)}\rangle,\quad 1\leq i,j\leq
N.$ (31)

We have then

$\displaystyle\langle\Phi^{(1)}|\Phi^{(2)}\rangle$ (32) $\displaystyle=$
$\displaystyle\frac{1}{N!}\sum_{P,Q\in
S_{N}}\langle\phi^{(1)}_{P_{1}}|\phi^{(2)}_{Q_{1}}\rangle\langle\phi^{(1)}_{P_{2}}|\phi^{(2)}_{Q_{2}}\rangle\ldots\langle\phi^{(1)}_{P_{N}}|\phi^{(2)}_{Q_{N}}\rangle$
$\displaystyle=$ $\displaystyle\sum_{R\in
S_{N}}\langle\phi^{(1)}_{1}|\phi^{(2)}_{R_{1}}\rangle\langle\phi^{(1)}_{2}|\phi^{(2)}_{R_{2}}\rangle\ldots\langle\phi^{(1)}_{N}|\phi^{(2)}_{R_{N}}\rangle$
$\displaystyle=$ $\displaystyle\operatorname{per}(A).$

Here $\operatorname{per}(A)$ denotes the permanent of the matrix $A$. As for
the matrix elements of $H_{1}$,

$\displaystyle\langle\Phi^{(1)}|H_{1}|\Phi^{(2)}\rangle$ $\displaystyle=$
$\displaystyle\frac{N}{N!}\sum_{P,Q\in
S_{N}}\langle\phi^{(1)}_{P_{1}}|K|\phi^{(2)}_{Q_{1}}\rangle\prod_{i=2}^{N}\langle\phi^{(1)}_{P_{i}}|\phi^{(2)}_{Q_{i}}\rangle$
(33) $\displaystyle=$
$\displaystyle\sum_{i_{1}=1}^{N}\sum_{j_{1}=1}^{N}\langle\phi_{i_{1}}^{(1)}|K|\phi_{j_{1}}^{(2)}\rangle\operatorname{per}(A;i_{1}|j_{1}),\quad\quad$

where $\operatorname{per}(A;i_{1}|j_{1})$ denotes the permanent of the
$(N-1)\times(N-1)$ minor obtained by deleting the $i_{1}$th row and the
$j_{1}$th column from $A$. Similarly, for the matrix elements of $H_{2}$,

$\displaystyle\langle\Phi^{(1)}|H_{2}|\Phi^{(2)}\rangle$ (34) $\displaystyle=$
$\displaystyle\frac{N(N-1)}{2(N!)}\sum_{P,Q\in
S_{N}}\langle\phi^{(1)}_{P_{1}}\phi^{(1)}_{P_{2}}|U|\phi^{(2)}_{Q_{1}}\phi^{(2)}_{Q_{2}}\rangle\prod_{i=3}^{N}\langle\phi^{(1)}_{P_{i}}|\phi^{(2)}_{Q_{i}}\rangle$
$\displaystyle=$ $\displaystyle\frac{1}{2}\sum_{i_{1}\neq
i_{2},1}^{N}\sum_{j_{1}\neq
j_{2},1}^{N}\langle\phi^{(1)}_{i_{1}}\phi^{(1)}_{i_{2}}|U|\phi^{(2)}_{j_{1}}\phi^{(2)}_{j_{2}}\rangle$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\times\operatorname{per}(A;i_{1},i_{2}|j_{1},j_{2}).$

Here in the last line, the summation $\sum_{i_{1}\neq i_{2},1}^{N}$ means that
$i_{1}$ and $i_{2}$ both run from $1$ to $N$, but they must take different
values. The second summation is interpreted similarly. By
$\operatorname{per}(A;i_{1},i_{2}|j_{1},j_{2})$ we mean the permanent of the
$(N-2)\times(N-2)$ minor of $A$ obtained by deleting row $i_{1}$, $i_{2}$ and
column $j_{1}$, $j_{2}$.

Above we see that to calculate the norm and physical expectation values of a
permanent state, we have to calculate the permanent of the overlap matrix $A$
and those of its minors. This is the price we have to pay for working with
non-orthonormal orbitals.

As is generally believed, unlike the determinant of a matrix, the permanent of
the matrix cannot be calculated in polynomial time. Currently, the best known
general exact algorithm is the Ryser algorithm ryser , which reduces the naive
$N\cdot N!$ evaluations to $O(N^{2}2^{N-1})$. By using the Gray code, a
further reduction by a factor of $N$ can be achieved fastryser . This is the
algorithm we use in this work code . With this algorithm, it takes about 1 sec
($0.01$ sec) to calculate the permanent of a $22\times 22$ ($15\times 15$,
respectively) real-valued matrix on a commercial laptop computer and with
MATLAB. We mention that in this paper, all calculation is done with MATLAB but
without invoking the parallel computing toolbox.

##  III Bose-Hubbard model at unit filling

To see whether a permanent wave function can be a good approximation of the
ground state of a bosonic system, we take the one-dimensional Bose-Hubbard
model with the periodic boundary condition and at unit filling as a case
study. The Hamiltonian, as is often written in the second-quantization
formalism, reads

$\displaystyle H$ $\displaystyle=$
$\displaystyle-\sum_{x=0}^{L-1}(a_{x}^{\dagger}a_{x+1}+a_{x+1}^{\dagger}a_{x})+\frac{g}{2}\sum_{x=0}^{L-1}a_{x}^{\dagger}a_{x}^{\dagger}a_{x}a_{x}.\quad$
(35)

Here the single-particle Hilbert space $\mathcal{H}$ is spanned by the
orthonormal site (Wannier) states $\\{|x\rangle,0\leq x\leq L-1\\}$. The
periodic boundary condition means $|x\rangle=|x+L\rangle$. Note that we have
taken the hopping strength as the unit of energy, and the Hamiltonian depends
only on the parameter $g\geq 0$, which characterizes the on-site interaction
strength.

For our purpose, it is often more convenient to work with the first-
quantization formalism of (28). The corresponding $K$ and $U$ operators have
matrix elements as

$\displaystyle K_{x,x^{\prime}}$ $\displaystyle=$
$\displaystyle-(\delta_{x,x^{\prime}+1}+\delta_{x,x^{\prime}-1}),$ (36)
$\displaystyle U_{x_{1}x_{2},x_{1}^{\prime}x_{2}^{\prime}}$ $\displaystyle=$
$\displaystyle
g\delta_{x_{1},x_{2}}\delta_{x_{1},x_{1}^{\prime}}\delta_{x_{1},x_{2}^{\prime}}.$
(37)

In this section, we shall confine ourself to the unit filling case, namely,
the case when the particle number $N=L$. The commensurate condition and the
translation symmetry suggest putting the $i$th particle in an orbital centered
at site $i$, with all the orbitals of the same shape and related to each other
by translations. Specifically,

$\displaystyle\phi_{i}(x)$ $\displaystyle=$ $\displaystyle\phi(x-i),\quad\quad
1\leq i\leq N,$ (38)

for some function $\phi$. Here the periodic boundary condition requires
$\phi(x)=\phi(x+L)$. The total variational wave function (VWF) is then
determined by the single-particle orbital $\phi$. To construct an $L$-periodic
function $\phi$, we can choose an arbitrary primitive function $\chi(x)$
defined on the whole axis $-\infty<x<+\infty$, and form the superposition

$\displaystyle\phi(x)$ $\displaystyle=$
$\displaystyle\sum_{j=-\infty}^{\infty}\chi(x-jL).$ (39)

We have tried two types of primitive functions, i.e., the Lorentz function and
the exponential function,

$\displaystyle\chi^{(l)}(x;\lambda)$ $\displaystyle=$
$\displaystyle(1+\lambda^{2}x^{2})^{-1},$ (40a)
$\displaystyle\chi^{(e)}(x;\lambda)$ $\displaystyle=$ $\displaystyle
e^{-\lambda|x|},$ (40b)

where $\lambda\geq 0$ is a parameter controlling the width of the functions.
For these two simple types of primitive functions, the summation in (39) can
be carried out analytically and yields

$\displaystyle\phi^{(l)}(x;\lambda)$ $\displaystyle=$
$\displaystyle\frac{1-e^{-4\pi/L\lambda}}{|1-e^{-2\pi/L\lambda}e^{i2\pi
x/L}|^{2}},$ (41a) $\displaystyle\phi^{(e)}(x;\lambda)$ $\displaystyle=$
$\displaystyle\frac{e^{-\lambda x}+e^{-\lambda(L-x)}}{1-e^{-\lambda L}},\quad
0\leq x\leq L.$ (41b)

Note that in the limit of $\lambda\rightarrow 0$, both $\phi^{(l)}$ and
$\phi^{(e)}$ reduce to the zero-momentum Bloch state on the periodic lattice,
while in the opposite limit of $\lambda\rightarrow\infty$, both of them reduce
to the Kronecker delta function $\delta_{x,0}$.

Figure 1: (Color online) Energy expectation value (per particle) of the
permanent variational state constructed with either the Lorentz-type (40a) or
the exponential-type (40b) primitive orbitals. The horizontal dashed line
indicates the exact ground state energy calculated by exact diagonalization.
The inset shows the optimal orbitals corresponding to the minima of the curves
$E_{var}(\lambda)$. The parameters are $N=L=10$ and $g=4$.  Figure 2: (Color
online) Ground state energy (per particle) estimated with the permanent
variational state constructed with either the Lorentz-type (40a) or the
exponential-type (40b) primitive orbitals. The solid line is the exact value
obtained by exact diagonalization (ED). In the region of $g\lesssim 8$, the
Lorenztian-type orbital yields a better estimate; while in the region of
$g\gtrsim 8$, the exponential-type orbital is better, although this latter
fact is hardly visible in the figure. The parameters are $N=L=10$. The
dimension of the Hilbert space is 92 378.

Figure 3: (Color online) Overlap between (a) the energy-minimizing permanent
variational state $\Phi_{e}$ and (b) the overlap-maximizing permanent
variational state $\Phi_{o}$ constructed with either the Lorentz-type (40a) or
the exponential-type (40b) primitive orbitals with the exact ground state
$|GS\rangle$. We see that for each lattice size $L$, in the region of
$g\lesssim 8$, the Lorent-type variational wave function has a higher overlap
with the exact ground state; while in the region of $g\gtrsim 10$, the
exponential-type variational wave function becomes better.

The strategy is then simply to vary the parameter $\lambda$ and calculate the
variational energy $E_{var}(\lambda;g)$ as a function of $\lambda$ by using
the formulae in Sec. II.2. Here and henceforth, the dependence of $E_{var}$ on
the primitive orbitals should be understood tacitly. By the variational
principle, $E_{var}$ is always above the exact ground state energy
$E_{gs}^{ext}$. The concern is whether the minimum of $E_{var}$ can be
sufficiently close to $E_{gs}^{ext}$. A case study with $N=L=10$ and $g=4$ is
shown in Fig. 1. We see that for both types of primitive orbitals, at some
value of $\lambda$, $E_{var}$ dips towards the horizontal line indicating
$E_{gs}^{ext}$. At the minima, the relative error is $2\%$ and $4\%$,
respectively, for the Lorentz-type and exponential-type orbital. This is very
encouraging and turns out to be typical.

By determining the minimum of the curve $E_{var}(\lambda;g)$ for a fixed value
of $g$, we can get a variational estimate (denoted as $E_{gs}^{var}$) of the
ground state energy and an optimal permanent variational wave function
$\Phi_{{e}}$. In Fig. 2, the variational energy $E_{gs}^{var}$ is compared
with the exact ground state energy $E_{gs}^{ext}$ obtained by exact
diagonalization (ED) ed . We see that over the full range of $0\leq g<\infty$,
the variational estimates agree with the exact values very well. The
discrepancy is apparent only for the exponential-type variational wave
function in the region around $g=4$, with a relative error about $5.6\%$.
However, in this region, the Lorentz-type wave function is a much better
approximation, reducing the relative error to about $2\%$. The general
observation is that for $g\lesssim 8$, the Lorentz-type wave function yields a
better upper bound for the ground state energy, while for $g\gtrsim 8$, the
exponential-type wave function is better, although the latter fact is hardly
visible in Fig. 2. Here in passing, we mention that if we take the Gross-
Pitaevskii approximation, the orbital occupied by all the particles should be
the zero-momentum Bloch state, as this choice minimizes the kinetic energy and
the interaction energy simultaneously. The energy per particle would be
$-2+g(L-1)/L$, i.e., linear in $g$, which is qualitatively wrong for large
values of $g$.

Figure 4: (Color online) Single-particle correlator $\langle
a_{0}^{\dagger}a_{x}\rangle$. In each panel, the $\ast$ markers are for the
exact ground state $|GS\rangle$, while the circles and squares are for the
energy-minimizing variational state $\Phi_{e}$ with Lorent-type or
exponential-type orbitals, respectively. The common parameters are $N=L=12$.
Note that because of the periodic boundary condition, the largest possible
distance between two sites is $6$.

Figure 5: (Color online) Single-particle correlator $\langle
a_{0}^{\dagger}a_{x}\rangle$. In each panel, the $\ast$ markers are for the
exact ground state $|GS\rangle$, while the circles and squares are for the
overlap-maximizing variational state $\Phi_{o}$ with Lorent-type or
exponential-type orbitals, respectively. The common parameters are $N=L=12$.
Note that because of the periodic boundary condition, the largest possible
distance between two sites is $6$.

Other than energy, a more stringent test for the accuracy of the variational
wave function is its overlap with the exact ground state. We have thus Fig.
3(a), in which the overlap between the exact ground state $|GS\rangle$, which
is obtained by ED, and the energy-minimizing (hence the subscript $e$) VWF
$|\Phi_{e}\rangle$, is shown as a function of $g$. We see that in the full
range, the overlap is at least $0.78$ for the exponential-type VWF for a
system as large as $N=L=12$, and this number is even higher (0.95) for the
Lorentz-type VWF. To appreciate these numbers, one should note that the
dimension of the many-body Hilbert space is as large as $1\,352\,078$. We also
note that Fig. 3(a) and Fig. 2 are consistent with each other. In Fig. 3(a),
all the curves show minima in the proximity of $g=4$. This is exactly where
the variational energies deviate most significantly from the exact one in Fig.
2. In this region, the Lorentz-type VWF has a much higher overlap than the
exponential-type VWF with the exact ground state and accordingly, in this
region the former has a lower energy as shown in Fig. 2. On the other hand,
when $g\gtrsim 8$, the exponential-type VWF becomes better by the overlap
criterion and accordingly, its energy is lower as shown in Fig. 2. The fact
that in the small-$g$ region, the Lorentz-type VWF wins over the exponential-
type VWF while in the large-$g$ region, the exponential-type VWF takes over
might be understandable in view of the superfluid-Mott insulator transition.
The large-$g$ region corresponds to the insulator phase, in which because of
the strong particle-particle repulsion, each particle tends to be localized in
its own site and the tunneling into neighboring sites should be exponentially
small. The small-$g$ region corresponds to the superfluid phase, in which the
particles are more mobile and a more extended orbital should be more
appropriate.

So far, we have been taking the energy minimizing state $\Phi_{e}$ among
either class of VWFs as an approximation of the exact ground state
$|GS\rangle$. Usually, this is the only thing one can do if the exact ground
state is unavailable. However, if the exact ground state is available in a
certain way, say, by exact diagonalization as we do here, a natural alternate
approximation of it should be the variational state having maximal overlap
with it. Let us denote it as $\Phi_{o}$, with the subscript meaning overlap.
Hence, we have two related by different optimization problems. One is energy
minimization and the other overlap maximization. There is no reason that the
two solutions $\Phi_{e}$ and $\Phi_{o}$ should be the same, and by definition,
we have $|\langle GS|\Phi_{e}\rangle|^{2}\leq|\langle GS|\Phi_{o}\rangle|^{2}$
necessarily. In Fig. 3(b), we show $|\langle GS|\Phi_{o}\rangle|^{2}$ as a
function of $g$. In comparison with Fig. 3(a), we see that all the curves
shift upwards as expected. For the exponential-type states, the increase of
the overlap is quite apparent. For instance, while in Fig. 3(a), the minimum
of $|\langle GS|\Phi_{e}\rangle|^{2}$ is about 0.78 for $N=L=12$, in Fig.
3(b), the minimum of $|\langle GS|\Phi_{o}\rangle|^{2}$ is about 0.85. This
strongly indicates that the two optimization problems are related but really
different. For the Lorentz-type states, the increase of the overlap is less
apparent but still visible. For $N=L=12$, the minimum of the overlap increases
from $0.95$ to $0.96$. From the curves in Fig. 3, by extrapolation one can
infer that even for a system as large as $N=L=20$, with a many-body Hilbert
space of dimension about $6.9\times 10^{10}$, the minimal values of the
overlaps $|\langle GS|\Phi_{e}\rangle|^{2}$ and $|\langle
GS|\Phi_{o}\rangle|^{2}$ would be about $0.9$ if we take the Lorentz-type
orbital. These numbers are very impressive.

We also note that in Fig. 3, all the curves, regardless of the primitive
orbital type or the criterion, show minima in the vicinity of $g=4$. This
should be anticipated in view of the superfluid-Mott insulator transition.
According to previous works dmrg ; elesin , the transition occurs at about
$g_{c}=3.61$. Close to the transition, the exact ground state should be most
complex and it is hardest to approximate it with some simple functions.

In hindsight, the large overlap between the variational states and the exact
ground state should be reasonable. There are at least three reasons that are
in favor of such a welcome result. First, both VWFs can reproduce the exact
ground state in either limit of $g=0$ and $g=\infty$. Second, it is well-known
that for an arbitrary value of $g$, by the Perron-Frobenius theorem horn , the
ground state is strictly positive everywhere in the Fock-state basis. This
property is shared by both VWFs by construction. Third, it is also known that
the ground state belongs to the trivial representation of the symmetry group
of the model (the dihedral group), or more specifically, it is invariant under
all translations and reflections, a property again shared by both VWFs by
construction.

Finally, as yet another check of the quality of the VWFs, we consider the
single-particle correlator $\langle a_{0}^{\dagger}a_{x}\rangle$. This
expression is convenient for the exact ground state, which is obtained by ED
in the Fock-state basis. For the VWFs, which are in the first-quantization
formalism, we note that the corresponding single-particle operator $\hat{C}$
has matrix elements $\langle m|\hat{C}|n\rangle=\delta_{m,0}\delta_{n,x}$, and
we can use (32) and (33) to calculate its expectation value

$\displaystyle\langle
a_{0}^{\dagger}a_{x}\rangle=\sum_{i,j=1}^{N}\phi_{i}^{*}(0)\phi_{j}(x)\frac{\operatorname{per}(A;i|j)}{\operatorname{per}(A)}.$
(42)

In Fig. 4, $\langle a_{0}^{\dagger}a_{x}\rangle$ is plotted against $x$ for
the energy-minimizing state $\Phi_{e}$ and the exact ground state
$|GS\rangle$. We see a picture consistent with Fig. 2 and Fig. 3. In Fig.
4(a)-(c), when $g$ is small and the lorent-type VWF is better, the correlator
predicted by the lorent-type VWF is very close to the exact one. In Fig.
4(e)-(f), when $g$ gets large, the exponential-type VWF is better, and
accordingly the correlator predicted by the exponential-type VWF is close to
the exact one. For any value of $g$, either the lorent-type or the
exponential-type VWF will be a good approximation by all the three criterions.
Here it is also interesting to note that while the exponential-type VWF always
underestimates the correlator, the lorent-type VWF slightly underestimates it
in the small-$g$ region, while overestimates it in the large-$g$ region. This
might be related to the superfluid-Mott insulator transition.

In Fig. 5, $\langle a_{0}^{\dagger}a_{x}\rangle$ is plotted against $x$ for
the overlap-maximizing states $\Phi_{o}$. It is normal to expect that
$\Phi_{o}$ reproduces the correlation function better than $\Phi_{e}$. This is
indeed the case. We see that for $g\leq 6$, when $|\langle
GS|\Phi_{o}\rangle|^{2}$ is significantly higher than $|\langle
GS|\Phi_{e}\rangle|^{2}$, the curves of $\langle a_{0}^{\dagger}a_{x}\rangle$
shift closer to the exact curve for both types of orbitals. In particular, the
curves with the Lorentz-type orbitals almost coincide with the exact curves.
For $g\geq 8$, $|\langle GS|\Phi_{e}\rangle|^{2}\simeq|\langle
GS|\Phi_{o}\rangle|^{2}$, indicating that $\Phi_{e}\simeq\Phi_{o}$, we do not
see any significant change of the curves.

##  IV Optimization algorithm

In the proceeding section, we have seen that for the one-dimensional Bose-
Hubbard model with the periodic boundary condition and at unit filling, it is
possible to construct some permanent state out of some simple orbitals to
approximate its exact ground state very well. This is checked by examining the
variational energy, the overlap with the exact ground state, and the single-
particle correlation function.

The close approximation is achieved with some _preassigned_ orbitals depending
on a single parameter $\lambda$. A natural question is whether the numbers can
be further improved by allowing more freedom of the orbitals. This leads to
two optimization problems. First, for a given bosonic system with Hamiltonian
(28), how can we find a set of orbitals $\\{\phi_{i},1\leq i\leq N\\}$, such
that the energy expectation value of the permanent state $\Phi$ constructed in
(2),

$\displaystyle E$ $\displaystyle=$
$\displaystyle\frac{\langle\Phi|H|\Phi\rangle}{\langle\Phi|\Phi\rangle},$ (43)

is minimized? Second, for a given normalized bosonic wave function
$|\Psi\rangle$, how can we find the permanent state which is the optimal
approximation of it? That is, how can we find the permanent state as in (2)
such that the overlap zhang1 ; zhang2

$\displaystyle O$ $\displaystyle=$
$\displaystyle\frac{\langle\Phi|\Psi\rangle\langle\Psi|\Phi\rangle}{\langle\Phi|\Phi\rangle}$
(44)

is maximized?

There exists a common simple strategy for both problems zhang1 ; zhang2 . For
clarity, let us focus on the first problem for the present. Let us fix $N-1$
orbitals, say, the orbitals $\phi_{2\leq i\leq N}$, and try to find an optimal
$\phi_{1}$. To this end, we note that with the orbitals $\phi_{2\leq i\leq N}$
fixed, the numerator and denominator in (43) are both hermitian forms of
$\phi_{1}$. That is, one can find operators $\hat{F}$ and $\hat{G}$ such that
$\langle\phi_{1}|\hat{F}|\phi_{1}\rangle=\langle\Phi|H|\Phi\rangle$ and
$\langle\phi_{1}|\hat{G}|\phi_{1}\rangle=\langle\Phi|\Phi\rangle$. We then can
rewrite the ratio as

$\displaystyle E(\phi_{1})$ $\displaystyle=$
$\displaystyle\frac{\langle\phi_{1}|\hat{F}|\phi_{1}\rangle}{\langle\phi_{1}|\hat{G}|\phi_{1}\rangle},$
(45)

By definition, $\hat{F}$ and $\hat{G}$ are hermitian single-particle operators
depending on the orbitals $\phi_{2\leq i\leq N}$. Importantly, $\hat{G}$ is
even positive definite as long as $\phi_{2\leq i\leq N}$ are all nonzero, as
by definition
$\langle\phi_{1}|\hat{G}|\phi_{1}\rangle=\langle\Phi|\Phi\rangle\geq 0$, with
the equality achieved only if $\phi_{1}=0$. Here we recall Proposition 1,
which asserts that $\Phi$ is necessarily nonzero if $\phi_{1\leq i\leq N}$ are
nonzero. It is a straightforward but lengthy calculation to derive the
explicit expressions of $\hat{F}$ and $\hat{G}$, so we defer it to the
Appendix. Suppose we have prepared the operators $\hat{F}$ and $\hat{G}$ (this
is the most time-consuming part of the iteration). The optimal $\phi_{1}$ that
will minimize the ratio in (45) is just the solution of the following
generalized eigenvalue problem

$\displaystyle\hat{F}\phi$ $\displaystyle=$
$\displaystyle\varepsilon\hat{G}\phi$ (46)

corresponding to the smallest eigenvalue $\varepsilon_{min}$, and the minimum
of the ratio is just $\varepsilon_{min}$. Note that we just need the smallest
eigenvalue. Therefore, we can resort to the Lanczos algorithm.

Once we have updated $\phi_{1}$, we can turn to $\phi_{2}$, and then to
$\phi_{3}$, and so on. However, for convenience of programming, we can just
make a circular shift of the orbitals $\phi_{i}\rightarrow\phi_{i-1}$, and
continue to update $\phi_{1}$. In this process, the variational energy
decreases monotonically, and as it is lower bounded by the exact ground state
energy, it will definitely converge.

We have to mention that (46) is essentially the self-consistency equation
derived by Heimsoth before by the method of performing variational
differentiation of abstract Hilbert space vectors martin1 ; martin2 . However,
hopefully here our different point of view has led to a more compact and
transparent formalism. As we shall see in the subsection below, this formalism
allows easy extension to the multiconfiguration case, which was not considered
previously.

Now it should be clear that the second problem can be treated similarly. The
ratio (44) can be written as

$\displaystyle O$ $\displaystyle=$
$\displaystyle\frac{\langle\phi_{1}|\gamma\rangle\langle\gamma|\phi_{1}\rangle}{\langle\phi_{1}|\hat{G}|\phi_{1}\rangle},$
(47)

where the single-particle orbital $\gamma$ is defined by the summation or
partial contraction

$\displaystyle\gamma(x)$ $\displaystyle=$
$\displaystyle\sqrt{N!}\sum_{x_{2},\ldots,x_{N}=1}^{L}\Psi(x,x_{2},\ldots,x_{N})\prod_{i=2}^{N}\phi_{i}^{*}(x_{i}).\quad$
(48)

Once $\gamma$ is calculated (again, this is the most time-consuming part), the
optimal $\phi_{1}$ can be obtained by solving a generalized eigenvalue
equation similar to (46), with $|\gamma\rangle\langle\gamma|$ replacing
$\hat{F}$.

Naively, the summation in (48) has the complexity of $L^{N-1}N$, as $x_{2\leq
i\leq N}$ run independently from $1$ to $L$. However, one should note that
$\Psi$ is invariant under permutations of $x_{2\leq i\leq N}$. Making use of
this fact and changing the dummy variables from $x_{i}$ to $y_{i-1}$, (48) can
be written as

$\displaystyle\gamma(x)$ $\displaystyle=$
$\displaystyle\sqrt{N!}\sum_{\textbf{y}}\Psi(x,\textbf{y})\frac{\operatorname{per}(\phi^{*}(\textbf{y}))}{\textbf{n}(\textbf{y})!},$
(49)

where the summation is over the ordered $(N-1)$-tuple
$\textbf{y}\equiv(y_{1},y_{2},\ldots,y_{N-1})$ with $1\leq y_{1}\leq
y_{2}\ldots\leq y_{N-1}\leq L$, and $\phi^{*}(\textbf{y})$ denotes the
$(N-1)\times(N-1)$ matrix with its $i$th row being the $y_{i}$th row of the
$L\times(N-1)$ matrix $(\phi_{2}^{*},\phi_{3}^{*},\ldots,\phi_{N}^{*})$. In
the denominator, $\textbf{n}(\textbf{y})!\equiv\prod_{i=1}^{L}n_{i}!$, with
$n_{i}$ denoting the times $i$ appears in y. The computational complexity is
now on the order of ${L+N-2\choose N-1}2^{N-2}N$.

###  IV.1 The multiconfiguration case

So far, we have assumed a single configuration. This is yet the only case
people have considered igor1 ; igor2 ; martin1 ; martin2 . For better
approximation, one can try $M>1$ sets of orbitals
$\\{\phi^{(\alpha)}_{i},1\leq\alpha\leq M,1\leq i\leq N\\}$, and let the
variational wave function be the sum of the permanent wave functions
constructed by each set of orbitals. Specifically,

$\displaystyle\Phi$ $\displaystyle=$
$\displaystyle\sum_{\alpha=1}^{M}\Phi^{(\alpha)}=\sum_{\alpha=1}^{M}\hat{\mathcal{S}}(\phi^{(\alpha)}_{1},\ldots,\phi^{(\alpha)}_{N}).$
(50)

Fixing the orbitals $\phi_{2\leq i\leq N}^{(\alpha)}$ in each set, the
variational energy (43) can be written as

$\displaystyle E$ $\displaystyle=$
$\displaystyle\frac{\sum_{\alpha,\beta=1}^{M}\langle\Phi^{(\alpha)}|H|\Phi^{(\beta)}\rangle}{\sum_{\alpha,\beta=1}^{M}\langle\Phi^{(\alpha)}|\Phi^{(\beta)}\rangle}$
(51) $\displaystyle=$
$\displaystyle\frac{\sum_{\alpha,\beta=1}^{M}\langle\phi^{(\alpha)}_{1}|\hat{F}^{(\alpha\beta)}|\phi_{1}^{(\beta)}\rangle}{\sum_{\alpha,\beta=1}^{M}\langle\phi^{(\alpha)}_{1}|\hat{G}^{(\alpha\beta)}|\phi_{1}^{(\beta)}\rangle},$

where $\hat{F}^{(\alpha\beta)}$ and $\hat{G}^{(\alpha\beta)}$ depend on the
fixed orbitals $\phi_{2\leq i\leq N}^{(\alpha)}$ and $\phi_{2\leq i\leq
N}^{(\beta)}$, and can be calculated with essentially the same formulae as in
the Appendix—Just add the superscript $\alpha$ to the orbitals in the bras and
$\beta$ to the orbitals in the kets. It is easily seen that

$\displaystyle(\hat{F}^{(\alpha\beta)})^{\dagger}=\hat{F}^{(\beta\alpha)},\quad(\hat{G}^{(\alpha\beta)})^{\dagger}=\hat{G}^{(\beta\alpha)}.$
(52)

Apparently, (51) can be cast in the same form as (45), if we identify
$\phi_{1}$ as the concatenated vector

$\displaystyle\phi_{1}$ $\displaystyle\equiv$
$\displaystyle(\phi_{1}^{(1)};\phi_{1}^{(2)};\ldots;\phi_{1}^{(M)}),$ (53)

which is of length $LM$, and define the block operators
$\hat{F}\equiv(\hat{F}^{(\alpha\beta)})$ and
$\hat{G}\equiv(\hat{G}^{(\alpha\beta)})$, which are of size $LM\times LM$. The
same update and iteration procedures can then be carried out.

Similarly, in the multiconfiguration case, the overlap (44) can be written as

$\displaystyle O$ $\displaystyle=$
$\displaystyle\frac{\sum_{\alpha,\beta=1}^{M}\langle\phi_{1}^{(\alpha)}|\gamma^{(\alpha)}\rangle\langle\gamma^{(\beta)}|\phi_{1}^{(\beta)}\rangle}{\sum_{\alpha,\beta=1}^{M}\langle\phi^{(\alpha)}_{1}|\hat{G}^{(\alpha\beta)}|\phi_{1}^{(\beta)}\rangle},$
(54)

where the single-particle orbital $\gamma^{(\alpha)}$, like $\gamma$ in (48),
is defined by the same summation with $\phi_{i}^{(\alpha)}$ replacing
$\phi_{i}$. Again, the same strategy as above can be applied.

###  IV.2 A pitfall with $N=2$

A pitfall is to be avoided in implementing the multi-configuration scheme. In
the single-configuration case, the hermitian operator $\hat{G}$ is strictly
positive definite as long as the orbitals $\phi_{2\leq i\leq N}$ are nonzero.
This is because by definition
$\langle\phi_{1}|\hat{G}|\phi_{1}\rangle=\langle\Phi|\Phi\rangle$, and by
Proposition 1, $\Phi$ is nonzero if $\phi_{1\leq i\leq N}$ are nonzero. In
contrast, in the multi-configuration case (50), we do not necessarily have
$\Phi\neq 0$ even if all the orbitals $\\{\phi^{(\alpha)}_{i},1\leq\alpha\leq
M,1\leq i\leq N\\}$ are nonzero—the configurations could cancel each other
out. This happens particularly in the two-particle case of $N=2$. For
illustration, let us consider the two-particle, two-configuration case. For
any value of $\\{\phi^{(1)}_{2},\phi^{(2)}_{2}\\}$, if we choose
$\\{\phi^{(1)}_{1},\phi^{(2)}_{1}\\}$ as
$\\{\phi^{(2)}_{2},-\phi^{(1)}_{2}\\}$,

$\displaystyle\Phi$ $\displaystyle=$
$\displaystyle\hat{\mathcal{S}}(\phi^{(1)}_{1},\phi^{(1)}_{2})+\hat{\mathcal{S}}(\phi^{(2)}_{1},\phi^{(2)}_{2})$
(55) $\displaystyle=$
$\displaystyle\hat{\mathcal{S}}(\phi^{(2)}_{2},\phi^{(1)}_{2})-\hat{\mathcal{S}}(\phi^{(1)}_{2},\phi^{(2)}_{2})=0.$

That the total wave function $\Phi$ could vanish means that $\hat{G}$ in the
multi-configuration case is just positive semi-definite—It could have zero
eigenvalues. Theoretically, this does not cause any problem because when the
denominator of (51) vanishes, its numerator vanishes too. In other words, the
eigenvectors of $\hat{G}$ with the zero eigenvalue are also eigenvectors of
$\hat{F}$ with the zero eigenvalue. In numerics, one has to restrict $\hat{F}$
and $\hat{G}$ to the subspace spanned by the eigenvectors of $\hat{G}$ with
nonzero (hence positive) eigenvalues. This can be easily implemented once
$\hat{G}$ is diagonalized (but a better approach is described below).

In practice, this cautious extra effort is necessary only for $N=2$. In our
extensive numerical simulations, we have never encountered a case of $\hat{G}$
becoming singular for $N\geq 3$. The reason is yet to be understood but it is
not surprising as when $N\geq 3$, for fixed $\phi_{2\leq i\leq N}^{(\alpha)}$,
unlike (55), it is hard to find
$\phi_{1}=(\phi_{1}^{(1)};\phi_{1}^{(2)};\ldots;\phi_{1}^{(M)})$ to make
$\Phi$ vanish; likely there is no solution. But for $N=2$, it occurs
necessarily. With $M$ configurations, generally $\hat{G}$ has $\binom{M}{2}$
eigenvectors with the eigenvalue zero. This number can be understood in view
of (55). For fixed orbitals $\phi_{2}^{(\alpha)}$, one can find $\binom{M}{2}$
linearly independent vectors $\phi_{1}$ such that $\Phi$ vanishes.

To handle the singularity of $\hat{G}$ for $N=2$ and $M>1$, we consider
instead of the ratio in (45), the modified ratio

$\displaystyle\tilde{E}(\phi_{1})$ $\displaystyle=$
$\displaystyle\frac{\langle\phi_{1}|\hat{F}+\mu_{1}\hat{P}|\phi_{1}\rangle}{\langle\phi_{1}|\hat{G}+\mu_{2}\hat{P}|\phi_{1}\rangle},$
(56)

where $\hat{P}$ is the projection operator onto $\ker{\hat{G}}$, the kernel of
$\hat{G}$, and $\mu_{1,2}$ are two real parameters which can be chosen
arbitrarily except that $\mu_{2}$ must be positive so that
$\hat{G}+\mu_{2}\hat{P}$ is positive definite. By construction, $\tilde{E}$
agrees with $E$ on the subspace orthogonal to $\ker{\hat{G}}$, and takes on
the value $\mu_{1}/\mu_{2}$ on $\ker{\hat{G}}$. Therefore, the modified
generalized eigenvalue problem

$\displaystyle(\hat{F}+\mu_{1}\hat{P})\phi$ $\displaystyle=$
$\displaystyle\varepsilon(\hat{G}+\mu_{2}\hat{P})\phi$ (57)

yields essentially the same generalized eigen-pairs as the initial one (46).
In particular, the minimal generalized eigenvalue (and its associated
eigenvector) is not altered as long as $\mu_{1}/\mu_{2}$ is not smaller. In
this paper, we find it okay to choose $\mu_{1}=0$ and $\mu_{2}=1$.

###  IV.3 Convergence of the algorithm

Figure 6: (Color online) Convergence of the energy of the permanent
variational state. Each curve corresponds to a different set of initial
orbitals. The solid (dotted) lines are with $M=1$ ($M=2$) configuration(s). In
each step, one orbital is updated. The energy after $n$ steps is denoted as
$E_{n}$. The limiting value $E_{\infty}$ is approximated by $E_{2400}$. The
setting is a one-dimensional Bose-Hubbard model in a harmonic trap with the
Hamiltonian (58). The parameters $(N,L,\kappa,g)=(3,13,0.2,5)$.

We take a concrete model to illustrate the convergence behavior of the
algorithm. Consider a one-dimensional Bose-Hubbard model in a harmonic trap.
The Hamiltonian is

$\displaystyle H_{trap}$ $\displaystyle=$
$\displaystyle-\sum_{x=-L_{0}}^{L_{0}-1}(a_{x}^{\dagger}a_{x+1}+\text{h.c.})+\frac{g}{2}\sum_{x=-L_{0}}^{L_{0}}a_{x}^{\dagger}a_{x}^{\dagger}a_{x}a_{x}$
(58) $\displaystyle+\kappa\sum_{x=-L_{0}}^{L_{0}}x^{2}a_{x}^{\dagger}a_{x}.$

Here $\kappa>0$ is the stiffness of the harmonic potential. For symmetry, we
have assume a lattice of size $L=2L_{0}+1$. We take the open boundary
condition.

We start from $MN$ random orbitals [see Eq. (50)], where $M$ is the
configuration number and $N$ is the particle number. Essentially, we just
generate an $ML\times N$ matrix with each element chosen randomly from the
interval $[0,1]$ according to the uniform distribution. The orbitals (the
columns) are then updated in a circular way. The variational state after $n$
updates (steps) will be denoted as $\Phi_{n}$, and its energy will be denoted
as $E_{n}$. Note that $E_{n}$ is obtained simultaneously in solving the
generalized eigenvalue problem (46).

By construction, $E_{n}$ decreases monotonically and will definitely converge.
The concern is in which way and with what rate it converges to its limit
$E_{\infty}$. This is studied in Fig. 6, where for a set of values of the
parameters $(N,L,\kappa,g)$ and with $M=1$ or $2$, some typical trajectories
of $E_{n}-E_{\infty}$ (here we just approximate $E_{\infty}$ by some $E_{n}$
with a large enough $n$) are displayed. For each value of $M$, we have three
trajectories corresponding to three different sets of initial orbitals.

We see that often the trajectory is not very regular—It is neither pure
exponential nor pure power law, but clearly divides into different parts. In
many cases, after some relaxation or transition stage, which can last for a
long time, the energy eventually enters an exponentially decreasing mode. In
our extensive numerical simulations, the observation is that the trajectory of
the energy depends not only on the model and the model parameters, but also on
the initial conditions and can differ significantly from run to run. We cannot
draw any definite rule for the convergence rate of the energy, but the feeling
is that the convergence tends to be slower with more configurations and weaker
interactions.

####  IV.3.1 The non-interacting case

Figure 7: (Color online) (a) and (b): Trajectories of the energy error
$E_{n}-E_{\infty}$ and the infidelity $1-|\langle GS|\Phi_{n}\rangle|^{2}$ in
a particular run in the non-interacting case ($g=0$). (c) and (d): The
orbitals $\phi_{i}(x)$, the particle density distribution $\rho(x,x)$, and the
one-particle correlator $\rho(0,x)$ at two snapshots. In (c1) and (d1), the
solid dots represent the single-particle ground state. The setting is the same
as in Fig. 6. The parameters are also the same except for $g$.

Possibly the most embarrassing and surprising thing is that the convergence is
slowest in the non-interacting limit. In this trivial case, the exact ground
state $|GS\rangle$ is simply a condensate-type (and hence a permanent) state
with all particles occupying the single-particle ground state. We do not need
to invoke the algorithm for its calculation, however, if we do, the
convergence is as slow as a power law. In Figs. 7(a) and 7(b), with $g=0$ but
the other parameters the same as in Fig. 6, trajectories of the energy error
$E_{n}-E_{\infty}$ and the infidelity $1-|\langle GS|\Phi_{n}\rangle|^{2}$ are
shown respectively for a particular run. Here we take $E_{\infty}$ to be the
exact ground state energy. In either figure, the curve drops down steeply at
about $n\sim N$, and afterwards it follows a straight line in the log-log
plot. Basically, the picture is that after the first round of update, i.e.,
when each orbital has been updated once, the variational wave function is
already very close to the exact state (the overlap is over $0.99$ in the
particular case). Afterwards, it improves slowly by a power law.

The exact reason behind the slow convergence is yet to be understood. Here we
just emphasize that the convergence is slow only in the asymptotic sense. In
the initial phase, the algorithm can already deliver the state $\Phi_{n}$
close enough to its limit $\Phi_{\infty}=|GS\rangle$.

In Figs. 7(c) and 7(d), we show two snapshots of the constituent orbitals, the
particle density distribution $\rho(x,x)$, and the correlation function
$\rho(0,x)$, which are respectively the diagonal and off-diagonal parts of the
single-particle reduced density matrix $\rho$ defined as

$\displaystyle\rho(x_{1},x_{2})$ $\displaystyle=$ $\displaystyle\langle
a_{x_{2}}^{\dagger}a_{x_{1}}\rangle.$ (59)

We see that even for an $n$ as small as $n=N=3$, the permanent variational
state $\Phi_{n}$ can already reproduce the exact values of $\rho(x,x)$ and
$\rho(0,x)$ to high precision. The intriguing thing is that while the total
wave function is already very close to the exact many-body ground state in
terms of energy, overlap, and some most relevant correlation functions, the
constituent orbitals are still far away or at least visibly different from the
exact single-particle orbital. In view of Proposition 2, which states that
different sets of orbitals necessarily result in different many-body wave
functions, the current observation implies that the latter is not necessarily
very sensitive to perturbations of the former in some cases. This in turn
implies that it might not be a good idea to use convergence of the orbitals as
a criterion in determining the termination of the algorithm martin2 .

####  IV.3.2 Local minima

Figure 8: (Color online) (a) Twenty trajectories of the variational energy
$E_{n}$. Each trajectory starts with a different set of randomly generated
initial orbitals. (b) Possible limiting values of the variational energy
$E_{\infty}$ and the overlap $|\langle GS|\Phi_{\infty}\rangle|^{2}$. Note
that for each value of $g$, we have 20 runs as in (a) and in each run the
number of iteration steps is 1500. The setting is a one-dimensional Bose-
Hubbard model in a harmonic trap with the Hamiltonian (58). The fixed
parameters are $(N,L,\kappa,M)=(3,13,0.2,1)$.

It should be no wonder that the algorithm can get stuck in a local minimum
like many other greedy algorithms. This is illustrated in Fig. 8(a) with a
model of (58) and some specific value of $g$. We see that many trajectories of
$E_{n}$ settle down on a secondary minimum. Note that for clarity, we have
shown only the first 100 steps, but actually the horizontal lines extend all
the way up to $n=1500$.

In Fig. 8(b), we plot the possible eventual values of the variational energy
$E_{n}$ and the overlap $|\langle GS|\Phi_{n}\rangle|^{2}$ as functions of
$g$. For each value of $g$, like in Fig. 8(a), we have run the algorithm 20
times, each time up to $n=1500$. We see that for $g$ smaller than some
critical value $g_{c}\simeq 8$, we get only a single value for either of
$E_{\infty}$ and $|\langle GS|\Phi_{\infty}\rangle|^{2}$, which indicates that
there is only a global minimum, however for $g$ larger than $g_{c}$, we get
two different values for either of $E_{\infty}$ and $|\langle
GS|\Phi_{\infty}\rangle|^{2}$, which indicates the presence of a second, local
minimum.

Our experience is that local minima are ubiquitous. Generally, their number
increases with the number $M$ of configurations. To avoid them and enhance the
probability of hitting the global minimum, we simply run the algorithm
multiple times, say $20$ times for $M=4$ configurations.

####  IV.3.3 Real versus complex

Figure 9: (Color online) Trajectories of the energy of the permanent
variational state with either (a) real or (b) complex orbitals. In each case,
we have five runs starting with five different sets of random orbitals. The
number of configurations is $M=1$. The horizontal dotted lines indicate the
exact ground state energy. The setting is a two-particle, two-site Bose-
Hubbard model with the Hamiltonian (60). The parameter $g=-1$.

So far, we have assumed the optimal orbitals to be real. For many systems with
the time reversal symmetry, the many-body ground state is real and the
assumption that the optimal orbitals should also be real seems very
reasonable. However, this is not the case.

We take a minimal model to illustrate the possibility that the optimal
orbitals could be complex although the total wave function is real. Consider a
two-particle, two-site Bose-Hubbard model, i.e., a model with $(N,L)=(2,2)$.
The Hamiltonian is

$\displaystyle
H_{ds}=-(a_{1}^{\dagger}a_{2}+a_{2}^{\dagger}a_{1})+\frac{g}{2}(a_{1}^{\dagger}a_{1}^{\dagger}a_{1}a_{1}+a_{2}^{\dagger}a_{2}^{\dagger}a_{2}a_{2}).\quad$
(60)

Consider a real wave function $\Psi(x_{1},x_{2})$, with $x_{1,2}=1,2$, of this
model. By Proposition 3 or Proposition 4, $\Psi$ can be written as a permanent
state

$\displaystyle\Psi$ $\displaystyle=$ $\displaystyle\hat{S}(u,v)$ (61)

with two orbitals $u$ and $v$. The question is whether $u$, $v$ can both be
real. Componentwise, (61) means

$\displaystyle\Psi(1,1)$ $\displaystyle=$
$\displaystyle\sqrt{2}u_{1}v_{1},\quad\Psi(2,2)=\sqrt{2}u_{2}v_{2},$
$\displaystyle\Psi(1,2)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}(u_{1}v_{2}+v_{1}u_{2}).$

We have then

$\displaystyle\Delta\equiv\Psi(1,2)^{2}-\Psi(1,1)\Psi(2,2)$ $\displaystyle=$
$\displaystyle\frac{1}{2}(u_{1}v_{2}-v_{1}u_{2})^{2}.$

We thus see that for (61) to have real solutions, a necessary condition is
that the quantity $\Delta$ be non-negative. It is easy to verify that this is
also sufficient.

Hence, when $\Delta<0$, in the single configuration approximation, the wave
function $\Psi$ can be exactly recovered with complex orbitals but not with
real orbitals. Note that such a condition is satisfied by a cat-type state,
which can be realized as the ground state of the model (60) with an attractive
on-site interaction $g<0$. In Fig. 9, we show the variational energies of the
ground state of such a model, calculated with either real or complex orbitals.
We see that while complex orbitals can deliver the exact value, real orbitals
miss it with some overestimation.

Therefore, we see that at least theoretically, complex orbitals are superior
to real orbitals for energy minimization. However, the observation is that for
all the models we consider in this paper, if the interaction is repulsive
($g\geq 0$) and if we take the single configuration approximation ($M=1$), the
complex approach leads to identical results with the real approach. More
specifically, in these circumstances, even if we start with complex orbitals,
the iteration will result in real orbitals. The reason is yet to be
understood. For multiple configurations ($M>1$), with complex orbitals, often
we do get lower energies, however, the improvement is not very significant. We
thus often confine ourselves to real orbitals in the following. Anyway, real
arithmetics are four times fast than complex arithmetics, and the disadvantage
in accuracy can be compensated by including more configurations.

Figure 10: (Color online) (a) Energy per particle of the energy-minimizing
permanent states $\Phi_{e}$ constructed with either the Lorentz-type (40a) or
the exponential-type (40b) primitive orbitals, or unrestricted orbitals. The
solid line indicates the exact ground state energy calculated by exact
diagonalization (ED). (b) Overlap between the energy-minimizing permanent
states $|\Phi_{e}\rangle$ and the exact ground state $|GS\rangle$. In this
figure, the parameters are $N=L=12$.

Figure 11: (Color online) Single-particle correlator $\langle
a_{0}^{\dagger}a_{x}\rangle$. In each panel, the $\ast$ markers are for the
exact ground state $|GS\rangle$, while the circles are for the unrestricted
energy-minimizing VWF $\Phi_{e}$. The common parameters are $N=L=12$. Note
that because of the periodic boundary condition, the largest possible distance
between two sites is $6$.

##  V Application of the algorithm

###  V.1 The Bose-Hubbard model revisited

With the numerical optimization algorithm above, we can handle more general
models. But let us start from the Bose-Hubbard model with the periodic
boundary condition and at unit filling, and see how much it can improve over
the results in Sec. III.

For this particular model, we take a single configuration ($M=1$). We never
encounter any local minimum, and the symmetry of the model is perfectly
preserved by the optimal orbitals. That is, although we always start from
random orbitals, the orbitals we eventually get are always of the same shape
and differ from each other just by translations, as described by (38).

In Fig. 10, we show the estimated ground state energy $E_{gs}$ and the overlap
$|\langle GS|\Phi_{e}\rangle|^{2}$ obtained with unrestricted orbitals. For
comparison, also shown are the results with Lorentz or exponential orbitals.
We see that in the region $g\leq 6$, the VWF with unrestricted orbitals does
not improve much over the VWF with Lorentz orbitals neither by energy nor by
overlap. Accordingly, the predicted correlation function $\langle
a_{0}^{\dagger}a_{x}\rangle$ is close to that predicted by the Lorentz VWF, as
can be seen by comparing Figs. 11(a)-(c) with Figs. 4(a)-(c). However, in the
region $g\geq 6$, unrestricted orbitals do lead to improvement over both the
Lorentz and the exponential orbitals, both by the criterions of energy and
overlap. For instance, at $g=8$ and with $N=L=12$, the overestimate in energy
(difference between the variational and exact ground state energy) reduces
from $0.0127$ (Lorentz) and $0.0106$ (exponential) to $0.0039$ (unrestricted),
and simultaneously the deficiency in overlap ($1-|\langle
GS|\Phi_{e}\rangle|^{2}$) reduces from 0.0213 (Lorentz) and 0.0203
(exponential) to 0.0026 (unrestricted). These numbers indicate that the
algorithm produces really good approximation of the exact ground state.
Indeed, as Figs. 11(d)-(f) show, now the VWF-predicted correlation function
$\langle a_{0}^{\dagger}a_{x}\rangle$ almost coincides with the exact values
for $g\geq 8$.

Overall, in Figs. 10 and 11, we see that for a system as large as $N=L=12$ and
in the whole range of $g$, the ground state can be very well approximated by a
permanent state. This is very impressive in view of the dimension of the many-
body Hilbert space, which is as large as $\mathcal{D}=1\,352\,078$. The exact
ground state is obtained by exact diagonalization and is a vector of this size
with the Fock states as a basis. In contrast, the permanent variational state
is constructed with $N=12$ orbitals, each of which is a vector of size $L=12$.
That the orbitals are not even independent but related to each other by
translations means that the permanent state is encoded with a $12\times 12$
circulant matrix with only 12 independent variables. From the data compression
point of view, with the permanent state as an approximation of the exact
ground state, the compression ratio is very high while the fidelity is still
very good.

###  V.2 More general models

Figure 12: (Color online) (a) Ground state energy (per particle) of a Bose-
Hubbard model at unit filling with the open boundary condition [see (62) for
the Hamiltonian]. The solid line is obtained by exact diagonalization, while
the dotted line by the permanent variational approach. (b) Overlap between the
exact ground state $|GS\rangle$ and the energy-minimizing permanent state
$\Phi_{e}$. The parameters are $N=L=11$.

Figure 13: (Color online) The orbitals $\phi_{i}(x)$, the particle density
distribution $\rho(x,x)$, and the one-particle correlator $\rho(0,x)$ for
three different values of $g$. The circles are for the permanent state, while
the $\ast$ markers are for the exact diagonalization results. The setting is a
one-dimensional Bose-Hubbard model with the open boundary condition as defined
in (62). The fixed parameters are $(N,L)=(11,11)$.

Figure 14: (Color online) (a1) and (b1): Ground state energy per particle of a
Bose-Hubbard model in a harmonic trap. The solid lines are obtained by exact
diagonalization, while the other lines are by the permanent variational wave
functions with different numbers of configurations. (a2) and (b2): Overlap
between the exact ground state $|GS\rangle$ and the energy-minimizing
variational state $\Phi_{e}$. The fixed parameters are
$(N,L,\kappa)=(3,13,0.2)$. The left (right) column is calculated with real
(complex) orbitals.

Figure 15: (Color online) The orbitals $\phi_{i}(x)$, the particle density
distribution $\rho(x,x)$, and the one-particle correlator $\rho(0,x)$ for
three different values of $g$. The circles are for the permanent state, while
the $\ast$ markers are for the exact diagonalization results. The setting is a
one-dimensional Bose-Hubbard model in a harmonic trap as defined in (58). The
fixed parameters are $(N,L,\kappa)=(5,13,0.25)$ as in Fig. 6. In the top
panels, the dotted line sketches the harmonic potential.

Figure 16: (Color online) Same as Fig. 15, but with $M=2$ configurations. Note
that the variationally calculated density distribution $\rho(x,x)$ and the
correlation function $\rho(0,x)$ are asymmetric.

Figure 17: (Color online) Same as Fig. 15 and Fig. 16, but with $M=3$
configurations.

We now turn to more general models.

First of all, let us break the periodic boundary condition of the Bose-Hubbard
model above and replace it with the open boundary condition. The Hamiltonian
is

$\displaystyle
H_{obc}=-\sum_{x=-L_{0}}^{L_{0}-1}(a_{x}^{\dagger}a_{x+1}+\text{h.c.})+\frac{g}{2}\sum_{x=-L_{0}}^{L_{0}}a_{x}^{\dagger}a_{x}^{\dagger}a_{x}a_{x}.\;\;\;\;\;$
(62)

We still assume unit filling, so the particle number $N=L=2L_{0}+1$. By the
mere change of the boundary condition, the translation symmetry is lost and
now the orbitals should differ in shape. It is then unclear what orbitals to
choose to construct the permanent state—We have to resort to the numerical
algorithm.

In Fig. 12(a), with $N=L=11$, the permanent estimated ground state energy is
compared with the exact diagonalization result, and in Fig. 12(b), the overlap
between the permanent variational state $\Phi_{e}$ and the exact ground state
$|GS\rangle$ is shown. We see that like the periodic boundary condition case,
across the full range of the on-site interaction $g$, the permanent state is a
very good approximation of the exact ground state. The overlap is as large as
$0.964$ even in the worst case, and the relative error in energy is at most
$2\%$.

In Fig. 13, for three different values of $g$, we show the constituent
orbitals, the density distribution function $\rho(x,x)$ and the correlation
function $\rho(0,x)$. For each value of $g$, we start from $N$ random
orbitals, and then update each orbital 100 times. In the top panels, we see
that besides the bulk orbitals which are similar to each other in shape, there
are two edge orbitals, which are maximal on the edges and decay into the bulk.
As the permanent state is very close to the exact ground state by overlap,
these orbitals provide a very good picture of the exact ground state. We also
see that the permanent state predicted values of the density distribution and
the correlation function agree with the exact results very well.

As a second example, let us consider the one-dimensional Bose-Hubbard model in
a harmonic trap, with the Hamiltonian of (58). In Fig. 14, we show the
variational ground state energy $E_{gs}$ and the overlap $|\langle
GS|\Phi_{e}\rangle|^{2}$, calculated with various numbers of configurations
and with either real or complex orbitals, as functions of $g$. We see that in
the single configuration case ($M=1$), the real and complex approaches agree
with each other exactly. The discrepancy between the variational energy and
the exact value is apparent, and the overlap drops to 0.76 at $g=15$. However,
the situation improves dramatically if we take $M=3$ configurations. With
three configurations, both the real and the complex approaches get the ground
state energy so accurate that the difference with the exact value is hardly
visible in the figures. Accordingly, the overlap $|\langle
GS|\Phi_{e}\rangle|^{2}$ is very close to unity throughout the range of $g$.
Actually, the minimal value of the overlap is as large as $0.995$ in the two
figures. This means that with three configurations, be the orbitals real or
complex, we can recover the exact ground state to very high precision. In the
intermediate case of $M=2$ configurations, the energy and the overlap are in-
between. A peculiarity is that with real orbitals, the curve of the overlap is
discontinuous at some point. This is due to the existence of local minima. At
the critical point, two local minima change order in energy, or more
precisely, the originally global minimum is surpassed by another minimum which
was originally just a local one.

In Figs. 15, 16, and 17, which correspond to $M=1$, $M=2$, and $M=3$,
respectively, we show the (real) orbitals $\phi_{i}^{(\alpha)}$, the density
distribution $\rho(x,x)$, and the correlation function $\rho(0,x)$ for three
different values of $g$. In Fig. 15, we see that the single configuration
approximation, in accord with Fig. 14, is quantitatively not very accurate for
large values of $g$. However, the orbitals are consistent with the
fermionization picture in the large-$g$ limit girardeau ; weiss ; yukalov . In
Fig. 16, with two configurations, the similarity between the variational
results and exact results improves, but the difference is still apparent for
$g=10$ and $15$. A notable feature of the two-configuration approximation is
that it breaks the parity symmetry of the model—the density distribution
$\rho(x,x)$ and the correlation function $\rho(0,x)$ are asymmetric. This kind
of phenomena is quite common in the conventional Hartree-Fock approximation,
and is equally so with us. Below we shall see more examples. In Fig. 17, we
have $M=3$ configurations, and now the variational results almost coincide
with the exact results, as is anticipated by the energy and overlap
information in Fig. 14.

Figure 18: (Color online) The orbitals $\phi_{i}(x)$, the particle density
distribution $\rho(x,x)$, and the one-particle correlator $\rho(0,x)$ for
three different values of $g$. In this figure, we have $N=4$ bosons on an open
chain of $L=13$ sites.  Figure 19: (Color online) (a) Ground state energy per
particle of a Bose-Hubbard model with the open boundary condition. The solid
line is obtained by exact diagonalization, while the other lines are by the
permanent variational approach with $M=1$ to $M=5$ configurations ($M$
increases in the direction of the arrow). The variational calculation is done
with real orbitals. (b) Overlap between the exact ground state and the energy-
minimizing permanent state. As in Fig. 18, the particle number $N=4$ and the
lattice size $L=13$.

Figure 20: (Color online) The orbitals $\phi_{i}^{(\alpha)}(x)$, the particle
density distribution $\rho(x,x)$, and the one-particle correlator $\rho(0,x)$
for three increasing values of $M$. The setting and parameters are the same as
in Fig. 18(c).

Figure 21: (Color online) Same as Fig. 20, but the value of $g$ is 15 instead
of 4.

In the three models above, we see that the permanent approach with a very
limited number of configurations can yield very accurate results. The
observation is that regardless of the strength of the interaction, this is
often the case if the system is not very dilute, or more precisely, if the
ratio $N/L_{eff}$ is not too small, where by $L_{eff}$ we mean the effective
volume of the system, i.e., the volume accessible to the particles. In the
Bose-Hubbard model at unit filling, regardless of the boundary condition, the
ratio is 1; in the Bose-Hubbard model in a harmonic trap above, although the
lattice is of size $L=13$, by the shapes of the orbitals or the density
distribution, we infer that the effective lattice size $L_{eff}\simeq 7$, so
the ratio is about $3/7$.

Our experience is that the most challenging situation for the permanent
variational approach is a dilute gas in the Tonks-Girardeau limit, in which
the system is dilute (i.e., $N\ll L_{eff}$) and the on-site interaction $g$ is
strong. We consider such a Bose-Hubbard model in Fig. 18, where we have $N=4$
particles on a flat, open lattice of size $L=13$. When $g$ is as small as $1$,
the single-configuration approximation can get the density distribution
$\rho(x,x)$ and the correlator $\rho(0,x)$ accurately. As $g$ increases to 2,
the error becomes visible but is still small. However, if $g$ further
increases to $4$, the discrepancy becomes very apparent.

The natural remedy is to take multiple configurations. In Fig. 19, we show how
the variational energy and the overlap improve as the number of configurations
increases. The improvement is steady but slow in comparison with Fig. 14. At
$g=4$, which is small by the scale of Fig. 14, even with $M=5$ configurations,
the deviation of the variational energy from the exact value and the deviation
of the overlap from unity are still visible. In Fig. 20, we show how the
situation in Fig. 18(c) improves by taking more and more configurations. We
see that unlike the situation in Fig. 17, with $M=3$ configurations, the
variational curves are still manifestly different from the exact ones. Only
with $M=5$, do the two almost coincide.

The trend of the curves in Fig. 19 suggests that for even larger values of
$g$, we would need even more configurations to get a good approximation of the
exact ground state. In Fig. 21, we examine the case of $g=15$. For this value
of $g$, the exact ground state is already very close to its fermionization
limit at $g=+\infty$. We see that the single configuration approximation fails
blatantly and it even breaks the parity symmetry of the model. With $M=5$
configurations, the situation improves but only with $M=10$ configurations do
the variational predicted density distribution and correlation function agree
with the exact values very well.

We thus see that generally taking multiple configurations can effectively
reduce the error. The concern is how much price we have to pay. From Sec.
IV.1, we see that with $M$ configurations, the number of blocks that we have
to calculate for preparation of $\hat{F}$ and $\hat{G}$ increases by a factor
of $\frac{1}{2}M(M+1)$. This polynomial growth is mild.

###  V.3 Symmetry breaking and restoration

Figure 22: (Color online) The orbitals $\phi_{i}(x)$, the particle density
distribution $\rho(x,x)$, and the one-particle correlator $\rho(0,x)$ for
$N=3$ bosons in a double-well potential [see Eqs. (63) and (64)]. The
parameters are $L=13$, $g=2$, $\kappa=0.5$, $h=4.5$, $\sigma=1$. In (a) and
(c), the variational state consists of $M=1$ or $M=2$ configuration(s),
respectively. In (b), the state is built out of the state in (a) by the
projection process in (65). In each column, the energy $E_{gs}$ of the
variational state and its overlap $O$ with the exact ground state are shown.
The exact value of the ground state energy per particle is $0.9583$.

Figure 23: (Color online) The orbitals $\phi_{i}(x)$, the particle density
distribution $\rho(x,x)$, and the one-particle correlator $\rho(0,x)$ for
$N=3$ bosons on a closed chain with $L=11$ sites. The on-site interaction
strength $g=5$. The top row corresponds to the pro-projection variational
states with $M=1$ to $M=3$ configurations. The bottom row corresponds to the
projected variational states. As in Fig. 22, for each state, its energy and
its overlap with the exact ground state are shown. The exact value of the
ground state energy per particle is $-1.8526$.

It is a common observation in the practice of Hartree-Fock approximation that
the solution often spontaneously breaks the symmetries of the Hamiltonian
symmetry1 ; symmetry . This could also happen with us, as we have seen in Fig.
16 and Fig. 21 above, where the permanent variational states do not respect
the $\mathbb{Z}_{2}$ symmetry of the models.

We have two options to restore symmetry, i.e., to construct a state sharing
the same symmetry with the exact ground state. The first approach is by brute
force, we can simply take more configurations. Hopefully, the accuracy of the
approximation will improve and the symmetry is restored alongside. This
happens in Fig. 17 and Fig. 21(c). The second approach is based on the group
representation theory. We can construct a projection operator corresponding to
the irreducible representation of the exact ground state, and let it act on
the variational state. The resultant state is also a permanent variational
state, but generally with more configurations.

Below we take two concrete models to illustrate and compare the two
approaches. In the first model, we have $N=3$ bosons in a double-well
potential. The Hamiltonian is

$\displaystyle H_{dw}$ $\displaystyle=$ $\displaystyle
H_{obc}+\sum_{x=-L_{0}}^{L_{0}}V(x)a_{x}^{\dagger}a_{x},$ (63)

where the double-well potential $V(x)$ is the superposition of a harmonic trap
and a Gaussian bump, i.e.,

$\displaystyle V(x)$ $\displaystyle=$ $\displaystyle\kappa
x^{2}+h\exp(-x^{2}/\sigma^{2}),$ (64)

where $\kappa$, $h$, and $\sigma$ are parameters. The system has a
$\mathbb{Z}_{2}$ symmetry, and it can be easily proven that the parity of the
ground state is even and consequently the particle density distribution and
the correlator are both even functions of $x$. However, in Fig. 22(a1), in the
single configuration approximation, we see that the orbitals are apparently
asymmetric, with one orbital located in the left well and the rest two in the
right well. Consequently, the particle density and the correlator are both
asymmetric, and the deviation from the exact values is significant.

To reinstall symmetry by the first approach, we can simply take one more
configuration as in Fig. 22(c), where we get symmetric orbitals, and the
variational results agree with the exact ones perfectly. We can also try the
second approach. Let the permanent state in Fig. 22(a) be
$\Phi_{e}=\hat{\mathcal{S}}(\phi_{1},\phi_{2},\phi_{3})$. An even-parity state
can then be easily constructed as

$\displaystyle\bar{\Phi}_{e}$ $\displaystyle=$
$\displaystyle\hat{\mathcal{S}}(\phi_{1},\phi_{2},\phi_{3})+\hat{\mathcal{S}}(\hat{P}\phi_{1},\hat{P}\phi_{2},\hat{P}\phi_{3}),$
(65)

where the inversion operator $\hat{P}$ is defined as
$(\hat{P}\phi)(x)=\phi(-x)$. The newly constructed variational state
$\bar{\Phi}_{e}$ consists of two configurations, with the new configuration
transformed from the old one by inversion. The physical quantities calculated
with the projected state $\bar{\Phi}_{c}$ are shown in Fig. 22(b). We see
significant improvement over the pro-projection state $\Phi_{e}$ in Fig.
22(a). Qualitatively, the symmetry is restored; quantitatively, the overlap
with the exact ground state has increased from 0.91786 to 0.99255, and the
predicted ground state energy per particle has decreased from 0.98983 to
0.97049 (the exact value is 0.95831). Of course, by construction, we do not
expect $\bar{\Phi}_{e}$ to be optimal in energy among all the two-
configuration variational states. Indeed, its energy is slightly higher than
that of the optimal state in Fig. 22(c). However, its advantage is that it is
obtained for free and is still a fairly good approximation.

The second model is simply the Bose-Hubbard model with the periodic boundary
condition. Because of the boundary condition, the lattice can be visualized as
a closed lattice ring, and the symmetry group of the model is recognized as
that of a regular polygon, i.e., the dihedral group $D_{L}$. The group
consists of $L$ translations (or rotations) and $L$ reflections

$D_{L}=\\{\hat{S}^{m}\hat{T}^{n},0\leq m\leq 1,0\leq n\leq L-1\\}.$ (66)

Here the generating operators $\hat{S}$ and $\hat{T}$ are defined as
$(\hat{S}\phi)(x)=\phi(-x)$ and $(\hat{T}\phi)(x)=\phi(x-1)$. By the Perron-
Frobenius theorem, the exact ground state $|GS\rangle$ is non-degenerate and
positive everywhere in the Fock basis. It then follows easily that
$|GS\rangle$ belongs to the trivial representation of the dihedral group. The
projection operator for this irreducible representation is simply

$\displaystyle\hat{\mathcal{P}}$ $\displaystyle=$
$\displaystyle\frac{1}{2L}\sum_{m=0}^{1}\sum_{n=0}^{L-1}\hat{S}^{m}\hat{T}^{n}.$
(67)

Therefore, if we obtain an $M$-configuration variational state in the form of
(50) by the algorithm, by projection we obtain immediately the following state
invariant under all the symmetry transforms of the model,

$\displaystyle\bar{\Phi}_{e}$ $\displaystyle=$
$\displaystyle\sum_{m=0}^{1}\sum_{n=0}^{L-1}\sum_{\alpha=1}^{M}\hat{\mathcal{S}}(\hat{S}^{m}\hat{T}^{n}\phi^{(\alpha)}_{1},\ldots,\hat{S}^{m}\hat{T}^{n}\phi^{(\alpha)}_{N}),$
(68)

which is a $2LM$-configuration state.

In Fig. 23, we consider such a model with $(N,L,g)=(3,11,5)$. In the top row,
we see that with $M=1$ to $M=3$, the variational state always breaks the
symmetry of the model completely, i.e., all the translation and reflection
symmetries are lost. This is most easily seen from the density distribution.
In the bottom row, with the projected states, the symmetry of the exact ground
state is recovered. We see that generally, the projection process not only
restores the expected symmetry, but also lowers the energy and increases the
overlap. Remarkably, it is most effective in the single configuration case.

From Fig. 22 and Fig. 23, we see that the symmetry restoration procedure is a
worthwhile post-processing to the time-consuming iteration procedure. It is
fairly cheap and reasonably effective.

###  V.4 When multiple configurations are perfect

In the proceeding sections, we have seen that taking multiple configurations
can get us very close to the target state. Here, we discuss the scenario that
a generic target state can be exactly recovered by multiple configurations.

By Proposition 4, if the single-particle Hilbert space is of dimension $L=2$,
then for arbitrary $N$, an $N$-boson state is a permanent state, i.e., it
consists of a single configuration. For higher values of $L$, a generic state
is not a permanent state and a natural question is, at least how many
configurations we need to recover it exactly. A quick lower bound is obtained
by mere dimension counting. The dimension $\mathcal{D}$ of the many-body
Hilbert space is given in (5). By (18), a permanent state has $d=N(L-1)+1$
degrees of freedom. Hence, we need at least

$\displaystyle M_{0}$ $\displaystyle=$
$\displaystyle\lceil\mathcal{D}/d\rceil$ (69)

configurations to recover a generic state, where $\lceil x\rceil$ denotes the
least integer no less than $x$. There is no reason that this lower bound can
be achieved. Indeed, it is an underestimate in the special case of $N=2$. By
Proposition 4, we need $\lfloor(L+1)/2\rfloor$ configurations to recover a
generic $2$-boson state. Here the number scales as $L/2$ for large $L$.
However, the estimate of (69) is $M_{0}\simeq L/4$ for large $L$.

Although (69) fails for $N=2$, there are evidences that for many pairs of
$(L,N)$, it does give the right answer. Specifically, extensive numerical
experiments indicate that that if the pair $(L,N)$ take values among the set
of $\\{(3,2),(3,3),(3,4),(3,5),(4,3)\\}$, the $N$-boson wave function can
always be written as the summation of $M_{0}=2$ permanent states. Similarly,
if $(L,N)$ take values among the set of $\\{(4,4),(5,3)\\}$, the $N$-boson
wave function can always be written as the sum of $M_{0}=3$ permanent states.

The numerical experiment is done in the following way. We first generate a
random $N$-boson state, then generate a set of $M_{0}N$ random single-particle
orbitals, with $M_{0}$ given by (69), and then use the overlap maximization
algorithm to update the orbitals. The concern is whether the overlap will
surpass the threshold $1-10^{-5}$ after $300$ rounds of update. If not, a new
set of random orbitals are generated and the optimization process is
restarted. This process is repeated until in some run the threshold is
surpassed. If so, we turn to a new random $N$-boson state and repeat the
check.

For the set of values of $(L,N)$ mentioned above, we have checked over
$10^{5}$ random $N$-boson states, and they all passed the check. This is
strong evidence that for such $(L,N)$, the naive lower bound of (69) is
achieved.

So far, we have failed to find a rigorous proof of the findings above. Here we
just reformulate the problem in pure mathematics so that it might be more
convenient for further study. In second quantization, a generic $N$-boson
state is

$\displaystyle\Psi$ $\displaystyle=$
$\displaystyle\sum_{\textbf{n}}C_{\textbf{n}}\prod_{j=1}^{L}\left(a_{j}^{\dagger}\right)^{n_{j}}|vac\rangle.$
(70)

Here the summation is over all occupation tuple
$\textbf{n}=(n_{1},n_{2},\ldots,n_{L})$ with $n_{j}\geq 0$ and
$\sum_{j=1}^{L}n_{j}=N$. By (8), that it can be written as the sum of $M$
permanent states means

$\displaystyle\Psi$ $\displaystyle=$
$\displaystyle\sum_{\alpha=1}^{M}\prod_{i=1}^{N}\left(\sum_{j=1}^{L}C_{ij}^{(\alpha)}a_{j}^{\dagger}\right)|vac\rangle,$
(71)

where $C_{ij}^{(\alpha)}$ are constants. Because the $a_{j}^{\dagger}$
operators commute, and the (unnormalized) Fock states
$\prod_{j=1}^{L}\left(a_{j}^{\dagger}\right)^{n_{j}}|vac\rangle$ are linearly
independent, this is equivalent to saying that the degree-$N$ homogeneous
polynomial in $L$ variables

$\displaystyle P(z_{1},z_{2},\ldots,z_{L})$ $\displaystyle=$
$\displaystyle\sum_{\textbf{n}}C_{\textbf{n}}\prod_{j=1}^{L}z_{j}^{n_{j}}$
(72)

is expressible as

$\displaystyle P(z_{1},z_{2},\ldots,z_{L})$ $\displaystyle=$
$\displaystyle\sum_{\alpha=1}^{M}\prod_{i=1}^{N}\left(\sum_{j=1}^{L}C_{ij}^{(\alpha)}z_{j}\right).$
(73)

Formulated in this way, we see the problem is very similar to the polynomial
Waring problem landsberg . The difference is just that while in the Waring
problem, one seeks to decompose a general degree-$N$ homogeneous polynomial in
$N$-th powers of linear forms (linear polynomials), here we seek a
decomposition in terms of $N$-th products of linear forms.

###  V.5 A stability problem

Figure 24: (Color online) (a) Ground state energy (per particle) of the two-
boson model (75) as calculated by exact diagonalization (ED), or the permanent
variational state with a single configuration, or under the Gross-Pitaevskii
approximation (GPA). The dashed line indicates the ionization threshold (76)
of the double-bound state. The parameters are $(g,U)=(1.5,4.5)$. (b) Overlap
of the permanent state and the Gross-Pitaevskii state with the exact ground
state $|GS\rangle$. Note that the very limited memory of our laptop (less than
8 GB) allows us to do ED on a cubic lattice only up to the size $L_{1}=16$.

Figure 25: (Color online) (a) Ground state energy (per particle) of the two-
boson model (75) as calculated by exact diagonalization (ED), or the permanent
approach with $M=1$ or $M=2$ configuration(s), or under the Gross-Pitaevskii
approximation (GPA). The parameters are $(g,U)=(-2.5,3.5)$. (b) Overlap
between the exact ground state $|GS\rangle$ and the various variational
states.

Having checked the accuracy, reliability, and flexibility of the algorithm, we
now apply it to a realistic problem.

In atomic physics, a famous problem is the stability of the negative ion of
hydrogen rau . The question is, can we add an extra electron to a hydrogen
atom, or can a proton bind two electrons simultaneously? The delicacy is that
the repulsion between the two electrons is as strong as the attraction of the
proton to either of them. Historically, it took people much imagination,
insights, and endeavors to devise appropriate variational wave functions to
establish the stability of the system bethe ; hylleraas ; chandrasekhar .
Among all the trial wave functions people came up with, the one by
Chandrasekhar is the simplest—It is a permanent state! Explicitly, it is of
the form chandrasekhar

$\displaystyle f(r_{1},r_{2})$ $\displaystyle=$ $\displaystyle e^{-\alpha
r_{1}-\beta r_{2}}+e^{-\beta r_{1}-\alpha r_{2}},$ (74)

where $r_{1,2}$ are the distances of the two electrons to the proton and
$\alpha$, $\beta$ are two variational parameters. The energy minimum is
achieved at $(\alpha,\beta)=(1.03925,0.28309)$. Note that $\alpha\gg\beta$,
which means that one orbital is “in”, while the other is far “out”. This
configuration is of course reasonable. The outer electron feels an almost
completely shielded potential, and therefore can only be loosely bound.

Here we consider a lattice variant of the problem bic1 . Suppose we have two
on-site interacting bosons on a cubic lattice with a single defect site. The
Hamiltonian is

$\displaystyle H=-\sum_{\langle
i,j\rangle}(a_{i}^{\dagger}a_{j}+\text{h.c.})+\frac{g}{2}\sum_{i}a_{i}^{\dagger}a_{i}^{\dagger}a_{i}a_{i}-Ua_{0}^{\dagger}a_{0}.\quad$
(75)

The defect site is taken as the origin of the lattice, which theoretically
should be infinite but in our numerical simulation will be assumed to be a
$L_{1}\times L_{1}\times L_{1}$ one and will be put as symmetric as possible
with respect to the defect site. Apparently, the defect site is introduced to
mimic the proton attraction and the on-site interaction to mimic the electron-
electron repulsion.

It is well known that if the defect potential $U$ is strong enough, i.e., if
$U$ is larger than some critical value $U_{c}$ which is about $3.957$, the
defect can induce a single-particle defect model localized around it green .
The two bosons can then be trapped in this defect mode simultaneously if they
are not interacting. The problem is whether they can still be bound by the
defect potential if the repulsion $g$ between them is turned on, and tuned to
a certain strength.

Denote the energy of the defect mode as $E_{df}$. At $U=U_{c}$, $E_{ef}=-6$.
As $U$ increases, $E_{df}$ decreases monotonically. Its exact value can be
easily calculated by exact diagonalization with a sufficiently large lattice,
say, $L_{1}=101$. The ionization threshold of the double-bound state is
apparently

$\displaystyle E_{th}$ $\displaystyle=$ $\displaystyle-6+E_{df},$ (76)

which corresponds to the state in which one boson is trapped in the defect
mode and the other boson free and at the bottom of the energy band.

The ground state of the two-boson system can then be shown to be a double-
bound state if we can show that the ground state energy is below the threshold
energy. In Fig. 24, we study a concrete case with $(g,U)=(1.5,4.5)$. Note that
for this value of $U$, the defect mode is very shallow, with the energy
$E_{df}=-6.136$, so the fate of the double-bound state is uncertain without
calculation. In Fig. 24(a), we see that the variational calculation with a
permanent state with a single configuration yields very accurate results—the
difference with those obtained by exact diagonalization is hardly visible even
by the minute scale of the plot. Limited by the memory capacity of our laptop,
we can do exact diagonalization only up to $L_{1}=16$, which as we see is
insufficient for proving binding. However, with the variational method, we can
work up to $L_{1}=21$ and succeed in pushing the ground state energy $E_{gs}$
below the threshold, thus proving binding. For comparison, we have also
checked the Gross-Pitaevskii approximation, which as a special case of the
permanent approximation, forces the two orbitals to be the same. We see that
the error is much larger—too large to prove binding actually. In Fig. 24(b),
we examine the quality of the permanent and the Gross-Pitaevskii variational
states by considering their overlap with the exact ground state $|GS\rangle$
obtained by exact diagonalization. Again, we see that the permanent state is
an exceedingly good variational state—its overlap with the exact state is as
high as $0.9994$ even for $L_{1}=16$. This vividly demonstrates the relevance
of the permanent state in bosonic systems. The Gross-Pitaevskii state is also
good, but not as good as the permanent state. Its errors in energy and overlap
are 15-50 times larger. We thus see how crucial it is to allow the orbitals to
vary independently.

As a second case, we consider $(g,U)=(-2.5,3.5)$ in Fig. 25. This value of $U$
is inadequate for the formation of a defect mode. However, as the interaction
between the two bosons is now attractive, a double-bound state is still
possible as the particle-particle attraction can reinforce the defect
potential. This is indeed the case. As we see in Fig. 25(a), both exact
diagonalization and the permanent variational state with $M=2$ configurations
can push the ground state energy below the threshold, which is now simply
$E_{th}=-12$. Unlike in Fig. 24, here the permanent approach with a single
configuration fails to prove binding. It yields identical results with the
Gross-Pitaevskii approach. This indicates that the permanent state actually
degenerates into the Gross-Pitaevskii state, a fact further confirmed by
studying the overlap of these states with the exact state in Fig. 25(b). In
both panels of Fig. 25, we see that in this particular problem and with this
particular set of parameters, by including one more configuration, the error
in energy and overlap can be reduced by one order of magnitude. This
exemplifies the necessity and effectiveness of taking multiple configurations.

Finally, we mention that for a two-body (or even a three-body) system like the
present one, the permanent calculation is never an issue. The preparation of
the operators $\hat{F}$ and $\hat{G}$ can be done swiftly, and the most time-
consuming part would be solving the generalized eigenvalue problem (46) or
(57) if the lattice size is in the order of thousands as here. But here are
some tricks. First, as we need just the smallest generalized eigenvalue, we
can invoke a Lanczos-type algorithm (in Matlab, just call `eigs`); second,
actually we do not need to construct $\hat{F}$ as a matrix explicitly. In a
Lanczos-type algorithm, we just need to know its action on a vector, which is
very simple and can be done swiftly as $\hat{F}$ is simply structured. With
these tricks, in Fig. 24, it takes only 1 (10, respectively) second(s) to
update an orbital on an $L_{1}=15$ ($L_{1}=21$, respectively) cubic lattice on
our laptop. The observation is that generally convergence is achieved after 10
rounds of update in this particular problem.

##  VI Conclusions and open problems

We have explored the potential of the permanent state as variational wave
functions for bosons. The result is very encouraging. First, we found that for
the one-dimensional Bose-Hubbard model with periodic boundary condition and at
unit filling, the exact ground state can be well approximated by a permanent
state with translation-related orbitals. The permanent state overlaps well
with the exact ground state, yields energy close to the exact value, and
produces correlation functions close to the exact ones. Then with an iteration
algorithm, we examined more general models. It is quite often that the a
single permanent state approximates the exact ground state very well, by all
the criterions of energy, overlap, and correlation functions. In case the
discrepancy is apparent, it can be remedied by including more configurations.

The algorithm has its advantages and disadvantages. Let us first address its
disadvantages. The primary drawback of a permanent-based approach is of course
the permanent computation, which scales unfavorably with the particle number
$N$. However, with current computational facilities, it is not prohibitively
expensive. On our laptop, in the single-configuration case, it takes about
$5.1$ seconds to update one orbital for $N=L=12$, and the time reduces to
$1.0$ seconds if $N=10$. Hence, it is totally feasible to study a large enough
few-boson system with the algorithm. Note that there is still room for
acceleration, as computing the permanents of the minors of the overlap matrix
can be easily parallelized. Note also that in practice many models of interest
have only $N=2$ or $N=3$ particles. For such small values of $N$, the
permanent computation is of course not an issue at all. While the scaling of
the complexity of the algorithm with respect to the particle number $N$ is not
that favorable, the scaling with respect to the system volume $L$ is quite
favorable. The observation is that for $L\lesssim N$, most time is spent on
the permanent calculation and the $L$-dependence is negligible. Only for $L\gg
N$, the time needed to update one orbital increases apparently with $L$, but
still it increases at most in a polynomial way. For example, for $N=12$, the
time is $5.5$, $7.6$, and $21.4$ seconds for $L=25,50,100$, respectively.
Roughly speaking, while we are indeed confined to a limited number of
particles, we have essentially no limitation on the system volume, nor the
dimensionality of the system. The memory needed by the algorithm is also
minimal. While in this paper we have only considered the $(N,L)$ pairs for
which the dimension of the many-body Hilbert space is at most on the order of
one million, so that exact diagonalization is possible and we have exact
results for comparison, the algorithm can handle other values of $(N,L)$
easily.

It might be a good idea to combine the current algorithm with the Lanczos
algorithm. That is, one can think of post-processing the single- or multi-
configurational variational state $\Phi$ obtained by the current algorithm by
a Lanczos-type process. Starting from $\Phi$, we apply the Hamiltonian $H$
repeatedly to it, and at each stage, we try to approximate the resultant wave
function with a combination of permanent states. In this way, we construct a
Krylov-type subspace spanned by a set of multi-configurational states. A
variational wave function better than $\Phi$ can then be sought in this
subspace by solving a generalized eigenvalue problem. The procedures can be
carried out most conveniently for the two-particle ($N=2$) case, thanks to the
Autonne-Takagi factorization (22). Some preliminary test yields encouraging
results.

In this tentative work, we have focused on some one-dimensional lattice models
in the Bose-Hubbard category. Generalization to continuum models zinner ;
pilati , higher dimensions, and multi-component systems should be
straightforward. Moreover, the simple strategy of the algorithm apparently is
also applicable to fermions—Actually, originally it was used for fermions
zhang1 ; zhang2 , although for some different problem. It should be worthwhile
to check how it works for fermions.

Below are some open problems.

The fact that for the one-dimensional Bose-Hubbard model with periodic
boundary condition and at unit filling, the exact ground state can be well
approximated by a permanent state with translation-related orbitals is very
impressive. This should arouse one’s interest in the permanent state in its
own right. In the field of cold atoms, the folklore is that Bose-Einstein
condensation occurs as the temperature lowers, the de Broglie wave lengths of
the particles increase, and the wave-packets overlap. Now, a permanent state
with translation-related orbitals is in accord with this picture. Given a
primitive wave packet $\phi(\vec{x})$, we can imagine a many-body permanent
state constructed with the wave packets $\phi(\vec{x}-\vec{R}_{m})$, where
$\vec{R}_{m}$ runs through an $n$-dimensional lattice. The overlap between
adjacent orbitals can be well adjusted by changing the length scale of the
primitive orbital $\phi$. The concern is, is it possible to realize a
transition by tuning its length scale? How does the correlation function
depend on the primitive orbital? What is the effect of the dimensionality of
the lattice? In short, can we use a permanent state with regularly distributed
identical wave-packets as a prototypical wave function to model the
condensation transition? Note that the problem does not refer to a
Hamiltonian.

That a single permanent state is often a very good approximation of the exact
ground state motivates two problems. First, is it possible to construct an
interacting Hamiltonian comment  whose ground state is exactly a permanent
state? Second, how _dense_ are the permanent states in the many-boson Hilbert
space? Quantitatively, is there a number $\delta>0$, such that for any many-
boson state there exists a permanent state whose overlap with it is at least
$\delta$? Preliminary study suggests that $\delta\geq 0.15$ for
$(N,L)=(3,13)$. This is not a small number in view of the fact that the
dimension of the few-body Hilbert space is 455. The ground state of a
realistic model is non-generic, so the largest possible value of the overlap
of a permanent state with it should be much higher.

In this paper, we have focused on the energy-minimization problem. The other
problem of overlap-maximization, i.e., the problem of finding the optimal
single- or multi-configurational permanent approximation of a given bosonic
wave function should also be a worthy one. It is about the structure of a
bosonic wave function. At least for fermions, it is now well-known that the
anti-symmetry condition entails deep structures of the fermionic wave
function, with consequences far beyond the commonplace of the Pauli exclusion
principle borland1 ; borland2 ; klyachko ; coleman ; schilling , and the
notion of optimal Slater approximation has proven to be useful in this study
zhang1 ; zhang2 ; mbl . It is fair to expect that for bosons, the symmetry
condition also has far-reaching consequences and hopefully, the notion of
optimal permanent approximation is a useful one too.

## Acknowledgments

The authors are grateful to J. Guo, K. Yang, Y. Xiang and K. Jin for their
helpful comments. This work is supported by the Science Challenge Project (NO.
TZ2018002) and the Foundation of LCP.

*

##  Appendix A Expressions of $\hat{F}$ and $\hat{G}$ in Sec. IV

Let us start with $\hat{G}$ which is simpler than $\hat{F}$. By an analogy of
the Laplace expansion for the determinant, we have

$\displaystyle\langle\Phi|\Phi\rangle$ $\displaystyle=$
$\displaystyle\operatorname{per}(A)=\langle\phi_{1}|\phi_{1}\rangle\operatorname{per}(A;1|1)+\sum_{j_{1}=2}^{N}\langle\phi_{1}|\phi_{j_{1}}\rangle\operatorname{per}(A;1|j_{1})$
(77) $\displaystyle=$
$\displaystyle\langle\phi_{1}|\phi_{1}\rangle\operatorname{per}(A;1|1)+\sum_{j_{1}=2}^{N}\sum_{i_{1}=2}^{N}\langle\phi_{1}|\phi_{j_{1}}\rangle\langle\phi_{i_{1}}|\phi_{1}\rangle\operatorname{per}(A;1,i_{1}|1,j_{1}).$

From this expression, we can read off the operator $\hat{G}$ defined by
$\langle\phi_{1}|\hat{G}|\phi_{1}\rangle=\langle\Phi|\Phi\rangle$. It is

$\displaystyle\hat{G}$ $\displaystyle=$
$\displaystyle\operatorname{per}(A;1|1)\hat{I}+\sum_{j_{1}=2}^{N}\sum_{i_{1}=2}^{N}\operatorname{per}(A;1,i_{1}|1,j_{1})|\phi_{j_{1}}\rangle\langle\phi_{i_{1}}|,$
(78)

where $\hat{I}$ is the identity operator. Note that $\hat{G}$ depends on the
orbitals $\phi_{2\leq j\leq N}$ but not on $\phi_{1}$.

We then turn to $\hat{F}$. For the single-particle part, we have

$\displaystyle\langle\Phi|H_{1}|\Phi\rangle$ $\displaystyle=$
$\displaystyle\sum_{i_{1}=1}^{N}\sum_{j_{1}=1}^{N}\langle\phi_{i_{1}}|K|\phi_{j_{1}}\rangle\operatorname{per}(A;i_{1}|j_{1})$
(79) $\displaystyle=$
$\displaystyle\langle\phi_{1}|K|\phi_{1}\rangle\operatorname{per}(A;1|1)+\sum_{j_{1}=2}^{N}\langle\phi_{1}|K|\phi_{j_{1}}\rangle\operatorname{per}(A;1|j_{1})+\sum_{i_{1}=2}^{N}\langle\phi_{i_{1}}|K|\phi_{1}\rangle\operatorname{per}(A;i_{1}|1)$
$\displaystyle+\sum_{i_{1}=2}^{N}\sum_{j_{1}=2}^{N}\langle\phi_{i_{1}}|K|\phi_{j_{1}}\rangle\operatorname{per}(A;i_{1}|j_{1})$
$\displaystyle=$
$\displaystyle\langle\phi_{1}|K|\phi_{1}\rangle\operatorname{per}(A;1|1)+\sum_{j_{1}=2}^{N}\sum_{i_{1}=2}^{N}\langle\phi_{1}|K|\phi_{j_{1}}\rangle\langle\phi_{i_{1}}|\phi_{1}\rangle\operatorname{per}(A;1,i_{1}|1,j_{1})$
$\displaystyle+\sum_{j_{1}=2}^{N}\sum_{i_{1}=2}^{N}\langle\phi_{1}|\phi_{j_{1}}\rangle\langle\phi_{i_{1}}|K|\phi_{1}\rangle\operatorname{per}(A;1,i_{1}|1,j_{1})+\sum_{i_{1}=2}^{N}\sum_{j_{1}=2}^{N}\langle\phi_{i_{1}}|K|\phi_{j_{1}}\rangle\langle\phi_{1}|\phi_{1}\rangle\operatorname{per}(A;1,i_{1}|1,j_{1})$
$\displaystyle+\sum_{i_{2}\neq i_{1},2}^{N}\sum_{j_{1}\neq
j_{2},2}^{N}\langle\phi_{i_{1}}|K|\phi_{j_{1}}\rangle\langle\phi_{1}|\phi_{j_{2}}\rangle\langle\phi_{i_{2}}|\phi_{1}\rangle\operatorname{per}(A;1,i_{1},i_{2}|1,j_{1},j_{2}).$

Here in the last line, the summation $\sum_{i_{2}\neq i_{1},2}^{N}$ means that
$i_{1}$ and $i_{2}$ both run from $2$ to $N$, but they must take different
values. Similar summation expressions below should be interpreted similarly.
We see that the contribution of $H_{1}$ to $\hat{F}$ is

$\displaystyle\operatorname{per}(A;1|1)K+\sum_{j_{1}=2}^{N}\sum_{i_{1}=2}^{N}K|\phi_{j_{1}}\rangle\langle\phi_{i_{1}}|*\operatorname{per}(A;1,i_{1}|1,j_{1})+\sum_{j_{1}=2}^{N}\sum_{i_{1}=2}^{N}|\phi_{j_{1}}\rangle\langle\phi_{i_{1}}|K*\operatorname{per}(A;1,i_{1}|1,j_{1})$
$\displaystyle+\sum_{i_{1}=2}^{N}\sum_{j_{1}=2}^{N}\langle\phi_{i_{1}}|K|\phi_{j_{1}}\rangle\hat{I}*\operatorname{per}(A;1,i_{1}|1,j_{1})+\sum_{i_{2}\neq
i_{1},2}^{N}\sum_{j_{1}\neq
j_{2},2}^{N}\langle\phi_{i_{1}}|K|\phi_{j_{1}}\rangle|\phi_{j_{2}}\rangle\langle\phi_{i_{2}}|*\operatorname{per}(A;1,i_{1},i_{2}|1,j_{1},j_{2}).\quad\quad$
(80)

Next we turn to the two-particle or the interaction term $H_{2}$. We have by
(34)

$\displaystyle\langle\Phi|H_{2}|\Phi\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{i_{1}\neq i_{2},1}^{N}\sum_{j_{1}\neq
j_{2},1}^{N}\langle\phi_{i_{1}}\phi_{i_{2}}|U|\phi_{j_{1}}\phi_{j_{2}}\rangle\operatorname{per}(A;i_{1},i_{2}|j_{1},j_{2}).$
(81)

For clarity, three cases will be considered separately. In the first case, two
$\phi_{1}$’s are associated with $U$. We have

$\displaystyle\frac{1}{2}\sum_{i_{2}=2}^{N}\sum_{j_{2}=2}^{N}\langle\phi_{1}\phi_{i_{2}}|U|\phi_{1}\phi_{j_{2}}\rangle\operatorname{per}(A;1,i_{2}|1,j_{2})+\frac{1}{2}\sum_{i_{1}=2}^{N}\sum_{j_{1}=2}^{N}\langle\phi_{i_{1}}\phi_{1}|U|\phi_{j_{1}}\phi_{1}\rangle\operatorname{per}(A;1,i_{1}|1,j_{1})$
(82)
$\displaystyle+\frac{1}{2}\sum_{i_{1}=2}^{N}\sum_{j_{2}=2}^{N}\langle\phi_{i_{1}}\phi_{1}|U|\phi_{1}\phi_{j_{2}}\rangle\operatorname{per}(A;1,i_{1}|1,j_{2})+\frac{1}{2}\sum_{i_{2}=2}^{N}\sum_{j_{1}=2}^{N}\langle\phi_{1}\phi_{i_{2}}|U|\phi_{j_{1}}\phi_{1}\rangle\operatorname{per}(A;1,i_{2}|1,j_{1})$
$\displaystyle=$
$\displaystyle\sum_{i_{2}=2}^{N}\sum_{j_{2}=2}^{N}\langle\phi_{1}\phi_{i_{2}}|U|\phi_{1}\phi_{j_{2}}\rangle\operatorname{per}(A;1,i_{2}|1,j_{2})+\sum_{i_{1}=2}^{N}\sum_{j_{2}=2}^{N}\langle\phi_{i_{1}}\phi_{1}|U|\phi_{1}\phi_{j_{2}}\rangle\operatorname{per}(A;1,i_{1}|1,j_{2})$
$\displaystyle=$
$\displaystyle\sum_{i_{2}=2}^{N}\sum_{j_{2}=2}^{N}\langle\phi_{1}\phi_{i_{2}}|U|\phi_{1}\phi_{j_{2}}\rangle\operatorname{per}(A;1,i_{2}|1,j_{2})+\sum_{i_{2}=2}^{N}\sum_{j_{2}=2}^{N}\langle\phi_{i_{2}}\phi_{1}|U|\phi_{1}\phi_{j_{2}}\rangle\operatorname{per}(A;1,i_{2}|1,j_{2}).$

Contributions of these expressions to the operator $\hat{F}$ can be easily
read off. For example, the operator $\mathcal{O}$ corresponding to the matrix
element $\langle\phi_{1}\phi_{i_{2}}|U|\phi_{1}\phi_{j_{2}}\rangle$ is defined
by
$\langle\phi_{1}|\mathcal{O}|\phi_{1}\rangle=\langle\phi_{1}\phi_{i_{2}}|U|\phi_{1}\phi_{j_{2}}\rangle$.
For the Bose-Hubbard model in which $U$ is an on-site interaction (37), the
operator $\mathcal{O}$ is actually a (generally complex) potential with the
explicit expression $\mathcal{O}(x)=g\phi_{i_{2}}^{*}(x)\phi_{j_{2}}(x)$.

In the second case, one $\phi_{1}$ is associated with $U$. We have

$\displaystyle\frac{1}{2}\sum_{i_{2}=2}^{N}\sum_{j_{1}\neq
j_{2},2}^{N}\langle\phi_{1}\phi_{i_{2}}|U|\phi_{j_{1}}\phi_{j_{2}}\rangle\operatorname{per}(A;1,i_{2}|j_{1},j_{2})+\frac{1}{2}\sum_{i_{1}=2}^{N}\sum_{j_{1}\neq
j_{2},2}^{N}\langle\phi_{i_{1}}\phi_{1}|U|\phi_{j_{1}}\phi_{j_{2}}\rangle\operatorname{per}(A;1,i_{1}|j_{1},j_{2})$
(83) $\displaystyle+\frac{1}{2}\sum_{i_{1}\neq
i_{2},2}^{N}\sum_{j_{2}=2}^{N}\langle\phi_{i_{1}}\phi_{i_{2}}|U|\phi_{1}\phi_{j_{2}}\rangle\operatorname{per}(A;i_{1},i_{2}|1,j_{2})+\frac{1}{2}\sum_{i_{1}\neq
i_{2},2}^{N}\sum_{j_{1}=2}^{N}\langle\phi_{i_{1}}\phi_{i_{2}}|U|\phi_{j_{1}}\phi_{1}\rangle\operatorname{per}(A;i_{1},i_{2}|1,j_{1})$
$\displaystyle=$ $\displaystyle\sum_{i_{2}=2}^{N}\sum_{j_{1}\neq
j_{2},2}^{N}\langle\phi_{1}\phi_{i_{2}}|U|\phi_{j_{1}}\phi_{j_{2}}\rangle\operatorname{per}(A;1,i_{2}|j_{1},j_{2})+\sum_{i_{1}\neq
i_{2},2}^{N}\sum_{j_{2}=2}^{N}\langle\phi_{i_{1}}\phi_{i_{2}}|U|\phi_{1}\phi_{j_{2}}\rangle\operatorname{per}(A;i_{1},i_{2}|1,j_{2})$
$\displaystyle=$ $\displaystyle\sum_{i_{1}\neq i_{2},2}^{N}\sum_{j_{1}\neq
j_{2},2}^{N}\langle\phi_{1}\phi_{i_{2}}|U|\phi_{j_{1}}\phi_{j_{2}}\rangle\langle\phi_{i_{1}}|\phi_{1}\rangle\operatorname{per}(A;1,i_{1},i_{2}|1,j_{1},j_{2})$
$\displaystyle+\sum_{i_{1}\neq i_{2},2}^{N}\sum_{j_{1}\neq
j_{2},2}^{N}\langle\phi_{1}|\phi_{j_{1}}\rangle\langle\phi_{i_{1}}\phi_{i_{2}}|U|\phi_{1}\phi_{j_{2}}\rangle\operatorname{per}(A;1,i_{1},i_{2}|1,j_{1},j_{2}).$

In the third case, none $\phi_{1}$ is associated with $U$. We have

$\displaystyle\frac{1}{2}\sum_{i_{1}\neq i_{2},2}^{N}\sum_{j_{1}\neq
j_{2},2}^{N}\langle\phi_{i_{1}}\phi_{i_{2}}|U|\phi_{j_{1}}\phi_{j_{2}}\rangle\operatorname{per}(A;i_{1},i_{2}|j_{1},j_{2})$
(84) $\displaystyle=$ $\displaystyle\frac{1}{2}\sum_{i_{1}\neq
i_{2},2}^{N}\sum_{j_{1}\neq
j_{2},2}^{N}\langle\phi_{i_{1}}\phi_{i_{2}}|U|\phi_{j_{1}}\phi_{j_{2}}\rangle\langle\phi_{1}|\phi_{1}\rangle\operatorname{per}(A;1,i_{1},i_{2}|1,j_{1},j_{2})$
$\displaystyle+\frac{1}{2}\sum_{i_{1}\neq i_{2}\neq
i_{3},2}^{N}\sum_{j_{1}\neq j_{2}\neq
j_{3},2}^{N}\langle\phi_{i_{1}}\phi_{i_{2}}|U|\phi_{j_{1}}\phi_{j_{2}}\rangle\langle\phi_{1}|\phi_{j_{3}}\rangle\langle\phi_{i_{3}}|\phi_{1}\rangle\operatorname{per}(A;1,i_{1},i_{2},i_{3}|1,j_{1},j_{2},j_{3}).$

In the last line, we see that we have to calculate the permanents of a series
of $(N-4)\times(N-4)$ matrices $\binom{N-1}{3}^{2}$ times. By the improved
Ryser algorithm, the total evaluation is on the order of $N^{7}2^{(N-5)}$.
This is the most time-consuming part in the preparation of $\hat{F}$ and
$\hat{G}$, if $N\geq 5$.

## References

  * (1) D. R. Hartree, Math. Proc. Camb. Philos. Soc. 24, 89 (1928).
  * (2) J. C. Slater, Phys. Rev. 32, 339 (1928).
  * (3) J. A. Gaunt, Math. Proc. Camb. Philos. Soc. 24, 328 (1928).
  * (4) V. A. Fock, Z. Phys. 61, 126 (1930).
  * (5) J. C. Slater, Phys. Rev. 35, 210 (1930).
  * (6) A. Szabo and N. S. Ostlund, _Modern Quantum Chemistry_ (Dover, New York, 1982).
  * (7) Actually, in the fermionic case, the orbitals $\phi_{1\leq i\leq N}$ do not need to be orthonormal neither. However, the resultant Slater state would be essentially the same if the orbitals are orthonormalized at the beginning.
  * (8) E. P. Gross, Nuovo Cimento 20, 454 (1961); E. P. Gross, J. Math. Phys. 4, 195 (1963).
  * (9) L. P. Pitaevskii, Sov. Phys. JETP 13, 451 (1961).
  * (10) F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 71, 463 (1999).
  * (11) A. J. Leggett, Rev. Mod. Phys. 73, 307 (2001).
  * (12) I. Romanovsky, C. Yannouleas and U. Landman, Phys. Rev. Lett. 93, 230405 (2004).
  * (13) I. Romanovsky, C. Yannouleas, L. O. Baksmaty, and U. Landman , Phys. Rev. Lett. 97, 090401 (2006).
  * (14) M. Heimsoth and M. Bonitz, Physica E 42, 420 (2010).
  * (15) M. Heimsoth, Diploma Thesis _Unrestricted Hartree-Fock Theory of Nonideal Bose Systems_ (Christian-Albrechts Universität zu Kiel, 2009).
  * (16) N. N. Bogoliubov and N. N. Bogoliubov Jr., _Introduction to Quantum Statistical Mechanics_ 2nd ed. (World Scientifc, Singapore, 2010).
  * (17) A. I. Streltsov, O. E. Alon, L. S. Cederbaum, Phys. Rev. A 73, 063626 (2006).
  * (18) B. Farid, arXiv:1912.00474.
  * (19) S. Rotman, _Advanced Linear Algebra_ 3rd ed. (Springer, Berlin, 2008).
  * (20) J. B. Fraleigh, _A First Course in Abstract Algebra_ 7th ed. (Addison-Wesley, Boston, 2003).
  * (21) Y. A. Aoto, M. F. da Silva, Phys. Rev. A 102, 052803 (2020).
  * (22) While it is known that the Slater determinant states form a manifold, namely the Grassmannian manifold, here it is conjectured that the set of permanent states is not a manifold.
  * (23) J. Links, A. Foerster, A. P. Tonel, and G. Santos, Ann. Henri Poincaré 7, 1591 (2006); A. P. Tonel, J. Links, and A. Foerster, J. Phys. A: Math. Gen. 38, 1235 (2005); M. T. Batchelor and A. Foerster, J. Phys. A: Math. Gen. 49, 173001 (2016).
  * (24) J. M. Zhang and M. Kollar, Phys. Rev. A 89, 012504 (2014).
  * (25) J. M. Zhang and N. J. Mauser, Phys. Rev. A 94, 032513 (2016).
  * (26) R. A. Horn and C. R. Johnson, _Matrix Analysis_ 2nd ed. (Cambridge University Press, 2013).
  * (27) R. Bhatia, _Matrix Analysis_ (Springer, New York, 1997).
  * (28) H. J. Ryser, _Combinatorial Mathematics_ (Wiley, New York, 1963).
  * (29) A. Nijenhuis and H. S. Wilf, _Combinatorial Algorithms_ (Academic Press, New York, 1978).
  * (30) Some MATLAB subroutines for permanent computation using the Ryser algorithm and taking advantage of the Gray code are in the Supplementary Material.
  * (31) We have actually also tried Gaussian functions. However, they lead to much poorer approximations, so we do not bother to mention them.