# Ground state of the hard-core Bose gas on lattice I. Energy estimates

## 1 Introduction

Bosons on a lattice interacting via an infinite on-site repulsion (hard-core bosons) represent a system of double interest. It is the simplest example of an interacting Bose gas and, thus, the most promising candidate for a rigorous treatment of Bose-Einstein condensation of interacting particles. On the other hand, the model is known to be equivalent to a system of $`\frac{1}{2}`$ spins  coupled via the $`XY`$\- and possibly the $`Z`$-components of neighboring spins and exposed to an external magnetic field in the $`Z`$ direction. Ordering of the planar component of the spins is equivalent to Bose-Einstein condensation or the appearance of off-diagonal long-range order (ODLRO) in the system of bosons. In spite of a long and extensive study the results about ordering are far from being complete. Apart from some exceptions, like bounds on the density of the condensate  or the discussion of the model on the full graph , the most interesting and difficult results were formulated in spin terminology . These works made use of a particular symmetry, the reflexion positivity. This introduced severe limitations as to the value of the external field (zero field) and the lattice type (essentially hypercubic lattices). Translated into the language of the boson gas, ODLRO was shown at half-filling on hypercubic lattices in the ground state in and above two dimensions, and for low enough temperatures above two dimensions.

In this paper we apply the boson terminology. Let $`𝐋`$ be an infinite lattice which, for the sake of simplicity, will be supposed to be regular with a constant coordination number (valency) $`k`$. Throughout the paper $`\mathrm{\Lambda }`$ denotes a finite part of $`𝐋`$ equipped with a periodic boundary condition so as to keep the valency constant (not really essential). The Hamiltonian we are going to study is

$$H=\underset{\{x,y\}E\mathrm{\Lambda }}{}(b_x^{}b_y+b_y^{}b_x)$$
(1)   
We write $`x,y,\mathrm{}`$ for the vertices of $`𝐋`$, and $`E\mathrm{\Lambda }`$ for the set of edges of $`\mathrm{\Lambda }`$; $`b_x^{}`$ and $`b_x`$ create, resp., annihilate a hard-core boson at $`x`$. Boson operators at different sites commute with each other while

$$b_x^{}b_x+b_xb_x^{}=1$$
(2)   
accounts for the hard-core condition. Correspondence with spin models is obtained by setting $`b_x=S_x^{}`$ and $`b_x^{}=S_x^+`$. The Hamiltonian conserves the number of bosons,

$$N=\underset{x\mathrm{\Lambda }}{}n_x=\underset{x\mathrm{\Lambda }}{}b_x^{}b_x$$
(3)   
and is also invariant under particle-hole transformation. We can, therefore, fix $`N`$ so that $`\rho =N/|\mathrm{\Lambda }|`$ is between 0 and $`\frac{1}{2}`$. (Here and below, if $`A`$ is a finite set, $`|A|`$ denotes the number of its elements.) Our concern in this paper is to provide nontrivial upper bounds to the ground state energy, $`E_0`$. That such an estimate may be useful in the study of qualitative properties of the ground state, was an interesting point of .

Let $`X,Y,\mathrm{}`$ denote $`N`$-point subsets of $`\mathrm{\Lambda }`$, called also configurations. A convenient basis is formed by the states

$$\varphi (X)=\underset{xX}{}b_x^{}\mathrm{\Psi }_{\mathrm{vac}}$$
(4)   
where $`\mathrm{\Psi }_{\mathrm{vac}}`$ is the vacuum state. Variational estimates of the ground state energy are of the form

$$E_0\psi |H|\psi /\psi |\psi .$$
(5)   
A trivial choice is

$$\psi =\varphi (X)$$
(6)   
where the summation goes over all $`N`$-point subsets of $`\mathrm{\Lambda }`$. It yields (cf. Section 2)

$$E_0k|\mathrm{\Lambda }|\left(\genfrac{}{}{0pt}{}{|\mathrm{\Lambda }|2}{N1}\right)\left(\genfrac{}{}{0pt}{}{|\mathrm{\Lambda }|}{N}\right)^1=k\rho (1\rho )|\mathrm{\Lambda }|+O(1).$$
(7)   
The bound is nothing else than minus the average size of the boundaries of the configurations: We call the boundary of $`X`$ the set of half-filled edges and denote it by $`X`$. Hence, (7) is equivalent to

$$|E_0|\overline{|X|}\left(\genfrac{}{}{0pt}{}{|\mathrm{\Lambda }|}{N}\right)^1|X|$$
(8)   
with summation over $`N`$-point configurations. If $`\mathrm{\Lambda }`$ is a full graph, i.e., any two sites are neighbors, (6) is an exact eigenvector and (8) holds with equality. In any other case we have a strict inequality.

The present work is about trying to improve on this bound, i.e., to make the right-hand side of (8) larger by an amount proportional to $`|\mathrm{\Lambda }|`$. It is important to know that this is possible for any $`k`$ and any $`\rho `$. In the opposite case, if $`k\rho (1\rho )`$ were to be the true ground state energy per site, one could easily conclude that as in the full graph, the Hamiltonian (1) has a product ground state in infinite volume,

$$\mathrm{\Psi }=\underset{x}{}(\sqrt{\rho }|n_x=1+\sqrt{1\rho }|n_x=0)$$
(9)   
with ODLRO and the value of the order parameter at its theoretical maximum ($`\rho (1\rho )`$, cf. ).

Since by putting equal weights on every configuration, as in (6), yields the average boundary size (8), we expect that a larger value can be obtained by giving larger weights to configurations with larger boundaries. Therefore, we calculate variational bounds by using trial functions of the form

$$\psi _v=v_X\varphi (X)v_X=v(|X|)$$
(10)   
i.e., $`v_X`$ depending on $`|X|`$ only, and $`v(n)`$ being concentrated on $`n>\overline{|X|}`$. Computability depends on our ability to estimate the number of configurations with a given size of boundary. The logarithm of this number is the entropy of the Ising model in a microcanonical ensemble with a fixed magnetization, $`_{x\mathrm{\Lambda }}\sigma _x=|\mathrm{\Lambda }|(12\rho )`$. Indeed, if we put $`\sigma _x=1`$ for $`x`$ in $`X`$ and $`\sigma _x=1`$ elsewhere, we get an Ising configuration with the corresponding union of contours $`X`$ whose total length, $`|X|`$, is the energy of the Ising configuration. In mathematical terms, the distribution of $`|X|`$ satisfies a large deviation principle whose rate function is, apart from a shift, minus the specific entropy of the Ising model, cf. . The exact form of the entropy is unknown for two- and higher-dimensional lattices. To circumvent this problem, we use an approximate formula for the probability of having a boundary of length $`n`$,

$$P_{\mathrm{\Lambda },N}(|X|=n)Z^1\mathrm{exp}\{\frac{(nM)^2}{2D^2}\}.$$
(11)   
Here $`Z`$ is for normalization and

$$M=\overline{|X|}=k\rho (1\rho )|\mathrm{\Lambda }|D^2=\overline{(|X|\overline{|X|})^2}.$$
(12)   
What is the approximation in the above expression? First, we replace the smooth and concave specific entropy $`s(ϵ)`$, having a maximum at the specific Ising energy $`ϵ_m=k\rho (1\rho )`$, by a parabola. Since the improved bound on $`|E_0|/|\mathrm{\Lambda }|`$ we are going to derive is $`(1+\delta )ϵ_m`$ with a $`\delta `$ never exceeding 0.2, cf. Table 4 below, this seems to be a consistent approximation. Second, we surmise that the second derivative of $`s(ϵ)`$ at the maximum is $`|\mathrm{\Lambda }|/D^2`$. An analogous statement holds true for the large deviations of a sum of identically distributed independent random variables. Now $`|X|`$ is the sum of identically distributed random variables, although not independent but finitely dependent ones, see Section 3. So this approximation is hopefully good. In any case, we use the formula (11) without looking for further justification. It may happen that energies of trial states depending on $`|X|`$ only can be quite close to $`E_0`$. Then the error of the approximation (11) may invalidate our estimates as rigorous upper bounds. They still can be useful as approximate formulas showing the $`k`$\- and $`\rho `$-dependence of the ground state energy. We apply the term ‘bound’ with this reservation.

In Section 2 we present the general setup for the estimates. In Section 3 we derive an expression for $`D^2`$ which seems to be valid for any vertex-transitive graph (whose all vertices are equivalent). This formula may also be of interest in site percolation problems and in the approximation of the entropy of the Ising model in nonvanishing magnetic fields. In Section 4 we compute variational bounds by using for $`v(n)`$ step-, exponential and Gaussian functions. The last two will be seen to yield improved bounds for any value of $`k`$ and $`\rho `$. Section 5 summarizes the results and indicates the way they can be extended to the grand-canonical ensemble.

In the true ground state $`v_X`$ is not a function of $`|X|`$ only. The detailed dependence on $`X`$ may not be relevant for computing the ground state energy but is absolutely crucial for qualitative properties, like ordering. Wave functions of the kind (10) trivially show an off-diagonal long-range order. In the true ground state fluctuations around a function of $`|X|`$ destroy ODLRO in one dimension and may decrease the order parameter in higher dimensions. A discussion of the ground state wave function will be given in a subsequent paper .

This work was supported by the Hungarian Scientific Research Fund (OTKA) under Grant No. T 030543.

## 2 Setup for the energy estimates

It seems advantageous to consider the eigenvalue problem of (1) as a problem of graph theory. The $`N`$-point configurations are vertices of a huge graph that we call the power graph of $`\mathrm{\Lambda }`$ of order $`N`$ and denote by $`G=G_{\mathrm{\Lambda },N}`$. Two configurations $`X`$ and $`Y`$ form an edge of $`G`$ if $`Y`$ can be obtained from $`X`$ by moving a single particle from a site of $`X`$ to a neighboring site, unoccupied in $`X`$. Thus, $`X`$ and $`Y`$ have $`N1`$ common sites and differ on an edge of $`\mathrm{\Lambda }`$. If $`VG`$ and $`EG`$ denote, respectively, the set of vertices and edges of $`G`$ then

$$|VG|=\left(\genfrac{}{}{0pt}{}{|\mathrm{\Lambda }|}{N}\right)|EG|=|E\mathrm{\Lambda }|\left(\genfrac{}{}{0pt}{}{|\mathrm{\Lambda }|2}{N1}\right)=\frac{1}{2}k|\mathrm{\Lambda }|\left(\genfrac{}{}{0pt}{}{|\mathrm{\Lambda }|2}{N1}\right).$$
(13)   
The boundary $`X`$ of $`X`$ can be given a new interpretation as the set of neighbors of $`X`$ in $`G`$. So if $`X`$ and $`Y`$ form an edge then $`YX`$ and $`XY`$. The action of $`H`$ on $`G`$ is that of the usual lattice Laplacian with the exception that there is no subtracted diagonal term. The matrix of $`H`$ in the basis (4) is the adjacency matrix $`A=A(G)`$ of $`G`$, that is, $`A_{XY}=1`$ if $`X`$ and $`Y`$ are neighbors and zero otherwise. We are interested in the largest eigenvalue of $`A`$, $`\lambda _1=|E_0|`$ and in the corresponding eigenvector, $`a=(a_X)`$. $`G`$ is connected if $`\mathrm{\Lambda }`$ is connected (that we suppose), therefore $`A`$ is irreducible (ergodic) and the Perron-Frobenius theorem applies: $`\lambda _1`$ is nondegenerate and largest to absolute value, and $`a_X>0`$ for all $`X`$. We note that $`G`$ is bipartite if and only if $`\mathrm{\Lambda }`$ is bipartite, and hence $`\lambda _1`$ is an eigenvalue of $`A`$ if and only if $`\mathrm{\Lambda }`$ is bipartite.

We call $`n_{\mathrm{min}}`$, resp., $`n_{\mathrm{max}}`$ the minimum, resp., maximum of $`|X|`$ among the $`N`$-point configurations. Clearly, $`n_{\mathrm{min}}=O(N^{(d1)/d})`$ if $`𝐋`$ is a d-dimensional lattice, and ($`\rho \frac{1}{2}`$)

$$n_{\mathrm{max}}kN=k\rho |\mathrm{\Lambda }|.$$
(14)   
We have the trivial inequalities

$$\overline{|X|}\lambda _1n_{\mathrm{max}}.$$
(15)   
The first is the variational bound (8): Setting $`v_X1`$ we find

$$\frac{|\psi _v|H|\psi _v|}{\psi _v|\psi _v}=\frac{\underset{X}{}\underset{YX}{}1}{_X1}=\overline{|X|}$$
(16)   
which is the same as

$$\frac{(v,Av)}{(v,v)}=\frac{1}{|VG|}\underset{X,Y:\{X,Y\}EG}{}1=\frac{2|EG|}{|VG|}.$$
(17)   
The upper bound in (15) follows from the eigenvalue equation: Let $`X_0`$ be a configuration on which $`a_X`$ reaches its maximum, then

$$\lambda _1a_{X_0}=\underset{YX_0}{}a_Y|X_0|a_{X_0}.$$
(18)   
We note that for $`\rho =\frac{1}{2}`$ there is a better upper bound, $`\lambda _1\frac{1}{4}(k+1)|\mathrm{\Lambda }|`$ (, Theorem C.1).

In Section 4 we shall see that $`n_{\mathrm{max}}`$ plays a role in optimizing the lower bound to $`\lambda _1`$. It is therefore important to know $`n_{\mathrm{max}}`$ exactly. In (14) there is equality if $`\rho `$ is small enough. In particular, $`n_{\mathrm{max}}=k\rho |\mathrm{\Lambda }|`$ for all $`\rho \frac{1}{2}`$ on bipartite lattices, and it is an easy graphical exercise to see that equality holds for $`\rho \frac{1}{3}`$ on the triangular ($`k=6`$) and on the Kagomé ($`k=4`$) lattices. Moreover, for both lattices $`n_{\mathrm{max}}`$ is constant between the densities $`\frac{1}{3}`$ and $`\frac{1}{2}`$: $`2|\mathrm{\Lambda }|`$ for the triangular and $`\frac{4}{3}|\mathrm{\Lambda }|`$ for the Kagomé lattice. This can be seen from the following argument. In general, $`n_{\mathrm{max}}`$ is the ground state energy of the antiferromagnetic Ising model under the restriction that the magnetization is fixed, $`_{x\mathrm{\Lambda }}\sigma _x=(12\rho )|\mathrm{\Lambda }|`$. However, we do not need to deal with the restriction. In both cases the (unrestricted) ground state is known to be highly degenerate. Among the exponentially large number of ground state configurations there are nonmagnetized ones, corresponding to $`\rho =\frac{1}{2}`$, others with concentration of down-spins $`\rho =\frac{1}{3}`$, and between these two limits $`\rho `$ can vary by steps of $`\frac{1}{|\mathrm{\Lambda }|}`$. The rule is to flip zero-energy spins one by one. The common energy of all these configurations is easy to compute from the fact that in each triangle there is precisely one unsatisfied bond. This fixes the value of $`n_{\mathrm{max}}`$ as given above.

The variational bound (5) reads $`\lambda _1B(v)`$ where

$$B(v)\frac{(v,Av)}{(v,v)}=\underset{X}{}v_X\underset{YX}{}v_Y/\underset{X}{}v_X^2.$$
(19)   
In the estimations below a crucial role is played by the inequality

$$||X||Y||2(k1)\text{if}YX.$$
(20)   
When passing from $`X`$ to $`Y`$ a neighboring particle-hole pair is interchanged. For both the particle and the hole the number of neighbors of the opposite kind can change by at most $`k1`$, whence (20) follows. If $`v_X=v(|X|)`$ and $`v(n)`$ is a nondecreasing sequence then

$$B(v)\frac{\underset{n=n_{\mathrm{min}}}{\overset{n_{\mathrm{max}}}{}}n|\mathrm{\Omega }_n|v(n2k+2)v(n)}{_{n=n_{\mathrm{min}}}^{n_{\mathrm{max}}}|\mathrm{\Omega }_n|v(n)^2}$$
(21)   
where $`\mathrm{\Omega }_n`$ denotes the set of configurations with boundary length $`n`$.

Above it is understood that $`v(n)=v(n_{\mathrm{min}})`$ if $`n<n_{\mathrm{min}}`$. We have computed bounds given by the right member of (21), using step-functions and functions with an exponential or a faster increase. In all cases we have found no improvement with respect to the trivial bound if $`\rho `$ was in a neighborhood of $`\frac{1}{2}`$. Apparently, we have lost too much in the inequality (21).

There is, however, a way to compute $`B(v)`$ in leading order ($`|\mathrm{\Lambda }|`$) by making a further hypothesis on some details of the large deviation principle (11). Let $`\nu _i(X)`$ be the number of those neighbors of $`X`$ having a boundary length $`|X|+i`$. (So $`_i\nu _i(X)=|X|`$.) For $`v_X=v(|X|)`$ we have

$$B(v)=\frac{\underset{n}{}v(n)\underset{i=2k+2}{\overset{2k2}{}}v(n+i)\underset{X\mathrm{\Omega }_n}{}\nu _i(X)}{_nv(n)^2|\mathrm{\Omega }_n|}.$$
(22)   
Equation (11) is equivalent to

$$|\mathrm{\Omega }_n|/|\mathrm{\Omega }_m|\mathrm{exp}\{[(nM)^2(mM)^2]/2D^2\}.$$
(23)   
An analogous statement for the ratio of $`\nu _i(X)`$ and $`\nu _j(X)`$ is certainly wrong for all $`X`$ separately but may be correct for the ratio of their sums over $`\mathrm{\Omega }_n`$. In this hope we formulate the hypothesis that

$$\frac{\underset{X\mathrm{\Omega }_n}{}\nu _i(X)}{_{X\mathrm{\Omega }_n}\nu _j(X)}\frac{e^{\frac{(n+iM)^2}{2D^2}}}{e^{\frac{(n+jM)^2}{2D^2}}}=e^{\frac{nM}{D^2}(ij)}2k+2i,j2k2.$$
(24)   
The third member of (24) is obtained by dropping a term of order $`|\mathrm{\Lambda }|^1`$. With (11) or (23) and (24) Eq. (22) reads

$`B(v)`$ $`=`$ $`\left({\displaystyle \underset{n}{}}e^{\frac{(nM)^2}{2D^2}}v(n)^2\right)^1{\displaystyle \underset{n}{}}ne^{\frac{(nM)^2}{2D^2}}v(n)v_n`$ (25)   
$`v_n`$ $`=`$ $`\left({\displaystyle \underset{j=2k+2}{\overset{2k2}{}}}e^{\frac{nM}{D^2}j}\right)^1{\displaystyle \underset{i=2k+2}{\overset{2k2}{}}}e^{\frac{nM}{D^2}i}v(ni).`$ (26)   
This is the starting point of the estimates of Section 4.

## 3 Statistics of the boundary lengths

In this section we show that the mean square deviation of $`|X|`$ is given by

$$D^2=[k2(2k1)\rho (1\rho )]k\rho (1\rho )|\mathrm{\Lambda }|=[k2(2k1)\rho (1\rho )]M.$$
(27)   
We have found this expression first for $`d`$-dimensional hypercubic lattices, and checked it later to hold also for the triangular, honeycomb and Kagomé lattices. This is somewhat surprising because our derivation below needs knowledge of the rather different local neighborhoods up to next-nearest neighbors. The only common feature of all these lattices seems to be that all sites are symmetry-related and thus equivalent. Therefore, it should be possible to prove Eq. (27) on the basis of vertex-transitivity alone.

Equation (27) is obtained by using the grand-canonical probabilities $`P_{\mathrm{\Lambda },\rho }`$, i.e., by filling the sites of $`\mathrm{\Lambda }`$ independently and with equal probability $`\rho `$. We expect smaller order corrections to appear if the canonical distribution is used. So in this section $`X`$ is a random subset of $`\mathrm{\Lambda }`$ whose probability to be selected is $`\rho ^{|X|}(1\rho )^{|\mathrm{\Lambda }||X|}`$, and $`n_x=n_x(X)`$ is a random variable taking the value 1 if $`x`$ is in $`X`$ and 0 otherwise. Then all $`n_x`$ are independent and take 1 with probability $`\rho `$ and 0 with probability $`1\rho `$. We define

$$f_x=n_x\underset{yx}{}(1n_y)$$
(28)   
where $`x`$ denotes the set of neighbors of $`x`$ in $`\mathrm{\Lambda }`$. Clearly,

$$\overline{f_x}=k\rho (1\rho ).$$
(29)   
The boundary length of $`X`$ is obtained as

$$|X|=\underset{x\mathrm{\Lambda }}{}f_x(X).$$
(30)   
Thus the mean value of (30) is

$$M=\underset{x\mathrm{\Lambda }}{}\overline{f_x}=k\rho (1\rho )|\mathrm{\Lambda }|$$
(31)   
as found earlier.

Let $`d(x,y)`$ denote the graph distance of $`x`$ and $`y`$ in $`\mathrm{\Lambda }`$, i.e., the length of the shortest walk between them. Since $`f_x`$ and $`f_y`$ are independent if $`d(x,y)>2`$, we find

$$D^2=\underset{x,y}{}\overline{(f_x\overline{f_x})(f_y\overline{f_y})}=\underset{x,y:d(x,y)2}{}r(x,y)$$
(32)   
$$r(x,y)=\overline{f_xf_y}\overline{f_x}^2.$$
(33)   
The computation of the different terms is straightforward by observing that $`n_x^2=n_x`$, $`(1n_x)^2=1n_x`$ and $`n_x(1n_x)=0`$. The contribution of the diagonal terms $`x=y`$ is the same for any $`k`$-regular lattice. Namely,

$$\overline{f_x^2}=k\rho (1\rho )+k(k1)\rho (1\rho )^2$$
(34)   
$$\underset{x\mathrm{\Lambda }}{}r(x,x)=|\mathrm{\Lambda }|r(x,x)=[k(2k1)\rho +k\rho ^2]M.$$
(35)   
The contribution of nearest neighbor pairs depends on the number of triangles containing a given edge. If there are $`\mathrm{}`$ such triangles then

$$\overline{f_xf_y}=\rho ^2[\mathrm{}(1\rho )+[(k1)^2\mathrm{}](1\rho )^2]$$
(36)   
$$\underset{x,y:d(x,y)=1}{}r(x,y)=k|\mathrm{\Lambda }|r(x,y)=\rho [\mathrm{}\rho (2k1)(1\rho )]M.$$
(37)   
If $`x`$ and $`y`$ are next-nearest neighbors to each other, they may have $`m`$ common nearest neighbors. Then

$`\overline{f_xf_y}`$ $`=`$ $`\rho ^2[m(1\rho )+(k^2m)(1\rho )^2]`$   
$`r(x,y)`$ $`=`$ $`m\rho ^3(1\rho ).`$ (38)   
In $`d`$-dimensional hypercubic lattices ($`k=2d`$) there are next-nearest neighbor pairs with $`m=1`$ and $`m=2`$. Their contribution to $`D^2`$ is

$`{\displaystyle \underset{x,y:d(x,y)=2}{}}r(x,y)`$ $`=`$ $`k|\mathrm{\Lambda }|r(x,y)_{m=1}+4\left({\displaystyle \genfrac{}{}{0pt}{}{d}{2}}\right)|\mathrm{\Lambda }|r(x,y)_{m=2}`$ (39)   
$`=`$ $`\rho ^2M+(k2)\rho ^2M=(k1)\rho ^2M.`$   
For the triangular lattice ($`k=6`$)

$$\underset{x,y:d(x,y)=2}{}r(x,y)=k|\mathrm{\Lambda }|[r(x,y)_{m=1}+r(x,y)_{m=2}]=(k3)\rho ^2M.$$
(40)   
In the honeycomb lattice ($`k=3`$) each site has 6 next-nearest neighbors, all of the type $`m=1`$. So

$$\underset{x,y:d(x,y)=2}{}r(x,y)=2k|\mathrm{\Lambda }|r(x,y)_{m=1}=(k1)\rho ^2M.$$
(41)   
In the Kagomé lattice ($`k=4`$) there are 8 next-nearest neighbors with $`m=1`$:

$$\underset{x,y:d(x,y)=2}{}r(x,y)=2k|\mathrm{\Lambda }|r(x,y)_{m=1}=(k2)\rho ^2M.$$
(42)   
Finally, we obtain $`D^2`$ by adding (35), (37) with $`\mathrm{}=0`$ and (39) for hypercubic lattices, (35), (37) with $`\mathrm{}=2`$ and (40) for the triangular lattice, (35), (37) with $`\mathrm{}=0`$ and (41) for the honeycomb lattice and (35), (37) with $`\mathrm{}=1`$ and (42) for the Kagomé lattice. All yield (27).

## 4 Energy estimates

By inspecting Eq. (25) it is clear that $`v(n)`$ has to be chosen in such a way that it shifts the expectation value of the Gaussian upwards. The appropriate choice can be either a rapidly – at least exponentially – increasing function or a function concentrated on values well above $`M`$. If in the numerator we had $`v(n)`$ also at the place of $`v_n`$, $`B(v)=n_{\mathrm{max}}`$ could be reached. However, typically $`v_n/v(n)<1`$ and decreases as we modify $`v`$ by putting increasing weights to larger boundary lengths. This limits the maximum of $`B(v)`$. We repeat here the DLS-bound , mentioned in Section 2,

$$\lambda _1M\left(1+\frac{1}{k}\right)(\rho =\frac{1}{2})$$
(43)   
and the trivial upper bound (15),

$$\lambda _1M\frac{1}{1\rho }.$$
(44)   
The right members of these inequalities are upper bounds to $`B(v)`$ as well.

Step functions

Let $`v(n)=1`$ if $`n_1nn_2`$ and 0 otherwise. Suppose that $`n_1M=m>c_1|\mathrm{\Lambda }|`$ and $`n_2n_1>c_2\sqrt{|\mathrm{\Lambda }|}`$ where $`c_1`$ and $`c_2`$ are positive constants. We have

$$B(v)=\left(\underset{n=n_1}{\overset{n_2}{}}e^{\frac{(nM)^2}{2D^2}}\right)^1\underset{n=n_1}{\overset{n_2}{}}e^{\frac{(nM)^2}{2D^2}}\left[n\frac{\underset{i}{}e^{\frac{nM}{D^2}i}v(ni)}{_ie^{\frac{nM}{D^2}i}}\right].$$
(45)   
Because of the lower cutoff, the Gaussian is sharply concentrated on $`n_1`$. Except for $`v(ni)`$, we can therefore replace $`n`$ by $`n_1`$ in the square bracket. Also, with the above choice of $`n_2`$, the dependence on $`n_2`$ is negligible. We find after some manipulations that apart from smaller order corrections

$$B(v)=n_1\left[1\left(\underset{i=2k+2}{\overset{2k2}{}}e^{\frac{m}{D^2}i}\right)^1\underset{i=1}{\overset{2k2}{}}e^{\frac{m}{D^2}i}\frac{\underset{n=n_1}{\overset{n_1+i1}{}}e^{\frac{(nM)^2}{2D^2}}}{_{nn_1}e^{\frac{(nM)^2}{2D^2}}}\right].$$
(46)   
We simplify the fraction with $`e^{m^2/2D^2}`$ and evaluate it to find $`1e^{mi/D^2}`$ plus a correction which vanishes as $`|\mathrm{\Lambda }|\mathrm{}`$. Inserting this into Eq. (46) we obtain

$`B(v)`$ $`=`$ $`M\left(1+{\displaystyle \frac{m}{M}}\right){\displaystyle \frac{2k2+\underset{i=0}{\overset{2k2}{}}e^{\frac{m}{D^2}i}}{_{i=2k+2}^{2k2}e^{\frac{m}{D^2}i}}}`$ (47)   
$`=`$ $`M\left\{1+\left[{\displaystyle \frac{D^2}{M}}{\displaystyle \frac{(k1)(2k1)}{4k3}}\right]{\displaystyle \frac{m}{D^2}}+O\left(\left({\displaystyle \frac{m}{D^2}}\right)^2\right)\right\}.`$   
From Eq. (27) we see that $`D^2/M`$ tends to $`k`$ as $`\rho `$ goes to zero. Hence, the coefficient of $`m/D^2`$ is positive if $`\rho `$ is small enough and we have the best (improved) bound for $`n_1`$ close to $`n_{\mathrm{max}}`$. Then $`m/D^2\rho /k`$ and we are not in conflict with (44). On the other hand, as $`\rho `$ goes to $`\frac{1}{2}`$, $`D^2/M`$ tends to $`\frac{1}{2}`$ and the coefficient of $`m/D^2`$ in (47) becomes negative. Thus, for any $`k2`$ there is a neighborhood of $`\rho =\frac{1}{2}`$ in which the maximum of $`B(v)`$ is reached for $`m=0`$, that is, we get no improved bound.

Exponential trial functions

Let $`v(n)=e^{xn}`$ with an $`x>0`$. Suppose first $`x<x_{\mathrm{max}}`$ where

$$x_{\mathrm{max}}=\frac{n_{\mathrm{max}}M}{2D^2}.$$
(48)   
Plugging $`v`$ into Eq. (25), due to the factor $`e^{2xn}`$ the expectation value of the Gaussian is shifted from $`M`$ to $`M+2xD^2<n_{\mathrm{max}}`$ both in the numerator and in the denominator. The new expectation value can replace $`n`$ in its other occurrences in the numerator. Apart from smaller order corrections we get

$$B(v)=M\left(1+2xD^2/M\right)F_k(x)MG_{k,\rho }(x)(x<x_{\mathrm{max}}),$$
(49)   
$$F_k(x)=\frac{\underset{i=2k+2}{\overset{2k2}{}}e^{xi}}{_{i=2k+2}^{2k2}e^{2xi}}=\frac{\mathrm{cosh}(2k\frac{1}{2})x\mathrm{cosh}(2k\frac{5}{2})x}{\mathrm{cosh}(4k\frac{5}{2})x\mathrm{cosh}(4k\frac{7}{2})x}.$$
(50)   
Expanding $`F_k(x)`$ around 0 we find

$$G_{k,\rho }(x)=\left(1+2xD^2/M\right)[1(k1)(2k1)x^2+O(x^4)]>1$$
(51)   
if $`x`$ is small enough, so we have an improved bound for any $`k`$ and $`\rho `$.

If $`x>x_{\mathrm{max}}`$, the maximum of the shifted Gaussian is outside the range of summation. Therefore the distribution is sharply concentrated on a small neighborhood of $`n_{\mathrm{max}}`$, and the computation of $`B(v)`$ can be done in analogy to the case of the step function. The result is

$$B(v)=\frac{n_{\mathrm{max}}}{_{i=2k+2}^{2k2}e^{2x_{\mathrm{max}}i}}\left[\underset{i=0}{\overset{2k2}{}}e^{(2x_{\mathrm{max}}x)i}+\underset{i=1}{\overset{2k2}{}}e^{xi}\right](x>x_{\mathrm{max}}).$$
(52)   
Notice that $`B(v)`$ varies continuously with $`x`$ and both (49) and (52) yield

$$B(v)=n_{\mathrm{max}}F_k\left(x_{\mathrm{max}}\right)\mathrm{if}x=x_{\mathrm{max}}.$$
(53)   
In order to obtain the best bound, $`B(v)`$ has to be maximized with respect to $`x`$. $`G_{k,\rho }(x)`$ has a unique maximum at some positive $`x`$. This can be the location of the maximum of $`B(v)`$ only when it is smaller than $`x_{\mathrm{max}}`$, and then the optimal $`B(v)`$ is computed with it from (49). Table 1 shows that this is the case for all $`k`$ if $`\rho `$ is not too close to 0. However, for small enough densities the maximum of $`G_{k,\rho }`$ is attained at an $`x`$ above $`x_{\mathrm{max}}`$. For these densities the first expression (49) increases with $`x`$ up to $`x_{\mathrm{max}}`$ while the second expression (52) decreases with $`x`$ for all densities. So the highest bound is provided by (53).

The conclusion is that for all densities we can obtain the best bound by maximizing (49) with respect to $`x`$ under the condition that $`xx_{\mathrm{max}}`$.

Numerical results on the best estimates are summarized in the tables below. They apply to bipartite lattices whenever $`D^2`$ is given by Eq. (27). According to Table 1, they are valid also to the Kagomé and triangular lattices, because the largest density for which the maximizing $`x`$ is $`x_{\mathrm{max}}`$ remains below $`\frac{1}{3}`$ for every $`k`$.

|  | k | | | | |
| --- | --- | --- | --- | --- | --- |
| $`\rho `$ | 2 | 3 | 4 | 6 | 8 |
| 0.1 | 0.038 | 0.026 | 0.020 | 0.014 | 0.010 |
| 0.2 | 0.120 | 0.089 | 0.071 | 0.040 | 0.027 |
| 0.3 | 0.206 | 0.082 | 0.047 | 0.024 | 0.016 |
| 0.4 | 0.167 | 0.057 | 0.030 | 0.013 | 0.008 |
| 0.5 | 0.152 | 0.048 | 0.023 | 0.009 | 0.005 |

Table 1. Rates of increase, $`x`$, of the exponential trial functions, maximizing $`G_{k,\rho }`$. Numbers with an asterix correspond to $`x_{\mathrm{max}}`$.

|  | k | | | | |
| --- | --- | --- | --- | --- | --- |
| $`\rho `$ | 2 | 3 | 4 | 6 | 8 |
| 0.1 | 1.106 | 1.103 | 1.102 | 1.099 | 1.098 |
| 0.2 | 1.199 | 1.157 | 1.128 | 1.100 | 1.088 |
| 0.3 | 1.155 | 1.075 | 1.051 | 1.033 | 1.027 |
| 0.4 | 1.094 | 1.035 | 1.019 | 1.009 | 1.006 |
| 0.5 | 1.077 | 1.024 | 1.012 | 1.004 | 1.002 |

Table 2. Maximal values of $`G_{k,\rho }(x)`$, computed with the entries of Table 1.

|  | k | | | | |
| --- | --- | --- | --- | --- | --- |
| $`\rho `$ | 2 | 3 | 4 | 6 | 8 |
| 0.1 | 0.199 | 0.298 | 0.397 | 0.594 | 0.791 |
| 0.2 | 0.383 | 0.555 | 0.722 | 1.056 | 1.393 |
| 0.3 | 0.485 | 0.677 | 0.882 | 1.302 | 1.725 |
| 0.4 | 0.525 | 0.745 | 0.978 | 1.453 | 1.932 |
| 0.5 | 0.538 | 0.768 | 1.012 | 1.507 | 2.005 |

Table 3. Best estimates of $`|E_0|/|\mathrm{\Lambda }|`$ obtained by multiplying $`k\rho (1\rho )`$ with the corresponding entry of Table 2.

| k | 2 | 3 | 4 | 6 | 8 |
| --- | --- | --- | --- | --- | --- |
| $`\rho `$ | 0.23 | 0.19 | 0.17 | 0.16 | 0.15 |
| $`x`$ | 0.158 | 0.080 | 0.051 | 0.031 | 0.021 |
| $`G_{k,\rho }(x)`$ | 1.206 | 1.159 | 1.143 | 1.129 | 1.123 |

Table 4. The maximum of $`G_{k,\rho }(x)`$ as a function of $`\rho `$ and $`x[0,x_{\mathrm{max}}]`$, together with the maximizing $`\rho `$ and $`x`$. In each case the latter is $`x_{\mathrm{max}}`$.

Gaussian trial functions

Let $`v(n)=e^{(nM_1)^2/4D_1^2}`$ where $`M_1>M`$ and $`D_1^2\mathrm{}`$ with increasing $`|\mathrm{\Lambda }|`$ but otherwise are free parameters. The product of the two Gaussians gives rise to a new Gaussian with mean value and variance

$$M_2=\frac{MD_1^2+M_1D^2}{D_1^2+D^2}D_2^2=\frac{D_1^2D^2}{D_1^2+D^2}.$$
(54)   
A straightforward computation yields the same $`B(v)`$, Eqs. (49), (52) and (53) as for the exponential trial function with

$$x=\frac{1}{2}\frac{M_1M}{D_1^2+D^2}.$$
(55)   
So Gaussian and exponential trial functions provide the same bound, at least in leading order. Moreover, for any fixed $`x>0`$ there is a one-parameter family of Gaussians giving the same result. An interesting choice is $`D_1^2=o(|\mathrm{\Lambda }|)`$, e.g. $`D_1^2=\mathrm{const}\sqrt{|\mathrm{\Lambda }|}`$. In this case $`M_2/M_1=1`$ and $`D_2^2/D_1^2=1`$ asymptotically.

## 5 Summary and concluding remarks

We have presented variational estimates of the ground state energy of a gas of hard-core bosons on regular lattices. The wave functions we have used depended only on the size of the boundary of the $`N`$-point configurations. Therefore, the estimates could be based on a large deviation principle governing the distribution of the boundary sizes. The corresponding rate function is related to the entropy of the (ferromagnetic) Ising model on the same lattice. We have derived a formula for the mean-square deviation of the boundary sizes and applied it in a quadratic approximation of the rate function.

The best estimates we have found are of the form (49) where for given $`k`$ and $`\rho `$ the parameter $`x`$ is to be determined numerically so as to maximize $`B(v)`$. This has been illustrated in Tables 1-3. The maximum is realized either by an exponential function or by a one-parameter family of Gaussian functions.

Extension to the grand-canonical ensemble is straightforward. Adding $`\mu _{x\mathrm{\Lambda }}n_x`$ to the Hamiltonian resumes in adding $`\mu \rho |\mathrm{\Lambda }|`$ to $`B(v)`$. After maximizing with respect to $`x`$ we have to maximize with respect to $`\rho `$. More precisely, if $`b(\rho )`$ is the maximum of $`B(v)/|\mathrm{\Lambda }|`$ with respect to $`x`$, first we extend $`b(\rho )`$ symmetrically to $`\rho >\frac{1}{2}`$ and then determine $`\widehat{b}(\mu )=\mathrm{max}_\rho [b(\rho )+\mu \rho ]`$ to get an estimate of the ground state energy per site in the full bosonic Fock space. It is more $`\widehat{b}(\mu )`$ than $`b(\rho )`$ which is useful for the equivalent spin model whose Hamiltonian now contains an external magnetic field in the $`Z`$-direction. The function $`\rho (\mu )`$ that we can find in determining $`\widehat{b}(\mu )`$ is only an approximation of the true relationship which exists between the chemical potential and the density in the ground state. Table 3 suggests that $`\rho (0)=\frac{1}{2}`$. This rigorously holds true for any positive temperature; even we have a much stronger ‘uniform density theorem’, as in the Hubbard model : Because of the particle-hole symmetry,

$$\frac{\mathrm{Tr}n_xe^{\beta H_{\mu =0}}}{\mathrm{Tr}e^{\beta H_{\mu =0}}}=\frac{\mathrm{Tr}(1n_x)e^{\beta H_{\mu =0}}}{\mathrm{Tr}e^{\beta H_{\mu =0}}}=\frac{1}{2}$$
(56)   
is valid in any $`\mathrm{\Lambda }`$ with free boundary condition.

In our estimates the lattice structure appears only through $`k`$, the coordination number. Therefore we obtain the same result for the Kagomé lattice as for the square lattice, and for the triangular lattice as for the simple cubic lattice. Although the energy corresponds to a nearest neighbor correlation, the details of the lattice geometry certainly influence the $`\rho `$-dependence of the exact ground state energy per site. The mark of the lattice could be recovered by using the exact Ising entropies instead of their quadratic approximants.