text:
# Antiferromagnetism in the Exact Ground State of the Half Filled Hubbard Model on the Complete-Bipartite Graph

## I Introduction

The Hubbard Hamiltonian is one of the most popular models to describe strongly correlated electron systems. In spite of its simple definition, it can be exactly solved only in few cases. An example is the milestone solution of the Hubbard ring by means of the Bethe-Ansatz extended to fermions. Other exact solutions were worked out for the one-dimensional Hubbard ring in the presence of an external magnetic field, for the one-dimensional Hubbard chain with open boundary conditions and thereafter for the $`SU(N)`$ one-dimensional Hubbard ring. However, when the space dimensionality becomes bigger than 1, few exact results are available and usually they concern the ground state properties. Among them, we mention the Lieb theorem on the ground state spin degeneracy that will be explicitly used in this work. Exact ground state wave functions are even more infrequent. To the best of our knowledge the only non-trivial results are the ferromagnetic ground-state solutions devised by Nagaoka (infinite repulsion and one hole over the half filled system), by Mielke and Tasaki (for which the kinetic energy spectrum is macroscopically degenerate or at least very flat) and by Wang (infinite-range hopping).

In the Hubbard Hamiltonian with infinite-range hopping a particle in any site can hop to any other site of the system; the associated graph is said to be complete. This model was numerically studied by Patterson on small clusters in 1972 and solved in the thermodynamic limit by van Dongen and Vollhardt only at the end of the eighties. Much more effort was needed to find the exact ground state(s) in the finite-size system. Verges et al. managed to accomplish this task for arbitrary numbers of particles and sites in the limit of infinite on-site repulsion $`U`$ by exploiting a scheme proposed by Brandt and Giesekus . Two years later Wang constructed explicitly the ground states of the system for arbitrary $`U`$ and number of particles above half filling. In the case of one particle added over the half-filled system the ferromagnetic ground-state solution follows as a special (and the easiest) case of general results by Mielke and Tasaki.

In this paper we find the exact ground state wave function of the half-filled Hubbard model on the Complete-Bipartite Graph (CBG) for arbitrary but finite $`U`$. The CBG $`\mathrm{\Lambda }=𝒜`$ has bonds connecting any element of $`𝒜`$ with all the elements of $``$ and can be considered as the natural further step (with respect to the complete graph previously described) towards the standard Hubbard Hamiltonian (defined on the hypercubic lattice and hopping between nearest-neighbours sites). Even if the CBG is still somewhat unrealistic, we feel that exact solutions should always be welcome, expecially because they lend themselves to be generalized. Furthermore, our solution is the first example of antiferromagnetic ground state in a model of itinerant electrons, in contrast with the ferromagnetic solutions mentioned above; it may provide useful hints about antiferromagnetism outside of the strong coupling regime (where the Hubbard model can be mapped onto the Heisenberg model).

The paper has been organized as follows. In Section II we define the model and discuss the physics of the non-interacting ($`U=0`$) Hamiltonian together with few relevant examples of finite-size realizations. In Section III we study the thermodynamic limit of a class of Hubbard-like models, which includes the Hubbard Hamiltonian on the complete graph and on the complete bipartite graph, having a non-extensive number of isolated one-particle energy levels plus a single level whose degeneracy is proportional to the size of the system. Following the reasoning of van Dongen and Vollhardt, we show that the kinetic term is totally decoupled from the interaction term and that the system behaves as a paramagnet: the spin-spin correlation length is zero. These results are a consequence of the trivial behaviour of such models any time the thermodynamic limit is taken first. On the contrary, much more difficult is to find exact properties in finite-size systems and different answers may be obtained if the thermodynamic limit is taken only at the end of the calculations (as we shall show in Section VI). In Section IV we introduce the essential tools to face the problem. Let $`N=N_A=N_B`$ be the number of sites of each sublattice $`𝒜`$ and $``$. First, we find a ($`2N2`$)-body determinantal eigenstate $`|\mathrm{\Phi }_{AF}`$ of the Hamiltonian with vanishing double occupation; then, we demonstrate that it is a key tool to build the ground state. We show that mapping the $`𝒜`$-sites onto the $``$-sites and viceversa, $`|\mathrm{\Phi }_{AF}`$ retains its form except for a spin-flip; we shall call this property the antiferromagnetic property. Several analogies with the properties of the half-filled standard $`N\times N`$ Hubbard model are pointed out at this stage. Here we also deal with the spin projection of $`|\mathrm{\Phi }_{AF}`$ onto the singlet and the triplet subspace and useful identities between the two spin-projected states are obtained. Then, in Section V we propose an Ansatz for the ground state wave function at half filling containing the singlet and triplet projections of $`|\mathrm{\Phi }_{AF}`$. We set up the Schrödinger equation and by exploiting the antiferromagnetic property we close the equations and get 3 exact eigenstates. The ground-state uniqueness proved by Lieb is used to show that the lowest-energy state of our Ansatz actually corresponds to the ground state of the system. Remarkably, the ground state energy is negative for any value of the repulsion $`U`$; qualitatively, we may say that the particles manage to avoid the double occupation very effectively. In Section VI we study the ground state energy $`E_0`$ as a function of $`U`$ and of the volume of the system and we discuss the implications of exchanging the thermodynamic limit with the limit of zero temperture. We find that $`E_0`$ is a monotonically increasing function of $`U`$ and $`N`$ due to the existence of non-trivial correlations even for large $`N`$. The nature of these correlations is investigated by computing the expectation value of the repulsion. We show that for any finite $`N`$ there is a critical value of $`U`$ yielding maximum repulsion. From the exact spin-spin correlation function we find that the ground state average of the staggered magnetization operator squared is exstensive. Hence, the ground state is antiferromagnetically ordered; we underline that this holds not only at strong coupling, but for any value of $`U`$. Finally, a summary of the main results and conclusions are drawn in Section VII.

## II Definition of the Model

Let $`\mathrm{\Lambda }=𝒜`$ with $`𝒜=\mathrm{}`$ be a collection of sites and

$$|𝒜|=N_A,||=N_B|\mathrm{\Lambda }|=N_A+N_B.$$
(1)
Here and in the following we shall denote by $`|𝒮|`$ the number of elements in the set $`𝒮`$. We consider the Hubbard Hamiltonian

$$H_{\mathrm{H}ubbard}=H_0+W,H_0=\underset{x,y\mathrm{\Lambda }}{}\underset{\sigma }{}t_{x,y}c_{x\sigma }^{}c_{y\sigma },W=U\underset{x\mathrm{\Lambda }}{}\widehat{n}_x\widehat{n}_x$$
(2)
where $`c_{x\sigma }`$ ($`c_{x\sigma }^{}`$) is the annihilation (creation) operator of a particle at site $`x`$ with spin $`\sigma =,`$ and $`\widehat{n}_{x\sigma }=c_{x\sigma }^{}c_{x\sigma }`$ is the corresponding particle number operator. The hopping matrix $`T`$ with elements $`T_{x,y}=t_{x,y}=t_{y,x}`$ is a real-symmetric matrix while $`U`$ is a positive constant. If

$$t_{x,y}=\{\begin{array}{cc}\stackrel{~}{t}\hfill & \mathrm{f}orx𝒜(x)\mathrm{a}ndy(y𝒜)\hfill \\ 0\hfill & \mathrm{o}therwise\hfill \end{array}$$
(3)
the graph $`\mathrm{\Lambda }`$ is said to be complete bipartite and we call the Hamiltonian in Eq.(2) the CBG-Hubbard Hamiltonian. The model is invariant under an arbitrary permutation of the $`𝒜`$-sites and/or of the $``$-sites. Therefore, the symmetry group includes $`S_{N_A}S_{N_B}`$, $`S_N`$ being the set of permutations of $`N`$ objects. As usual, the full group is much bigger: the presence of spin and pseudospin symmetries leads to an $`SO(4)`$ internal symmetry Group and in the case $`N_A=N_B`$ there is a $`Z_2`$ symmetry because of the $`𝒜`$ exchange.

In Fig.(1) we have drawn a few examples of finite-size systems. For $`N=1`$ and $`N=2`$ the model is equivalent to a one dimensional ring of length $`L=2,\mathrm{\hspace{0.33em}4}`$ respectively. For $`N=3`$ we have what can be understood as a prototype, (1,1) nanotube model, the one of smallest length $`L=1`$, with periodic boundary conditions. For general $`N`$, one can conceive a gedanken device, like the one illustrated pictorially for $`N=6`$ in Fig.(1.$`d`$). (The whole device should actually have the topology of a torus, and the two horizontal faces should coincide, but this is not shown in the Figure for the sake of clarity.) The $`N`$ vertical lines represent a realization of the $`𝒜`$ sublattice while the $``$ sublattice is mimicked by the central object. The radial tracks in the Figure represent conducting paths linking each $`𝒜`$ site to each $``$ site according to the topology of our model. The $`𝒜`$ sites are represented by one-dimensional electron-in-a-box systems of length $`L`$. Each $`𝒜`$ site has the one-body energy spectrum $`\epsilon _{\mathrm{f}ree}n^2/mL^2`$ where $`n`$ is an integer and $`m`$ is the electron mass. We assume that $`L`$ is so small that the excited states are at high energy and can be disregarded in the low-energy sector; this requires $`mL^21/U`$. Here, of course, $`U`$ denotes the Coulomb self-energy of a box with two electrons. The $``$ sites are hosted by $`N`$ quantum dots, that can be represented by $`\delta `$-function-like attractive potential wells of depth $`V`$; if $`|V|`$ is large and the radius of the dots is $`L`$, these are practically independent of each other. The one-body energy of each $``$ site is $`\epsilon _{\mathrm{b}ound}mV^2`$; the unbound states can be neglected if $`mV^2U`$. We are assuming for simplicity that the $`U`$ of the $``$ sites is the same as for the $`𝒜`$ ones. Turning on a constant potential $`V_0`$ on the central object, one can arrange that the Hubbard Hamiltonian on the CBG $`H_{\mathrm{H}ubbard}`$ gives a good description of the system, with filling $`2`$ per site.

It is convenient to label each site with an integer in such a way that $`x=1,\mathrm{},N_A`$ corresponds to sites in $`𝒜`$ and $`x=N_A+1,\mathrm{},N_A+N_B`$ corresponds to sites in $``$. In the special case $`N_A=N_B=N`$ the hopping matrix can be written as

$$T=\stackrel{~}{t}\left(\begin{array}{cc}Z\hfill & J\hfill \\ J\hfill & Z\hfill \end{array}\right)$$
(4)
where $`Z`$ is the $`N\times N`$ null matrix and $`J`$ is the $`N\times N`$ matrix whose generic element $`J_{x,y}=1`$. The one-body spectrum has three different eigenvalues

$`\epsilon _g`$ $`=`$ $`N\stackrel{~}{t}t`$ (5)
$`\epsilon _0`$ $`=`$ $`0`$ (6)
$`\epsilon _{\overline{g}}`$ $`=`$ $`N\stackrel{~}{t}t`$ (7)
with degeneracy $`d_g=1`$, $`d_0=2N2`$ and $`d_{\overline{g}}=1`$ respectively. We use the convention $`t>0`$ so that $`\epsilon _g`$ is the lowest level and we shall call $`𝒮_{hf}`$ the set of zero-energy one-body eigenstates.

The orthogonal matrix $`O`$ that diagonalizes $`T`$ can be written in the form

$$O=\left(\begin{array}{cc}\begin{array}{c}\frac{1}{\sqrt{2N}}\frac{1}{\sqrt{2N}}\mathrm{}\\ R\\ Z_r\\ \frac{1}{\sqrt{2N}}\frac{1}{\sqrt{2N}}\mathrm{}\end{array}\hfill & \begin{array}{c}\frac{1}{\sqrt{2N}}\frac{1}{\sqrt{2N}}\mathrm{}\\ Z_r\\ R\\ \frac{1}{\sqrt{2N}}\frac{1}{\sqrt{2N}}\mathrm{}\end{array}\hfill \end{array}\right)$$
(8)
where $`Z_r`$ is the $`(N1)\times N`$ rectangular null matrix, while $`R`$ is an $`(N1)\times N`$ rectangular matrix whose rows are orthonormal vectors which are orthogonal to the $`N`$dimensional vector $`(1,1,\mathrm{},1,1)`$. The zero-energy one-body orbitals can be visualized by a simple argument. Consider any pair $`x,y`$, with $`xy,`$ of sites belonging to the same sublattice, say $`𝒜`$, and a wavefunction $`\psi _{x,y}`$ taking the values 1 and -1 on the pair and 0 elsewhere in $`𝒜`$ and in $``$. It is evident that $`\psi _{x,y}`$ belongs to $`𝒮_{hf}`$. Operating on $`\psi _{x,y}`$ by the permutations of $`S_N`$ we can generate a (non-orthogonal) basis of $`N1`$ eigenfunctions*** There is a simple trick to see that the number of independent wavefunctions obtained in this way is $`N1`$. Let $`PS_N`$ be a permutation such that $`P(x)=x`$ and $`P(y)x`$. We get the full basis since $`P(y)`$ takes $`N1`$ values with independent $`\psi _{x,P(y)}`$. vanishing in $``$; further, by means of the $`Z_2`$ symmetry, we obtain the remaining orbitals of $`𝒮_{hf}`$, which vanish on $`𝒜`$. This exercise shows that the group considered above justifies the ($`2N2`$)-fold degeneracy of the one-body spectrum.

Having obtained the one-body basis, one can form a suitable $`R`$ matrix by orthogonalization. Let us define the transformed operators:

$`g{\displaystyle \underset{x=1}{\overset{2N}{}}}O_{1,x}c_x={\displaystyle \frac{1}{\sqrt{2N}}}{\displaystyle \underset{x=1}{\overset{N}{}}}(c_xc_{x+N})`$ (9)
$`a_i{\displaystyle \underset{x=1}{\overset{2N}{}}}O_{i+1,x}c_x={\displaystyle \underset{x=1}{\overset{N}{}}}R_{i,x}c_x,i=1,\mathrm{},N1`$ (10)
$`b_i{\displaystyle \underset{x=1}{\overset{2N}{}}}O_{i+N,x}c_x={\displaystyle \underset{x=1}{\overset{N}{}}}R_{i,x}c_{x+N},i=1,\mathrm{},N1`$ (11)
$`\overline{g}{\displaystyle \underset{x=1}{\overset{2N}{}}}O_{2N,x}c_x={\displaystyle \frac{1}{\sqrt{2N}}}{\displaystyle \underset{x=1}{\overset{N}{}}}(c_x+c_{x+N}).`$ (12)
The inverse transformation reads

$`c_x={\displaystyle \frac{1}{\sqrt{2N}}}(g+\overline{g})+{\displaystyle \underset{i}{}}R_{i,x}a_i,x𝒜`$ (13)
$`c_x={\displaystyle \frac{1}{\sqrt{2N}}}(g+\overline{g})+{\displaystyle \underset{i}{}}R_{i,xN}b_i,x`$ (14)
and the kinetic term $`H_0`$ becomes

$$H_0=t\underset{\sigma }{}(g_\sigma ^{}g_\sigma \overline{g}_\sigma ^{}\overline{g}_\sigma ).$$
(15)
Hence, if we do not rescale the hopping constant, the average kinetic energy remains an extensive quantity proportional to $`N\stackrel{~}{t}`$: the kinetic energy of the two particles in the lowest level coincides with the kinetic energy of the whole system. On the other hand, if $`\stackrel{~}{t}1/N`$ the two energy gaps $`\epsilon _0\epsilon _g`$ and $`\epsilon _{\overline{g}}\epsilon _0`$ remain constant as $`N`$ increases. As we shall show in the next Section, the model can be exactly solved in the thermodynamic limit whatever is the dependence of $`\epsilon _g`$ and $`\epsilon _{\overline{g}}`$ on $`|\mathrm{\Lambda }|`$.

## III Thermodynamic Limit

The Hubbard Hamiltonian on the CBG (as on the complete graph) belongs to a class of Hubbard-like models that can be exactly solved in the thermodynamic limit. Here we generalize the results by van Dongen and Vollhardt to the case of a non-extensive number of non-vanishing one-body eigenenergies and in the presence of a local external magnetic field coupled to the spin of the particles. Even if the reasoning is similar to that of Ref., a detailed presentation of the main results is needed to clarify more subtle points and make our work self-contained. Furthermore, some of the results are absent in the original paper and we believe that their derivation will facilitate the reader in the comprehension of what follows.

Let us consider the Hamiltonian

$$H=H_{\mathrm{H}ubbard}\mu \underset{x\mathrm{\Lambda }}{}\underset{\sigma }{}\widehat{n}_{x\sigma }\underset{x\mathrm{\Lambda }}{}h_x(\widehat{n}_x\widehat{n}_x)H_{\mathrm{H}ubbard}+H_\mu +H_h$$
(16)
where $`H_{\mathrm{H}ubbard}=H_0+W`$ is defined in Eq.(2), $`\mu `$ is the chemical potential and $`H_h`$ represents the coupling with an external magnetic field. Let us assume that only a finite number f of eigenvalues of the hopping matrix $`T`$ are non-vanishing, $`\epsilon _i0`$ for $`i=1,\mathrm{},`$f, while the remaining $`|\mathrm{\Lambda }|`$f are identically zero. Denoting with $`\varphi _i`$ and $`\varphi _i^{}`$ the annihilation and creation operators that diagonalize $`H_0`$ we have

$$H_0=\underset{i=1}{\overset{\mathrm{f}}{}}\underset{\sigma }{}\epsilon _i\varphi _{i\sigma }^{}\varphi _{i\sigma }$$
(17)
and the CBG corresponds to the case f=2, see Eq.(15). The Gran-Partition Function $`𝒵_\beta [T,U,\{h_x\}]`$ can always be written as

$$𝒵_\beta [T,U,\{h_x\}]=𝒵_\beta [T,0,\{h_x=0\}]e^{𝒲_\beta [T,U,\{h_x\}]},𝒲_\beta [T,U,\{h_x\}]=e^{_0^\beta H_I(\tau )𝑑\tau }_0^c$$
(18)
where $`\beta `$ is the inverse Temperature, $`H_I(\tau )=W(\tau )+H_h(\tau )`$ is the interaction Hamiltonian in the imaginary-time Heisenberg representation and $`\mathrm{}_0^c`$ is the thermal average in terms of $`𝒵_\beta [T,0,\{h_x=0\}]𝒵_\beta ^0[T]`$ where, by the linked-cluster theorem, one has to retain only those contributions represented by connected diagrams. As far as f remains fixed, we can safely substitute $`𝒲_\beta [T,U,\{h_x\}]`$ with $`𝒲_\beta [T=0,U,\{h_x\}]`$ in Eq.(18). Indeed,

$$𝒲_\beta [T,U,\{h_x\}]=𝒲_\beta [T=0,U,\{h_x\}]+O(\mathrm{f}/|\mathrm{\Lambda }|)$$
(19)
as long as the unperturbed one-body Green’s functions $`G_{i\sigma }(\tau )=𝒯\varphi _{i\sigma }(\tau )\varphi _{i\sigma }^{}(0)_0`$, $`i=1,\mathrm{},`$f, do not diverge in the thermodynamic limit; here $`𝒯`$ is the Wick time-ordering operator and $`\mathrm{}_0`$ is the thermal average with $`U=0`$ and $`\{h_x=0\}`$. From the explicit expression of $`G_{i\sigma }(\tau )`$,

$$G_{i\sigma }(\tau )=\{\begin{array}{cc}e^{(\epsilon _i\mu )\tau }(1f_i)\hfill & \tau >0\hfill \\ e^{(\epsilon _i\mu )\tau }f_i\hfill & \tau 0\hfill \end{array},$$
(20)
with $`f_i=[1+e^{\beta (\epsilon _i\mu )}]^1`$ the Fermi distribution function, it follows that $`G_{i\sigma }(\tau )`$ converges to a finite value whenever $`lim_{|\mathrm{\Lambda }|\mathrm{}}\epsilon _i`$ is well defined. One can use the result in Eq.(19) to write the exact Gran-Partition Function in the thermodynamic limit

$`𝒵_\beta [T,U,\{h_x\}]`$ $`=`$ $`{\displaystyle \frac{𝒵_\beta ^0[T]}{𝒵_\beta ^0[T=0]}}𝒵_\beta [T=0,U,\{h_x\}]=`$ (21)
$`=`$ $`{\displaystyle \underset{i=1}{\overset{\mathrm{f}}{}}}\left[{\displaystyle \frac{1+2ze^{\beta \epsilon _i}+z^2e^{2\beta \epsilon _i}}{(1+z)^2}}\right]{\displaystyle \underset{x\mathrm{\Lambda }}{}}\left[1+2z\mathrm{cosh}\beta h_x+z^2e^{\beta U}\right],`$ (22)
where $`z=e^{\beta \mu }`$ is the fugacity, and the thermodynamic potential $`\mathrm{\Omega }_\beta [T,U,\{h_x\}]=\beta ^1\mathrm{ln}𝒵_\beta [T,U,\{h_x\}]`$

$$\mathrm{\Omega }_\beta [T,U,\{h_x\}]=\frac{1}{\beta }\left\{2\mathrm{f}\mathrm{ln}(1+z)\underset{i=1}{\overset{\mathrm{f}}{}}\mathrm{ln}[1+2ze^{\beta \epsilon _i}+z^2e^{2\beta \epsilon _i}]\underset{x\mathrm{\Lambda }}{}\mathrm{ln}\left[1+2z\mathrm{cosh}\beta h_x+z^2e^{\beta U}\right]\right\}.$$
(23)
Eq.(23) reduces to the result of van Dongen and Vollhardt in the case f=1, $`h_x=0`$ $`x\mathrm{\Lambda }`$ and $`\epsilon _1|\mathrm{\Lambda }|<0`$. All the thermodynamic quantities can be derived from Eq.(23) and we defer the interested reader to the original paper by van Dongen and Vollhardt for further details. Here we calculate the magnetization at site $`x`$, $`m_x=\frac{1}{2}(\mathrm{\Omega }_\beta /h_x)=\frac{1}{2}\widehat{n}_x\widehat{n}_xS_x^z`$ and the connected spin-spin correlation function $`S_x^zS_y^z^cS_x^zS_y^zS_x^zS_y^z=(1/4\beta )(^2\mathrm{\Omega }_\beta /h_xh_y)`$, where $`\mathrm{}`$ means the average in the gran-canonical ensemble. Due to the exact decoupling of the kinetic energy, we expect from the begining a paramagnetic behaviour. We have

$$m_x=\frac{1}{2}\frac{\mathrm{\Omega }_\beta }{h_x}=\frac{z\mathrm{sinh}\beta h_x}{1+2z\mathrm{cosh}\beta h_x+z^2e^{\beta U}}$$
(24)
where $`z`$ can be expressed in terms of the number density by means of the relation

$$n=\frac{1}{|\mathrm{\Lambda }|}\frac{\mathrm{\Omega }_\beta }{\mu }=\frac{2z}{|\mathrm{\Lambda }|}\underset{x\mathrm{\Lambda }}{}\frac{\mathrm{cosh}\beta h_x+ze^{\beta U}}{1+2z\mathrm{cosh}\beta h_x+z^2e^{\beta U}}.$$
(25)
Since $`z`$ is finite for any non-zero Temperature, $`m_x0`$ when $`h_x0`$, i.e.

$$\underset{\beta \mathrm{}}{lim}\underset{h_x0}{lim}\underset{|\mathrm{\Lambda }|\mathrm{}}{lim}m_x=0.$$
(26)
On the other hand, $`m_x0`$ if we exchange the last two limits in Eq.(26). Indeed, taking $`h_x=h`$=const$`>0`$ for the sake of clarity, Eq.(25) can be exactly solved for $`z`$:

$$z=\frac{\sqrt{(n1)^2\mathrm{cosh}^2\beta h+n(2n)e^{\beta U}}(1n)\mathrm{cosh}\beta h}{(2n)e^{\beta U}}_\beta \mathrm{}\{\begin{array}{cc}e^{\beta (U+h)}(n1)/(2n)\hfill & n>1\hfill \\ e^{\beta U/2}[n/(2n)]^{1/2}\hfill & n=1\hfill \\ e^{\beta h}n/(1n)\hfill & n<1\hfill \end{array}.$$
(27)
Substitution of these asymptotic behaviours in Eq.(24) yields

$$\underset{h_x0^+}{lim}\underset{\beta \mathrm{}}{lim}\underset{|\mathrm{\Lambda }|\mathrm{}}{lim}m_x=\{\begin{array}{cc}(2n)/2\hfill & n>1\hfill \\ 1/2\hfill & n=1\hfill \\ n/2\hfill & n<1\hfill \end{array}.$$
(28)
The case $`h<0`$ is similar and one can show that Eq.(28) remains true in the limit $`h_x0^{}`$ if $`m_xm_x`$. Therefore, the magnetization does not depend on $`U`$ and the zero-Temperature limit of the model is the same of the paramagnetic Hamiltonian $`H_{\mathrm{p}ara}=H_\mu +H_h`$, where $`H_\mu `$ and $`H_h`$ are defined in Eq.(16). From Eq.(24) we also conclude that

$$S_x^zS_y^z^c\frac{1}{4\beta }\left(\frac{^2\mathrm{\Omega }_\beta }{h_xh_y}\right)=0S_x^zS_y^z=S_x^zS_y^z=m_xm_y,$$
(29)
a result which is independent of the local external field configuration: in the thermodinamic limit two localized spins do not interact and the spin-spin correlation length is zero. Denoting with $`G_{\mathrm{s}pin}(x,y;\beta ,h,|\mathrm{\Lambda }|)`$ the thermal average of $`S_x^zS_y^z`$ at inverse Temperature $`\beta `$, external field $`h`$ and size of the system $`|\mathrm{\Lambda }|`$, from Eq.(26) we get

$$\underset{\beta \mathrm{}}{lim}\underset{h0}{lim}\underset{|\mathrm{\Lambda }|\mathrm{}}{lim}S_x^zS_y^z=\underset{\beta \mathrm{}}{lim}\underset{h0}{lim}\underset{|\mathrm{\Lambda }|\mathrm{}}{lim}G_{\mathrm{s}pin}(x,y;\beta ,h,|\mathrm{\Lambda }|)=0,$$
(30)
that is, the spin-spin correlation function vanishes if we first take the limit $`h0`$ and then the limit of $`\beta \mathrm{}`$. We emphasize that the above results are correct only if we first take the limit $`|\mathrm{\Lambda }|\mathrm{}`$. The thermodynamic limit makes the model trivial and different graph structures, like the complete graph and the complete bipartite graph, all yield the same thermodynamic behaviour. On the contrary, a much more difficult task consists in finding exact results in finite-size systems. To face such problems is not only a mathematical exercise. After that in Section V the exact ground state of the half-filled Hubbard model on the CBG is explicitly written down for arbitrary value of $`U`$ and $`|\mathrm{\Lambda }|`$ and in the absence of an external field $`h`$, we shall see in Section VI that the thermodynamic limit, $`|\mathrm{\Lambda }|\mathrm{}`$, and the limits of zero Temperature and external field, $`\beta \mathrm{},h0`$, do not commute.

## IV The $`W=0`$ States

In this Section we introduce the essential tools to deal with the antiferromagnetic ground-state solution. Let us consider the one-body eigenstate $`a_i^{}|0`$ with vanishing amplitudes on the $``$ sublattice, that is $`0|c_xa_i^{}|0=0`$ if $`x`$. Similarly, $`b_i^{}|0`$ has vanishing amplitude on the $`𝒜`$ sublattice and therefore the $`(2N2)`$-body state

$$|\mathrm{\Phi }_{AF}^{(\sigma )}=a_{1\sigma }^{}\mathrm{}a_{N1\sigma }^{}b_{1\overline{\sigma }}^{}\mathrm{}b_{N1\overline{\sigma }}^{}|0,\overline{\sigma }=\sigma $$
(31)
is an eigenstate of $`H_0`$ and of $`W`$ with vanishing eigenvalue. In the following we shall use the wording $`W=0`$ state to denote any eigenstate of $`H_{\mathrm{H}ubbard}`$ in the kernel of $`W`$. It is worth to observe that by mapping the $`𝒜`$-sites onto the $``$-sites and viceversa, $`|\mathrm{\Phi }_{AF}`$ retains its form except for a spin-flip; we call this property the antiferromagnetic property for obvious reasons.

In the non-interacting ($`U=0`$) half filled system, the structure of the ground state is trivial: two particles of opposite spin sit in the lowest energy level $`g`$, while $`|𝒮_{hf}|=(2N2)`$ particles lie on the shell $`𝒮_{hf}`$ of zero energy. In the spin $`S^z=0`$ subspace this ground state is $`\left(\begin{array}{c}2N2\\ N1\end{array}\right)^2`$ times degenerate. To first order in $`W`$, the degeneracy is only partially removedAt weak coupling it makes sense to speak about particles in filled shells, which behave much as core electrons in atomic physics, and particles in partially filled, or valence, shells. In first order, we must diagonalize the perturbation $`W`$ in the degenerate subspace, dropping the filled shells. Since $`W`$ is a positive semi-definite operator, with vanishing eigenvalue, 0 is necessarily the lowest one. Typically, the first-order solution is a crude approximation to the exact one, but below, when discussing the exact results, we shall demonstrate that many of the qualitative features hold for all $`U`$.. Indeed, if $`P_{S,M_S}`$ is the projection operator onto the subspace of total spin $`S`$ with $`z`$-component $`M_S`$, the structure of the determinantal state in Eq.(31) yields

$$P_{S,0}|\mathrm{\Phi }_{AF}^{(\sigma )}|\mathrm{\Phi }_{AF}^{S,0}0,S=0,\mathrm{},N1;$$
(32)
therefore, the states $`\{g_{}^{}g_{}^{}|\mathrm{\Phi }_{AF}^{S,0},S=0,\mathrm{},N1\}`$ belong to the ground state multiplet in first order perturbation theory (being $`W=0`$ states). Since the Lieb’s theorem ensures that the interacting ground state $`|\mathrm{\Psi }_0(N,U)`$ is a singlet, only $`g_{}^{}g_{}^{}|\mathrm{\Phi }_{AF}^{0,0}`$ can have a non-vanishing overlap with $`|\mathrm{\Psi }_0(N,U)`$ in the limit $`U0^+`$. In Section V we shall prove that $`g_{}^{}g_{}^{}|\mathrm{\Phi }_{AF}^{0,0}`$ is the ground state for $`U=0^+`$, that is $`lim_{U0}\mathrm{\Psi }_0(N,U)|g_{}^{}g_{}^{}|\mathrm{\Phi }_{AF}^{0,0}=1`$. In this model, the first-order solution has a peculiar significance. Usually, the exact ground state has no overlap with the non-interacting one, in the thermodynamic limit: this is because the $`U=0`$ ground states define a proper subspace of the full Hilbert space. Below, we prove that $`\mathrm{\Psi }_0(N,U)|g_{}^{}g_{}^{}|\mathrm{\Phi }_{AF}^{0,0}0U`$; that is, the lowest approximation keeps a finite weight.

At this stage, we note the analogy of the above results with those relevant to the standard Hubbard model defined on a $`N\times N`$ square lattice with periodic boundary conditions (hopping matrix elements only between nearest neighbor sites). The determinantal state $`|\mathrm{\Phi }_{AF}^{(\sigma )}`$ resembles the $`W=0`$ state that we obtained for even $`N`$ in Ref.. The degeneracy of the zero-energy one-body eigenspace was $`2N2`$ and it was shown that $`N1`$ zero-energy eigenfunctions have vanishing amplitudes on a sublattice while the remaining $`N1`$ vanish on the other.

As we shall see in the next Section, we need the explicit expression for the singlet $`|\mathrm{\Phi }_{AF}^{0,0}`$ and the triplet $`|\mathrm{\Phi }_{AF}^{1,m}`$, $`m=0,\pm 1`$, to calculate the ground state of the CBG-Hubbard model. Therefore, we now briefly review how to get the projection of the determinantal state $`|\mathrm{\Phi }_{AF}^{(\sigma )}`$ onto the singlet and the triplet spin subspaces.

To obtain the singlet $`|\mathrm{\Phi }_{AF}^{0,0}`$ one has to antisymmetrize each product $`a_i^{}b_j^{}`$, getting a two-body spin singlet operator, and subsequently antisymmetrize with respect to the $`N1`$ indices of the $`b^{}`$’s; one sees easily that this entails the antisymmetrization of the $`a^{}`$’s. Hence

$$|\mathrm{\Phi }_{AF}^{0,0}=\underset{P}{}()^P\underset{i=1}{\overset{N1}{}}\sigma _{i,P(i)}^{}|0$$
(33)
where $`\sigma _{i,j}^{}`$ creates a two-body singlet state:

$$\sigma _{i,j}^{}=\frac{1}{\sqrt{2}}(a_i^{}b_j^{}a_i^{}b_j^{}),$$
(34)
and $`P`$ is a permutation of the indices $`1,\mathrm{},N1`$. Since $`W`$ commutes with the total spin operators and the determinant $`|\mathrm{\Phi }_{AF}^{(\sigma )}`$ is a $`W=0`$ state, $`W|\mathrm{\Phi }_{AF}^{0,0}=0`$ as already noted. Applying $`W`$ on the state in Eq.(33), one can verify this property by direct inspection.

In a similar way one obtains the triplet projection. We define $`\tau _{i,j}^{(0)^{}}`$ as the two-body triplet creation operator with vanishing $`z`$-component $`m=0`$ in the orbitals $`a_i`$ and $`b_j`$:

$$\tau _{i,j}^{(0)^{}}=\frac{1}{\sqrt{2}}(a_i^{}b_j^{}+a_i^{}b_j^{}).$$
(35)
Then, one has

$$|\mathrm{\Phi }_{AF}^{1,0}=\underset{P}{}()^P\underset{j=1}{\overset{N1}{}}\tau _{j,P(j)}^{(0)^{}}\underset{ij}{}\sigma _{i,P(i)}^{}|0.$$
(36)
The $`m=\pm 1`$ components of the triplet in Eq.(36) can be otained by means of the total raising and lowering spin operators. Let us introduce the following notations

$`S_i^{+a}=a_i^{}a_i,S_i^a=[S_i^{+a}]^{},S^{+a}={\displaystyle \underset{i=1}{\overset{N1}{}}}S_i^{+a},S^a=[S^{+a}]^{}`$ (37)
$`S_i^{+b}=b_i^{}b_i,S_i^b=[S_i^{+b}]^{},S^{+b}={\displaystyle \underset{i=1}{\overset{N1}{}}}S_i^{+b},S^b=[S^{+b}]^{}`$ (38)
$`\widehat{n}_{i\sigma }^a=a_{i\sigma }^{}a_{i\sigma },S_i^{za}={\displaystyle \frac{1}{2}}(\widehat{n}_i^a\widehat{n}_i^a),\widehat{n}_\sigma ^a={\displaystyle \underset{i=1}{\overset{N1}{}}}\widehat{n}_{i\sigma }^a,\widehat{n}^a=\widehat{n}_{}^a+\widehat{n}_{}^a,S^{za}={\displaystyle \frac{1}{2}}(\widehat{n}_{}^a\widehat{n}_{}^a)`$ (39)
$`\widehat{n}_{i\sigma }^b=b_{i\sigma }^{}b_{i\sigma },S_i^{zb}={\displaystyle \frac{1}{2}}(\widehat{n}_i^b\widehat{n}_i^b),\widehat{n}_\sigma ^b={\displaystyle \underset{i=1}{\overset{N1}{}}}\widehat{n}_{i\sigma }^b,\widehat{n}^b=\widehat{n}_{}^a+\widehat{n}_{}^b,S^{zb}={\displaystyle \frac{1}{2}}(\widehat{n}_{}^b\widehat{n}_{}^b)`$ (40)
$`.`$ (41)
The states $`|\mathrm{\Phi }_{AF}^{1,\pm 1}`$ can be written as

$`|\mathrm{\Phi }_{AF}^{1,1}={\displaystyle \frac{1}{\sqrt{2}}}(S^{+a}+S^{+b})|\mathrm{\Phi }_{AF}^{1,0}={\displaystyle \underset{P}{}}()^P{\displaystyle \underset{j=1}{\overset{N1}{}}}\tau _{j,P(j)}^{(+1)^{}}{\displaystyle \underset{ij}{}}\sigma _{i,P(i)}^{}|0`$ (42)
$`|\mathrm{\Phi }_{AF}^{1,1}={\displaystyle \frac{1}{\sqrt{2}}}(S^a+S^b)|\mathrm{\Phi }_{AF}^{1,0}={\displaystyle \underset{P}{}}()^P{\displaystyle \underset{j=1}{\overset{N1}{}}}\tau _{j,P(j)}^{(1)^{}}{\displaystyle \underset{ij}{}}\sigma _{i,P(i)}^{}|0`$ (43)
with

$$\tau _{i,j}^{(+1)^{}}=a_i^{}b_j^{},\tau _{i,j}^{(1)^{}}=a_i^{}b_j^{}.$$
(44)
Equivalently, the triplet state $`|\mathrm{\Phi }_{AF}^{1,m}`$, $`m=0,\pm 1`$, can also be expressed in terms of the singlet $`|\mathrm{\Phi }_{AF}^{0,0}`$. It is a simple exercise to prove the following identities

$`{\displaystyle \underset{i=1}{\overset{N1}{}}}(S_i^{+a}S_i^{+b})|\mathrm{\Phi }_{AF}^{0,0}=(S^{+a}S^{+b})|\mathrm{\Phi }_{AF}^{0,0}=\sqrt{2}|\mathrm{\Phi }_{AF}^{1,1}`$ (45)
$`{\displaystyle \underset{i=1}{\overset{N1}{}}}(\widehat{n}_i^a\widehat{n}_i^b)|\mathrm{\Phi }_{AF}^{0,0}={\displaystyle \underset{i=1}{\overset{N1}{}}}(\widehat{n}_i^a+\widehat{n}_i^b)|\mathrm{\Phi }_{AF}^{0,0}=(S^{za}+S^{zb})|\mathrm{\Phi }_{AF}^{0,0}=|\mathrm{\Phi }_{AF}^{1,0}`$ (46)
$`{\displaystyle \underset{i=1}{\overset{N1}{}}}(S_i^a+S_i^b)|\mathrm{\Phi }_{AF}^{0,0}=(S^a+S^b)|\mathrm{\Phi }_{AF}^{0,0}=\sqrt{2}|\mathrm{\Phi }_{AF}^{1,1}.`$ (47)
As for the singlet, one can verify that $`W|\mathrm{\Phi }_{AF}^{1,m}=0`$, $`m=0,\pm 1`$, using the definitions in Eqs.(36-42-43).

## V The Ground State at Half Filling

We are now ready to calculate the ground state $`|\mathrm{\Psi }_0(N,U)`$ of the half filled CBG-Hubbard model described in Section II. As already observed, $`g_{}^{}g_{}^{}|\mathrm{\Phi }_{AF}^{0,0}`$ is a good candidate for the non-interacting ground state; its quantum numbers are the same as those of $`|\mathrm{\Psi }_0(N,U)`$ and moreover it belongs to the first-order ground state multiplet. Using the short-hand notations

$$|g^0g_{}^{}g_{}^{}|0,|\overline{g}^0\overline{g}_{}^{}\overline{g}_{}^{}|0,|[g\overline{g}]^{0,0}\frac{1}{\sqrt{2}}(g_{}^{}\overline{g}_{}^{}g_{}^{}\overline{g}_{}^{})|0$$
(48)
for the three different two-body singlets that one gets from the lowest and the highest energy orbitals $`g`$ and $`\overline{g}`$ and

$$|[g\overline{g}]^{1,1}g_{}^{}\overline{g}_{}^{}|0,|[g\overline{g}]^{1,0}\frac{1}{\sqrt{2}}(g_{}^{}\overline{g}_{}^{}+g_{}^{}\overline{g}_{}^{})|0,|[g\overline{g}]^{1,1}g_{}^{}\overline{g}_{}^{}|0$$
(49)
for the triplet, we propose the following Ansatz for the interacting ground state $`|\mathrm{\Psi }_0(N,U)`$:

$$|\mathrm{\Psi }_0(N,U)=\left[\gamma _g|g^0+\gamma _{\overline{g}}|\overline{g}^0\right]|\mathrm{\Phi }_{AF}^{0,0}+\gamma _0\underset{m=1}{\overset{1}{}}()^m|[g\overline{g}]^{1,m}|\mathrm{\Phi }_{AF}^{1,m},$$
(50)
where the $`\gamma `$’s are $`c`$-numbers. It is worth to note that in $`|\mathrm{\Psi }_0(N,U)`$ the number of particles in the shell $`𝒮_{hf}`$ is a constant given by $`2N2`$. A priori, there is no reason for this choice since the total number operator $`\widehat{n}^a+\widehat{n}^b`$ of particles in the shell $`𝒮_{hf}`$ does not commute with the Hamiltonian. Nevertheless, we shall see that the scatterings which do not preserve this number cancel out provided in $`𝒮_{hf}`$ the $`(2N2)`$-body state is a $`W=0`$ state. We shall refer to this very remarkable property as to the Shell Population Rigidity. We emphasize that even constraining the particle-number in $`𝒮_{hf}`$ to be $`2N2`$ with vanishing total spin $`z`$-component, there are $`\left(\begin{array}{c}2N2\\ N1\end{array}\right)^2`$ configurations that can contribute to the ground state expansion in the interacting case, while the Ansatz (50) contains only $`|\mathrm{\Phi }_{AF}^{0,0}`$ and $`|\mathrm{\Phi }_{AF}^{1,m}`$, $`m=0,\pm 1`$. The reason why the $`W=0`$ states $`|\mathrm{\Phi }_{AF}^{S,M_S}`$ with $`S>1`$ do not enter into the ground state expansion (50) comes from the Lieb’s theorem: the ground state must be a singlet and with only two particles outside $`𝒮_{hf}`$ (in the $`g`$ and/or $`\overline{g}`$ states) the angular momentum composition law forbids $`W=0`$ states with $`S>1`$. However, we observe that the state $`|[g\overline{g}]^{0,0}|\mathrm{\Phi }_{AF}^{0,0}`$ is a singlet and $`|\mathrm{\Phi }_{AF}^{0,0}`$ is a $`W=0`$ state. It has the right quantum numbers and, in principle, it could have a non-zero overlap with $`|\mathrm{\Psi }_0(N,U)`$. Nevertheless, the matrix elements of $`W`$ between this state and the ones in the Ansatz are zero. This is the reason why we have dropped it in the expansion (50).

The three states in Eq.(50) are eigenstates of the kinetic-energy operator $`H_0`$:

$`H_0`$ $`|g^0|\mathrm{\Phi }_{AF}^{0,0}=2t|g^0|\mathrm{\Phi }_{AF}^{0,0}`$ (51)
$`H_0`$ $`|\overline{g}^0|\mathrm{\Phi }_{AF}^{0,0}=2t|\overline{g}^0|\mathrm{\Phi }_{AF}^{0,0}`$ (53)
$`H_0`$ $`{\displaystyle \underset{m=1}{\overset{1}{}}}()^m|[g\overline{g}]^{1,m}|\mathrm{\Phi }_{AF}^{1,m}=0.`$ (54)
Expanding $`c_xc_x`$ occurring in $`W`$ in $`a,b,g,\overline{g}`$ operators of Eq.(12), we may take advantage from the fact that $`|\mathrm{\Phi }_{AF}^{0,0}`$ and $`|\mathrm{\Phi }_{AF}^{1,m}`$, $`m=0,\pm 1`$, are $`W=0`$ states. Indeed, since one cannot annihilate $`g`$ or $`\overline{g}`$ over them, taking for instance $`x𝒜`$,

$$0=c_xc_x|\mathrm{\Phi }_{AF}^{S,M_S}=\underset{i,j}{}R_{i,x}R_{j,x}a_ia_j|\mathrm{\Phi }_{AF}^{S,M_S};$$
(55)
then, multiplying by $`g_{}^{}g_{}^{}`$, and the like, one obtains that the contribution proportional to

$$\underset{i,j}{}R_{i,x}R_{j,x}a_ia_j|\mathrm{\Psi }_0(N,U)$$
(56)
yield nothing. Of course, if $`x`$, by the same reasoning we get

$$\underset{i,j}{}R_{i,x}R_{j,x}b_ib_j|\mathrm{\Psi }_0(N,U)=0.$$
(57)
Hence, we can write

$$W=W_𝒜+W_{}$$
(58)
with

$$W_𝒜:=U\underset{x=1}{\overset{N}{}}c_x^{}c_x^{}\left[\frac{1}{2N}(g_{}+\overline{g}_{})(g_{}+\overline{g}_{})+\frac{1}{\sqrt{2N}}\underset{i=1}{\overset{N1}{}}R_{i,x}[(g_{}+\overline{g}_{})a_i+a_i(g_{}+\overline{g}_{})]\right]$$
(59)
and

$$W_{}:=U\underset{x=N+1}{\overset{2N}{}}c_x^{}c_x^{}[\frac{1}{2N}(g_{}+\overline{g}_{})(g_{}+\overline{g}_{})+\frac{1}{\sqrt{2N}}\underset{i=1}{\overset{N1}{}}R_{i,xN}[(g_{}+\overline{g}_{})b_i+b_i(g_{}+\overline{g}_{})],$$
(60)
where $`:=`$ means that the two sides are equivalent when acting on $`|\mathrm{\Psi }_0(N,U)`$. Next, we transform the remaining two $`c^{}`$ operators, using Eq.(14). We note that the terms containing the sequence $`g^{}g^{}ga`$ and similar terms with $`b`$ instead of $`a`$ and/or some $`g`$ replaced by $`\overline{g}`$ vanish by symmetry; this is evident since $`_{x=1}^NR_{i,x}=0`$ for any $`i=1,\mathrm{},N1`$. Similarly, the sequences like $`a^{}g^{}gg`$ also do not contribute. The sequences like $`a^{}a^{}gg`$ are diagonal in the orbital indices of the two $`a^{}`$ operators, since $`_{x=1}^NR_{i,x}R_{j,x}=\delta _{ij}`$; hence, they annihilate $`|\mathrm{\Psi }_0(N,U)`$ (recall that in this state all the $`a`$ orbitals are singly occupied, as we can see by direct inspection of the expressions for $`|\mathrm{\Phi }_{AF}^{0,0}`$ and $`|\mathrm{\Phi }_{AF}^{1,m}`$, $`m=0,\pm 1`$, of Section IV). By a particle-hole transformation one can show that the remaining terms that create pairs in the shell $`𝒮_{hf}`$, like $`a^{}a^{}ga`$, etc., have no effect, as in Eqs.(56-57); we postpone the proof of this last spectacular cancellation until Appendix A. To sum up, no term in $`W`$ can change the number of particles in $`𝒮_{hf}`$, that is, that number is independent of $`t`$ and $`U`$! This proves the remarkable property that we have named Shell Population Rigidity: it holds for any finite $`N`$ and it is not specific of the thermodynamic limit. We shall see below that the Shell Population Rigidity characterizes the ground state and a suitable subspace of the total Hilbert space. In Fig.2 we have drawn a physical picture of the cancellation among diagrams which do not preserve the number $`n_a+n_b`$. We can say that the $`W=0`$ states entering in our Ansatz are stable with respect to the Hubbard interaction $`W`$ also in the presence of other particles in the system. This fact allows to exactly solve the Schrödinger equation.

Taking these cancellations into account one obtains

$`W_𝒜`$ $`:=`$ $`{\displaystyle \frac{U}{4N}}(g_{}^{}+\overline{g}_{}^{})(g_{}^{}+\overline{g}_{}^{})(g_{}+\overline{g}_{})(g_{}+\overline{g}_{})+`$ (61)
$`+`$ $`{\displaystyle \frac{U}{2N}}{\displaystyle \underset{i=1}{\overset{N1}{}}}\left[a_i^{}(g_{}^{}+\overline{g}_{}^{})+(g_{}^{}+\overline{g}_{}^{})a_i^{}\right]\left[(g_{}+\overline{g}_{})a_i+a_i(g_{}+\overline{g}_{})\right]`$ (62)
and

$`W_{}`$ $`:=`$ $`{\displaystyle \frac{U}{4N}}(g_{}^{}+\overline{g}_{}^{})(g_{}^{}+\overline{g}_{}^{})(g_{}+\overline{g}_{})(g_{}+\overline{g}_{})+`$ (63)
$`+`$ $`{\displaystyle \frac{U}{2N}}{\displaystyle \underset{i=1}{\overset{N1}{}}}\left[b_i^{}(g_{}^{}+\overline{g}_{}^{})+(g_{}^{}+\overline{g}_{}^{})b_i^{}\right]\left[(g_{}+\overline{g}_{})b_i+b_i(g_{}+\overline{g}_{})\right].`$ (64)
With $`W`$ written in the transformed picture one can calculate $`W|\mathrm{\Psi }_0(N,U)`$. After very long, but standard algebra one finds

$`W|g^0|\mathrm{\Phi }_{AF}^{0,0}`$ $`=`$ $`U\{{\displaystyle \frac{2N1}{2N}}|g^0|\mathrm{\Phi }_{AF}^{0,0}+{\displaystyle \frac{1}{2N}}|\overline{g}^0|\mathrm{\Phi }_{AF}^{0,0}+`$ (65)
$`+`$ $`{\displaystyle \frac{1}{2N}}{\displaystyle \underset{i=1}{\overset{N1}{}}}|[g\overline{g}]^{1,1}(S_i^a+S_i^b)|\mathrm{\Phi }_{AF}^{0,0}+`$ (66)
$`+`$ $`{\displaystyle \frac{1}{2N}}{\displaystyle \underset{i=1}{\overset{N1}{}}}\left[\overline{g}_{}^{}g_{}^{}(\widehat{n}_i^an_i^b)+\overline{g}_{}^{}g_{}^{}(\widehat{n}_i^a+n_i^b)\right]|\mathrm{\Phi }_{AF}^{0,0}+`$ (67)
$`+`$ $`{\displaystyle \frac{1}{2N}}{\displaystyle \underset{i=1}{\overset{N1}{}}}|[g\overline{g}]^{1,1}(S_i^{+a}S_i^{+b})|\mathrm{\Phi }_{AF}^{0,0}\}`$ (68)
and

$`W|\overline{g}^0|\mathrm{\Phi }_{AF}^{0,0}`$ $`=`$ $`U\{{\displaystyle \frac{2N1}{2N}}|\overline{g}^0|\mathrm{\Phi }_{AF}^{0,0}+{\displaystyle \frac{1}{2N}}|g^0|\mathrm{\Phi }_{AF}^{0,0}+`$ (69)
$`+`$ $`{\displaystyle \frac{1}{2N}}{\displaystyle \underset{i=1}{\overset{N1}{}}}|[g\overline{g}]^{1,1}(S_i^aS_i^b)|\mathrm{\Phi }_{AF}^{0,0}+`$ (70)
$`+`$ $`{\displaystyle \frac{1}{2N}}{\displaystyle \underset{i=1}{\overset{N1}{}}}\left[g_{}^{}\overline{g}_{}^{}(\widehat{n}_i^a+n_i^b)+g_{}^{}\overline{g}_{}^{}(\widehat{n}_i^an_i^b)\right]|\mathrm{\Phi }_{AF}^{0,0}+`$ (71)
$`+`$ $`{\displaystyle \frac{1}{2N}}{\displaystyle \underset{i=1}{\overset{N1}{}}}|[g\overline{g}]^{1,1}(S_i^{+a}+S_i^{+b})|\mathrm{\Phi }_{AF}^{0,0}\},`$ (72)
while for the singlet $`_{m=1}^1()^m|[g\overline{g}]^{1,m}|\mathrm{\Phi }_{AF}^{1,m}`$ one obtains

$`W{\displaystyle \underset{m=1}{\overset{1}{}}}`$ $`()^m`$ $`|[g\overline{g}]^{1,m}|\mathrm{\Phi }_{AF}^{1,m}=U\{{\displaystyle \frac{N+1}{N}}{\displaystyle \underset{m=1}{\overset{1}{}}}()^m|[g\overline{g}]^{1,m}|\mathrm{\Phi }_{AF}^{1,m}+{\displaystyle \frac{1}{2N}}(|g^0|\overline{g}^0)`$ (73)
$``$ $`{\displaystyle \underset{i=1}{\overset{N1}{}}}[(S_i^{+a}S_i^{+b})|\mathrm{\Phi }_{AF}^{1,1}+{\displaystyle \frac{1}{\sqrt{2}}}(n_i^an_i^an_i^b+n_i^b)|\mathrm{\Phi }_{AF}^{1,0}+(S_i^a+S_i^b)|\mathrm{\Phi }_{AF}^{1,1}]\}.`$ (74)
By means of the identities in Eqs.(45-46-47) one can write Eqs.(68-72) in a compact form

$$W|g^0|\mathrm{\Phi }_{AF}^{0,0}=U\left\{\frac{2N1}{2N}|g^0|\mathrm{\Phi }_{AF}^{0,0}+\frac{1}{2N}|\overline{g}^0|\mathrm{\Phi }_{AF}^{0,0}+\frac{\sqrt{2}}{2N}\underset{m=1}{\overset{1}{}}()^m|[g\overline{g}]^{1,m}|\mathrm{\Phi }_{AF}^{1,m}\right\}$$
(75)
and

$$W|\overline{g}^0|\mathrm{\Phi }_{AF}^{0,0}=U\left\{\frac{2N1}{2N}|\overline{g}^0|\mathrm{\Phi }_{AF}^{0,0}+\frac{1}{2N}|g^0|\mathrm{\Phi }_{AF}^{0,0}\frac{\sqrt{2}}{2N}\underset{m=1}{\overset{1}{}}()^m|[g\overline{g}]^{1,m}|\mathrm{\Phi }_{AF}^{1,m}\right\}.$$
(76)
Let us now consider the second row in Eq.(74). If our Ansatz is correct, the state in the square brackets must be proportional to $`|\mathrm{\Phi }_{AF}^{0,0}`$. In this way one could close the equations and find an exact eigenstate. In Appendix B we prove that this is the case and in particular that

$$\underset{i=1}{\overset{N1}{}}\left[(S_i^{+a}S_i^{+b})|\mathrm{\Phi }_{AF}^{1,1}+\frac{1}{\sqrt{2}}(n_i^an_i^an_i^b+n_i^b)|\mathrm{\Phi }_{AF}^{1,0}+(S_i^a+S_i^b)|\mathrm{\Phi }_{AF}^{1,1}\right]=\sqrt{2}(N^21)|\mathrm{\Phi }_{AF}^{0,0}.$$
(77)
Hence, Eq.(74) yields

$`W{\displaystyle \underset{m=1}{\overset{1}{}}}()^m|[g\overline{g}]^{1,m}|\mathrm{\Phi }_{AF}^{1,m}=U\{{\displaystyle \frac{N+1}{N}}{\displaystyle \underset{m=1}{\overset{1}{}}}()^m|[g\overline{g}]^{1,m}|\mathrm{\Phi }_{AF}^{1,m}+`$ (78)
$`+\sqrt{2}{\displaystyle \frac{N^21}{2N}}(|g^0|\overline{g}^0)|\mathrm{\Phi }_{AF}^{0,0}\}.`$ (79)
This result, together with Eqs.(75-76) and Eqs.(54), allows us to reduce the Schrodinger equation $`(H_{\mathrm{H}ubbard}E)|\mathrm{\Psi }_0(N,U)=0`$ to the diagonalization of the matrix

$$H(N,U)=\left(\begin{array}{ccc}2t+\frac{2N1}{2N}U& \frac{U}{2N}& \frac{\sqrt{2(N^21)}}{2N}U\\ & & \\ \frac{U}{2N}& 2t+\frac{2N1}{2N}U& \frac{\sqrt{2(N^21)}}{2N}U\\ & & \\ \frac{\sqrt{2(N^21)}}{2N}U& \frac{\sqrt{2(N^21)}}{2N}U& \frac{N+1}{N}U\end{array}\right),$$
(80)
where the relations

$$\left[g^0|\mathrm{\Phi }_{AF}^{0,0}|\right]\left[|g^0|\mathrm{\Phi }_{AF}^{0,0}\right]=\left[\overline{g}^0|\mathrm{\Phi }_{AF}^{0,0}|\right]\left[|\overline{g}^0|\mathrm{\Phi }_{AF}^{0,0}\right]$$
(81)
and

$`g^0|\mathrm{\Phi }_{AF}^{0,0}|W{\displaystyle \underset{m=1}{\overset{1}{}}}()^m|[g\overline{g}]^{1,m}|\mathrm{\Phi }_{AF}^{1,m}={\displaystyle \underset{m=1}{\overset{1}{}}}()^m\mathrm{\Phi }_{AF}^{1,m}|[g\overline{g}]^{1,m}|W|g^0|\mathrm{\Phi }_{AF}^{0,0},`$ (82)
$`\overline{g}^0|\mathrm{\Phi }_{AF}^{0,0}|W{\displaystyle \underset{m=1}{\overset{1}{}}}()^m|[g\overline{g}]^{1,m}|\mathrm{\Phi }_{AF}^{1,m}={\displaystyle \underset{m=1}{\overset{1}{}}}()^m\mathrm{\Phi }_{AF}^{1,m}|[g\overline{g}]^{1,m}|W|\overline{g}^0|\mathrm{\Phi }_{AF}^{0,0}`$ (83)
have been used to express $`_{m=1}^1\mathrm{\Phi }_{AF}^{1,m}|\mathrm{\Phi }_{AF}^{1,m}`$ in terms of $`\mathrm{\Phi }_{AF}^{0,0}|\mathrm{\Phi }_{AF}^{0,0}`$.

The thermodynamic limit $`N\mathrm{}`$ is well defined and non-trivial if we rescale $`\stackrel{~}{t}`$ in such a way that $`N\stackrel{~}{t}=t=`$const (which corresponds to the case of a fixed energy gap between the first two energy levels). The Hamiltonian matrix in Eq.(80) becomes

$$H(\mathrm{},U)=\left(\begin{array}{ccc}2t+U& 0& U/\sqrt{2}\\ 0& 2t+U& U/\sqrt{2}\\ U/\sqrt{2}& U/\sqrt{2}& U\end{array}\right)$$
(84)
with eigenvalues

$$_0=U,_\pm =U\pm \mathrm{\Delta },\mathrm{\Delta }\sqrt{U^2+4t^2}.$$
(85)
The (unnormalized) eigenvector corresponding to the lowest eigenvalue $`_{}`$ gives

$$|\mathrm{\Psi }_0(\mathrm{},U)=\left[\left(\frac{U}{t}\frac{2\mathrm{\Delta }}{U}(2+\frac{\mathrm{\Delta }}{t})\right)|g^0+\frac{U}{t}|\overline{g}^0\right]|\mathrm{\Phi }_{AF}^{0,0}+\sqrt{2}(2+\frac{\mathrm{\Delta }}{t})\underset{m=1}{\overset{1}{}}()^m|[g\overline{g}]^{1,m}|\mathrm{\Phi }_{AF}^{1,m}.$$
(86)
So far we have proved that the Ansatz in Eq.(50) gives three exact eigenstates of the CBG-Hubbard model at half filling. The eigenvalue of lowest energy $`_{}`$ reduces to the non-interacting ground state energy for $`U0`$, while the associated eigenstate (86) to the state $`|g^0|\mathrm{\Phi }_{AF}^{0,0}`$. On the other hand, when $`U\mathrm{}`$ the original Hamiltonian can be mapped onto what we can define the CBG-Heisenberg Hamiltonian:

$$H_{\mathrm{H}eisenberg}=\frac{4\stackrel{~}{t}^2}{U}\underset{x𝒜}{}\underset{y}{}(𝐒_x𝐒_y\frac{1}{4}).$$
(87)
This model is exactly solvable and the ground state is obtained by projecting onto the singlet the state where $`N`$ particles lie on the $`𝒜`$-sites with spin up and the remaining $`N`$ on the $``$-sites with spin down (Neel state): one finds

$$|\mathrm{\Psi }_0(N,\mathrm{})=\underset{P}{}()^P\underset{x=1}{\overset{N}{}}\sigma _{x,P(x+N)}^{}|0$$
(88)
with $`\sigma _{x,y}^{}=\frac{1}{\sqrt{2}}(c_x^{}c_y^{}c_x^{}c_y^{})`$. The ground state energy is $`2t^2/U`$ for $`N\mathrm{}`$ and is equal to the first order approximation of $`_{}`$ in the small parameter $`(t/U)^2`$ (the same conclusion holds for any finite $`N`$, although it is more tedious to prove). Furthermore, by direct inspection of the lowest energy eigenvector of $`H(N,U)`$ one can show that $`|\mathrm{\Psi }_0(N,U)`$ becomes the state in Eq.(88) when $`U\mathrm{}`$. Therefore, $`|\mathrm{\Psi }_0(N,U)`$ reduces to the ground state of the CBG-Hubbard model in the two opposite limits $`U0`$ and $`U\mathrm{}`$.

To prove that the lowest energy eigenstate of $`H(N,U)`$ is the unique ground state we are looking for, one can exploit the ground state uniqueness of the Heisenberg Hamiltonian proved in Ref.. Since $`|\mathrm{\Psi }_0(N,\mathrm{})`$ is the ground state of the Heisenberg model and since the ground state of the half-filled Hubbard model is unique, a level crossing for some value of $`U`$ would be required if $`|\mathrm{\Psi }_0(N,U)`$ were an excited state, contradicting the uniqueness.

In conclusion we have proved that the Ansatz in Eq.(50) with $`(\gamma _g,\gamma _{\overline{g}},\gamma _0)`$ given by the lowest energy eigenvector of the matrix in Eq.(80) is the exact ground state of the half-filled Hubbard model defined on the CBG of size $`|\mathrm{\Lambda }|=2N`$ and repulsion $`U`$. In the next Section we shall calculate some physical quantities, as the spin-spin correlation function, and we shall prove that the particles are antiferromagnetically ordered.

## VI Results and Discussion

The exact ground state solution for arbitrary but finite $`N`$ allows to study the ground state energy $`E_0(N,U)`$ as a function of $`N`$ and $`U`$. Taking $`\stackrel{~}{t}N=t=1`$, in Fig.(3) we have plotted $`E_0(N,U)`$ in the range $`N=1,..,10`$ and $`0<U<20`$.

The figure shows that $`E_0(N,U)`$ is a monotonically increasing function of $`N`$ and $`U`$, that is $`E_0(N+1,U)E_0(N)>0`$ and $`E_0(N,U)/U>0`$. The limit $`N\mathrm{}`$ at fixed $`U`$ is given by $`_{}`$ in Eq.(85):

$$\underset{|\mathrm{\Lambda }|\mathrm{}}{lim}\underset{\beta \mathrm{}}{lim}E^{(\beta )}(|\mathrm{\Lambda }|,t,U)=U\sqrt{U^2+4t^2}<0$$
(89)
where $`E^{(\beta )}(|\mathrm{\Lambda }|,t,U)`$ is the internal energy of the system at inverse Temperature $`\beta `$ and number of lattice sites $`|\mathrm{\Lambda }|=2N`$. On the other hand, one can calculate the internal energy at zero Temperature in thermodynamic limit: $`E^{(\beta )}(|\mathrm{\Lambda }|\mathrm{},t,U)=\mathrm{\Omega }_\beta +k_BS/\beta +\mu |\mathrm{\Lambda }|n`$ with $`\mathrm{\Omega }_\beta `$ from Eq.(23) and $`S`$ the entropy. Following the procedure of van Dongen and Vollhardt one gets

$$\underset{\beta \mathrm{}}{lim}\underset{|\mathrm{\Lambda }|\mathrm{}}{lim}E^{(\beta )}(|\mathrm{\Lambda }|,t,U)=\{\begin{array}{cc}2t=2N\stackrel{~}{t}<0\hfill & \mathrm{i}f\stackrel{~}{t}\mathrm{i}s\mathrm{n}ot\mathrm{r}escaled\hfill \\ O(1)\hfill & \mathrm{i}f\stackrel{~}{t}1/|\mathrm{\Lambda }|\hfill \end{array}.$$
(90)
In the case $`\stackrel{~}{t}`$=const the two limits commute and the trivial thermodynamic behavior of the model originates from the infinite gap between the lowest level and the zero-energy levels: the two particles in the $`g`$ orbitals are completely decoupled from the dynamics of the system. On the other hand, if $`\stackrel{~}{t}`$ is rescaled the gap remains frozen and taking $`|\mathrm{\Lambda }|\mathrm{}`$ first, we can only say that the ground state energy is $`O(1)`$, but we cannot predict the exact amount within the scheme proposed by van Dongen and Vollhardt. This is a consequence of the fact that in the limit $`|\mathrm{\Lambda }|\mathrm{}`$, we retain only the extensive contributions to the internal energy.

The ground state energy Eq.(89) is always higher than that of the non-interacting case and hence non-trivial correlations survive in the thermodynamic limit. We have calculated the ground state average of the number of doubly occupied sites:

$`\overline{D}(N,U){\displaystyle \frac{\overline{W}(N,U)}{U}}{\displaystyle \frac{1}{U}}\mathrm{\Psi }_0(N,U)|W|\mathrm{\Psi }_0(N,U)`$ $`=`$ $`{\displaystyle \frac{1}{U}}\mathrm{\Psi }_0(N,U)|H_{\mathrm{H}ubbard}H_0|\mathrm{\Psi }_0(N,U)=`$ (91)
$`=`$ $`{\displaystyle \frac{1}{U}}E_0(N,U)+{\displaystyle \frac{2t}{U}}[\gamma _g^2(N,U)\gamma _{\overline{g}}^2(N,U)]`$ (92)
where $`\gamma _g`$ and $`\gamma _{\overline{g}}`$ are the first two components of the normalized ground state vector, see Eq.(50). In Fig.(4$`a`$) we have plotted the trend of $`\overline{D}(N,U)`$ as $`N`$ increases for different values of $`U`$. As expected $`\overline{D}(N,U)`$ is a monotonically decreasing function of $`U`$, $`\overline{D}(N,U)/U<0`$, approaching zero for $`U\mathrm{}`$ \[where the exact ground state reduces to the antiferromagnetic Neel state, see Eq.(88)\]. Nevertheless, $`\overline{D}(N,U)`$ shows a non-trivial behaviour as $`N`$ increases at fixed $`U`$ values. In the weak coupling regime, $`Ut`$, the number of doubly occupied sites grows as $`N`$ becomes larger and larger converging to a finite value when $`N\mathrm{}`$. The opposite trend is observed in the strong coupling regime, $`Ut`$, where $`\overline{D}(N,U)`$ decreases as $`N`$ increases. In the intermediate regime $`Ut`$, $`\overline{D}(N,U)`$ is an increasing function of $`N`$ for small $`N`$, but becomes a decreasing function for large $`N`$. Hence, there is a critical value $`U_c(N)`$ where the analytic continuation of $`\overline{D}(N,U)`$ to real $`N`$ verifies

$$\left(\frac{^2\overline{D}}{UN}\right)_{U_c(N)}=0$$
(93)
In Fig.(4$`b`$) we report the ground state average of the interacting term $`W`$ as a function of $`U`$ for three different values of $`N`$. There is always a critical value for $`U`$ where the repulsion is maximized. Moreover, $`lim_N\mathrm{}\overline{W}(N,U)\overline{W}_{\mathrm{}}(U)0`$ for any finite $`U`$ and the system cannot avoid double occupation neither in the thermodynamic limit. In the limit $`U\mathrm{}`$, $`\overline{W}(N,U)\frac{1}{U}`$, that is, $`\overline{W}(N,U)0`$ in agreement with the Gutzwiller Ansatz.

Next, we have calculated the spin-spin correlation function

$$G_{\mathrm{s}pin}(x,y)\mathrm{\Psi }_0(N,U)|S_x^zS_y^z|\mathrm{\Psi }_0(N,U)=\underset{\beta \mathrm{}}{lim}G_{\mathrm{s}pin}(x,y;\beta ,h=0,|\mathrm{\Lambda }|)=\underset{\beta \mathrm{}}{lim}\underset{h0}{lim}G_{\mathrm{s}pin}(x,y;\beta ,h,|\mathrm{\Lambda }|)$$
(94)
where $`G_{\mathrm{s}pin}(x,y;\beta ,h,|\mathrm{\Lambda }|)`$ was defined in Section III, $`S_x^z=\frac{1}{2}(\widehat{n}_x\widehat{n}_x)`$ and $`|\mathrm{\Psi }_0`$ is the normalized ground state

$$|\mathrm{\Psi }_0=\left[\frac{\gamma _g}{𝒩_0}|g^0+\frac{\gamma _{\overline{g}}}{𝒩_0}|\overline{g}^0\right]|\mathrm{\Phi }_{AF}^{0,0}+\frac{1}{\sqrt{3}}\frac{\gamma _0}{𝒩_1}\underset{m=1}{\overset{1}{}}()^m|[g\overline{g}]^{1,m}|\mathrm{\Phi }_{AF}^{1,m},$$
(95)
with $`𝒩_0^2=\mathrm{\Phi }_{AF}^{0,0}|\mathrm{\Phi }_{AF}^{0,0}`$, $`𝒩_1^2=\mathrm{\Phi }_{AF}^{1,m}|\mathrm{\Phi }_{AF}^{1,m}`$, $`m=0,\pm 1`$, and $`\gamma _g^2+\gamma _{\overline{g}}^2+\gamma _0^2=1`$. Due to the $`S_NS_NZ_2`$ symmetry, $`G_{\mathrm{s}pin}(x,y)`$ can be written as

$$G_{\mathrm{s}pin}(x,y)=\{\begin{array}{cc}G_0\hfill & x=y\hfill \\ G_{\mathrm{o}n}\hfill & x𝒜(x)\mathrm{a}ndy𝒜(y)\hfill \\ G_{\mathrm{o}ff}\hfill & x𝒜(x)\mathrm{a}ndy(x𝒜)\hfill \end{array}$$
(96)
and by exploiting the sum rule $`_{y\mathrm{\Lambda }}G_{\mathrm{s}pin}(x,y)=G_0+(N1)G_{\mathrm{o}n}+NG_{\mathrm{o}ff}=0`$, one can express $`G_{\mathrm{o}n}`$ in terms of $`G_{\mathrm{o}ff}`$ and $`G_0`$:

$$G_{\mathrm{o}n}=\frac{N}{N1}G_{\mathrm{o}ff}\frac{G_0}{N1}.$$
(97)
The problem is then reduced to the evaluation of $`G_0`$ and $`G_{\mathrm{o}ff}`$. We have

$`G_0`$ $`=`$ $`G_{\mathrm{s}pin}(x,x)=\mathrm{\Psi }_0|(S_x^z)^2|\mathrm{\Psi }_0={\displaystyle \frac{1}{2N}}{\displaystyle \underset{x\mathrm{\Lambda }}{}}\mathrm{\Psi }_0|(S_x^z)^2|\mathrm{\Psi }_0=`$ (98)
$`=`$ $`{\displaystyle \frac{1}{8N}}{\displaystyle \underset{x\mathrm{\Lambda }}{}}\mathrm{\Psi }_0|\widehat{n}_x+\widehat{n}_x2\widehat{n}_x\widehat{n}_x|\mathrm{\Psi }_0={\displaystyle \frac{1}{4}}(1{\displaystyle \frac{\overline{D}}{N}}).`$ (99)
More effort is needed to calculate $`G_{\mathrm{o}ff}`$. Expanding the $`c`$ operators in $`g`$, $`\overline{g}`$, $`a`$ and $`b`$ operators of Eq.(12) and taking into account that $`(_{x\mathrm{\Lambda }}S_x^z)^2|\mathrm{\Psi }_0=0`$, we get

$`G_{\mathrm{o}ff}={\displaystyle \frac{1}{N^2}}{\displaystyle \underset{x𝒜}{}}{\displaystyle \underset{y𝒜}{}}\mathrm{\Psi }_0|S_x^zS_y^z|\mathrm{\Psi }_0={\displaystyle \frac{1}{N^2}}[{\displaystyle \frac{1}{16}}\mathrm{\Psi }_0|\left({\displaystyle \underset{\sigma }{}}()^\sigma (g_\sigma ^{}+\overline{g}_\sigma ^{})(g_\sigma +\overline{g}_\sigma )\right)^2|\mathrm{\Psi }_0+`$ (100)
$`+{\displaystyle \frac{1}{2}}\mathrm{\Psi }_0|S^{za}{\displaystyle \underset{\sigma }{}}()^\sigma (g_\sigma ^{}+\overline{g}_\sigma ^{})(g_\sigma +\overline{g}_\sigma )|\mathrm{\Psi }_0+\mathrm{\Psi }_0|S^{za}S^{za}|\mathrm{\Psi }_0]`$ (101)
with $`()^\sigma =+,`$ for $`\sigma =,`$. After long but standard algebra one obtains

$`G_{\mathrm{o}ff}={\displaystyle \frac{1}{N^2}}\{{\displaystyle \frac{1}{8}}[(\gamma _g\gamma _{\overline{g}})^2+2\gamma _0^2]\sqrt{{\displaystyle \frac{2}{3}}}{\displaystyle \frac{\gamma _0(\gamma _g\gamma _{\overline{g}})}{𝒩_0𝒩_1}}\mathrm{\Phi }_{AF}^{1,0}|S^{za}|\mathrm{\Phi }_{AF}^{0,0}+{\displaystyle \frac{\gamma _g^2+\gamma _{\overline{g}}^2}{𝒩_0^2}}\mathrm{\Phi }_{AF}^{0,0}|S^{za}S^{za}|\mathrm{\Phi }_{AF}^{0,0}+`$ (102)
$`+{\displaystyle \frac{\gamma _0^2}{3𝒩_1^2}}[\mathrm{\Phi }_{AF}^{1,1}|S^{za}|\mathrm{\Phi }_{AF}^{1,1}\mathrm{\Phi }_{AF}^{1,1}|S^{za}|\mathrm{\Phi }_{AF}^{1,1}+{\displaystyle \underset{m=1}{\overset{1}{}}}\mathrm{\Phi }_{AF}^{1,m}|S^{za}S^{za}|\mathrm{\Phi }_{AF}^{1,m}]\}.`$ (103)
and exploiting the $`Z_2`$ symmetry $`a_ib_i`$ and $`b_ia_i`$, which implies $`|\mathrm{\Phi }_{AF}^{S,M_S}()^S|\mathrm{\Phi }_{AF}^{S,M_S}`$,

$$\mathrm{\Phi }_{AF}^{1,0}|S^{za}|\mathrm{\Phi }_{AF}^{0,0}=\mathrm{\Phi }_{AF}^{1,0}|S^{zb}|\mathrm{\Phi }_{AF}^{0,0}=\frac{1}{2}\mathrm{\Phi }_{AF}^{1,0}|S^{za}S^{zb}|\mathrm{\Phi }_{AF}^{0,0}=\frac{1}{2}𝒩_1^2$$
(104)
$$\mathrm{\Phi }_{AF}^{1,1}|S^{za}|\mathrm{\Phi }_{AF}^{1,1}=\frac{1}{2}\mathrm{\Phi }_{AF}^{1,1}|S^{za}+S^{zb}|\mathrm{\Phi }_{AF}^{1,1}=\frac{1}{2}𝒩_1^2$$
(105)
$$\mathrm{\Phi }_{AF}^{1,1}|S^{za}|\mathrm{\Phi }_{AF}^{1,1}=\frac{1}{2}\mathrm{\Phi }_{AF}^{1,1}|S^{za}+S^{zb}|\mathrm{\Phi }_{AF}^{1,1}=\frac{1}{2}𝒩_1^2$$
(106)
where in Eq.(104) we have used the identity Eq.(46). Finally we recall that the $`W=0`$ states $`|\mathrm{\Phi }_{AF}^{S,M_S}`$ are eigenstates of the square of the total spin operators of each sublattice $`S_𝒜^2(S^{za})^2+\frac{1}{2}(S^{+a}S^a+S^aS^{+a})`$ and $`S_{}^2(S^{zb})^2+\frac{1}{2}(S^{+b}S^b+S^bS^{+b})`$ with eigenvalue $`(\frac{N1}{2})(\frac{N1}{2}+1)=\frac{1}{4}(N^21)`$; hence

$$\underset{M_S=S}{\overset{S}{}}\mathrm{\Phi }_{AF}^{S,M_S}|(S^{za})^2|\mathrm{\Phi }_{AF}^{S,M_S}=\frac{N^21}{12}(2S+1)𝒩_S^2;𝒩_S^2=\mathrm{\Phi }_{AF}^{S,M_S}|\mathrm{\Phi }_{AF}^{S,M_S},M_S=S,\mathrm{},S.$$
(107)
Substitution of Eqs.(104-105-106-107) in Eq.(103) yields

$$G_{\mathrm{o}ff}=\frac{1}{N^2}\left\{\frac{N^21}{12}+\frac{1}{8}[(\gamma _g\gamma _{\overline{g}})^2+2\gamma _0^2]\frac{1}{3}\gamma _0^2\sqrt{\frac{N^21}{18}}\gamma _0(\gamma _g\gamma _{\overline{g}})\right\}$$
(108)
where we have taken into account Eq.(83) to evaluate the ratio $`𝒩_1/𝒩_0=\sqrt{(N^21)/3}`$.

In Fig.(5) we report the trend of $`G_{\mathrm{o}n}`$ and $`G_{\mathrm{o}ff}`$ versus $`U`$ for three different values of $`N=2,\mathrm{\hspace{0.33em}3},\mathrm{\hspace{0.33em}10}`$. According to the Shen-Qiu-Tian theorem, $`G_{\mathrm{o}n}`$ is always larger than zero while $`G_{\mathrm{o}ff}`$ is always negative. Now consider the ground state average of the square of the staggered magnetization operator

$$m_{AF}^2\frac{1}{|\mathrm{\Lambda }|}\mathrm{\Psi }_0|[\underset{x\mathrm{\Lambda }}{}ϵ(x)S_x^z]^2|\mathrm{\Psi }_0,ϵ(x)=1,1\mathrm{f}orx𝒜,.$$
(109)
The Shen-Qiu-Tian theorem implies that each term in the expansion of Eq.(109) is non-negative. We emphasize, however, that for $`|𝒜|=||`$ this does not imply that $`m_{AF}^2`$ is an extensive quantity! Remarkably, our results show that $`m_{AF}^2=G_0+(N1)G_{\mathrm{o}n}NG_{\mathrm{o}ff}`$ is extensive for any value of the on-site repulsion parameter $`U`$ and constitute the first example of antiferromagnetism in the ground state of an itinerant-electrons model of arbitrary size.

We also observe that the thermodynamic limit yields a non-vanishing result:

$$\underset{|\mathrm{\Lambda }|\mathrm{}}{lim}G_0=\frac{1}{4},\underset{|\mathrm{\Lambda }|\mathrm{}}{lim}G_{\mathrm{o}n}=\underset{|\mathrm{\Lambda }|\mathrm{}}{lim}G_{\mathrm{o}ff}=\frac{1}{12},$$
(110)
or

$$\underset{|\mathrm{\Lambda }|\mathrm{}}{lim}\underset{\beta \mathrm{}}{lim}\underset{h0}{lim}|G_{\mathrm{s}pin}(x,y;\beta ,h,|\mathrm{\Lambda }|)|=\frac{1}{12}\mathrm{f}orxy.$$
(111)
By comparing Eq.(111) with Eq.(30), we see that the $`lim_\beta \mathrm{}lim_{h0}`$ do not commute with the $`lim_{|\mathrm{\Lambda }|\mathrm{}}`$.

## VII Summary and Conclusions

We have introduced the CBG-Hubbard model and explicitly written the ground state wave function at half filling for any repulsion parameter $`U`$. The topology of the graph $`\mathrm{\Lambda }`$ allows $`(2N2)`$-body eigenstates $`|\mathrm{\Phi }_{AF}^{S,M_S}`$ of the Hamiltonian which are free of double occupation. They can be obtained projecting the determinantal state $`|\mathrm{\Phi }_{AF}`$ on a given spin-$`S`$ subspace. $`|\mathrm{\Phi }_{AF}`$ has the antiferromagnetic property, that is the map $`𝒜`$ is equivalent to a spin-flip.

We have shown that the ground state has a fixed number of particles in the zero-energy one-body shell $`𝒮_{hf}`$, independent of $`t`$ and $`U`$. This very remarkable Shell Population Rigidity holds in any finite-size system and it provides the closure of the equations. This implies extra-conservation laws in a suitable subspace of the full Hilbert space including the ground state. Only the singlet $`|\mathrm{\Phi }_{AF}^{0,0}`$ and the triplet $`|\mathrm{\Phi }_{AF}^{1,m}`$, $`m=0,\pm 1`$, are involved in the ground state expansion. Qualitatively, we may say that according to Eq.(86) the particles in the shell $`𝒮_{hf}`$ manage to avoid the double occupation.

We have calculated the spin-spin correlation function and shown that the ground state exhibits an antiferromagnetic order for any non-zero $`U`$ even in the thermodynamic limit. The kinetic term induces non-trivial correlations among the particles and an antiparallel spin configuration in the two sublattices turns out to be energetically favoured. Therefore, the model contains the basic ingredients to understand how the subtle competition between the delocalization induced by $`H_0`$ and the localization induced by $`W`$ may give rise to a magnetically ordered state even outside the strong coupling regime (where the Hubbard model is equivalent to the Heisenberg model). On the other hand, within the scheme proposed by van Dongen and Vollhardt, the model is well described by a paramagnetic Hamiltonian; the difference stems from taking the thermodynamic limit before the limit of zero Temperature, because then the kinetic term completly decouples from the interaction and loses any r$`\widehat{o}`$le in the physical response functions.

The present formalism lends itself to solve more realistic Hubbard Hamiltonians with an increasing number of negative and positive energy levels. Higher-spin projections $`|\mathrm{\Phi }_{AF}^{S,M_S}`$ are involved in the ground state expansions; the results will be published elsewhere. We are currently investigating group-theory aspects of this model and its extensions.

Finally, we recall that the standard Hubbard model on a $`N\times N`$ square lattice and periodic boundary conditions also has $`|\mathrm{\Phi }_{AF}^{S,M_S}`$-like eigenstates. It could be that the Shell Population Rigidity holds in this case too, but the proof is lacking.

## A Shell Population Rigidity

We have seen in Section V that the most part of the contributions which do not preserve the number of particles in the shell $`𝒮_{hf}`$ yield nothing. To prove the Shell Popoulation Rigidity property we need to show that the sequences like $`a^{}a^{}ga`$ and similar terms with $`b`$ instead of $`a`$ and/or $`g`$ replaced by $`\overline{g}`$ vanish or in formulæ

$$\underset{x𝒜}{}\underset{ijm}{}R_{i,x}R_{j,x}R_{m,x}a_i^{}a_j^{}a_{m\sigma }|\mathrm{\Phi }_{AF}^{S,M_S}=0$$
(A1)
and the like with $`a`$ replaced by $`b`$. This last spectacular cancellation can be seen by particle-hole transforming Eq.(A1). Under a particle-hole transformation $`c_{x\sigma }c_{x\sigma }^{}`$ and the $`W=0`$ states $`|\mathrm{\Phi }_{AF}^{S,M_S}|g^0\overline{g}^0|\mathrm{\Phi }_{AF}^{S,M_S}`$, modulo a phase factor. Hence, Eq.(A1) is equivalent to

$$\underset{x𝒜}{}\underset{ijm}{}R_{i,x}R_{j,x}R_{m,x}a_ia_ja_{m\sigma }^{}|\mathrm{\Phi }_{AF}^{S,M_S}=0.$$
(A2)
Let the spin of the operator $`a_{m\sigma }^{}`$ be up for the sake of definiteness (the same reasoning holds in the case of down spin). The l.h.s. of Eq.(A2) can be rewritten as

$$\underset{x𝒜}{}\underset{ijm}{}R_{i,x}R_{j,x}R_{m,x}[\delta _{i,m}a_m^{}a_i]a_j|\mathrm{\Phi }_{AF}^{S,M_S}=\underset{j}{}\left[\underset{x}{}\underset{i}{}R_{i,x}^2R_{j,x}\right]a_j|\mathrm{\Phi }_{AF}^{S,M_S}$$
(A3)
where we have exploited the lack of double occupation of the $`W=0`$ state, see Eq.(55). Therefore, we are left to prove that

$$\underset{x=1}{\overset{N}{}}\underset{i=1}{\overset{N1}{}}R_{j,x}R_{i,x}^2=0,j.$$
(A4)
We recall that $`R`$ is an $`(N1)\times N`$ rectangular matrix whose rows are orthonormal vectors which are orthogonal to the $`N`$dimensional vector $`(1,1,\mathrm{},1,1)`$. It is a simple exercise to verify that

$$R=\left(\begin{array}{cccccc}\frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}& 0& 0& \mathrm{}& 0\\ \frac{1}{\sqrt{6}}& \frac{1}{\sqrt{6}}& \frac{2}{\sqrt{3}}& 0& \mathrm{}& 0\\ :& :& :& :& :& :\\ \frac{1}{\sqrt{N(N+1)}}& \frac{1}{\sqrt{N(N+1)}}& \frac{1}{\sqrt{N(N+1)}}& \frac{1}{\sqrt{N(N+1)}}& \mathrm{}& \sqrt{\frac{N}{(N+1)}}\end{array}\right)R_{m,x}=\{\begin{array}{cc}0\hfill & x>m+1\hfill \\ \sqrt{\frac{m}{m+1}}\hfill & x=m+1\hfill \\ \frac{1}{m}\sqrt{\frac{m}{m+1}}\hfill & xm\hfill \end{array}$$
(A5)
is a correct choice and that Eq.(A4) is identically verified.

## B Proof of Eq.(77)

To prove Eq.(77) we shall use the definitions (36-42-43) for the triplet $`W=0`$ state $`|\mathrm{\Phi }_{AF}^{1,m}`$, $`m=0,\pm 1`$. The first term yields

$`{\displaystyle \underset{k}{}}(S_k^{+a}S_k^{+b})|\mathrm{\Phi }_{AF}^{1,1}`$ $`=`$ $`{\displaystyle \underset{P}{}}()^P{\displaystyle \underset{k}{}}(S_k^{+a}S_{P(k)}^{+b}){\displaystyle \underset{j}{}}\tau _{j,P(j)}^{(1)^{}}{\displaystyle \underset{ij}{}}\sigma _{i,P(i)}^{}|0=`$ (B1)
$`=`$ $`{\displaystyle \underset{P}{}}()^P{\displaystyle \underset{j}{}}(S_j^{+a}S_{P(j)}^{+b})\tau _{j,P(j)}^{(1)^{}}{\displaystyle \underset{ij}{}}\sigma _{i,P(i)}^{}|0+`$ (B2)
$`+`$ $`{\displaystyle \underset{P}{}}()^P{\displaystyle \underset{j}{}}\tau _{j,P(j)}^{(1)^{}}{\displaystyle \underset{kj}{}}(S_k^{+a}S_{P(k)}^{+b})\sigma _{k,P(k)}^{}{\displaystyle \underset{ij,k}{}}\sigma _{i,P(i)}^{}|0=`$ (B3)
$`=`$ $`\sqrt{2}{\displaystyle \underset{P}{}}()^P{\displaystyle \underset{j}{}}\sigma _{j,P(j)}^{}{\displaystyle \underset{ij}{}}\sigma _{i,P(i)}^{}|0+`$ (B4)
$``$ $`\sqrt{2}{\displaystyle \underset{P}{}}()^P{\displaystyle \underset{j}{}}{\displaystyle \underset{kj}{}}\tau _{j,P(j)}^{(1)^{}}\tau _{k,P(k)}^{(+1)^{}}{\displaystyle \underset{ij,k}{}}\sigma _{i,P(i)}^{}|0=`$ (B5)
$`=`$ $`\sqrt{2}(N1)|\mathrm{\Phi }_{AF}^{0,0}\sqrt{2}{\displaystyle \underset{P}{}}()^P{\displaystyle \underset{j}{}}{\displaystyle \underset{kj}{}}\tau _{j,P(j)}^{(1)^{}}\tau _{k,P(k)}^{(+1)^{}}{\displaystyle \underset{ij,k}{}}\sigma _{i,P(i)}^{}|0,`$ (B6)
while second term yields

$`{\displaystyle \underset{k}{}}(\widehat{n}_k^a\widehat{n}_k^a`$ $``$ $`\widehat{n}_k^b+\widehat{n}_k^b)|\mathrm{\Phi }_{AF}^{1,0}={\displaystyle \underset{P}{}}()^P{\displaystyle \underset{k}{}}(\widehat{n}_k^a\widehat{n}_k^a\widehat{n}_{P(k)}^b+\widehat{n}_{P(k)}^b){\displaystyle \underset{j}{}}\tau _{j,P(j)}^{(0)^{}}{\displaystyle \underset{ij}{}}\sigma _{i,P(i)}^{}|0=`$ (B7)
$`=`$ $`{\displaystyle \underset{P}{}}()^P{\displaystyle \underset{j}{}}(\widehat{n}_j^a\widehat{n}_j^a\widehat{n}_{P(j)}^b+\widehat{n}_{P(j)}^b)\tau _{j,P(j)}^{(0)^{}}{\displaystyle \underset{ij}{}}\sigma _{i,P(i)}^{}|0+`$ (B8)
$`+`$ $`{\displaystyle \underset{P}{}}()^P{\displaystyle \underset{j}{}}\tau _{j,P(j)}^{(0)^{}}{\displaystyle \underset{kj}{}}(\widehat{n}_k^a\widehat{n}_k^a\widehat{n}_{P(k)}^b+\widehat{n}_{P(k)}^b)\sigma _{k,P(k)}^{}{\displaystyle \underset{ij,k}{}}\sigma _{i,P(i)}^{}|0=`$ (B9)
$`=`$ $`2{\displaystyle \underset{P}{}}()^P{\displaystyle \underset{j}{}}\sigma _{j,P(j)}^{}{\displaystyle \underset{ij}{}}\sigma _{i,P(i)}^{}|0+`$ (B10)
$`+`$ $`2{\displaystyle \underset{P}{}}()^P{\displaystyle \underset{j}{}}{\displaystyle \underset{kj}{}}\tau _{j,P(j)}^{(0)^{}}\tau _{k,P(k)}^{(0)^{}}{\displaystyle \underset{ij,k}{}}\sigma _{i,P(i)}^{}|0=`$ (B11)
$`=`$ $`2(N1)|\mathrm{\Phi }_{AF}^{0,0}+2{\displaystyle \underset{P}{}}()^P{\displaystyle \underset{j}{}}{\displaystyle \underset{kj}{}}\tau _{j,P(j)}^{(0)^{}}\tau _{k,P(k)}^{(0)^{}}{\displaystyle \underset{ij,k}{}}\sigma _{i,P(i)}^{}|0.`$ (B12)
The third term can be computed following the same steps of Eq.(B6) and the final result is

$$\underset{k}{}(S_k^a+S_k^b)|\mathrm{\Phi }_{AF}^{1,1}=\sqrt{2}(N1)|\mathrm{\Phi }_{AF}^{0,0}\sqrt{2}\underset{P}{}()^P\underset{j}{}\underset{kj}{}\tau _{j,P(j)}^{(+1)^{}}\tau _{k,P(k)}^{(1)^{}}\underset{ij,k}{}\sigma _{i,P(i)}^{}|0.$$
(B13)
Denoting by $`[..]`$ the left hand side of Eq.(77), from Eqs.(B6-B12-B13) one obtains

$`[..]=3\sqrt{2}(N1)|\mathrm{\Phi }_{AF}^{0,0}\sqrt{2}{\displaystyle \underset{P}{}}()^P{\displaystyle \underset{j}{}}{\displaystyle \underset{kj}{}}[\tau _{j,P(j)}^{(1)^{}}\tau _{k,P(k)}^{(+1)^{}}\tau _{j,P(j)}^{(0)^{}}\tau _{k,P(k)}^{(0)^{}}+\tau _{j,P(j)}^{(+1)^{}}\tau _{k,P(k)}^{(1)^{}}]{\displaystyle \underset{ij,k}{}}\sigma _{i,P(i)}^{}|0.`$ (B14)
(B15)
Let us consider the first and the third term in the square bracket of the right hand side of Eq.(B15). We have

$`{\displaystyle \underset{P^{}}{}}()^P^{}\left[\tau _{j,P^{}(j)}^{(1)^{}}\tau _{k,P^{}(k)}^{(+1)^{}}+\tau _{j,P^{}(j)}^{(+1)^{}}\tau _{k,P^{}(k)}^{(1)^{}}\right]{\displaystyle \underset{ij,k}{}}\sigma _{i,P^{}(i)}^{}|0=`$ (B16)
$`={\displaystyle \underset{P}{}}()^P\left[\tau _{j,P(k)}^{(1)^{}}\tau _{k,P(j)}^{(+1)^{}}+\tau _{j,P(k)}^{(+1)^{}}\tau _{k,P(j)}^{(1)^{}}\right]{\displaystyle \underset{ij,k}{}}\sigma _{i,P(i)}^{}|0`$ (B17)
where the permutations $`P`$ and $`P^{}`$ satisfy

$$P^{}(k)=P(j),P^{}(j)=P(k),P^{}(i)=P(i)ij,k;$$
(B18)
hence $`()^P^{}=()^P`$. Substituting Eq.(B17) into Eq.(B15) and taking into account that

$$\tau _{j,P(k)}^{(1)^{}}\tau _{k,P(j)}^{(+1)^{}}+\tau _{j,P(j)}^{(0)^{}}\tau _{k,P(k)}^{(0)^{}}+\tau _{j,P(k)}^{(+1)^{}}\tau _{k,P(j)}^{(1)^{}}=\sigma _{j,P(j)}^{}\sigma _{k,P(k)}^{}$$
(B19)
one obtains

$$[..]=3\sqrt{2}(N1)|\mathrm{\Phi }_{AF}^{0,0}+\sqrt{2}\underset{P}{}()^P\underset{j}{}\underset{kj}{}\sigma _{j,P(j)}^{}\sigma _{k,P(k)}^{}\underset{ij,k}{}\sigma _{i,P(i)}^{}|0=\sqrt{2}(N^21)|\mathrm{\Phi }_{AF}^{0,0},$$
(B20)
that is, Eq.(77).

REFERENCES


format: text

text:
# The decoupling problem at RHIC

## 1 Introduction

Long-standing theoretical interests of the RMKI KFKI research group have been hydrodynamical, decoupling and transport phenomena. A recent theoretical puzzle is why the pion source extracted via HBT interferometry in Au+Au at $`\sqrt{s}=130A`$ GeV at RHIC has approximately identical ’out’ and ’side’ radii for average pair momenta $`K_{}<0.7`$ GeV. Most naively, this indicates very short pion emission duration, corresponding to almost instantaneous decoupling.

The observed ratio $`R_{out}/R_{side}0.91.1`$ contradicts theoretical expectations based on cylindrically symmetric and longitudinally boost invariant ideal hydrodynamics, indicating $`R_{out}/R_{side}`$ significantly larger than one . On the other hand, a hydrodynamical approach to the decoupling problem is inherently questionable because in the decoupling region the mean free path exceeds the side of the system, while hydrodynamics is only valid if the mean free path is zero. As a workaround, hydrodynamics is usually supplemented with an arbitrary freezeout prescription postulating sudden freezeout on a three-dimensional hypersurface, which is typically parametrized by a constant freezeout temperature or energy density. The unsatisfactory nature of this procedure has been realized recently and several improvements have been proposed .

A dynamical framework with a self-consistent freezeout is provided by classical transport theory, where the rate of interactions is given by the product of the density and the cross section $`\lambda ^1=n\sigma `$. As the system expands, it rarifies, and the rate of particle interactions drops accordingly. Because the theory is not limited to equilibrium dynamics, it also allows for systematic studies of dissipative effects. Dissipation is crucial in order to reproduce the observed saturation of elliptic flow at $`p_{}>2`$ GeV.

The equilibrium assumption was partly relaxed in an earlier study, where hydrodynamical evolution was only assumed until hadronization at which point particles were fed into a hadronic cascade. The calculation gave $`R_{out}/R_{side}1.21.4`$, still above one. On the other hand, it is not clear whether the hadronic densities were large enough to support equilibrium on the freezeout hypersurface, which would be necessary to justify the assumed equilibrium evolution until the end of hadronization. In fact, antiprotons are known to scatter only once on average in this approach, indicating the need for a nonequilibrium approach before hadronization.

In this study we investigate how nonequilibrium evolution in the partonic phase influences the pion distributions at freezeout. Utilizing the MPC numerical technique, we compute the covariant freezeout distributions for a variety of initial conditions for Au+Au at $`\sqrt{s}=130A`$ GeV based on perturbative QCD, gluon saturation models, or a combination of both. Varying the transport opacity, which is an effective measure of the degree of collectivity, we study under what conditions $`R_{out}<R_{side}`$ is possible.

## 2 Two-particle HBT interferometry

Hanbury-Brown and Twiss interferometry is a truly remarkable technique to deduce space-time geometrical information purely via measuring momentum correlations. We review here only a few essential features and refer the reader to Refs.  and references therein.

Detailed studies of pion production by an ensemble of classical currents showed that for a chaotic boson source $`\rho (\stackrel{}{x},t)`$, in the limit of a large particle numbers, and also neglecting final state interactions, the two-particle correlation function is given by the space-time Fourier transform of the source as

$$C(q,K)C(p_1,p_2)\frac{N(p_1,p_2)}{N(p_1)N(p_2)}=1+\left|d^4x\rho (\stackrel{}{x},t)e^{iqx}\right|^2.$$
(1)
Here $`q^\mu p_1^\mu p_2^\mu `$ is momentum difference, while $`K^\mu (p_1^\mu +p_2^\mu )/2`$ is the average momentum of the boson pair.

Unfortunately, Eq. (1) cannot be uniquely inverted to obtain the source distribution, because the phase information is lost and, in addition, only a three-dimensional projection of the four-dimensional correlation function is measurable because $`q^\mu `$ is constrained by

$$q^\mu K_\mu =(m_1^2m_2^2)/2=0(\mathrm{for}\mathrm{identical}\mathrm{particles}).$$
(2)
Despite the inversion problem, HBT measurements provide strong constraints on the freezeout distributions of particles and therefore on the possible dynamical scenarios in heavy-ion collisions.

Conventionally, $`C(\stackrel{}{q},K)`$ is expressed in terms of the ’out-side-long’ variables, $`q_O`$, $`q_S`$, and $`q_L`$. These are defined in the longitudinally boosted LCMS reference frame, where $`K_z=0`$ with $`z`$ or ’long’ being the beam direction, ’out’ the direction of $`\stackrel{}{K}_{}`$, and ’side’ the one orthogonal to both. In this frame, $`K_{LCMS}^\mu =(\stackrel{~}{K}^0,K_{},0,0)`$, $`\stackrel{}{q}_{LCMS}=(q_O,q_S,q_L)`$, $`x_{LCMS}^\mu (\stackrel{~}{t},x_O,x_S,x_L)`$ and Eq. (2) gives

$$q^\mu x_\mu =q_O(x_O\stackrel{~}{t}K_{}/\stackrel{~}{K}^0)q_Sx_Sq_Lx_L.$$
(3)
Clearly, the measurable projected correlation function is completely blind to the dependence of the source function on the fourth orthogonal coordinate $`\stackrel{~}{t}+x_OK_{}/\stackrel{~}{K}^0`$.

Experimentally, the measured correlation function (after correcting for final state interactions) is fitted with a Gaussian, which for central collisions ($`b=0`$) and at midrapidity is constrained by symmetry to be of the form

$$C(\stackrel{}{q},K)=\lambda (K)\mathrm{exp}\left[q_O^2R_O^2(K)q_S^2R_S^2(K)q_L^2R_L^2(K)\right].$$
(4)
Here $`R_O`$, $`R_S`$ and $`R_L`$ are the ’out’, ’side’, and ’long’ HBT radii. For a perfectly Gaussian source, the correlation function is Gaussian and

$`R_O^2(K)`$ $`=`$ $`\mathrm{\Delta }x_O^2_K+v_{}^2\mathrm{\Delta }\stackrel{~}{t}^2_K2v_{}\mathrm{\Delta }x_O\mathrm{\Delta }\stackrel{~}{t}_K`$
$`R_S^2(K)`$ $`=`$ $`\mathrm{\Delta }x_S^2_K,\mathrm{and}R_L^2(K)=\mathrm{\Delta }x_L^2_K,`$ (5)
where $`v_{}K_{}/\stackrel{~}{K}^0`$. Thus $`R_S`$ and $`R_L`$ have simple geometric interpretation as the ’side’ and ’long’ widths of the distribution function, while $`R_O`$ is a mixture of the ’out’ width, time spread, and the $`x_O\stackrel{~}{t}`$ correlation.

## 3 Covariant parton transport theory

### 3.1 Transport equation

We consider here, as in Refs. , the simplest but nonlinear form of Lorentz-covariant Boltzmann transport theory in which the on-shell phase space density $`f(x,\stackrel{}{p})`$, evolves with an elastic $`22`$ rate as

$$p_1^\mu _\mu f_1=S(x,\stackrel{}{p}_1)+\underset{2}{}\underset{3}{}\underset{4}{}\left(f_3f_4f_1f_2\right)W_{1234}\delta ^4(p_1+p_2p_3p_4).$$
(6)
Here $`W`$ is the square of the scattering matrix element, the integrals are shorthands for $`\underset{i}{}\frac{gd^3p_i}{(2\pi )^3E_i}`$, where $`g`$ is the number of internal degrees of freedom, while $`f_jf(x,\stackrel{}{p}_j)`$. The initial conditions are specified by the source function $`S(x,\stackrel{}{p})`$. For our applications below, we neglect quark degrees of freedom and interpret $`f(x,\stackrel{}{p})`$ as describing an ultrarelativistic massless gluon gas with $`g=16`$ (8 colors, 2 helicities).

The elastic gluon scattering matrix elements in dense parton systems were modeled with the isotropic form $`d\sigma _{el}/dt=\sigma _0(s)/s`$. The simplified angular dependence is justified by our previous study, where we showed that the covariant transport solutions do not depend explicitly on the differential cross section but only on the transport opacity

$$\chi \frac{\sigma _{tr}}{\sigma _{el}}n\sigma _{tr}𝑑z\rho \left(𝐱_0+z\widehat{𝐧},\tau =\frac{z}{c}\right).$$
(7)
Here

$$\sigma _{tr}(s)𝑑\sigma _{el}\mathrm{sin}^2\theta _{cm}=𝑑t\frac{d\sigma _{el}}{dt}\frac{4t}{s}\left(1\frac{t}{s}\right)$$
(8)
is the transport cross section (in our case, $`\sigma _{tr}=2\sigma _0/3`$), while $`n`$ is the average number of scatterings a parton undergoes. We neglected the weak logarithmic energy dependence of the total cross section and took a constant $`\sigma _0`$ for simplicity. For a fixed nuclear geometry, a given transport opacity $`\chi `$ represents a whole class of initial conditions and partonic matrix elements, which is demonstrated by the approximate proportionality $`\chi \sigma _{tr}dN_g(\tau _0)/d\eta `$.

We solved Eq. (6) numerically via the MPC parton cascade algorithm , which utilizes the parton subdivision technique to obtain the correct, covariant solutions to the transport equation. Parton subdivision is essential to eliminate numerical artifacts caused by acausal (superluminal) propagation due to action at a distance. For initial partonic densities expected at RHIC, the severe violation of Lorentz covariance in the naive cascade algorithm that employs no subdivision artificially reduces elliptic flow and heats up the $`p_{}`$ spectra.

### 3.2 Initial conditions

We modeled central Au+Au collisions at RHIC with three different classes of initial conditions, corresponding to different assumptions on the initial transverse momentum distribution and density distribution of partons. In all three cases, the evolution started from a longitudinally boost invariant Bjorken tube at proper time $`\tau _0=0.1`$ fm/$`c`$, with isotropic momentum distribution in the local rest frame and with uniform pseudorapidity $`\eta 1/2\mathrm{log}[(t+z)/(tz)]`$ distribution between $`|\eta |<5`$.

The first class was essentially the same as the minijet initial conditions in Ref. . The initial transverse density distribution was proportional to the binary collision distribution for two Woods-Saxon distributions, while the $`p_{}`$ distribution was a fit to the minijet distribution predicted by HIJING (without shadowing and jet quenching). To include the transport effect of the soft partonic component present in HIJING in addition to the minijets, in the present study we increased the parton density by a factor five, to yield $`dN_g(\tau _0)/d\eta =1000`$.

The second class was based on gluon saturation models. Here we assumed a uniform Bjorken cylinder with a radius $`R_0=6`$ fm, $`dN_g/d\eta =1000`$ and a constant $`d^2N/p_{}dp_{}`$ distribution that vanishes for $`p_{}>Q_{sat}=1.2`$ GeV.

The third class was a combination of the first two. We took the minijet part with the original HIJING normalization ($`dN/d\eta =210`$) and added four times as many saturated soft gluons.

## 4 Results

Using MPC, we computed the freezeout distribution numerically for central Au+Au at 130$`A`$ GeV as a function of the average pair momentum and transport opacity, for all three classes of initial conditions.

The freezeout distribution $`d^4N/d^4x`$ was defined as the distribution of space-time coordinates for the point of last interaction of the test particles, with the strong assumption that this point is not affected by hadronization. We furthermore neglected resonance contributions to the pion yield. The same two hadronization schemes, parton-hadron duality and independent fragmentation, were applied as in Ref. .

Of our primary interest were the integrated distributions $`d^2N/rdrd\stackrel{~}{t}`$, the equivalent of the usual hydrodynamical freezeout curves in the $`r\tau `$ plane, and $`d^2N/dx_Odx_S`$, which is the source projected onto the ’out-side’ plane. We worked in the approximation $`p_{1,}p_{2,}K_{}`$, which is valid up to corrections of $`O(|\stackrel{}{q}|/|\stackrel{}{K}|)`$. To get a simple estimate of the radii, we used Eq. (5) even though our sources are non-Gaussian. A more detailed analysis of the radii based on the momentum-space correlation function is in progress.

Figure 1 shows the freezeout distribution in $`r\stackrel{~}{t}`$ coordinates ($`r\sqrt{x_O^2+x_S^2}`$) as a function of the pair momentum and transport opacity for the minijet initial condition with hadronization via parton-hadron duality. Unlike the sharp freezeout imposed in hydrodynamical models, the transport theory freezeout is a continuous, evaporation-like process. For a given (nonzero) opacity, the larger the $`p_{}`$ of the particle, the earlier it decouples and the closer it is to the surface of the nuclei. Thus, low-$`p_{}`$ partons freeze out in the center at late times, while high-$`p_{}`$ ones escape from the surface at early times. Furthermore, a larger opacity increases decoupling times, especially for low-$`p_{}`$ particles. We observed the same features for the other two classes of initial conditions.

Figure 2 shows the shape of the source in $`x_Ox_S`$ coordinates as a function of the pair momentum and transport opacity for the minijet initial condition with hadronization via parton-hadron duality. Results for the saturation and the ’combined’ initial condition were very similar. The black crosses are positioned at the first moment of the distributions $`(x_O,x_S)`$, while their extension corresponds to the square root of the second moments $`(\sqrt{\mathrm{\Delta }x_O^2},\sqrt{\mathrm{\Delta }x_S^2})`$. As the opacity increases, the source shrinks in the ’out’ direction and the higher the $`p_{}`$, the more pronounced the reduction of $`\mathrm{\Delta }x_O^2`$ is. From Eq. (5) it is clear that this will help reduce $`R_O`$, at least if the other two terms do not change. On the other hand, we find that the source size in the ’side’ direction is approximately independent of both the opacity and $`p_{}`$. Thus, surprisingly, $`R_S`$ is unaffected by the strong collective dynamics.

The mean $`x_O`$ does not coincide with the maximum of the distributions showing that the distributions are asymmetric in the ’out’ direction, i.e., they are clearly non-Gaussian. In fact, the change in the shape of contour lines with increasing opacity is much more pronounced than the slow decrease of the second moment $`\mathrm{\Delta }x_O^2`$ because there is a long tail in the negative ’out’ direction. We expect that $`R_O`$ extracted via a least-squared Gaussian fit would turn out to be somewhat smaller than from Eq. (5).

The left column in Fig. 3 shows the dependence of the $`R_O/R_S`$ and $`R_L/R_O`$ ratios on transport opacity and $`K_{}`$ for the minijet initial condition. For $`\chi =0`$, both $`R_O`$ and $`R_S`$ are determined by the nuclear size to be $`3`$ fm, therefore $`R_O/R_S=1`$. However, as the opacity increases, the ratio increases above one for $`K_{}<K_c`$, while decreases below one for $`K_{}>K_c`$, where the location of the turning point $`K_c1.32`$ GeV depends on the hadronization scheme.

The naive $`R_O^2=R_S^2+\mathrm{\Delta }t^2`$ expectation does not hold now because in the LCMS frame the source is not cylindrical (see Fig. 2) and there are dynamical correlations. Figure 4 illustrates how the different terms in Eq. (5) contribute to $`R_O`$ as a function of $`K_{}`$. For all three classes of initial conditions, the time spread $`\mathrm{\Delta }\stackrel{~}{t}^2`$ is comparable to the spatial width of the source $`3`$ fm, which has the effect to increase $`R_O`$. However, this effect is compensated by the large positive $`x_O\stackrel{~}{t}`$ correlation term, reducing $`R_O`$ below $`R_S`$.

The right column in Fig. 3 shows a comparison to the HBT radii measured at RHIC. In the transverse opacity range $`\chi 08`$ we studied, the transport theory results are smaller than the observed $`R_O`$ and $`R_L`$. This is in sharp contrast to ideal hydrodynamics, which overpredicts both radii. The monotonic dependence of $`R_O`$ and $`R_L`$ on opacity indicates that a much better agreement with the data is possible at a larger but finite opacity $`\chi >10`$.

In Fig. 3, the reason for the anomalous small $`R_L`$ for zero opacity is instantaneous freezeout at the formation time $`\tau =\tau _0`$. Since $`R_L^2\tau ^2[\mathrm{\Delta }(\eta y)]^2`$, $`R_L`$ depends on the decoupling time and the strength of the $`\eta y`$ correlation. For our thermally correlated initial condition $`[\mathrm{\Delta }(\eta y)]^2T/m_{}`$ and $`\tau _0=0.1`$ fm$`/c`$, which for $`\chi =0`$ gives $`R_L<0.5`$ fm for all $`K_{}`$ bins we studied.

However, as the transport opacity increases, $`R_L`$ becomes much larger because the decoupling time increases as is evident from Fig. 1. The largest increase $`\tau /\tau _0R/\tau _050`$ is for low-$`p_{}`$ partons, which freeze out latest. Thus the observed $`R_L(K_{})`$ carries important constraints on the dynamics of freezeout. For example, perfect inside-outside correlation, i.e., $`\eta =y`$ as assumed in classical Yang-Mills approaches cannot be reconciled with the RHIC data because it gives $`R_L=0`$.

Finally, the biggest puzzle in Fig. 3 is why $`R_S(K_{})const3`$ fm contrary to the data and furthermore independent of the transport opacity. This suggests that $`R_S`$ is independent of the degree of collective dynamics, i.e., it is determined solely by the initial nuclear geometry. The same problem with $`R_S`$ being too small and insensitive to the dynamical assumptions has been observed in hydrodynamical calculations as well.

## Acknowledgements

D. M. would like to acknowledge enlightening discussions with Tamás Csörgő and Ulrich Heinz.

This work was supported by the Director, Office of Energy Research, Division of Nuclear Physics of the Office of High Energy and Nuclear Physics of the U.S. Department of Energy under contract No. DE-FG-02-93ER-40764.