text:
# Thermodynamics of a 4-site Hubbard model by analytical diagonalization

## I Introduction

The Hubbard model  with the restraint of electron hopping to nearest neighbours only is certainly one of the simplest many particle models. It was utilized for the description of correlated electrons in the intermediate and strong coupling regimes. First it was used in chemistry to describe unsaturated hydrocarbons . In solid state physics it was introduced to describe the behaviour of narrow-band transition metal compounds, with the main emphasis given to the explanation of the metal-insulator transition  and the itinerant ferromagnetism . Later on it played a central role in the heavy fermion context, and especially the discovery of high Tc superconductivity in the cuprates renewed the interest, since its involved spin dynamics is widely considered as a good candidate to understand at least qualitatively the pairing mechanism in the copper oxide planes, being common to all high Tc superconductors. It is clear that after about four decades of research there exist a vast literature and meanwhile also some good reviews dealing with different aspects of the model , where good reasons to deal with that model are given, so I refer to them. In contrast to the large amount of work spent to that model, there is astonishing few exact knowledge. The initial hope that by reducing the description of the correlated electrons to a one-band model with hopping restricted to next nearest neighbours and correlation restricted to the on-site term a model will result, which can be solved analytically, has not yet been fulfilled in general. Nevertheless, Lieb and Wu  calculated the ground-state-energy of the one dimensional model with a filling of one electron per lattice site by help of the Bethe ansatz. Starting from this paper there were later on attempts to get the complete spectrum and thus the thermodynamics for the half filled model and during the last decade also for arbitrary electron filling, at least approximately. The state of art for the one dimensional model was recently summarized by Deguchi et. al. . All theories starting from the Bethe ansatz result in a set of coupled nonlinear equations, the Lieb-Wu equations, which have to be solved numerically for every parameter set. Thus all studies of physical quantities like susceptibilities or excitation spectra resting upon the solution of the Lieb-Wu equations have to be done numerically, which often limits the application. The transfer matrix method introduced by Shastry  supplies an alternate solution, which is in principle exact. The drawback for calculating physical quantities is that the lowest eigenvalue of the transfer matrix has to be calculated, what requires again a considerable numerical effort. For the cases where results exist, they were shown to be equivalent to the Bethe ansatz method . Another way to achieve exact information is the numerical treatment of finite systems. The techniques by far most intensively used for this task are the Lanczos algorithm  and Monte Carlo methods . Lanczos calculations. Albeit giving essentially exact results for ground-state and low lying eigenvectors and single particle excitation spectra both methods suffer from the relative small cluster size treatable with the contemporary computer facilities, since the basis set of vectors grows exponentially with the system size. The main drawback in the MC calculations for fermions are the relatively large uncertainties due to the so called sign problem, although there exist a enormous number of derivate methods developed to overcome this problem. Direct exact diagonalization, e.g. with the householder method, is restricted to much smaller problems and exact is then always meant numerically. Along this way it is in principle impossible to get closed analytical expressions in dependence of the model parameters. In most cases a numerically exact calculation may be enough, but there are problems which need the explicit analytical dependence of the eigenvalues and the coefficients of the eigenvectors with respect to a suitable basis set. The problem of determining the crossing temperature of the high temperature specific heat, a problem we address in the third section, is such a case, since it invokes the third derivative of the thermodynamical potential. Another problem, where the analytical knowledge of the spectrum is inevitable, is the term-by-term comparison of series expansions with exact solutions, and, last but not least, the study of the influence of a multitude of external parameters, (e.g. in the model under consideration the ratio t/U between the hopping energy and the on-site correlation energy, the electron density n, the temperature, and the external applied fields) is considerably simplified, if the eigensystem of an Hamiltonian is known analytically, especially, if the calculation aims at the grand canonical partition sum, where the chemical potential has to be determined from the particle number. To get the explicit dependencies on the model parameters it is necessary to diagonalize the Hamiltonian in an analytical manner, what restricts the calculation to very small systems of course. We have chosen a four-point square for this task for the following reasons: (i) It is the earliest  and most often studied non-trivial object in numerical diagonalization. (ii) There were some attempts in the nineties which remained incomplete, since they discussed the half filled case only, where the eigenvalues were given. The authors speculate about hidden symmetries, since they were not able to resolve all degeneracies . (iii) It is a special case of the N-site rings, which are subject of the above cited Bethe ansatz and quantum transfer matrix methods, and might be useful as a proving tool. (iv) In the following sense the four point Hubbard square may serve as a crude model to the square lattice: If you neglect every second hopping in x- and y-direction respectively you will get a system of N/4 independent 4-site-clusters being the topic of this paper, see Fig. 1. The model is thus a little bit better than the so-called atomic limit, were every hopping is neglected. The main fault of this cluster-gas model is that there is no direct electron transfer between the clusters. This is in some sense corrected by applying the grand canonical distribution, since it reflects the situation where the subsystems are in contact with a particle reservoir allowing for particle number fluctuations, thus establishing an indirect particle transfer via the bath. Very recently the spectral function of the square-lattice Hubbard model was calculated by a cluster-perturbation method, which used the square-clusters as unperturbed system and the inter-square-hopping as perturbation , albeit they diagonalized the cluster numerically. (iv) The rigorous solution of any strong interacting many particle system is a value independent of its applicability to realistic systems, at least for pedagogical reasons. (v) The limited resources of the available computers, we used for symbolic formula handling.

The paper is organized as follows. In the next two sections the diagonalization procedure is explained and the many-particle spectrum is discussed. In the fourth section we give the complete thermodynamics, based on the grand canonical potential. The main focus lies on the dependence of the physical quantities on the electron density, whereas the comparison of our results to the literature for the half-filled case serves as independent proof. In the closing section we discuss, what we could learn from our analytical solution.

## II The analytical diagonalization procedure

We are interested in the eigenvalue sequence for the operator

$`𝐇`$ $`=`$ $`t{\displaystyle \underset{ij\sigma }{}}𝐜_{i\sigma }^+𝐜_{j\sigma }+{\displaystyle \underset{i\sigma }{}}\left({\displaystyle \frac{U}{2}}𝐧_{i\sigma }𝐧_{i\sigma }(\mu +\sigma h)𝐧_{i\sigma }\right)`$ (1)
Here $`𝐜_{i\sigma }^+`$ and $`𝐜_{i\sigma }`$ are the creation and destruction operators in Wannier representation. The chemical potential $`\mu `$ and the magnetic field h in z-direction, which we choose to be normal to the square, are introduced to take into account the effects of doping and applying external magnetic fields. The related spectrum we will call grand-canonical, due to its relation to the weights within the grand canonical distribution. The lowest level of this operator determines the state with the highest weight, becoming one for $`T0`$, consequently, we will call it the grand canonical groundstate, which of course may be different for different values of $`\mu `$ and h. The many particle Hilbert space is the direct product of the four single-site Hilbert spaces, each containing four states, i.e. the empty state, the two single occupied states with spin up and spin down, and the double occupied state. Thus we have 256 states. For example, $`|u,0,2,d`$ indicates the state which has a spin-up electron at the first point, no electron at the second position, two electrons at the third, and one spin-down electron at the last position. Of course it is impossible to diagonalize the resulting 256x256 matrix in an analytic manner straight forwardly, instead one has to use all the known symmetries or the related conservation laws to reduce the matrix to block diagonal form with all blocks being at least 4x4.

The Hubbard model is very good suited for the task of analytical diagonalization by help of the conservations laws, since it contains several conserved quantities which are easy to handle. First of all, there are the two commuting SU algebras of the total spin and the total pseudospin  providing us with four operators commuting with each other and with the hamiltonian, i.e.

$`[𝐒^2,𝐇]=[𝐒_z,𝐇]=[𝐑^2,𝐇]=[𝐑_z,𝐇]=0`$ (2)
with

$`𝐒^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(𝐒^+𝐒^{}+𝐒^{}𝐒^+\right)+𝐒_z^2`$ (3)
$`𝐒^+`$ $`=`$ $`{\displaystyle \underset{i}{}}𝐜_i^+𝐜_i,`$ (4)
$`𝐒^{}`$ $`=`$ $`{\displaystyle \underset{i}{}}𝐜_i^+𝐜_i,`$ (5)
$`𝐒_z`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{1}{2}}\left(𝐧_i𝐧_i\right)`$ (6)
$`𝐑^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(𝐑^+𝐑^{}+𝐑^{}𝐑^+\right)+𝐑_z^2`$ (7)
$`𝐑^+`$ $`=`$ $`{\displaystyle \underset{i}{}}(1)^i𝐜_i^+𝐜_i^+,`$ (8)
$`𝐑^{}`$ $`=`$ $`{\displaystyle \underset{i}{}}(1)^i𝐜_i𝐜_i,`$ (9)
$`𝐑_z`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{1}{2}}\left(𝐧_i+𝐧_i1\right)`$ (10)
We mention here that the total pseudospin is conserved only if the hopping matrix elements fulfill the condition

$`t_{ij}\left((1)^i+(1)^j\right)=0`$ (11)
what is the case if one restricts e.g. the hopping to next neighbours as we supposed in eq. (1). The expectation value of $`𝐑_z`$ is nothing but the difference of the model from the half filled case, or in the context of high-$`T_c`$ superconductivity it is the doping of the model. It can be controled by adjusting the chemical potential $`\mu `$ of course. Dealing with a special cluster one can use the symmetries of it. In our case the actual symmetry is C4v. Indeed there are two works on the four-site square, however restricted to the half filled case, which tried to use these symmetries. Within the first  the authors start from a parametric expression of the Hubbard interaction, which they took obviously from Ref. . But instead of using both spin and pseudospin symmetry they used the spin symmetry first and then the point symmetry of the four site cluster. This way the biggest block in the hamiltonian remains 7x7, what in general prevents an analytical solution. Nevertheless they found the roots. More elegantly the same task was done in Ref. , where both spin and pseudospin were used first and the translation invariance afterwards. This is enough to bring the hamiltonian matrix being 70x70 in the case of the half filled model into the wanted block form with none of the remaining blocks being bigger than 4. In both papers the authors mentioned a degeneration of the energy levels, which they ascribed to some hidden symmetries, a point which will be clarified within the present work.
In this paper we use the same method for arbitrary filling. Due to the commutation of $`𝐑_z`$ and $`𝐒_z`$ with the hamiltonian, we started with the basis of the Hilbert subspace with a given value of $`𝐑_z`$ or, equivalently, with a given number of electrons, and a given value of $`𝐒_z`$. Afterwards we calculated the matrix of $`𝐑^2`$ and diagonalized it resulting in a common eigensystem of $`𝐑_z,𝐒_z,𝐑^2`$. The next step was the diagonalization of $`𝐒^2`$ in every subspace for given eigenvalues $`𝐑_z,𝐒_z,𝐑^2`$ thus getting common eigenvectors to $`𝐑_z,𝐒_z,𝐑^2,𝐒^2`$. The last step, if necessary, was the calculation of the matrix of the translation operator, which fulfills $`𝐔^4|\mu =|\mu `$ for a four-site ring and therefore has the eigenvalues 1, i, -1, and -i. We took the operator $`𝐔`$ in the form given in Ref. .

$`𝐔`$ $`=`$ $`𝐔_{}𝐔_{}`$ (12)
$`𝐔_\sigma `$ $`=`$ $`𝐔_{12\sigma }𝐔_{23\sigma }𝐔_{34\sigma }`$ (13)
$`𝐔_{ij\sigma }`$ $`=`$ $`1(𝐜_{i\sigma }^+𝐜_{j\sigma }^+)(𝐜_{i\sigma }𝐜_{j\sigma })`$ (14)
The operators $`𝐔`$, $`𝐔^\mathrm{𝟐}`$, $`𝐔^3`$ are operator representations of the group elements C4z, C2z, and C$`{}_{}{}^{1}{}_{4z}{}^{}`$, which form together with the 1 standing for E an Abelian subgroup. The relation to the conserved momentum is obvious. The operator $`𝐔_{ij\sigma }`$ permutes two electrons and is therefore well-suited to construct the missing group elements. We find for

$`𝐏(IC_{2x})`$ $`:=`$ $`𝐔_{14}𝐔_{23}𝐔_{12}𝐔_{34}`$ (15)
$`𝐏(IC_{2y})`$ $`:=`$ $`𝐔_{14}𝐔_{24}𝐔_{12}𝐔_{34}`$ (16)
$`𝐏(IC_{2a})`$ $`:=`$ $`𝐔_{13}𝐔_{13}`$ (17)
$`𝐏(IC_{2b})`$ $`:=`$ $`𝐔_{24}𝐔_{24}`$ (18)
Since these operators neither commute with each other nor with the $`𝐔`$-operators we have to construct new operators commuting with all the other operators. The result is

$`𝐏^{++}`$ $`=`$ $`𝐏(IC_{2x})𝐏(IC_{2y})`$ (20)
$`+𝐏(IC_{2a})𝐏(IC_{2b})`$
$`𝐏^+`$ $`=`$ $`𝐏(IC_{2x})𝐏(IC_{2y})`$ (22)
$`𝐏(IC_{2a})+𝐏(IC_{2b})`$
$`𝐏^+`$ $`=`$ $`𝐏(IC_{2x})+𝐏(IC_{2y})`$ (24)
$`+𝐏(IC_{2a})𝐏(IC_{2b})`$
$`𝐏^{}`$ $`=`$ $`𝐏(IC_{2x})+𝐏(IC_{2y})`$ (26)
$`𝐏(IC_{2a})+𝐏(IC_{2b})`$
Grosse  constructed further conserved quantities for the one-dimensional system with periodic boundary conditions, which are in general highly nontrivial due to their dependence on the interaction strength. The operator given in the first theorem of Ref.  reads

$`𝐊^1`$ $`=`$ $`𝐊_0^1+{\displaystyle \frac{U}{t}}𝐊_1^1`$ (27)
$`𝐊_0^1`$ $`=`$ $`{\displaystyle \underset{i,\sigma }{}}\left(𝐜_{i+2\sigma }^+𝐜_{i\sigma }𝐜_{i\sigma }^+𝐜_{i+2\sigma }\right)`$ (28)
$`𝐊_1^1`$ $`=`$ $`{\displaystyle \underset{i,\sigma }{}}(𝐜_{i+1\sigma }^+𝐜_{i\sigma }𝐜_{i\sigma }^+𝐜_{i+1\sigma })\times `$ (30)
$`\times \left(𝐜_{i,\sigma }^+𝐜_{i,\sigma }+𝐜_{i+1,\sigma }^+𝐜_{i+1,\sigma }\right)`$
For our four-site ring $`𝐊_0^1`$ is zero. In Ref.  there another operator $`𝐊^2`$ is given, which does not give new results due to our short periodicity length. As mentioned above, for our small system it is not necessary to use $`𝐊^1`$, but, since it commutes with all the other operators, it must have the same eigensystem like the hamiltonian. Since it is an antihermitian operator the eigenvalues have to be zero or pure imaginary. Due to its rather complicated character we used it among others as a powerful tool for checking the eigensystem. All the calculations were done by help of algebraic programming. The resulting eigenvectors and the related eigenvalues were proved by help of the eigenvalue equation. Also the matrices for the special case $`N_e=4`$ are compared to Ref. , where they are given. The matrices $`𝐇_{0,0,1}`$ and $`𝐇_{1,0,1}`$  differ in the sign for two off-diagonal matrix elements, but the eigenvalues are the same and the matrices $`𝐇_{0,0,1}`$ and $`𝐇_{0,1,1}`$ are interchanged. Since the resulting spectrum for the half filled case is not altered by these differences, the conclusions drawn in Ref.  remain true. The complete set of the eigenvalues and eigenvectors $`|m_r,r(r+1),m_s,s(s+1),u,k_{11},p_{++},p_+`$, with $`m_r,m_s,s,r,u,k_{11},p_{++},p_+`$ indicating the eigenvalues of the operator-set $`𝐑_z,𝐑,𝐒_z,𝐒^2,𝐔,𝐊_{11},𝐏^{++},𝐏^+`$, in its analytical dependence on the model parameters are given in the appendix. Inspecting the spectrum carefully we found all the states classified by the eigenvalues of the conserved quantities in a unique manner. Thus we clarified the reason of the mysterious degeneracies mentioned in Ref. . Also in Ref.  hidden degeneracies remained, although not mentioned by the authors. For instance the two degenerate states No. 95 and 96 have identical eigenvalues of the conserved quantities $`𝐑_z,𝐑,𝐒_z,𝐒^2,𝐔`$, which was the set exploited in Ref. . This degeneracy is due to the conserved Grosse operator $`𝐊_{11}`$. Heilmann and Lieb  showed for the benzene that there are a multitude of level crossings. Since they have constructed all correlation-independent symmetries they ascribe this fact to the existence of unknown U-dependent symmetries. We believe, that the Grosse operators  are candidates for the missing symmetry, since they were left out in Ref. . For instance, it is obvious that for a six-site ring the operator $`𝐊^1`$ is U-dependent. Nevertheless, in our special case it is not, since the operator $`𝐊_0^1`$ vanishes. Thus within the four-site model $`𝐊_1^1`$ represents a further U-independent symmetry. By this way we remark that there is no contradiction to the non-crossing rule, since all the eigenstates differ in the symmetry quantum numbers of mutual commuting operators, a question raised in Ref. .

## III The grand canonical eigensystem

### A The low lying states

In the following discussion of the spectrum we always mean the grand canonical levels, which are numbered for referencing according to the appendix, exceptions are explicitly indicated. The well-known particle-hole symmetry allows to fix the electron number to four if we choose $`\mu `$=U/2. Looking at the lower part of the spectrum of the half filled model, as plotted in Fig. 2, we find the state No. 111 being the groundstate. As already mentioned in Ref.  this is both a spin and a pseudospin singlet (mr=0 and ms=0) in accordance with the statements of Lieb’s positive-U theorem  and the lemma of Shen and Qiu . In Fig. 3 we show the results for different temperatures. As the first exited states we find a spin triplet (states No. 102, 128, 154), regarding the pseudospin it remains a singlet. The next group consists again of a singlet-triplet type sequence, where the singlet differs from the groundstate by its eigenvalue of the translation operator only (it is 1 in comparison with U=-1 of the groundstate), the triplet states have pure imaginary U eigenvalues and are therefore twice degenerated, in agreement with a statement given in Ref. . It is interesting, that the next exited states deviate in the charge. Being both spin and pseudospin doublets they are eightfold degenerated, if we do not apply neither an external magnetic field nor a deviation of the chemical potential from U/2. Since we calculated Fig. 2 with a small magnetic field, the spin degeneracy is lifted. We also find a pseudospin triplett (not shown in Fig. 2) among the low lying states, i.e. the states with No. 22 (n=2), No.131 (n=4) and No. 232 (n=6). A small deviation of $`\mu `$ from U/2 will lift the degeneracy of these states, since $`\mu `$ couples to $`𝐑_z`$ in the same manner as the magnetic field couples to $`𝐒_z`$. But in between the pseudospin singlet and triplet states we find a further pseudospin-singlet-spin-quintet-structure and pseudospin-doublet-spin-quartet structure. It was quite surprising for us, that the low lying spectrum of the half-filled model contains states with electron numbers differing from four by one and two. For finite temperatures and deviation from half-filling we have to determine the chemical potential from the $`N(\mu )`$ curves. Due to the analytically given expression, this is in principle possible in any accuracy. For example, the chemical potential used in the Fig. 5, where we depicted the lowest part of the spectrum for $`N_e`$=3, results in a deviation of the electron number from 3 which is less than $`10^{12}`$ for $`T=10^7t`$. If the temperature reaches zero, we have a step curve. Surprisingly the model steps from Ne=2 to Ne=4 with no plateau at Ne=3. This corresponds to the interesting fact that there is an extreme small $`\mu `$-window, where a n=3 state is groundstate. This is demonstrated by help of the inset of Fig. 4 where we plotted the groundstate in dependence of the chemical potential. Thus it is interesting to look at the spectrum with Ne=3 also. Since we work with the grand-canonical potential we have to determine first the chemical potential. Fig. 5 makes clear that we find again states with different electron numbers among the lowest levels. The degeneration of the pseudospin-triplet is now lifted due to the large deviation of the chemical potential from U/2, and its lowest state (No. 22) is the first exited state above the groundstate, which is No. 46 and No. 50. Both states belong to the Ne=3 pseudospin doublets in the half-filled model. As mentioned above the degeneracy of the groundstate is due to the imaginary eigenvalue of the translation operator. Among the low lying states the No. 111, being the groundstate of the half-filled case, remained. If we do not apply a small magnetic field, the states No. 46 and No. 50 are degenerated with the states No. 70 and 74, and so do the states No. 22 and No. 111, which become the groundstates of the model for lower correlation, although their occupation eigenvalue is not 3. This situation is depicted in Fig. 6. Due to the interplay of the four s=1/2 states, belonging to Ne=3, and the two s=0 states belonging to Ne=2 and Ne=4 respectively one has to expect a complex magnetic behaviour of the model. For example a small magnetic field lifts the degeneracy of the Ne=3 states whereby the states No 46 and 50 are lowered below the states No. 22 and No. 111. This process is shown in Fig. 7. By further increasing the correlation energy U we found a spin change of the groundstate, since the S=3/2 states no. 38, 58, 82, and 90 become the groundstates. This was first discovered by Callaway et.al.  via numerical diagonalization. We add that the groundstate is fourfold degenerated, and the exact U value, where the transition happens is U/t=4(2+$`\sqrt{7}`$), what differs slightly from the numerical value 1/0.055 given in Ref. . The difference is in the second digit and may serve as an estimation of the accuracy of their calculation. For Ne=5 they did not find such a change in the groundstate spin, what is contradictory to the electron-hole symmetry, which was otherwise explicitly mentioned in their paper - and of course: the same transition happens at the same U value.

### B The contents of double occupied states

A further question of interest is whether projecting out the double-occupied states from beginning is a good approximation. Since we have all the eigenstates in analytical form, it was easy to calculate the contents of double occupied states in each state. The U-dependence of the admixing of double occupied states to the lowest lying states is depicted in Fig. 8. Among the considered states the groundstate has the highest content for bigger U, being about ten per cent for U=10t and falling to about two per cent for U=30t. If the correlation energy U is lower than the bandwidth the double occupied states admix considerably. Thus one has to be careful with the projection, especially if one aims at a description of the Hubbard model in the region were correlation energy and bandwidth are of the same order.

### C The spectral function

By help of the eigenstates and the spectrum it is possible to calculate the spectral function for arbitrary operators. Since the $`\delta `$-like peaks result in discontinuous curves being inconvenient for the graphics we substituted the $`\delta `$-functions by Lorentz distributions, i.e. we defined

$`J_{BA}(\omega )`$ $`=`$ $`{\displaystyle \frac{2\pi }{Z}}{\displaystyle \underset{m}{}}{\displaystyle \underset{n}{}}e^{\beta E_m}\mathrm{\Phi }_m|𝐁|\mathrm{\Phi }_n\mathrm{\Phi }_n|𝐀|\mathrm{\Phi }_m\times `$ (32)
$`\times {\displaystyle \frac{1}{\pi }}{\displaystyle \frac{\delta }{\delta ^2+(\omega E_m+E_n)^2}}`$
with Z being the grand canonical partition sum, and $`E_m`$ being the grand canonical levels. The most interesting function due to its relation to spectroscopy we get with $`𝐁=𝐜_{i\sigma }`$ and $`𝐀=𝐜_{i\sigma }^+`$. In Fig. 9 we show the results in dependence of U for the case of half-filling. The most ocular feature of the spectral function is the opening of the gap of width U with increasing correlation strength. This is in qualitative accordance with the common picture of the Hubbard-Mott transition.

## IV Thermodynamics

Once calculated the spectrum in an analytical way it is straight forward to get the grand canonical potential and the complete thermodynamics. For the four point cluster there exist several papers, which got the thermodynamics by numerical diagonalization. Heinig and Monecke calculated for this model the specific heat , the susceptibility and local magnetic moments  for the half filled case. Shiba and Pincus did the same within a comprehensive work devoted to the one dimensional half-filled-band Hubbard model. The discrepancies of these two works are discussed in a paper of Cabib and Kaplan, who calculated in contrast both the canonical and the grand canonical ensemble . They especially discussed the specific heat and the spin-spin correlation functions, where the main focus was again the half-filled case. In the following we give our results for both the half-filled case (for comparison) and the doped case.

### A Thermodynamic density of states

In their paper  Cabib and Kaplan argued that the grand canonical ensemble for a finite N system reflects the infinite N case much better than the canonical ensemble, especially for the large U case. We agree with that point of view, especially regarding the thermodynamical density of states (TDOS). We calculated this quantity, i.e. the derivation of the particle number with respect to the chemical potential.

$`D_T(\mu )`$ $`=`$ $`{\displaystyle \frac{N(\mu )}{\mu }}.`$ (33)
In Fig. 10 we depicted the results for different values of U at the relatively high temperature of $`T=0.1t`$ for graphical reasons. It is obvious, that the TDOS changes its character qualitatively with increasing U starting from a three-peak structure for small on site correlation changing to an eight-peak structure at large correlation. Around $`U=1t`$ there opens a gap around $`\mu =U/2`$. This means, adding a further electron to the half-filled model needs an energy of order U. We interpret this, in the same sense as above for the spectral function, as the variant of the “Hubbard-Mott-transition” for our small model. Of course, the “critical U”, where the gap opens, depends on the temperature. For $`T=0`$ the gap opens for infinitesimal small values of the correlation. This is a special case for the same statement, proved for the one-dimensional Hubbard model by help of the exact solution for the groundstate given by Lieb and Wu .

### B Specific heat

The temperature dependence was controversial in Refs. . For the half-filled case we find exactly the same features as given in Ref. . Cabib et. al. stated the qualitative picture given in Ref.  to be right and that of Ref.  to be wrong for the case U$`<`$6. Unfortunately a quantitative comparison with Ref.  was impossible due to inconsistencies in that paper. In Fig. 11 we show the specific heat curves in dependence on the temperature. The curves with U=0.7t and U=8t respectively are exactly the same as given in Fig. 1 of Ref.  and thus confirms this paper and serves as independent proof of our calculation. Furthermore, the change from the double peak shape to the triple peak structure is due to the behaviour of state 119, what is evident from Fig. 12. Variation of the electron density causes dramatic changes of the specific heat. In Fig. 13 we calculated for a constant electron density and constant magnetic field (for simplicity set to zero). The shape of the specific heat curves is qualitatively affected if the electron density is changed. In the right column of Fig. 13 the ratio of correlation energy U and band width is chosen to 1. This is the most interesting case, since it is not accessible by means of perturbation theory. We see that with deviation from half-filling the small peak at kB=0.1t vanishes, whereas a low energy peak develops at kBT=0.002t, being maximum at Ne=3, i.e. ne=0.75. The reason is the small distance between the states 46, 50, 70, 74 belonging to Ne=3, and the state 111, belonging to Ne=4, and state 22, belonging to Ne=2, as may be seen from Fig. 5. Since the Ne=3 groundstate is magnetic and the groundstates for Ne=2 and Ne=4 respectively are not, these excitations are both spin excitations and charge excitations. In the right column of Fig. 13 the correlation strength is twice the bandwidth. In that case we find no small energy peak for Ne deviating from 4. Instead the double peak structure smears out due to the increasing influence of magnetic excitations which are maximum if the electron number is fixed to three. Further decreasing the electron number results in vanishing of the double peak structure. The reason is that the number of accessible states with double occupation is reduced drastically.

### C The crossing points

A further subject of interest are the so called crossing points, which attracted much interest recently, since they are present in several strong correlated systems. Vollhardt clarified this point . In a second paper  he and the coworkers specified to Hubbard models and found a (nearly) universal crossing point of the C(T,U)/kB-curves. Since the arguments in Ref.  never use the thermodynamic limit, we feel that they should work for our small model also. In Fig. 14 we depicted the specific heat curves nearly in the same fashion as was done in Fig. 1 (b) of Ref. . The similarity is ocularly. The most striking point the authors made in Ref.  is the discovery of a nearly universal value of about 0.34 of the specific heat at the high temperature crossing point for small U. We find for our model a value of the crossing temperature T(U$``$0)/t=1.334820204t and C(U$``$0,T)=0.3351877115 NkB. All the Hubbard systems considered in Ref.  are in the thermodynamic limit at half-filling. Whether our small cluster is comparable to infinite systems is a difficult question. Usually the physical quantities calculated from the canonical and grand-canonical potentials differ in small systems. Therefore this difference may serve as a rough indicator how relevant the results are for the thermodynamic limit. In Fig. 15 we plotted the specific heat for the canonical ensemble. We found again a high temperature crossing point, whereas there is no sharp crossing point around kBT/t=0.8, in contrast to Fig. 14. The related specific heat value is C(U$``$0,T)=0.3744085787 NkB with T=0.940347332t. This is only a little bit bigger than the grand-canonical value. Nevertheless the big difference in T shows that there are differences at all. The crossing points at different U values were determined according to Ref.  by help of the equation

$`{\displaystyle \frac{C(T,U,\mu )}{U}}`$ $`=`$ $`0\text{with}\mu =U/2`$ (34)
for the grand canonical ensemble and by

$`{\displaystyle \frac{C(T,U)}{U}}`$ $`=`$ $`0`$ (35)
for the canonical ensemble with Ne=4 respectively. These functions significantly differ from each other, as one can see from Fig. 16.

### D Susceptibility

The isothermal susceptibility shows the interplay between the low-lying nonmagnetic states and the low-lying magnetic states. As shown above, the fact whether the groundstate is magnetic or not is dependent on the occupation. That is why we find an increasing peak due to the low-energy magnetic excitations if we decrease the electron number from four to three, which vanishes, if we subsequently decrease the number of electrons towards two. This is illustrated in Fig. 17. In the strong correlated case (U/t=8) we find a complete different behaviour. Now the distance between the states with different electron occupation is large. Consequently, we do not find any low-energy peaks as visible in Fig. 18. Instead we see a divergent behaviour due to the Sz degenerated groundstate if the electron number deviates from four. It is well believed, that the Hubbard model in three dimensions exhibits antiferromagnetism near the half-filled limit and for low temperatures. In our small model long range order does not exist, nevertheless we can ask how the correlation functions depend upon change of the correlation strength, of temperature and electron filling. We adopt the notation of Cabib and Kaplan  for the spin-spin correlation functions

$`L_n`$ $`=`$ $`4𝐒_i^z𝐒_{i+n}^z`$ (36)
If not otherwise mentioned we always take the expectation with respect to the grand canonical ensemble. By help of the functions $`L_n`$ one can calculate the susceptibility for our small cluster according to

$`\chi ^T`$ $`=`$ $`4\beta (L_0+2L_1+L_2)4\beta 𝐒^z.`$ (37)
The susceptibility calculated from eq. 37 has to agree with the expression we got from the differentiation of the thermodynamical potential. We used this fact to check the calculated correlation functions. In Fig. 19 the three different spin-spin-correlation functions are shown in dependence on the temperature for the half-filled case. The model shows the well known antiferromagnetic coupling of adjacent spins. Of course this breaks down if the temperature is increased.The breakdown shifts to lower temperatures when the strong coupling regime is reached. This is exactly what is expected, due to the fact that for large values of U the model can be mapped onto an antiferromagnetic Heisenberg model with the coupling constant being $``$ 1/U. Obviously the U=4t case and the U=8t case have nearly the same breakdown temperature, fortifying the above stated fact, that U=8t, i.e. U is twice the unperturbed bandwidth, belongs not to the strong coupling limit. It is exactly this parameter regime where we find an increase of the local magnetic moment with increasing temperature before it breaks down to its high temperature limit. The strong correlation behaviour of the model is shown in Fig. 20. The most striking fact is the change from antiferromagnetic correlation to ferromagnetic correlation if the electron density is lowered from one to three quarters. Also the next nearest neighbour correlation changes its sign when the electron occupation changes from three towards two. While for the neutral model the nearest neighbours are antiferromagnetically correlated, the next nearest neighbours are ferromagnetic. For three electrons both correlation functions are ferromagnetic and for two electrons they are both antiferromagnetic. The temperature dependence of the density-density correlations is depicted for the neutral model in Fig. 21. Regarding $`𝐧_1𝐧_1=𝐧_1+2𝐝`$ with $`𝐝`$ being the double occupancy operator $`𝐝=𝐧_{}𝐧_{}`$ it is obvious how the onsite electron correlation suppresses the probability to find two electrons at the same lattice site as long, as the correlation energy is small with respect to the temperature. Especially the correlation function $`𝐧_1𝐧_2`$ makes evident that the tendency to charge ordering is lowered. We remark, that the TDOS may be calculated via the density-density correlation functions, giving one more possibility to prove the inner consistency of the derived formulae. In Fig. 22 we depicted the nearest-neighbour density correlation for strong correlation in dependence of the electron filling of the model. With decreasing electron number the electrons get more space. Due to the repulsive electron-electron correlation the electrons tend to increase their distance. Therefore the nearest-neighbour density-density correlation is reduced considerably.

## V Discussion

The benefit one gets from an analytical solution of a model is the easy access to its complete physics. Due to the limited space, we had to restrict the presentation of the results to a few points. We focused here mainly to two points. The first one was the study of eigenvalues and related eigenstates in dependence on the correlation strength. This way we clarified the open questions about possibly hidden symmetries. It was shown, that every state is distinguished from each other by using the two SU symmetries of spin and pseudospin, the translation symmetry, and the eigenvalues of the Grosse operator $`𝐊_1`$. Since therefore all the eigenstates are off different symmetries, there remains no problem with the multitude of level crossings present within the model. The lowest states in the half-filled model yield the singlet-triplet picture, widely used in the context of high-Tc superconductivity. In the half-filled model the charge excitations are a little bit higher in energy, but, these excitations are lowered considerably in energy, if the electron number is driven towards three by appropriate tuning of the chemical potential, thus enlightening the complex interplay of spin and charge (more exactly: pseudospin) degrees of freedom. We were really surprised that the states with adjacent electron numbers are still lower in energy if three is chosen for the mean occupation, even if the correlation equals the bandwidth. This situation is completely changed if one applies a (relatively) small magnetic field or if the ratio U/t is increased slightly ( what may be done by help of pressure in real systems). Since our model is a special case of the n-site Hubbard rings which are the pet of the Bethe ansatzers, there must be a correspondence between the Bethe ansatz eigenstates and the eigenstates given in the Appendix. However, it is not obvious, how the quantum numbers characterizing the states in our solution correspond to the quantum numbers of the theories based on the Bethe ansatz. Two of them, the z-component of spin and pseudospin respectively are used within the Bethe ansatz, but regarding the other quantum numbers of the Bethe ansatz, i.e the charge momenta and the spin rapidities, we have no idea how they are related to the remaining quantum numbers of our model, i.e. the modules of spin and pseudospin, the eigenvalue of the translation operator, and the eigenvalue of the Grosse operator. The complicated numerical solution of the Lieb-Wu equations delivering the charge momenta and spin rapidities prevents a deeper insight, so it remains an open question. The interesting problem whether the eigenstates may be classified as k strings and k-$`\lambda `$ strings as Deguchi et. al.  did with respect to solutions of the Lieb-Wu equations following the string hypothesis of Takahashi , or in the way as was done in a very recent paper  remains a unsolved question for the time being. A second point we were interested in, is whether the model gives some hints towards the metal-insulator transition. Both the spectral density and the thermodynamical density of states show the expected behaviour, i.e. the opening of a gap of width U for T$`>`$0, so there is nothing spectacular in this direction. The general theorems for the Hubbard model are fulfilled of course. Our analytical solution allowed a detailed study of the crossing points of the specific heat. Our results confirm the explanation given by Vollhardt et. al. completely. At the first glance this seems astonishing, since we have a level spectrum instead of the various continuous DOS’ used in Ref. . On the other hand, if one deals with the high-temperature crossing point the temperature smoothes the level structure in some sense. Furthermore, the argument, given in Ref. , that the energy scale for T is essentially t if we are at higher temperatures and of the order of the singlet-triplet excitation at lower temperatures holds for our model also. The comparison of the canonical and grand canonical ensemble at half filling exhibit differences regarding the temperature dependence of the derivatives of specific heat with respect to the correlation strength.The reason is that a “sampling” of the DOS is much poorer with 72 states instead of 256. Thus, it is not unexpected that the values of C+ differ slightly from the value calculated in Ref.  for an infinite linear chain by help of second order perturbation theory in U/t. Considering the specific heat for deviations of the mean occupation number from half-filling we find another much lower energy scale, due to the fact that states with different electron numbers are nearly degenerated if we fix the mean electron number to three. Regarding the magnetic behaviour we see the interplay between the even-occupied clusters, which are nonmagnetic and exhibit low-lying singlet-triplet excitations and the odd occupied clusters which exhibit a magnetic degenerated groundstate with doublet-quartet excitations if we have Ne=3 or Ne=5. By fixing the mean occupation number to non-integers we enforce a mixture of both cluster groups yielding a diverging susceptibility due to the contents of magnetic clusters. The spin-spin correlations confirm the common accepted picture, that the strong on-site correlation favours antiferromagnetism instead ferromagnetism. Otherwise, next-nearest neighbours are coupled ferromagnetic. This holds for the half filled model. This situation is changed completely if the mean occupation number is driven towards three. There we find ferromagnetic coupling between adjacent spins for low temperature. A completely unexpected behaviour is the change to an antiferromagnetic coupling (albeit small) when the temperature is increased and becomes of the order of the hopping matrix element. For higher temperatures the off-site spin-spin correlation vanishes. Since we have a small cluster long range magnetic order is impossible, nevertheless, if one regards the nearest neighbour correlation as an hint to the magnetic phases realized in an infinite system, our findings support the common accepted phase diagram for the (simple cubic) three-dimensional Hubbard model. The suppression of electron motion with increasing correlation strength is evident from the density-density-correlations. If U is large enough we find exact one electron at every site. Motion seems forbidden, what hints towards an insulator. Regarding the square lattice we believe that it is more appropriate, with reference to the conserved quantities of the Hubbard model, to speak about spin-pseudospin separation instead of spin-charge separation, since the “charge” is the “z-component” of the pseudospin only. From models in magnetism it is common knowledge that quantum fluctuations reduce the magnetic moment, due to the fact that the total spin is conserved, while the local spin is usually not. In the Hubbard model the same holds for the pseudospin also. If the hopping between the clusters is treated by a perturbation theory as was done in Ref. , the calculated cluster levels will form bands, and the related cluster-spins and -pseudospins will not be conserved furthermore, thus “pseudospin-waves” will arise as well as spin-waves. The undoped model has a vanishing total pseudospin. Reducing the number of electrons from one per lattice cite creates a certain number of empty places. At these sites a local spin is missing now and a local pseudospin is created. Further doping dilutes the spins and increases the distance, what reduces the spin-spin interaction. Otherwise the number of pseudospins increases; their mean distance is reduced and the pseudospin-pseudospin interaction is increased, thus favouring long-range order of the pseudospins. As we have shown, among the lowest eigenstates of our cluster we find states with different spins and pseudospins. Since the related states are many-particle states the spins and the pseudospins are not coupled to individual electrons or holes. If these many-particle states become itinerant by switching on the inter-cluster-hopping, there is no reason, that the resulting bands must exhibit identical shapes, or in other words that the resulting quasiparticles have the same mass. On the contrary, one would expect, that a quasiparticle made from a cluster-state with three electrons and spin 1/2 will differ significantly from a quasiparticle made from a spinless cluster-state with four (or two) electrons, both in the effective mass and in the velocity of propagation. This is exactly, what was calculated in Ref.  and what is called spin-charge separation. Finally, we mention that the given solution is applicable to the four-site Hubbard model with attractive interaction as well, simply by applying the Shiba map  $`𝐜_i^+𝐚_i^+`$ and $`𝐜_i𝐚_i^+`$, which exchanges the spin with the pseudospin. The tables given in the appendix apply as well for negative U hamiltonian. The only task is to replace $`UU`$, $`h+U/2\mu `$, and $`\mu +U/2h`$ for the energy eigenvalues and to exchange the spin and pseudospin quantum numbers in the assignment of the eigenvectors. All the presented figures, may be calculated for the attractive model as well. There is a multitude of other interesting features, e.g. the transport properties within the cluster, which are nevertheless easy to calculate by help of the analytically given eigenstates and eigenvectors. We abstain from presenting the related figures to prevent overloading of the present paper. In summary we may say that the analytical solution of the four-site Hubbard model for arbitrary electron filling, arbitrary interaction strength and arbitrary temperature is highly non-trivial. It allows easy access to all interesting features of the model and may serve as a reference object for various numerical or perturbation methods dealing with more complex models of strong electron correlation.

###### Acknowledgements.

I would like to thank S.-L. Drechsler and R. Hayn for usefull discussion, J. Monecke for some helpful information on their early work, and P. Vollhardt for encouraging remarks on the crossing point calculation.

## Eigenvectors and Eigenvalues

Here we present some eigenstates and the complete spectrum. Since by far the most papers deal with the groundstate of the Hubbard model we restrict to presentation of the groundstate of the neutral model with Ne=4, i.e. state 111, and one of the groundstates for N$`{}_{e}{}^{}=3`$, i.e. state 46. The analytical form of the other 254 eigenvectors the interested reader may find at the web-page www.physik.tu-dresden.de/itp/members/schumann/research.html.

$`\mathrm{\Psi }_{111}`$ $`=`$ $`C_1\left(|0022|0220|2002+|2200\right)`$ (44)
$`+C_2\left(|02du+|02ud+|0du2|0ud2\right)`$
$`+C_2\left(|20du+|20ud+|2du0|2ud0\right)`$
$`+C_2\left(+|d02u+|d20u|du02|du20\right)`$
$`+C_2\left(|u02d|u20d+|ud02+|ud20\right)`$
$`+C_3\left(|dudu+|udud\right)`$
$`+C_4\left(|dduu+|duud+|uddu+|uudd\right)`$
$`C_1`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{3}}}+{\displaystyle \frac{4i\sqrt{3}UY^{\frac{1}{3}}}{N_1}}`$ (45)
$`C_1`$ $`=`$ $`\mathrm{96\hspace{0.17em}2}^{\frac{1}{3}}(\sqrt{3}i)t^2+\mathrm{6\hspace{0.17em}2}^{\frac{1}{3}}(\sqrt{3}i)U^2`$ (47)
$`12iUY^{\frac{1}{3}}2^{\frac{2}{3}}(\sqrt{3}+i)Y^{\frac{2}{3}}`$
$`C_2`$ $`=`$ $`{\displaystyle \frac{U}{8\sqrt{3}t}}+{\displaystyle \frac{(1+i\sqrt{3})\mathrm{\hspace{0.17em}2}^{\frac{1}{3}}t}{\sqrt{3}Y^{\frac{1}{3}}}}`$ (49)
$`+{\displaystyle \frac{(1+i\sqrt{3})U^2}{\mathrm{8\hspace{0.17em}2}^{\frac{2}{3}}\sqrt{3}tY^{\frac{1}{3}}}}+{\displaystyle \frac{(1i\sqrt{3})Y^{\frac{1}{3}}}{\mathrm{48\hspace{0.17em}2}^{\frac{1}{3}}\sqrt{3}t}}`$
$`C_3`$ $`=`$ $`1+{\displaystyle \frac{1}{\sqrt{3}}}`$ (50)
$`C_4`$ $`=`$ $`{\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{\sqrt{3}}}`$ (51)
$`Y`$ $`=`$ $`216t^2U+6\sqrt{3}X`$ (52)
$`X`$ $`=`$ $`\sqrt{4096t^6336t^4U^248t^2U^4U^6}`$ (53)
$`\mathrm{\Psi }_{46}`$ $`=`$ $`C_1\left(i|0duu+|duu0|u0du+i|uu0d\right)`$ (59)
$`+C_2\left(|0uud+i|d0uu+|ud0ui|uud0\right)`$
$`+C_3\left(i|0udu+i|du0u|u0ud+|udu0\right)`$
$`+C_4\left(i|002u|02u0i|2u00+|u002\right)`$
$`+C_5\left(i|020ui|0u02|20u0+|u020\right)`$
$`+C_6\left(|00u2+i|0u20i|200u|u200\right)`$
$`C_1`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{2}}}+{\displaystyle \frac{\frac{3i}{4}U}{\sqrt{2}t}}{\displaystyle \frac{X}{4\sqrt{2}t^2}}+{\displaystyle \frac{\frac{3i}{4}Y}{\sqrt{2}t}}`$ (62)
$`+{\displaystyle \frac{8i\sqrt{2}tY}{U^2}}{\displaystyle \frac{i\sqrt{2}XY}{tU^2}}+{\displaystyle \frac{6i\sqrt{2}t}{U}}`$
$`{\displaystyle \frac{\frac{3i}{2}X}{\sqrt{2}tU}}+{\displaystyle \frac{\sqrt{2}Y}{U}}{\displaystyle \frac{XY}{4\sqrt{2}t^2U}}`$
$`C_2`$ $`=`$ $`{\displaystyle \frac{i}{3\sqrt{2}}}{\displaystyle \frac{6\sqrt{2}t}{U}}{\displaystyle \frac{3U}{4\sqrt{2}t}}+{\displaystyle \frac{iX}{6\sqrt{2}t^2}}`$ (65)
$`+{\displaystyle \frac{3X}{2\sqrt{2}tU}}{\displaystyle \frac{Y}{2\sqrt{2}t}}{\displaystyle \frac{i\sqrt{2}Y}{U}}`$
$`+{\displaystyle \frac{iXY}{4\sqrt{2}t^2U}}+{\displaystyle \frac{8\sqrt{2}tY}{U^2}}+{\displaystyle \frac{\sqrt{2}XY}{tU^2}}`$
$`C_3`$ $`=`$ $`{\displaystyle \frac{1}{3\sqrt{2}}}+{\displaystyle \frac{X}{6\sqrt{2}t^2}}{\displaystyle \frac{2\sqrt{2}Y}{3U}}+{\displaystyle \frac{XY}{6\sqrt{2}t^2U}}`$ (66)
$`C_4`$ $`=`$ $`{\displaystyle \frac{i}{\sqrt{2}}}{\displaystyle \frac{4\sqrt{2}t}{3U}}+{\displaystyle \frac{X}{3\sqrt{2}tU}}{\displaystyle \frac{8\sqrt{2}tY}{3U^2}}+{\displaystyle \frac{\sqrt{2}XY}{3tU^2}}`$ (67)
$`C_5`$ $`=`$ $`{\displaystyle \frac{4\sqrt{2}t}{3U}}+{\displaystyle \frac{X}{3\sqrt{2}tU}}{\displaystyle \frac{8\sqrt{2}tY}{3U^2}}+{\displaystyle \frac{\sqrt{2}XY}{3tU^2}}`$ (68)
$`C_6`$ $`=`$ $`{\displaystyle \frac{i}{\sqrt{2}}}+{\displaystyle \frac{\sqrt{2}}{3}}\left({\displaystyle \frac{4t}{U}}{\displaystyle \frac{X}{2tU}}+{\displaystyle \frac{8tY}{U^2}}{\displaystyle \frac{XY}{tU^2}}\right)`$ (69)
$`X`$ $`=`$ $`\sqrt{64t^4+3t^2U^2}`$ (70)
$`Y`$ $`=`$ $`\sqrt{32t^2+U^2+4X}`$ (71)
In following Tables I-XIII we give the results for the whole spectrum. The eigenvalues of the operators $`𝐑_z`$ ,$`𝐑`$,$`𝐒_z`$,$`𝐒^2`$ ,$`𝐔`$,$`𝐊_{11}`$, $`𝐏^{++}`$,$`𝐏^+`$, $`𝐇`$ are assigned as $`m_r`$, $`m_s`$, $`r`$, $`s`$, $`u`$, $`k_{11}`$, $`p_{++}`$, $`p_+`$, $`\epsilon `$. In the first column we number the eigenstates. In the second column the energy eigenstates are given in the form $`|m_r,m_s,r(r+1),s(s+1),u,k_{11},p_{++},p+`$. The related energy eigenvalues in dependence on the parameters t and U are given in the third column. For the grand canonical eigenvalues these energy eigenvalues have to be completed by $`\mu N_ehm_s`$, therefore we give the related values in the subheads of the tables. In the fourth column we calculated the numbers for the energy for t=1, U=5, $`\mu =U/2`$, and $`h=0.01`$, i.e. the half filled case. The small magnetic field was introduced to lift spin-degeneracy.

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8

Figure 9

Figure 10

Figure 11

Figure 12

Figure 13

Figure 14

Figure 15

Figure 16

Figure 17

Figure 18

Figure 19

Figure 20

Figure 21

Figure 22