text:
# The 2-site Hubbard and 𝑡-𝐽 models

## I Introduction

Two are the aspects that gave so much popularity to the Hubbard model: the richness of its dynamics that is thought to permit a description of many puzzling issues like metal-insulator transition, itinerant magnetism, electronic superconductivity, and the simplicity of the Hamiltonian structure that let one speculate about the possibility of finding the exact and complete solution for any realization of the underlying lattice. Anyway, although the model has been studied more than any other one in the last fifty years, very few exact results are available and what we have mainly regards either finite clusters or the infinite chain (i.e., the 1D case). For finite clusters of 2Harris and Lange (1967); Shiba and Pincus (1972) or 4Shiba and Pincus (1972); Heinig and Monecke (1972a, b); Cabib and Kaplan (1973); Schumann (2002) sites it is possible to find the complete set of eigenstates and eigenvalues of the Hamiltonian and compute any quantity by means of the thermal averages. However, it is not easy at all, although possible in principle, to extract valuable and scalable (i.e., which can be used to find the solutions of bigger and bigger clusters and, ultimately, of the infinite lattice cases) information regarding the effective excitations present in the system, the operators describing them and their dynamics. For the infinite chain neither, we have all the information we wish; the Bethe Ansatz is a very powerful tool, but is severely limited as regards the range of applicability of the self-consistent equations it supplies and the quantities for which it gives an answer.

In this manuscript, in order to overcome the limitations discussed above, we have exactly solved the Hubbard model, on a 2-site cluster, completely within the equations of motion and the Green’s function formalism. By using this approach, we have had the possibility to find the complete set of eigenoperators of the Hamiltonian and the corresponding eigenenergies. This information has been really fundamental as it permitted a deeper comprehension of the features shown by the properties we have analyzed. It is worth noting that the Hubbard model on a 2-site cluster is the smallest system where both terms of the Hamiltonian (i.e., kinetic and electrostatic) are effective and contributes to the dynamics.

By properly tuning the value of the temperature, we have found that the 2-site system can almost perfectly mime, as regards relevant properties such as the kinetic energy, the double occupancy, the energy, the specific heat, the entropy and fundamental issues such as the metal insulator transition, the behavior of larger clusters and of the infinite chain. The tuning of the temperature is necessary in order to get a comparable effective level spacing (bigger the cluster, lower the spacing), i.e. to excite the correct levels: the relevant energy scales are present although the relative positions of the levels are affected by the size of the system (only two k points!). The very positive comparisons with exact results (Bethe Ansatz, exact diagonalization) and numerical data (quantum Monte Carlo, Lanczos) support the idea that a lot of physics can be described and understood within this very small system, for which there is the possibility to know the analytic expressions for all the quantities under study. As regards the relevant scales of energy, this 2-site system has demonstrated to contain all the necessary ingredients to describe many features coming from the strong electronic correlations and also appearing in the lattice case. Two of the three relevant energy scales, which are thought to be present in the Hubbard model, naturally emerge: the local Coulomb interaction $`U`$ and the spin-exchange one $`J=\frac{4t^2}{U}`$ that is, in principle, extraneous to the original purely electrostatic Hamiltonian and is dynamically generated by the combined actions of the two terms of the Hamiltonian. As useful guide to better understand the low-energy dynamics we have also solved the $`t`$-$`J`$ model and presented the solution in parallel with the one found for the Hubbard model.

This analysis, which have resulted to be really relevant by itself as we have got a much better understanding of some energy scales and internal parameter dynamics, has been worth to be performed also as a prelude of the lattice analysis. In fact, the eigenoperators we have found, both in the fermionic and the bosonic sectors, can be used in the lattice case as a basis for the Green’s functions. In the strongly correlated systems, the interactions can alter so radically the dynamics of the original particles that these latter lose completely their own identitiesUmezawa (1993). Actually, some new objects are generated by the interactions and dictate the physical response. They are not so easy to be identified: their number, exact expression and relevance can only be suggested by the experience and, when available, by exact and/or numerical results. For instance, one can choose: the higher order fields emerging from the equations of motion, the eigenoperators of some relevant interacting terms, the eigenoperators of the problem reduced to a small cluster, etc. In the last years, we have been focusing our activity on the study of strongly correlated electronic models like Kondo, $`t`$-$`J`$, $`p`$-$`d`$, Hubbard by means of the *Composite Operator Method*COM (a, b); Mancini (1998); Fiorentino:00a; Mancini and Avella (2003) that is based on two main ideas: one is the use of composite fields as basis for our Green’s functions, in accordance to what has been discussed above, and the other one is the exploitation of algebra constraints (e.g., the Pauli principle, the particle-hole symmetry, the Ward-Takahashi identities, …) to fix the correct representation of the Green’s functions and to recover the links among the spin and charge configurations dictated by the symmetries. It is worth noticing that the Composite Operator Method is exact in itself. An additional approximation treatment is needed when we deal with large or infinite degree-of-freedom systems; in this case we have to treat in an approximate way the otherwise intractable hierarchy of the equations of motion generated by the projection procedure. If no approximation is necessary (finite and reasonably small degree-of-freedom systems), the COM cannot do else than give the exact solution. According to this, the COM gives the exact solution also for the two systems under analysis in this manuscript: the two-site Hubbard and t-J models. Whenever, instead, we should resort to an approximate treatment to close the hierarchy of the equations of motion generated by the projection procedure, we expect some limitations connected with the chosen approximation. For instance, we get only the first moments correct if we truncate the equations of motion hierarchy Mancini (1998). On the other hand, it is really worth noting that we properly take into account: the interaction term of the Hamiltonian by using as basis operators its eigenoperators Avella et al. (2003a); the short-range correlations by using as basic fields the eigenoperators of the problem reduced to a small cluster Avella et al. (2003b); the presence of a Kondo-like singlet at low-energy by properly closing the equation of motion of an *ad hoc* chosen composite operator Villani et al. (2000).

Obviously, we are aware that the exact diagonalization of this very small system takes less than one afternoon to any graduate student. Then, the reader could wonder why we decided to study such a system so in detail. Well, the reasons are many and some have been already pointed out above:

* Any graduate student can surely compute eigenstates and eigenvalues of these Hamiltonians in one afternoon, but the analytic computation of Green’s functions and correlation functions, in terms of the former eigen*stuff*, is not that straightforward as one can think. At the end of the day, the time saved in computing eigenstates and eigenvalues instead of eigenoperators is almost fully recovered if you also take into account the time needed by the computation of the physical properties (Green’s functions and correlations functions). Then, the possibility to have *scalable* information putting, altogether, almost the same effort become quite tempting for anyone. At any rate and for the sake of completeness, we report in Appendix the expressions of Green’s functions and correlations functions in terms of eigenstates and eigenvalues.
* The knowledge of the exact eigenoperators of a system is invaluable as they could be used as correct starting point for the application of the projection methodsMor ; Rowe (1968); Roth (1969); Tserkovnikov (1981); Fulde (1995) to strongly correlated systems whose minimal modelMancini and Avella (2003) is the exactly solved one.
* The Green’s function formalism for composite operators is extremely complicatedMancini and Avella (2003). The 2-site Hubbard model can be used as toy model to fully explore this formalism and evidence the difficulties connected with the treatment of composite operators with non-canonical commutation relations. Within this system, the appearance of zero-frequency functions can be safely handled and resolved. The links among the different channels (fermion, charge, spin, pair) can be studied in detail.
* The 2-site Hubbard model, according to its status of minimal modelMancini and Avella (2003), contains the main scales of energy related with the interactions present in the Hamiltonian. The exact solution in terms of eigenoperators permits to individuate which are the composite fields responsible for the relevant transitions. These latter can be then used to efficiently study larger clusters.
* The absence of the three-site term, for obvious geometrical reasons, in the derivation of the 2-site $`t`$-$`J`$ model from the 2-site Hubbard one permits to push further the study of the relations between the two models. It is possible to exactly individuate the low-energy contributions and study them separately.
* Last, but not least, the capability of miming the numerical results for larger clusters (sometimes tuning the temperature and, consequently, the effective level spacing) opens the possibility to provide a testing-ground for the numerical techniques.

In the following, we define the models, give the self-consistent solutions in the fermionic and bosonic sectors and study the local, single-particle, thermodynamic and response properties of the systems; the eigenoperators, eigenstates and eigenvalues of the system are also given and analyzed in detail.

## II The 2-site Hubbard and $`t`$-$`J`$ models

The Hamiltonian of the Hubbard modelHub  for a $`N`$-site chain reads as

$$H=\underset{ij}{}\left(t_{ij}\mu \delta _{ij}\right)c^{}(i)c(j)+U\underset{i}{}n_{}(i)n_{}(i)$$
(1)
$`\mu `$ is the chemical potential, $`c^{}(i)=(c_{}^{}(i),c_{}^{}(i))`$ is the electronic creation operator at the site $`i`$ in the spinorial notation, $`U`$ is the on-site Coulomb interaction strength, $`n_\sigma (i)=c_\sigma ^{}(i)c_\sigma (i)`$ is the charge density operator for spin $`\sigma `$ at the site $`i`$ and

$$t_{ij}=2t\alpha _{ij}=t\underset{k}{}\mathrm{cos}\left[k\left(ij\right)\right]\alpha (k)$$
(2)
where $`t`$ is the hopping integral, $`a`$ is the lattice constant, $`\alpha _{ij}`$ is the projection operator on the nearest-neighbor sites and $`\alpha (k)=\mathrm{cos}(ka)`$. In the momentum representation, the kinetic term $`H_t`$ of the Hamiltonian (1) reads as

$$H_t=2t\underset{k}{}\alpha (k)c^{}(k)c(k)$$
(3)
where $`k`$ assumes the values $`\frac{2\pi }{Na}l`$ with $`l=0,\mathrm{},N1`$ for periodic boundary conditions and $`\frac{\pi }{(N+1)a}l`$ with $`l=1,\mathrm{},N`$ for open boundary conditions. We will study a 2-site cluster within periodic boundary conditions: $`N=2`$ and $`k=0,\frac{\pi }{a}`$.

The Hamiltonian of the $`t`$-$`J`$ modelChao et al. (1977); Anderson (1987); Zhang and Rice (1988) for the same cluster and boundary conditions reads as

$$\begin{array}{c}H=\underset{k}{}\left[2t\alpha (k)\mu \right]\xi ^{}(k)\xi (k)\hfill \\ \hfill +\frac{1}{2}J\underset{ij}{}\alpha _{ij}\nu ^\mu (i)\nu _\mu (j)\end{array}$$
(4)
where $`\xi (i)=\left[1n(i)\right]c(i)`$ is the fermionic composite operator describing the transitions $`n=0n=1`$, $`\nu _\mu (i)=\xi ^{}(i)\sigma _\mu \xi (i)`$ is the total charge ($`\mu =0`$) and spin ($`\mu =1,`$ $`2,`$ $`3`$) density operator at the site $`i`$, $`\sigma _\mu =(1,\stackrel{}{\sigma })`$, $`\sigma ^\mu =(1,\stackrel{}{\sigma })`$ and $`\stackrel{}{\sigma }`$ are the Pauli matrices.

We will extensively use the following definition for any field operator $`\mathrm{\Psi }(i)`$

$$\mathrm{\Psi }^\alpha (i)=\underset{j}{}\alpha _{ij}\mathrm{\Psi }(j)$$
(5)
In particular, for the 2-site system we have $`\mathrm{\Psi }^\alpha \left(0\right)=\mathrm{\Psi }\left(a\right)`$, $`\mathrm{\Psi }^\alpha \left(a\right)=\mathrm{\Psi }\left(0\right)`$ and $`\mathrm{\Psi }^{\alpha ^2}(i)=\left(\mathrm{\Psi }^\alpha (i)\right)^\alpha =\mathrm{\Psi }(i)`$.

## III The fermionic sector

### III.1 The equations of motion and the basis

#### III.1.1 The Hubbard model

After the Hubbard Hamiltonian \[Eq. (1)\], the electronic field $`c(i)`$ satisfies the following equation of motion

$$\mathrm{i}\frac{}{t}c(i)=\mu c(i)2tc^\alpha (i)+U\eta (i)$$
(6)
with $`\eta (i)=n(i)c(i)=\frac{1}{3}\sigma _kn_k(i)c(i)`$. According to this, we can decompose $`c(i)`$ as

$$c(i)=\xi (i)+\eta (i)$$
(7)
where $`\xi (i)`$ and $`\eta (i)`$ are the Hubbard operators and describe the transitions $`n=0n=1`$ and $`n=1n=2`$, respectively. Moreover, they are the local eigenoperators of the local term of the Hubbard Hamiltonian and describe the original electronic field dressed by the on-site charge and spin excitations. They satisfy the following equations of motion

$$\begin{array}{cc}& \mathrm{i}\frac{}{t}\xi (i)=\mu \xi (i)2tc^\alpha (i)2t\pi (i)\hfill \\ & \mathrm{i}\frac{}{t}\eta (i)=(U\mu )\eta (i)+2t\pi (i)\hfill \end{array}$$
(8)
with

$$\pi (i)=\frac{1}{2}\sigma ^\mu n_\mu (i)c^\alpha (i)+\xi (i)c^\alpha (i)\eta (i)$$
(9)
where $`n_\mu (i)=c^{}(i)\sigma _\mu c(i)`$ is the total charge ($`\mu =0`$) and spin ($`\mu =1,`$ $`2,`$ $`3`$) density operator at the site $`i`$ in the Hubbard model. We use a different symbol to distinguish it from the analogous operator, $`\nu _\mu `$, we have defined in the $`t`$-$`J`$ model. The field $`\pi (i)`$ contains the nearest neighbor charge, spin and pair excitations dressing the electronic field $`c(i)`$.

The field $`\pi (i)`$ can be also decomposed as

$$\pi (i)=\xi _s(i)+\eta _s(i)$$
(10)
where

$$\begin{array}{cc}& \xi _s(i)=\frac{1}{2}\sigma ^\mu n_\mu (i)\xi ^\alpha (i)+\xi (i)\eta _{}^{\alpha }{}_{}{}^{}(i)\eta (i)\hfill \\ & \eta _s(i)=\frac{1}{2}\sigma ^\mu n_\mu (i)\eta ^\alpha (i)+\xi (i)\xi _{}^{\alpha }{}_{}{}^{}(i)\eta (i)\hfill \end{array}$$
(11)
These latter, which are non-local eigenoperators of the local term of the complete Hubbard Hamiltonian, satisfy the following equations of motion

$$\begin{array}{cc}& \mathrm{i}\frac{}{t}\xi _s(i)=\mu \xi _s(i)+4t\eta (i)+2t\xi _s^\alpha (i)+4t\eta _s^\alpha (i)\hfill \\ & \mathrm{i}\frac{}{t}\eta _s(i)=(U\mu )\eta _s(i)+2t\eta (i)+2t\xi _s^\alpha (i)\hfill \end{array}$$
(12)
By choosing these four operators as components of the basic field

$$\psi (i)=\left(\begin{array}{c}\xi (i)\\ \eta (i)\\ \xi _s(i)\\ \eta _s(i)\end{array}\right)$$
(13)
we obtain a closed set of equations of motion. In the momentum space, we have

$$\mathrm{i}\frac{}{t}\psi (k)=\epsilon (k)\psi (k)$$
(14)
where the energy matrix $`\epsilon (k)`$ is

$$\epsilon (k)=\left(\begin{array}{cccc}\mu 2t\alpha (k)& 2t\alpha (k)& 2t& 2t\\ 0& U\mu & 2t& 2t\\ 0& 4t& \mu +2t\alpha (k)& 4t\alpha (k)\\ 0& 2t& 2t\alpha (k)& U\mu \end{array}\right)$$
(15)
It is worth noting that $`\psi (i)`$ is an eigenoperator of the Coulomb term of the Hubbard Hamiltonian for any lattice structure according to the local nature of the interactionAvella et al. (2003a).

#### III.1.2 The $`t`$-$`J`$ model

For the $`t`$-$`J`$ model we can obtain a closed set of equations of motion by choosing as basic field

$$\psi (i)=\left(\begin{array}{c}\xi (i)\\ \zeta (i)\end{array}\right)$$
(16)
where

$$\zeta (i)=\frac{1}{2}\sigma ^\mu \nu _\mu (i)\xi ^\alpha (i)$$
(17)
We have

$$\mathrm{i}\frac{}{t}\psi (k)=\epsilon (k)\psi (k)$$
(18)
with

$$\epsilon (k)=\left(\begin{array}{cc}\mu 2t\alpha (k)& 2t+2J\alpha (k)\\ 0& \mu +2t\alpha (k)4J\end{array}\right)$$
(19)
It should be noted that the $`t`$-$`J`$ model exactly reproduces the Hubbard model in the regime of very strong coupling (i.e., $`Ut`$). The three-site term, appearing in the derivation of the $`t`$-$`J`$ model from the Hubbard one, is completely absent in the 2-site system for obvious geometric reasons and the given Hamiltonians are exactly equivalent in the strong coupling limit. In particular, we note the following limits: $`lim_{Ut}\xi _s(i)=\zeta (i)`$ and $`lim_{Ut}\eta (i)=lim_{Ut}\eta _s(i)=0`$.

### III.2 The Green’s function

#### III.2.1 The Hubbard model

Let us now compute the thermal retarded Green’s function $`G(k,\omega )=\left[\psi (i)\psi ^{}(j)\right]`$ that, after the equation of motion of $`\psi (k)`$ (14), satisfies the following equation

$$\left[\omega \epsilon (k)\right]G(k,\omega )=I(k)$$
(20)
where $`I(k)=\{\psi (i),\psi ^{}(j)\}_{E.T.}`$ is the normalization matrix. The subscript $`E.T.`$ means that the anticommutator $`\{\mathrm{}\}`$ is evaluated at equal times. $``$ and $`[\mathrm{}]`$ stand for the Fourier transformation and the usual retarded operator, respectively. $`\mathrm{}`$ indicates the thermal average in the grand canonical ensemble. The solution of Eq. 20, by taking into account the retarded boundary conditions, is the following

$$G(k,\omega )=\underset{n=1}{\overset{4}{}}\frac{\sigma ^{(n)}(k)}{\omega E_n(k)+\mathrm{i}\delta }$$
(21)
The spectral weights $`\sigma ^{(n)}(k)`$ can be computed by means of the following expression

$$\sigma _{ab}^{(n)}(k)=\mathrm{\Lambda }_{an}\underset{c}{}\mathrm{\Lambda }_{nc}^1I_{cb}(k)$$
(22)
where $`\mathrm{\Lambda }`$ is a matrix whose columns are the eigenvectors of the energy matrix $`\epsilon (k)`$.

The energy spectra $`E_n(k)`$ are the eigenvalues of the energy matrix $`\epsilon (k)`$

$$\begin{array}{cc}& E_1(k)=\mu 2t\alpha (k)\hfill \\ & E_2(k)=\mu 2t\alpha (k)+U\hfill \\ & E_3(k)=\mu +2t\alpha (k)4J_U\hfill \\ & E_4(k)=\mu +2t\alpha (k)+4J_U+U\hfill \end{array}$$
(23)
where $`J_U=\frac{1}{8}\left(\sqrt{U^2+64t^2}U\right)`$. It is worth noting that $`lim_{Ut}J_U=\frac{4t^2}{U}`$ which is the value of $`J`$ we get from the derivation of the $`t`$-$`J`$ model from the Hubbard one in the strong coupling regime. One can easily check that the energy spectra $`E_n(k)`$ represent the energy of the single-particle transitions among the eigenstates given in App. A.

The normalization matrix $`I(k)`$ has the structure

$$I(k)=\left(\begin{array}{cccc}I_{11}& 0& I_{13}(k)& 0\\ 0& I_{22}& 0& I_{24}(k)\\ I_{13}(k)& 0& I_{33}(k)& 0\\ 0& I_{24}(k)& 0& I_{44}(k)\end{array}\right)$$
(24)
The explicit form of its entries are

$$\begin{array}{cc}& I_{11}=1\frac{n}{2}\hfill \\ & I_{22}=\frac{n}{2}\hfill \\ & I_{13}(k)=\mathrm{\Delta }+\alpha (k)(pI_{22})\hfill \\ & I_{24}(k)=\mathrm{\Delta }\alpha (k)p\hfill \\ & I_{33}(k)=2(pI_{22})2\alpha (k)\mathrm{\Delta }\hfill \\ & I_{44}(k)=I_{22}\hfill \end{array}$$
(25)
with

$`n=n(i)`$ (26a)
$`\mathrm{\Delta }=\xi ^\alpha (i)\xi ^{}(i)\eta ^\alpha (i)\eta ^{}(i)`$ (26b)
$`p={\displaystyle \frac{1}{4}}\chi ^\alpha d`$ (26c)
$`\chi ^\alpha ={\displaystyle \underset{\mu =0}{\overset{3}{}}}n_\mu ^\alpha (i)n_\mu (i)`$ (26d)
$`d=\xi _{}(i)\eta _{}(i)[\eta _{}^{}(i)\xi _{}^{}(i)]^\alpha `$ (26e)
$`\mathrm{\Delta }`$ gives a measure of the difference in mobility between the two Hubbard subbands. $`p`$ contains charge, spin and pair correlation functions. In order to compute the electronic Green’s function $`G_{cc}=\left[cc^{}\right]=G_{11}+2G_{12}+G_{22}`$ we only need the following quantities

$$\begin{array}{cc}& \sigma _{11}^{(1)}(k)=I_{11}+\frac{1}{2}\alpha (k)I_{13}(k)\hfill \\ & \sigma _{11}^{(2)}(k)=0\hfill \\ & \sigma _{11}^{(3)}(k)=\frac{1}{2}\alpha (k)I_{13}(k)\left(1\frac{4J_U}{U+8J_U}\right)\hfill \\ & \sigma _{11}^{(4)}(k)=\frac{1}{2}\alpha (k)I_{13}(k)\frac{4J_U}{U+8J_U}\hfill \end{array}$$
(27)
$$\begin{array}{cc}& \sigma _{12}^{(1)}(k)=\sigma _{12}^{(2)}(k)=0\hfill \\ & \sigma _{12}^{(3)}(k)=\sigma _{12}^{(4)}(k)=I_{13}(k)\frac{2t}{U+8J_U}\hfill \end{array}$$
(28)
$$\begin{array}{cc}& \sigma _{22}^{(1)}(k)=0\hfill \\ & \sigma _{22}^{(2)}(k)=I_{22}+\frac{1}{2}\alpha (k)I_{13}(k)\hfill \\ & \sigma _{22}^{(3)}(k)=\frac{1}{2}\alpha (k)I_{13}(k)\frac{4J_U}{U+8J_U}\hfill \\ & \sigma _{22}^{(4)}(k)=\frac{1}{2}\alpha (k)I_{13}(k)\left(1\frac{4J_U}{U+8J_U}\right)\hfill \end{array}$$
(29)
Equation (21) does not define uniquely the Green’s function, but only its functional dependenceMancini and Avella (2003). The knowledge of the spectral weights $`\sigma ^{(n)}(k)`$ and the spectra $`E_n(k)`$ requires the determination of the parameters $`\mathrm{\Delta }`$ and $`p`$, together with that of the chemical potential $`\mu `$. The connection between the parameter $`\mathrm{\Delta }`$ and the elements of the Green’s function and the equation fixing the filling allow us to determine these two parameters as functions of the parameter $`p`$

$$\{\begin{array}{c}n=2(1C_{11}C_{22})\hfill \\ \mathrm{\Delta }=C_{11}^\alpha C_{22}^\alpha \hfill \end{array}$$
(30)
where the correlation functions $`C_{ab}=\psi _a(i)\psi _b^{}(i)`$ and $`C_{ab}^\alpha =\psi _a^\alpha (i)\psi _b^{}(i)`$ can be computed, after the spectral theorem, by means of the following equations

$$\begin{array}{cc}& C_{ab}=\frac{1}{4}\underset{k}{}\underset{n=1}{\overset{4}{}}[1+T_n(k)]\sigma _{ab}^{(n)}(k)\hfill \\ & C_{ab}^\alpha =\frac{1}{4}\underset{k}{}\underset{n=1}{\overset{4}{}}\alpha (k)[1+T_n(k)]\sigma _{ab}^{(n)}(k)\hfill \end{array}$$
(31)
with

$$T_n(k)=\mathrm{tanh}\left(\frac{E_n(k)}{2T}\right)$$
(32)
The presence of parameters related to bosonic correlation functions within the fermionic dynamics (e.g., the parameter $`p`$) is characteristic of strongly correlated systems. In these systems, the elementary fermionic excitations are described by composite operators whose non-canonical anticommutation relations contain bosonic operators. Then, according to the general projection procedure, which is approximate for an *infinite* (i.e., with infinite degrees of freedom) system and coincident with the exact diagonalization for a *finite* system, the energy and normalization matrices generally contain correlation functions of these bosonic operators.

The parameter $`p`$ could be computed through its definition (see Eq. 26c). In this case we need to open the bosonic sector (i.e., the charge, spin and pair channels), to determine the eigenoperators and the normalization and energy matrices and to require the complete self-consistency between the fermionic and the bosonic sectors. Moreover, we run into the problem of fixing the value of the zero-frequency constantsZFC . This procedure will be widely discussed in the next section where the bosonic sector and the relative channels will be opened and completely solved. However, we can use another procedure to fix the parameter $`p`$ without resorting to the bosonic sector. We have to think over the reason why the parameters $`\mu `$, $`\mathrm{\Delta }`$ and $`p`$ appear into our equations: we have not fixed yet the representation of the Hilbert space where the Green’s functions are realized. The proper representation is the one where all the relations among operators coming from the algebra (e.g., $`n_\sigma ^2=n_\sigma `$) and the symmetries (e.g., particle-hole, spin rotation, …) are verified as relations among matrix elements. In the case of usual electronic operators the representation is fixed by simply determining the value of the chemical potential. In the case of composite fields, which do not satisfy canonical anticommutation relations, fixing the representation is more involved and the presence of internal parameters (e.g., $`\mu `$, $`\mathrm{\Delta }`$ and $`p`$) is essential to the process of determining the proper representation. The requirement that the algebra is satisfied also at macroscopic level generates a set of self-consistent equations, *the local algebra constraints*, that will fix the value of the internal parametersMancini and Avella (2003).

In particular, for the 2-site system the equation

$$\xi (i)\eta ^{}(i)=0$$
(33)
coming from the algebra constraint $`\xi (i)\eta ^{}(i)=0`$, together with Eqs. 30 constitute a complete set of self-consistent equations which exactly solves the fermionic dynamics allowing to compute the internal parameters (i.e., $`\mu `$, $`\mathrm{\Delta }`$ and $`p`$) for any value of the model ($`t`$, $`U`$) and thermodynamical ($`n`$, $`T`$) parameters. This procedure is clearly extremely simpler than that requiring the opening of the bosonic sector. It is worth mentioning that the system of self-consistent equations can be analytically solved as regards the parameters $`\mathrm{\Delta }`$ and $`p`$ as functions of the chemical potential $`\mu `$. It is possible to show that the equation for this latter parameter, as function of the internal (i.e., model and thermodynamical) parameters, exactly agrees with that coming from the thermal averages (cfr. Eq. A.3.1) as it should be since the model is exactly solved. This further confirms the validity, the effectiveness and the power of the method used to fix the representation.

In the next section, we will see that also in the case of the bosonic sector, the determination of the proper representation will be obtained by means of the local algebra constraints which will fix the value of the zero-frequency constants. It will be also shown that the two procedures for computing the parameter $`p`$ are equivalent, as it should be, and give exactly the same results for the fermionic dynamics although with remarkably different effort. It is worth noticing that the use of the local algebra constraints, which is unavoidable to fix the representation both in the fermionic and in the bosonic sectorsMancini and Avella (2003), permits to close the fermionic sector on itself without resorting to the bosonic one also in the lattice caseCOM (b) where a fully self-consistent solution of both sectors, although approximate, is very difficult to obtain.

#### III.2.2 The $`t`$-$`J`$ model

In the $`t`$-$`J`$ model, the energy spectra (i.e., the eigenvalues of $`\epsilon (k)`$) read as

$$\begin{array}{cc}& E_1(k)=\mu 2t\alpha (k)\hfill \\ & E_2(k)=\mu +2t\alpha (k)4J\hfill \end{array}$$
(34)
In particular, $`E_2(k)`$ corresponds to the single-particle transitions between the single occupied states and the 2-site singlet. The latter is the ground state at half filling. $`\zeta (k)`$ drives the corresponding fermionic excitation between the single occupied states and the 2-site singlet and the appearance of the most relevant scale of energy at low temperatures $`J`$.

The normalization matrix $`I(k)`$ has the following entries

$$\begin{array}{cc}& I_{11}=1\frac{n}{2}\hfill \\ & I_{12}(k)=C_{11}^\alpha +\alpha (k)\left(\frac{1}{4}\chi ^\alpha \frac{n}{2}\right)\hfill \\ & I_{22}(k)=2I_{12}(k)\hfill \end{array}$$
(35)
In the $`t`$-$`J`$ model the spin and charge correlator $`\chi ^\alpha `$ is defined as $`\chi ^\alpha =_{\mu =0}^3\nu _\mu ^\alpha (i)\nu _\mu (i)`$ in agreement with the strong coupling nature of the model. The spectral weights $`\sigma ^{(n)}(k)`$ have the following expressions

$$\begin{array}{cc}& \sigma _{11}^{(1)}(k)=I_{11}+\frac{1}{2}I_{12}(k)\hfill \\ & \sigma _{11}^{(2)}(k)=\frac{1}{2}I_{12}(k)\hfill \end{array}$$
(36)
$$\begin{array}{cc}& \sigma _{12}^{(1)}(k)=0\hfill \\ & \sigma _{12}^{(2)}(k)=I_{12}(k)\hfill \end{array}$$
(37)
$$\begin{array}{cc}& \sigma _{22}^{(1)}(k)=0\hfill \\ & \sigma _{22}^{(2)}(k)=I_{22}(k)\hfill \end{array}$$
(38)
Also in this case, as for the Hubbard model, we could compute $`\chi ^\alpha `$ opening the bosonic sector for the $`t`$-$`J`$ model. Again, following the reasoning given in the previous section, a simpler and completely equivalent procedure relies on the exploitation of the algebra constraint $`\zeta (i)\zeta ^{}(i)=2\zeta (i)\xi ^{}(i)`$. According to this, the parameters $`\mu `$ and $`\chi ^\alpha `$ can be computed by means of the following set of self-consistent equations

$$\{\begin{array}{c}n=1C_{11}\hfill \\ C_{22}=2C_{12}\hfill \end{array}$$
(39)
## IV The bosonic sector

### IV.1 The equations of motion and the basis

#### IV.1.1 The Hubbard model

##### The spin and charge sectors.

After the Hubbard Hamiltonian \[Eq. (1)\], the charge ($`\mu =0`$) and spin ($`\mu =1`$, $`2`$, $`3`$) density operator $`n_\mu (i)=c^{}(i)\sigma _\mu c(i)`$ satisfies a closed set of equations of motion, which describes the spin and charge dynamics in the system under study, once we choose as basic field

$$\varphi _\mu (i)=\left(\begin{array}{c}\varphi _\mu ^{(1)}(i)\\ \varphi _\mu ^{(2)}(i)\\ \varphi _\mu ^{(3)}(i)\\ \varphi _\mu ^{(4)}(i)\\ \varphi _\mu ^{(5)}(i)\\ \varphi _\mu ^{(6)}(i)\end{array}\right)=\left(\begin{array}{c}n_\mu (i)\\ g_\mu (i)\\ w_\mu (i)+w_\mu ^{}(i)\\ w_\mu (i)w_\mu ^{}(i)\\ h_\mu (i)h_\mu ^{}(i)\\ h_\mu (i)+h_\mu ^{}(i)\end{array}\right)$$
(40)
where

$`g_\mu (i)`$ $`=c^{}(i)\sigma _\mu c^\alpha (i)c^\alpha (i)\sigma _\mu c(i)`$ (41)
$`w_\mu (i)`$ $`=d_\mu (i)d_\mu ^\alpha (i)`$ (42)
$`h_\mu (i)`$ $`=f_\mu (i)f_\mu ^\alpha (i)`$ (43)
$`d_\mu (i)`$ $`=\xi ^{}(i)\sigma _\mu \eta ^\alpha (i)`$ (44)
$`f_0(i)`$ $`=\eta ^{}(i)\eta (i)d_0^{}(i)d_0^\alpha (i)`$
$`+\eta ^{}(i)\eta (i)\xi ^\alpha (i)\xi ^\alpha (i)`$ (45)
$`f_k(i)`$ $`=\xi ^{}(i)\xi (i)n_k^\alpha (i){\displaystyle \frac{1}{2}}\mathrm{i}\epsilon ^{kpq}n_q(i)n_q^\alpha (i)`$ (46)
The components of $`\varphi _\mu (i)`$, suggested by the hierarchy of the equations of motion, are either hermitian or anti-hermitian as densities or currents should be at any order in time differentiation. $`\varphi _\mu (i)`$ satisfies the following equation of motion in momentum space

$$\mathrm{i}\frac{}{t}\varphi _\mu (k)=\omega _b(k)\varphi _\mu (k)$$
(47)
where

$$\omega _b(k)=\left(\begin{array}{cccccc}0& 2t& 0& 0& 0& 0\\ 4t[1\alpha (k)]& 0& U& 0& 0& 0\\ 0& 0& 0& U& 2t& 0\\ 0& 0& U& 0& 0& 2t\\ 0& 0& 8t& 0& 0& 0\\ 0& 0& 0& 8t& 0& 0\end{array}\right)$$
(48)
##### The pair channel.

Within the analysis of the dynamics of the Hubbard model, another relevant bosonic operator is the pair operator $`p(i)=c_{}(i)c_{}(i)`$. The set of composite fields

$$P(i)=\left(\begin{array}{c}P^{(1)}(i)\\ P^{(2)}(i)\\ P^{(3)}(i)\\ P^{(4)}(i)\\ P^{(5)}(i)\\ P^{(6)}(i)\end{array}\right)$$
(49)
where

$`P^{(1)}(i)`$ $`=p(i)`$ (50)
$`P^{(2)}(i)`$ $`=c_{}(i)c_{}^\alpha (i)`$ (51)
$`P^{(3)}(i)`$ $`=c_{}(i)\eta _{}^\alpha (i)+\eta _{}(i)c_{}^\alpha (i)`$ (52)
$`P^{(4)}(i)`$ $`=p(i)n^\alpha (i)`$ (53)
$`P^{(5)}(i)`$ $`=2\eta _{}(i)\eta _{}^\alpha (i)`$ (54)
$`P^{(6)}(i)`$ $`=p(i)\eta ^\alpha (i)\eta ^\alpha (i)`$ (55)
satisfies the following closed set of equations of motion

$$\mathrm{i}\frac{}{t}P(k)=\omega _p(k)P(k)$$
(56)
where

$$\omega _p(k)=\left(\begin{array}{cccccc}U2\mu & 2t\left[1+\alpha (k)\right]& 0& 0& 0& 0\\ 2t\left[1+\alpha (k)\right]& 2\mu & U& 0& 0& 0\\ 0& 0& U2\mu & 2t\left[1+\alpha (k)\right]& U& 0\\ 0& 0& 2t\left[1+\alpha (k)\right]& U2\mu & 0& 0\\ 0& 0& 0& 0& 2(U\mu )& 2t\left[1+\alpha (k)\right]\\ 0& 0& 0& 0& 2t\left[1+\alpha (k)\right]& U2\mu \end{array}\right)$$
(57)
#### IV.1.2 The $`t`$-$`J`$ model

The basis for the bosonic sector in the $`t`$-$`J`$ model is given by

$$\phi (i)=\left(\begin{array}{c}\phi _0^{(1)}(i)\\ \phi _0^{(2)}(i)\end{array}\right)$$
(58)
for the charge channel, and by

$$\phi _k(i)=\left(\begin{array}{c}\phi _k^{(1)}(i)\\ \phi _k^{(2)}(i)\\ \phi _k^{(3)}(i)\\ \phi _k^{(4)}(i)\end{array}\right)$$
(59)
for the spin channel. We have defined

$`\phi _\mu ^{(1)}(i)`$ $`=\nu _\mu (i)`$ (60)
$`\phi _\mu ^{(2)}(i)`$ $`=\xi ^{}(i)\sigma _\mu \xi ^\alpha (i)\xi ^\alpha (i)\sigma _\mu \xi (i)`$ (61)
$`\phi _k^{(3)}(i)`$ $`=\nu (i)\nu _k^\alpha (i)`$ (62)
$`\phi _k^{(4)}(i)`$ $`=\mathrm{i}\epsilon ^{kpq}\nu _p(i)\nu _q^\alpha (i)`$ (63)
In the momentum space we have the following closed sets of equations of motion

$`\mathrm{i}{\displaystyle \frac{}{t}}\phi (q)`$ $`=\omega _c(q)\phi (q)`$ (64)
$`\mathrm{i}{\displaystyle \frac{}{t}}\phi _k(q)`$ $`=\omega _s(q)\phi _k(q)`$ (65)
with

$$\omega _c(q)=\left(\begin{array}{cc}0& 2t\\ 4t[1\alpha (q)]& 0\end{array}\right)$$
(66)
and

$$\omega _s(q)=\left(\begin{array}{cccc}0& 2t& 0& 2J\\ 4t[1\alpha (q)]& 0& 4t[1\alpha (q)]& 0\\ 0& 0& 0& 2J\\ 0& 0& 4J[1\alpha (q)]& 0\end{array}\right)$$
(67)
### IV.2 The Green’s function

In Ref. Mancini and Avella, 2003, we have shown that the retarded and causal Green’s functions contain substantially different information according to the unavoidable presence of the zero frequency constants (ZFC). In particular, we have reported on the relations between different types of Green’s functions and on the correct order in which they should be computed. According to this, in the bosonic sector, we have to start from the causal Green’s function and not from the retarded one, as we correctly did in the fermionic sector.

#### IV.2.1 The Hubbard model

##### The spin and charge sectors.

After Eq. 47, the thermal causal Green’s function $`G_\mu (k,\omega )=𝒯\left[\varphi _\mu (i)\varphi _\mu ^{}(j)\right]`$ satisfies the following equation

$$\left[\omega \omega _b\left(k\right)\right]G_\mu (k,\omega )=I_\mu (k)$$
(68)
where the relevant entries of the normalization matrix $`I_\mu (k)=[\varphi _\mu (i),\varphi _\mu ^{}(j)]_{E.T.}`$ are

$`I_{11\mu }(k)`$ $`=0`$ (69)
$`I_{12\mu }(k)`$ $`=4[1\alpha (k)]C_{cc}^\alpha `$ (70)
$`I_{13\mu }(k)`$ $`=0`$ (71)
$`I_{14\mu }(k)`$ $`=8[1\alpha (k)]C_{12}^\alpha `$ (72)
$`I_{15\mu }(k)`$ $`=[1\alpha (k)]q^{(\mu )}`$ (73)
$`I_{16\mu }(k)`$ $`=0`$ (74)
with

$`C_{cc}^\alpha =c^\alpha (i)c^{}(i)`$ (75)
$`q^{(0)}=16d`$ (76)
$`q^{(k)}={\displaystyle \frac{8}{3}}\chi _s^\alpha `$ (77)
$`\chi _s^\alpha ={\displaystyle \underset{k=1}{\overset{3}{}}}n_k^\alpha (i)n_k(i)`$ (78)
$`𝒯[\mathrm{}]`$ stands for the usual time ordering operator.

The solution of Eq. 68 is the following oneMancini and Avella (2003)

$$\begin{array}{c}G_\mu (k,\omega )=2\mathrm{i}\pi \mathrm{\Gamma }_\mu (k)\delta (\omega )\hfill \\ \hfill +\underset{i=1}{\overset{6}{^{}}}\frac{\sigma _\mu ^{(i)}(k)}{1\mathrm{e}^{\beta \omega }}\left[\frac{1}{\omega E_{ib}(k)+\mathrm{i}\delta }\frac{\mathrm{e}^{\beta \omega }}{\omega E_{ib}(k)\mathrm{i}\delta }\right]\end{array}$$
(79)
where

$$\begin{array}{c}\mathrm{\Gamma }_\mu (k)=\frac{1}{2}\underset{\omega 0}{lim}\omega \hfill \\ \hfill \times \left[\theta (t_it_j)\varphi _\mu (i)\varphi _\mu ^{}(j)\theta (t_jt_i)\varphi _\mu ^{}(j)\varphi _\mu (i)\right]\end{array}$$
(80)
is the zero frequency function and is undetermined at this level unless to compute the quite *anomalous* Green’s function appearing in its definition, which involves anticommutators of bosonic operators Mancini and Avella (2003). The general definition of a zero frequency function in terms of eigenvectors and eigenvalues is given in Eq. 205. The primed sum is restricted to values of $`i`$ for which $`E_{ib}(k)0`$.

The $`E_{ib}(k)`$ are the eigenvalues of the energy matrix $`\omega _b\left(k\right)`$

$`E_{1b}(k)`$ $`=2t\sqrt{2[1\alpha (k)]}`$ (81)
$`E_{2b}(k)`$ $`=2t\sqrt{2[1\alpha (k)]}`$ (82)
$`E_{3b}(k)`$ $`=U4J_U`$ (83)
$`E_{4b}(k)`$ $`=4J_U`$ (84)
$`E_{5b}(k)`$ $`=4J_U`$ (85)
$`E_{6b}(k)`$ $`=U+4J_U`$ (86)
According to this, the primed sum in Eq. 79 does not contain the $`i=1`$ and $`i=2`$ elements for $`k=0`$ ($`E_{1,2b}(0)=0`$) and the zero frequency function $`\mathrm{\Gamma }_\mu (k)`$ reduces to the constant $`\mathrm{\Gamma }_\mu =\mathrm{\Gamma }_\mu (0)`$. The spectral weights $`\sigma _{ab\mu }^{(i)}(k)`$ with $`a`$ and $`b=1`$, $`2`$, computed through Eq. 22 have the following expressions

$`\sigma _{11\mu }^{(1)}(k)`$ $`=\sigma _{11\mu }^{(2)}(k)={\displaystyle \frac{\sqrt{1\alpha (k)}\left[U\left\{tq^{(\mu )}\left[1+\alpha (k)\right]+4U\left[1\alpha (k)\right]C_{12}^\alpha \right\}+2C^\alpha \left\{8\left[1+\alpha (k)\right]^2t^2U^2\left[1\alpha (k)\right]\right\}\right]}{\sqrt{2}\left\{8\left[1+\alpha (k)\right]^2t^2U^2\left[1\alpha (k)\right]\right\}}}`$ (87)
$`\sigma _{11\mu }^{(3)}(k)`$ $`=\sigma _{11\mu }^{(6)}(k)={\displaystyle \frac{4tU\left[1\alpha (k)\right]\left[4(U+4J_U)C_{12}^\alpha tq^{(\mu )}\right]}{\left(U+8J_U\right)\left\{16t^2\left[1+\alpha (k)\right]+2U(U+4J_U)\right\}}}`$ (88)
$`\sigma _{11\mu }^{(4)}(k)`$ $`=\sigma _{11\mu }^{(5)}(k)={\displaystyle \frac{4tU\left[1\alpha (k)\right]\left[16J_UC_{12}^\alpha +tq^{(\mu )}\right]}{\left(U+8J_U\right)\left\{16t^2\left[1+\alpha (k)\right]8UJ_U\right\}}}`$ (89)
and

$`\sigma _{12\mu }^{(i)}(k)`$ $`={\displaystyle \frac{E_{ib}(k)}{2t}}\sigma _{11\mu }^{(i)}(k)`$ (90)
$`\sigma _{22\mu }^{(i)}(k)`$ $`={\displaystyle \frac{E_{ib}(k)}{2t}}\sigma _{12\mu }^{(i)}(k)=\left({\displaystyle \frac{E_{ib}(k)}{2t}}\right)^2\sigma _{11\mu }^{(i)}(k)`$ (91)
We note the sum rules

$$\underset{i=1}{\overset{6}{}}\sigma _{ab\mu }^{(i)}(k)=I_{ab\mu }(k)$$
(92)
In order to finally compute the Green’s function $`G_{11\mu }(k,\omega )`$ we should fix the internal parameters $`\chi _s^\alpha `$, $`d`$ and the zero frequency constant $`\mathrm{\Gamma }_{11\mu }`$.

The parameter $`\chi _s^\alpha `$ is directly connected to the Green’s function. Let us consider the correlation function $`D_\mu (k,\omega )=\varphi _\mu (i)\varphi _\mu ^{}(j)`$ which is linked to the causal Green’s function through the spectral theorem

$$\begin{array}{c}D_\mu (k,\omega )=\left(1+\mathrm{tanh}\frac{\omega }{2T}\right)\mathrm{}\left[G_\mu (k,\omega )\right]\hfill \\ \hfill =2\pi \mathrm{\Gamma }_\mu \delta (\omega )+2\pi \underset{i=1}{\overset{6}{^{}}}\delta \left[\omega E_{ib}(k)\right]\frac{\sigma _\mu ^{(i)}(k)}{1\mathrm{e}^{\beta E_{ib}(k)}}\end{array}$$
(93)
Then,

$$\chi _s^\alpha =3\varphi _3^{(1)}(i)\varphi _3^{(1)\alpha }(i)$$
(94)
The computation of the parameter $`d`$ requires instead the opening of the pair sector. Otherwise, it could be computed through the following relation with the parameter $`p`$

$$d=\frac{1}{4}\chi ^\alpha p$$
(95)
already given in the fermionic sector (see Eq. 26c), once we note that $`\chi ^\alpha =\chi _s^\alpha +\varphi _0^{(1)}(i)\varphi _0^{(1)\alpha }(i)`$.

The zero frequency constant $`\mathrm{\Gamma }_{11\mu }`$ cannot be directly connected to any correlation function at this level and its determination requires the use of local algebra constraints. In particular, we can use the following relation

$$\varphi _\mu ^{(1)}(i)\varphi _\mu ^{(1)}(i)=n_\mu (i)n_\mu (i)=n+2D(2\delta _{\mu 0}1)$$
(96)
where the double occupancy $`D=n_{}(i)n_{}(i)`$ is given by $`D=I_{22}C_{22}`$.

Equations 94, 95 and 96 constitute a complete set of self-consistent equations which allow to compute the parameters $`\chi _s^\alpha `$, $`d`$ and $`\mathrm{\Gamma }_{11\mu }`$ and then to determine the Green’s function $`G_{11\mu }(k,\omega )`$.

It is worth noting that, once $`C_{cc}^\alpha `$, $`C_{12}^\alpha `$, $`p`$ and $`D`$ are computed in the fermionic sector, the bosonic correlation functions can be easily obtained. We have the following expressions for the relevant correlators and zero-frequency constants

$`\chi _s^\alpha =3{\displaystyle \frac{n2D2a}{1+2b}}`$ (97)
$`\chi _c^\alpha =\varphi _0^{(1)}(i)\varphi _0^{(1)\alpha }(i)=n+2D2a`$
$`{\displaystyle \frac{b}{1+b}}\left[{\displaystyle \frac{3(n2D2a)}{1+2b}}+n+2D2a4p\right]`$ (98)
$`d={\displaystyle \frac{1}{4(1+b)}}\left[{\displaystyle \frac{3(n2D2a)}{1+2b}}+n+2D2a4p\right]`$ (99)
$`\mathrm{\Gamma }_{110}={\displaystyle \frac{1}{1+b}}[2n+4D2a`$
$`+b(n+2D3{\displaystyle \frac{n2D2a}{1+2b}}+4p]`$ (100)
$`\mathrm{\Gamma }_{113}={\displaystyle \frac{2(n2D)(1+b)2a}{1+2b}}`$ (101)
where

$`a`$ $`=(C_{cc}^\alpha 2C_{12}^\alpha )\mathrm{coth}{\displaystyle \frac{2t}{T}}`$
$`+{\displaystyle \frac{8tC_{12}^\alpha }{U+8J_U}}\left[\mathrm{coth}{\displaystyle \frac{E_{3b}(\pi )}{2T}}+\mathrm{coth}{\displaystyle \frac{E_{4b}(\pi )}{2T}}\right]`$ (102)
$`b`$ $`={\displaystyle \frac{4J_U}{U+8J_U}}\mathrm{coth}{\displaystyle \frac{E_{3b}(\pi )}{2T}}+{\displaystyle \frac{U+4J_U}{U+8J_U}}\mathrm{coth}{\displaystyle \frac{E_{4b}(\pi )}{2T}}`$ (103)
Taking into account the following local algebra constraints

$`D_{120}(i,i)`$ $`=n(i)g^{}(i)=2C_{cc}^\alpha `$ (104)
$`D_{220}(i,i)`$ $`=g(i)g^{}(i)=2n\chi ^\alpha 4d`$
$`=2(n4d2p)=2(n\chi ^\alpha +2p)`$ (105)
$`D_{123}(i,i)`$ $`=n_3(i)g_3^{}(i)=2C_{cc}^\alpha `$ (106)
$`D_{223}(i,i)`$ $`=g_3(i)g_3^{}(i)=2n+{\displaystyle \frac{1}{3}}\chi _s^\alpha \chi _c^\alpha +4d`$
$`=2(n{\displaystyle \frac{2}{3}}\chi _s^\alpha 2p)`$ (107)
and performing similar calculations we can also obtain (we omit the expressions of $`\mathrm{\Gamma }_{22\mu }`$ for the sake of brevity)

$`n^\alpha (i)g^{}(i)=2C_{cc}^\alpha `$ (108)
$`g^\alpha (i)g^{}(i)=2(n4d2p)`$ (109)
$`\mathrm{\Gamma }_{12\mu }=0`$ (110)
##### The pair channel.

After Eq. 56, the thermal causal Green’s function $`G_p(k,\omega )=𝒯\left[P(i)P^{}(j)\right]`$ satisfies the following equation

$$\left[\omega \omega _p\left(k\right)\right]G_p(k,\omega )=I_p(k)$$
(111)
where the relevant entries of the normalization matrix $`I_p(k)=[P(i),P^{}(j)]_{E.T.}`$ are

$`I_{11p}(k)`$ $`=1n`$ (112)
$`I_{12p}(k)`$ $`=[1+\alpha (k)]C_{cc}^\alpha `$ (113)
$`I_{13p}(k)`$ $`=I_{12p}(k)`$ (114)
$`I_{14p}(k)`$ $`=n\chi _c^\alpha +2d\alpha (k)`$ (115)
$`I_{15p}(k)`$ $`=2[1+\alpha (k)]C_{12}^\alpha `$ (116)
$`I_{16p}(k)`$ $`=2D2\gamma +2d\alpha (k)`$ (117)
where $`\gamma =n^\alpha (i)n_{}(i)n_{}(i)`$.

The solution of Eq. 111 is the following oneMancini and Avella (2003)

$$\begin{array}{c}G_p(k,\omega )=2\mathrm{i}\pi \mathrm{\Gamma }_p(k)\delta (\omega )\hfill \\ \hfill +\underset{i=1}{\overset{6}{^{}}}\frac{\sigma _p^{(i)}(k)}{1\mathrm{e}^{\beta \omega }}\left[\frac{1}{\omega E_{ip}(k)+\mathrm{i}\delta }\frac{\mathrm{e}^{\beta \omega }}{\omega E_{ip}(k)\mathrm{i}\delta }\right]\end{array}$$
(118)
where

$$\begin{array}{c}\mathrm{\Gamma }_p(k)=\frac{1}{2}\underset{\omega 0}{lim}\omega \hfill \\ \hfill \times \left[\theta (tt^{})P(i)P^{}(j)\theta (t^{}t)P^{}(j)P(i)\right]\end{array}$$
(119)
The primed sum is again restricted to values of $`i`$ for which $`E_{ip}(k)0`$.

The $`E_{ip}(k)`$ are the eigenvalues of the energy matrix $`\omega _p\left(k\right)`$

$`E_{1p}(k)`$ $`=2\mu +U2t[1+\alpha (k)]`$ (120)
$`E_{2p}(k)`$ $`=2\mu +U+2t[1+\alpha (k)]`$ (121)
$`E_{3p}(k)`$ $`=2\mu +{\displaystyle \frac{1}{2}}[UQ(k)]`$ (122)
$`E_{4p}(k)`$ $`=2\mu +{\displaystyle \frac{1}{2}}[3UQ(k)]`$ (123)
$`E_{5p}(k)`$ $`=2\mu +{\displaystyle \frac{1}{2}}[U+Q(k)]`$ (124)
$`E_{6p}(k)`$ $`=2\mu +{\displaystyle \frac{1}{2}}[3U+Q(k)]`$ (125)
with $`Q(k)=\sqrt{U^2+16t^2[1+\alpha (k)]^2}`$. The $`E_{ip}(k)`$ are zero only in isolated points of the parameter space ($`n`$, $`T`$, $`U`$) (see Eqs. 120 \- 125). According to this, the zero frequency function $`\mathrm{\Gamma }_p(k)`$ is identically zero except in these points, where it could be finite. In this treatment, for the sake of simplicity, we neglect these isolated points.

The spectral weights $`\sigma _{11p}^{(i)}(k)`$, computed through Eq. 22, have the following expressions

$`\sigma _{11p}^{(1)}(k)`$ $`={\displaystyle \frac{1}{2}}[I_{13p}(k)+I_{14p}(k)I_{15p}(k)I_{16p}(k)]`$ (126)
$`\sigma _{11p}^{(2)}(k)`$ $`={\displaystyle \frac{1}{2}}[I_{13p}(k)+I_{14p}(k)+I_{15p}(k)I_{16p}(k)]`$ (127)
$`\sigma _{11p}^{(3)}(k)`$ $`={\displaystyle \frac{t[1+\alpha (k)]}{Q(k)}}I_{15p}(k)`$
$`+{\displaystyle \frac{Q(k)U}{4Q(k)}}[2I_{11p}(k)2I_{14p}(k)+I_{16p}(k)]`$ (128)
$`\sigma _{11p}^{(4)}(k)`$ $`={\displaystyle \frac{t[1+\alpha (k)]}{Q(k)}}I_{15p}(k)+{\displaystyle \frac{Q(k)+U}{4Q(k)}}I_{16p}(k)`$ (129)
$`\sigma _{11p}^{(5)}(k)`$ $`={\displaystyle \frac{t[1+\alpha (k)]}{Q(k)}}I_{15p}(k)`$
$`+{\displaystyle \frac{Q(k)+U}{4Q(k)}}[2I_{11p}(k)2I_{14p}(k)+I_{16p}(k)]`$ (130)
$`\sigma _{11p}^{(6)}(k)`$ $`={\displaystyle \frac{t[1+\alpha (k)]}{Q(k)}}I_{15p}(k)+{\displaystyle \frac{Q(k)U}{4Q(k)}}I_{16p}(k)`$ (131)
In order to finally compute the Green’s function $`G_{11p}(k,\omega )`$ we should fix the internal parameter $`\gamma `$. The determination of the parameter $`\gamma `$ requires the computation of another Green’s function $`G_D(k,\omega )=𝒯[n(i)D(j)]`$, where $`D(i)=n_{}(i)n_{}(i)`$. Otherwise, we could resort to the following local algebra constraint

$$p(i)p^{}(i)=1n+D$$
(132)
It is worth noting that within the pair channel we can also compute the parameter $`d`$ directly from its definition: $`d=p^\alpha (i)p^{}(i)`$. This finally opens the possibility to compute fully self-consistently the fermionic and bosonic sectors at once. It can be shown that the two procedures (the fully self-consistent one and the one presented at length in the previous sections) give exactly the same results, as it should be according to the exact nature of the proposed treatment.

#### IV.2.2 The $`t`$-$`J`$ model

After Eqs. 64, the thermal causal Green’s functions $`G_c(q,\omega )=𝒯\left[\phi (i)\phi ^{}(j)\right]`$ and $`G_{sk}(q,\omega )=𝒯\left[\phi _k(i)\phi _k^{}(j)\right]`$ satisfy the following equations

$`\left[\omega \omega _c\left(q\right)\right]G_c(q,\omega )=I_c(q)`$ (133)
$`\left[\omega \omega _s\left(q\right)\right]G_{sk}(q,\omega )=I_s(q)`$ (134)
where the relevant entries of the normalization matrices $`I_c(q)=[\phi (i),\phi ^{}(j)]_{E.T.}`$ and $`I_s(q)=[\phi _k(i),\phi _k^{}(j)]_{E.T.}`$ are

$`I_{11c}(q)`$ $`=0`$ (135)
$`I_{12c}(q)`$ $`=4\left[1\alpha (q)\right]C_{11}^\alpha `$ (136)
$`I_{22c}(q)`$ $`=0`$ (137)
and

$`I_{11s}(q)`$ $`=0`$ (138)
$`I_{12s}(q)`$ $`=4\left[1\alpha (q)\right]C_{11}^\alpha `$ (139)
$`I_{13s}(q)`$ $`=0`$ (140)
$`I_{14s}(q)`$ $`={\displaystyle \frac{4}{3}}\left[1\alpha (q)\right]\chi _s^\alpha `$ (141)
with

$$\chi _s^\alpha =\underset{k=1}{\overset{3}{}}\nu _k^\alpha (i)\nu _k(i)$$
(142)
The spin rotational invariance makes $`G_{sk}`$ independent from the index $`k`$. According to this and for the sake of simplicity, we have omitted in the text the index $`k`$ in the expressions of the related quantities: $`I_s`$, $`E_s`$, $`\sigma _s`$ and $`\mathrm{\Gamma }_s`$.

The general solutions of Eqs. 133 are the following onesMancini and Avella (2003)

$`G_c(q,\omega )=2\mathrm{i}\pi \mathrm{\Gamma }_c(q)\delta (\omega )`$
$`+\underset{i=1}{\overset{2}{{\displaystyle ^{}}}}{\displaystyle \frac{\sigma _c^{(i)}(q)}{1\mathrm{e}^{\beta \omega }}}\left[{\displaystyle \frac{1}{\omega E_{ic}(q)+\mathrm{i}\delta }}{\displaystyle \frac{\mathrm{e}^{\beta \omega }}{\omega E_{ic}(q)\mathrm{i}\delta }}\right]`$ (143)
$`G_{sk}(q,\omega )=2\mathrm{i}\pi \mathrm{\Gamma }_s(q)\delta (\omega )`$
$`+\underset{i=1}{\overset{4}{{\displaystyle ^{}}}}{\displaystyle \frac{\sigma _s^{(i)}(q)}{1\mathrm{e}^{\beta \omega }}}\left[{\displaystyle \frac{1}{\omega E_{is}(q)+\mathrm{i}\delta }}{\displaystyle \frac{\mathrm{e}^{\beta \omega }}{\omega E_{is}(q)\mathrm{i}\delta }}\right]`$ (144)
where

$`\mathrm{\Gamma }_c(q)={\displaystyle \frac{1}{2}}\underset{\omega 0}{lim}\omega `$
$`\times \left[\theta (tt^{})\phi (i)\phi ^{}(j)\theta (t^{}t)\phi ^{}(j)\phi (i)\right]`$ (145)
$`\mathrm{\Gamma }_s(q)={\displaystyle \frac{1}{2}}\underset{\omega 0}{lim}\omega `$
$`\times \left[\theta (tt^{})\phi _k(i)\phi _k^{}(j)\theta (t^{}t)\phi _k^{}(j)\phi _k(i)\right]`$ (146)
The primed sum is again restricted to values of $`i`$ for which $`E_{ic,s}(k)0`$.

The $`E_{ic}(q)`$ and $`E_{is}(q)`$ are the eigenvalues of the energy matrices $`\omega _c(q)`$ and $`\omega _s(q)`$, respectively

$`E_{1c}(q)`$ $`=2t\sqrt{2[1\alpha (q)]}`$ (147)
$`E_{2c}(q)`$ $`=2t\sqrt{2[1\alpha (q)]}`$ (148)
and

$`E_{1s}(q)`$ $`=2t\sqrt{2[1\alpha (q)]}`$ (149)
$`E_{2s}(q)`$ $`=2t\sqrt{2[1\alpha (q)]}`$ (150)
$`E_{3s}(q)`$ $`=2J\sqrt{2[1\alpha (q)]}`$ (151)
$`E_{4s}(q)`$ $`=2J\sqrt{2[1\alpha (q)]}`$ (152)
According to this, the primed sums in Eqs. IV.2.2 do not contain the elements for $`q=0`$ and the zero frequency functions $`\mathrm{\Gamma }_{c,s}(q)`$ reduce to the constants $`\mathrm{\Gamma }_c=\mathrm{\Gamma }_c(0)`$ and $`\mathrm{\Gamma }_s=\mathrm{\Gamma }_s(0)`$, respectively. The spectral weights $`\sigma _{c,s}^{(i)}(k)`$, computed through Eq. 22 have the following expressions

$`\sigma _{11c}^{(1)}(q)`$ $`=\sigma _{11c}^{(2)}(q)=C_{11}^\alpha \sqrt{2}\sqrt{1\alpha (q)}`$ (153)
$`\sigma _{12c}^{(1)}(q)`$ $`=\sigma _{12c}^{(2)}(q)=2C_{11}^\alpha \left[1\alpha (q)\right]`$ (154)
$`\sigma _{22c}^{(1)}(q)`$ $`=\sigma _{22c}^{(2)}(q)=2\sqrt{2}C_{11}^\alpha \left[1\alpha (q)\right]^{\frac{3}{2}}`$ (155)
$`\sigma _{11s}^{(1)}(q)`$ $`=\sigma _{11s}^{(2)}(q)=C_{11}^\alpha \sqrt{2}\sqrt{1\alpha (q)}`$ (156)
$`\sigma _{11s}^{(3)}(q)`$ $`=\sigma _{11s}^{(4)}(q)={\displaystyle \frac{1}{3}}\chi _s^\alpha \sqrt{2}\sqrt{1\alpha (q)}`$ (157)
In order to finally compute the Green’s function $`G_c(q,\omega )`$ we should fix the zero frequency constant $`\mathrm{\Gamma }_c`$. This latter cannot be directly connected to any correlation function at this level and its determination requires the use of local algebra constraints. Let us consider the correlation function $`D_c(q,\omega )=\phi (i)\phi ^{}(j)`$ which is linked to the causal Green’s function through the spectral theorem

$$\begin{array}{c}D_c(q,\omega )=\left(1+\mathrm{tanh}\frac{\omega }{2T}\right)\mathrm{}\left[G_c(q,\omega )\right]\hfill \\ \hfill =2\pi \mathrm{\Gamma }_c\delta (\omega )+2\pi \underset{i=1}{\overset{2}{^{}}}\delta \left[\omega E_{ic}(q)\right]\frac{\sigma _c^{(i)}(q)}{1\mathrm{e}^{\beta E_{ic}(q)}}\end{array}$$
(158)
Then, we can use the following relations in order to compute $`\mathrm{\Gamma }_c`$

$`D_{11c}(i,i)`$ $`=\nu (i)\nu (i)=n`$ (159)
$`D_{12c}(i,i)`$ $`=\nu (i)\mathrm{\Pi }^{}(i)=2C_{11}^\alpha `$ (160)
$`D_{22c}(i,i)`$ $`=\mathrm{\Pi }(i)\mathrm{\Pi }^{}(i)=2\left(n\chi _c^\alpha \right)`$ (161)
$`\chi _c^\alpha `$ is directly connected to the Green’s function

$$\chi _c^\alpha =\nu (i)\nu ^\alpha (i)$$
(162)
In order to finally compute the Green’s function $`G_{11s}(k,\omega )`$ we should fix the internal parameter $`\chi _s^\alpha `$ and the zero frequency constant $`\mathrm{\Gamma }_{11s}`$.

The parameter $`\chi _s^\alpha `$ is directly connected to the Green’s function. Let us consider the correlation function $`D_{sk}(q,\omega )=\phi _k(i)\phi _k^{}(j)`$ which is linked to the causal Green’s function through the spectral theorem

$$\begin{array}{c}D_{sk}(q,\omega )=\left(1+\mathrm{tanh}\frac{\omega }{2T}\right)\mathrm{}\left[G_{sk}(q,\omega )\right]\hfill \\ \hfill =2\pi \mathrm{\Gamma }_s\delta (\omega )+2\pi \underset{i=1}{\overset{4}{^{}}}\delta \left[\omega E_{is}(q)\right]\frac{\sigma _s^{(i)}(q)}{1\mathrm{e}^{\beta E_{is}(q)}}\end{array}$$
(163)
Then,

$$\chi _s^\alpha =3\nu _3(i)\nu _3^\alpha (i)$$
(164)
The zero frequency constant $`\mathrm{\Gamma }_s`$ cannot be directly connected to any correlation function at this level and its determination requires the use of local algebra constraints. In particular, we can use the following relation

$$D_{11sk}(i,i)=\nu _3(i)\nu _3(i)=n$$
(165)
Equations 164 and 165 constitute a complete set of self-consistent equations which allow to compute the Green’s function $`G_{11sk}(q,\omega )`$.

We have the following expressions for the relevant zero frequency constants and correlators

$`\mathrm{\Gamma }_{11c}=n+\mathrm{coth}{\displaystyle \frac{2t}{T}}C_{11}^\alpha `$ (166)
$`\mathrm{\Gamma }_{12c}=0`$ (167)
$`\mathrm{\Gamma }_{22c}=0`$ (168)
$`\mathrm{\Gamma }_{11s}=3n+3\mathrm{coth}{\displaystyle \frac{2t}{T}}C_{11}^\alpha +\mathrm{coth}{\displaystyle \frac{2J}{T}}\chi _s^\alpha `$ (169)
$`\chi _c^\alpha =n+2\mathrm{coth}{\displaystyle \frac{2t}{T}}C_{11}^\alpha `$ (170)
$`\nu ^\alpha (i)\mathrm{\Pi }^{}(i)=2C_{11}^\alpha `$ (171)
$`\mathrm{\Pi }^\alpha (i)\mathrm{\Pi }^{}(i)=4\mathrm{coth}{\displaystyle \frac{2t}{T}}C_{11}^\alpha `$ (172)
$`\chi _s^\alpha ={\displaystyle \frac{3n+6\mathrm{coth}\frac{2t}{T}C_{11}^\alpha }{12\mathrm{coth}\frac{2J}{T}}}={\displaystyle \frac{3\chi _c^\alpha }{12\mathrm{coth}\frac{2J}{T}}}`$ (173)
## V The Results

In the previous two sections, we have given a detailed summary of the analytical calculations which lead to computable expressions for internal parameters, zero frequency constants, correlation and response functions for the fermionic and bosonic (charge, spin and pair) sectors of the Hubbard and $`t`$-$`J`$ models. In this section, we analyze the related results for the local, single-particle, thermodynamical and response properties. Hereafter, $`t`$ will be used as energy scale.

### V.1 The $`\mathrm{\Delta }`$ and $`p`$ parameters

The $`\mathrm{\Delta }`$ and $`p`$ parameters are computed by solving the system of self-consistent Equations 30 and 33. The behavior of $`\mathrm{\Delta }`$, as function of the filling $`n`$, is shown in Fig. 1 (top panel). $`\mathrm{\Delta }`$ is defined as the difference between the hopping amplitudes computed between empty and singly occupied sites (i.e., $`\xi (i)\xi ^{}(i)`$) and between singly and doubly occupied sites (i.e., $`\eta (i)\eta ^{}(i)`$). In this sense, $`\mathrm{\Delta }`$ gives a measure of the electron mobility. The following features have been observed: the hopping amplitude coming from $`\xi `$ ($`\eta `$) prevails below (above) half-filling and is always negative \[see inset Fig. 1 (top panel)\]; the absolute value of $`\mathrm{\Delta }`$ diminishes on increasing the temperature $`T`$ and on decreasing the Coulomb repulsion $`U`$. Below half filling, the electrons move preferably among empty and single occupied sites; the same happens to the holes above half filling. This explains the prevalence of one hopping amplitude per each region of filling. In addition, by decreasing the temperature or by increasing the Coulomb repulsion we can effectively reduce the number of doubly occupied sites present in the system and, consequently, favor the mobility of the electrons and increase the absolute value of $`\mathrm{\Delta }`$. These are also the reasons behind the filling dependence: its absolute value first increases with the number of available almost free moving particles ($`n0.5`$), then decreases when the number of particles is such to allow some double occupancies and to force the reduction of the number of empty sites ($`0.5n1`$) and finally vanishes when, in average, the sites are all singly occupied ($`n=1`$); the behavior above half filling is equivalent and is related to the holes.

The behavior of $`C_{11}^\alpha `$ within the $`t`$-$`J`$ model, as function of the filling $`n`$, is shown in Fig 1 (bottom panel). Below half filling and for rather high values of the Coulomb repulsion $`U`$, $`\mathrm{\Delta }`$ and $`C_{11}^\alpha `$ exactly coincide as the $`\eta `$ component of $`\mathrm{\Delta }`$ completely vanishes. According to this, all the features discussed above for the parameter $`\mathrm{\Delta }`$ are also observed for $`C_{11}^\alpha `$ ($`J`$ acts like $`1/U`$). Anyway, we have to report a really small dependence on the exchange interaction $`J`$ that is fully effective only at $`n=1`$ where the mobility is already zero as no empty site at all is left. It is worth noting that the exact result for the 2-site system at $`T=2`$ almost exactly reproduces the Exact Diagonalization (ED)Dagotto et al. (1992) data for a $`4\times 4`$ cluster at $`T=0`$. This shows that, by increasing the temperature in a 2-site system it is possible to mime a cluster of bigger size at a lower temperature, at least as regards some of its properties. This can be understood by thinking to the level spacing in the two cases (i.e., bigger the cluster lower the spacing) and to the value of temperature needed to excite those levels. Clearly, not all the properties are in a so strong relation with the relative level positions, but depend on the absolute positions of them. Anyway, this also shows that the relative energy scales/levels are already present in the 2-site system. We have used a value of $`J`$ twice smaller than the one used within the numerical analysis according to the difference, between the 2-site system and any larger system, regarding the value of the exchange energy $`J`$ appearing in the derivation of the $`t`$-$`J`$ model from the Hubbard oneCastellani et al. (1979). Actually, the derivation is pathological just for the 2-site system and gives a $`J`$ larger of a factor $`2`$.

As regards the $`p`$ parameter, in Fig 2 we report its dependence on filling $`n`$ for two values of the on-site Coulomb repulsion $`U=0`$ and $`U=4`$ and for two temperatures $`T=0`$ and $`T=2`$. At low temperatures and for $`n1.26`$, the value of the $`p`$ parameter is mainly negative according to the antiferromagnetic nature of the spin correlations in the ground state; at higher temperatures the antiferromagnetic spin correlations get weaker and the sign of the parameter changes. The filling $`n`$, the temperature $`T`$ and the Coulomb interaction $`U`$ rule the balance between the spin, the charge and the pair correlations.

### V.2 The chemical potential and the density of states

The chemical potential $`\mu `$ is one of the self-consistent parameters and has been computed by solving the system of self-consistent Equations 30 and 33. As a function of $`U`$ the chemical potential increases for all studied values of temperature $`T`$ and filling $`n`$. In particular, we have the maximum increment for $`U20`$; for higher values of $`U`$ the chemical potential saturates. A similar behavior has been reported by numerical data on larger clustersDagotto et al. (1991, 1992); Furukawa and Imada (1992, 1993). In Fig. 3, we show the filling dependence of the chemical potential for some values of the temperature $`T`$ at $`U=4`$ (top panel) and for $`U=4`$, $`12`$ and $`T=1`$, $`2`$ (bottom panel). The Finite Temperature Lanczos (FTL) data ($`4\times 4`$) are taken from Ref. Prelovsek, . At low temperatures, the chemical potential has a step-like behavior which recalls the discreteness of the energy levels. At higher temperatures, the step-like feature is smeared out due to the thermal hybridization of the energy levels. It is again really noteworthy the agreement between the 2-site results and the numerical data. No temperature adjustment has shown to be necessary in this case as the temperature is already quite high \[cfr. Fig. 12 (bottom)\]. Generally, the agreement is better for high temperatures and high values of the Coulomb repulsion. Anyway, there is a competition between the temperature smearing effects on the level spacing and the possibility to access the states with energy $`J_U`$ at lower temperatures for higher values of $`U`$ ($`J_U1/U`$). Actually, the number of accessible states is comparable for medium temperatures and too large for very high temperatures. This should explain the deviations for $`U=4`$ and $`T=1`$, and $`U=12`$ and $`T=2`$ at low doping.

By looking at these plots we could be induced to consider the possibility of a metal-insulator transition driven by the Coulomb repulsion $`U`$ and controlled by the temperature $`T`$. Obviously, no such transition is possible in a finite system that is always in a metal paramagnetic state for any finite or zero value of the Coulomb repulsion $`U`$; we will show, in Sec. V.8, that the Drude weight is always finite. Such tricky behavior of the chemical potential can be understood by thinking at the nature of the grand canonical ensemble we are using for our thermal averages. As we are very far from the thermodynamic limit, there is no equivalence at all among the micro-, the grand- and the canonical ensembles. In particular, whenever we speak about temperature and chemical potential in a finite system we have just to think in terms of mixtures of quantum mechanical states with different energies and numbers of particles. To get a deeper comprehension of this issue we can define a thermodynamic density of states $`N_T(\mu )`$ through the following equation

$$n=_{\mathrm{}}^\mu 𝑑\omega N_T(\omega )$$
(174)
which is a generalization to finite temperatures of the identical relation existing, at zero temperature, between the filling $`n`$ and the usual density of states $`N(\omega )`$. In general, the two densities of states coincide only at zero temperature and for systems, like the non-interacting ones, whose density of states is independent on the chemical potential $`\mu `$. From the definition (174), we can simply compute the thermodynamic density of states $`N_T(\mu )`$ by differentiating the filling $`n`$ with respect to the chemical potential $`\mu `$ \[$`N_T(\mu )=\text{d}n/\text{d}\mu `$\]. It is worth mentioning that the thermodynamic density of states $`N_T(\mu )`$, apart from a factor $`\frac{1}{n^2}`$, is simply the compressibility of the system. A vanishing value of this quantity denotes the impossibility for the system to accept more particles and, obviously, the presence of a gap. The thermodynamic density of states $`N_T(\mu )`$ enormously facilitates the comprehension of the chemical potential features that are singled out by the differentiation procedure. The usual density of states $`N(\omega )`$ is computed as follows

$$N(\omega )=\frac{1}{2}\underset{n=1}{\overset{4}{}}\underset{k}{}\delta \left[\omega E_n(k)\right]\sigma _{cc}^{(n)}(k)$$
(175)
Both densities of states are reported in Figs. 4 and 5 as function of the chemical potential ($`\mu `$) and the frequency ($`\omega `$), respectively. We have computed $`N(\omega )`$ by substituting, for obvious graphical reasons, the Dirac deltas ($`\delta (\omega )`$) with Lorentzian functions ($`\frac{1}{\pi }\frac{\delta }{\omega ^2+\delta ^2}`$) with $`\delta =0.1`$. In reporting the thermodynamic density of states $`N_T(\mu )`$ we have scaled the chemical potential $`\mu `$ by its value at half filling ($`U/2`$) as we want to make the comparison for this particular value of $`n`$. While $`N(\omega )`$ always presents a gap at $`\omega =0`$ (the overlapping tails are due to the finite $`\delta `$ and the non-zero kinetic energy ensures a metallic behavior), $`N_T(\mu )`$ shows a gap only above a certain value of the Coulomb repulsion for a fixed $`T`$ (see Fig. 4) and below a certain temperature for a fixed $`U`$ (see Fig. 5). We want to emphasize that, at high temperature, $`N_T(\mu )`$ is capable to mime the behavior expected for the lattice system. In particular, $`N_T(\mu )`$ presents two well defined structures (i.e., the two Hubbard subbands) that continuously separate on increasing $`U`$ (see Fig. 4). This behavior is the one expected, after many analytical and numerical results, for the bulk at dimension greater than one and it is a realization of the metal-insulator transition according to the Mott-Hubbard mechanism. It is worth noting that, for this small system, the critical Coulomb repulsion, at which the gap appears, is a function of the temperature $`T`$.

The density of states $`N(\omega )`$ is also reported in order to give an idea of the positions of the poles and of the relative intensity of the spectral weights, although we had to cut the highest peaks in drawing the pictures. We can observe two relevant features. In the non-interacting case ($`U=0`$) only two poles/peaks are present; those coming from the band $`E_1(k)`$ and separated by the *bandwidth* $`W=4t`$. On increasing the Coulomb potential $`U`$ the three scales of energy present in the system clearly manifest themselves: the exchange interaction $`J_U`$, the *bandwidth* $`W`$ and the Coulomb repulsion $`U`$. The eight poles group in two main structures separated by $`U`$. Within any structure the four poles are separated according to the combination of $`\pm 2t`$ (*bandwidth* separation) and the presence or absence of the exchange interaction $`J_U`$. The *bandwidth* $`W`$ is obviously $`U`$ independent and generates a rigid two peak structure for any band $`E_n(k)`$. On the contrary, the exchange energy $`J_U`$ decreases on increasing $`U`$ with a consequent reduction of the resolution of the peaks within the Hubbard subbands. The other relevant feature is the redistribution of the total spectral weight on increasing the temperature. At low temperatures the poles/peaks that do not contain the exchange interaction $`J_U`$ have negligible spectral weights; the 2-site singlet is the ground state at half-filling. On increasing temperature the total spectral weight redistributes and the other poles get more and more spectral weight due to the thermal hybridization of the energy levels.

### V.3 The double occupancy $`D`$

In the non-interacting case ($`U=0`$) we have $`D=n^2/4`$. At zero temperature, the double occupancy $`D`$ vanishes for $`n0.5`$; below this value of filling we have, in average, less than one electron in the system and only for a finite temperature we can get some contributions by states with a finite double occupancy. At half filling we have the following exact formula for the double occupancy at zero temperature

$$D=\frac{2J_U}{U+8J_U}$$
(176)
The filling dependence of the double occupancy $`D`$ (top panel) is reported in Fig. 6 for several values of the Coulomb repulsion $`U`$ at temperature $`T=2`$. There are also reported some quantum Monte Carlo (qMC) dataMoreo et al. (1990) for a bigger cluster ($`12\times 12`$). As in the case of $`C_{11}^\alpha `$ within the $`t`$-$`J`$ model, we note that the exact results for the 2-site system at temperature $`T=1`$ very well reproduce the quantum Monte Carlo (qMC) data for the bigger cluster at $`T=1/6`$. The same explanation obviously holds. A comparison with the exact results from Bethe Ansatz (BA) are also reported in Fig. 6 (bottom panel). Also in this case the 2-site system manage to reproduce the *bigger* cluster data (actually, the Bethe Ansatz (BA) system is the 1D bulk) by increasing the temperature. The discrepancy at half filling can be understood as a consequence of the difference in the definition of the exchange energy $`J`$ discussed in the Section regarding the internal energy. Obviously, the discrepancy is larger where the exchange interaction is mainly effective (i.e., at half filling and for intermediate-strong values of the Coulomb interaction).

### V.4 Thermodynamics

#### V.4.1 The energy $`E`$, the specific heat $`C`$ and the entropy $`S`$

The energy per site $`E`$ is computed as thermal average of the Hamiltonian divided by the number of sites. Its filling dependence is reported in Fig. 7 (top panel) together with some Exact Diagonalization (ED) data for bigger clustersDagotto et al. (1992). As expected, the energy $`E`$ increases as the Coulomb repulsion $`U`$ increases. For $`n<1`$ and high values of the Coulomb repulsion the kinetic term prevails on the Coulomb one. The opposite behavior is observed for $`n>1`$. The behavior as a function of the temperature $`T`$ and Coulomb repulsion $`U`$ are reported in Fig. 7 (middle and bottom panels, respectively) together with some data coming from numerical analysis for bigger clustersMoreo et al. (1990); Dagotto et al. (1992); Duffy and Moreo (1997). By using a higher temperature for the 2-site data is again possible to get an extremely good agreement with the numerical results \[see Fig. 7 (top and bottom panels)\]. In the case of the temperature behavior instead, the agreement is obtained without tuning any parameter \[see Fig. 7 (middle panel)\].

The discrepancy at half filling and low temperatures \[see Fig. 7 (middle panel)\] is a consequence of the different value of the exchange energy $`J`$ appearing in the derivation of the 2-site $`t`$-$`J`$ model from the Hubbard one (see detailed discussion about the behavior of $`C_{11}^\alpha `$ reported in Fig. 3 (bottom panel), Sec. V.1). According to this, at half-filling and low temperatures we have also reported the 2-site results for $`U=24`$ and $`J=1/6`$ ($`t`$-$`J`$ model) and obtained the expected agreement. The reason why no effect is evident for lower fillings and high temperatures is that the exchange interaction is really effective only at half-filling (in average one spin per site is necessary) and low temperatures (the fluctuations should be small).

In Fig. 8 we show the temperature dependence of the specific heat $`C=dE/dT`$ for several values of $`U`$ at $`n=0.75`$ and $`n=1`$. For $`n1`$ we have three peaks at temperatures of the order half ($`0.3÷0.6`$) the scales of energies $`J_U`$, $`t`$ and $`U`$. This is in agreement with the behavior of a *pure* two- or three- gap system in the canonical ensemble. The effect of increasing the Coulomb repulsion $`U`$ is a better resolution of the three peaks. In fact, the first one moves towards lower temperatures, the second one is stable and the third moves to higher temperatures; this is in perfect agreement with their origins: $`1/U`$, $`t`$ and $`U`$, respectively. At half filling nothing changes except for the absence of the middle peak: the single occupied states do not contribute. It is worth mentioning that the kinetic peak appears as the energy levels corresponding to different values of the momentum are here discrete as in any finite system. No such a peak is present in infinite systems where the kinetic energy just spreads over a band the energy levels coming from the interactions.

The specific heat $`C`$ is a really valuable property to quantify how good is the $`t`$-$`J`$ model to describe the low energy dynamics of the Hubbard model. We have studied $`C`$ as a function of both $`4t^2/U`$ and $`J`$ in the Hubbard and $`t`$-$`J`$ model, respectively. We need a Coulomb repulsion $`U24t`$ ($`4t^2/U\frac{1}{6}t`$) to have the possibility of a faithful mapping at low temperatures; only for such values of $`U`$ we manage to sufficiently resolve the three scales of energy in the Hubbard model. For instance, at $`U=32J=0.125`$ the mapping can be absolutely trusted as shown in Fig. 9 where we report the specific heat for $`n=0.75`$. As regards the exchange and kinetic peaks the agreement is perfect; obviously the $`t`$-$`J`$ model cannot reproduce the Coulomb peak, which is extraneous to its dynamics. The absence of the kinetic peak at half filling is also perfectly reproduced (see Fig. 11).

In Fig. 10 (top panel) we report the behavior, in the $`t`$-$`J`$ model, of the specific heat as function of the temperature $`T`$ for different values of the exchange constant $`J`$. It is worth mentioning that in the $`t`$-$`J`$ model the two peaks in the specific heat (i.e., the exchange and kinetic peaks) can be easily studied separately as they just come from the corresponding terms of the Hamiltonian. This simple analysis cannot be performed in the Hubbard model where the exchange peak gets contribution from both terms of the Hamiltonian.

In Fig. 10 (bottom panel) we report the entropy $`S`$ of $`t`$-$`J`$ model computed by means of the usual thermodynamic relations from the specific heat $`C`$

$$S(T)=_0^T\frac{d\stackrel{~}{T}}{\stackrel{~}{T}}C(\stackrel{~}{T})$$
(177)
We can easily put in correspondence the peaks of the specific heat \[see Fig. 10 (middle panel)\] and the change in the slope of the entropy. Again, a lowering of the exchange energy $`J`$ helps resolving the scales of energy and makes much more visible the difference in the slopes.

In Fig. 11, it is reported the specific heat as a function of the temperature for the Hubbard model (top panel) and the $`t`$-$`J`$ model (bottom panel). For this latter, it has been also reported the entropy. The Finite Temperature Lanczos (FTL) data for the Hubbard model ($`4\times 4`$) and for the $`t`$-$`J`$ model ($`N=20`$) are taken from Ref. Bonca and Prelovsek, 2003 and Ref. Jaklic and Prelovsek, 2000, respectively. The agreement is absolutely noteworthy on the whole temperature range. Again, the discrepancy at low temperatures in the Hubbard model case is a consequence of the difference regarding the value of the exchange energy $`J`$. According to this, we have also reported the 2-site results for $`U=24`$ and $`J=1/6`$ ($`t`$-$`J`$ model) and obtained the expected agreement. In the $`t`$-$`J`$ model case, we have used a value of $`J`$ twice smaller than the one used within the numerical analysis according to the boundary conditions we applied to the 2-site system. These results show, once more, that the correct and necessary scales of energies are already present in the 2-site system.

#### V.4.2 The temperature dependence of <br>the chemical potential $`\mu `$ and the double occupancy $`D`$

The temperature dependence of the chemical potential $`\mu `$ is given in Fig. 12 (top panel). $`\mu `$ has only one maximum for $`U8`$ and two maxima and a minimum for higher values of $`U`$. In the $`t`$-$`J`$ model, the chemical potential $`\mu `$ shows a behavior at low temperatures in agreement with the one found for the Hubbard model (i.e., $`lim_{Ut}J_U=J`$). On the contrary, for high temperatures the divergence is upward instead of downward \[see Fig. 12 (bottom panel)\].

The temperature dependence of the chemical potential can be explained as follows. In the Hubbard model and at zero temperature, there exists a critical value of the Coulomb interaction, function of the filling, above which the potential energy becomes larger, in absolute value, than the kinetic one. This critical value is a decreasing function of the filling: lower is the doping higher the double occupancy, bigger the potential energy and smaller, in absolute value, the kinetic energy. Then, by increasing the filling we increase the energy and the chemical potential increases. At higher and higher temperatures, the system behaves like a free system at zero temperature and has a decreasing chemical potential. In particular, for a diverging temperature we have a negatively diverging chemical potential (i.e., $`\mu =dE/dnTdS/dn`$) as the entropy is an increasing function of the filling at very high temperatures (i.e., in the almost free case) and $`n<1`$. For intermediate temperatures, the chemical potential has maxima and minima in coincidence with the peaks in the specific heat. Any peak in the specific heat marks the temperature at which the system gets the freedom to occupy a new state. The spin (i.e., $`TJ_U`$) and the charge (i.e., $`TU`$) peaks reflect in two maxima as they signal the availability of triplet states and of doubly occupied states, respectively. Both states has zero kinetic energy. On the contrary, the kinetic (i.e., $`Tt`$) peak signals the availability of further singly occupied states which maximize the absolute value of the kinetic energy. This reasoning strongly relies on the region of filling the figure refers to (i.e., $`n0.75`$). For other doping regions the situation can be quite different owing to the behavior in filling of the double occupancy, the entropy and the positions and existence of the peaks in the specific heat. These considerations can also explain the quite different behavior we have in the $`t`$-$`J`$ model. At $`n=0.75`$ the $`t`$-$`J`$ model is approaching the full-filled system (i.e., $`n=1`$), on the contrary the Hubbard model is approaching the half-filled system. At zero temperature for $`n>0.5`$, increasing the filling we lowers, in absolute value, the kinetic energy as the single electron states, which are the only ones with a finite kinetic energy, are replaced by the singlet and triplet states: the chemical potential is positive. Anyway, the exchange energy can effectively lower its value on increasing $`J`$. At very high temperatures (i.e., in the almost free case), the entropy is now a decreasing function of the filling and leads to the positive divergence. At intermediate temperatures, the reasoning is identical to the one given for the Hubbard model except for the obvious absence of the charge peak in the $`t`$-$`J`$ model.

In Fig. 12 (bottom), it is reported the chemical potential as function of the temperature $`T`$ for $`U=12`$ and $`n=0.5`$, $`0.6`$, $`0.7`$, $`0.8`$ and $`0.9`$. The FTL data ($`4\times 4`$) are from Ref. Prelovsek, . The agreement is really excellent except at low temperatures and for the higher values of the filling ($`n=0.8`$ and $`n=0.9`$). At low temperatures, owing to the finite level spacing the chemical potential has a step behavior and, in particular, between $`n=0.5`$ and $`n=0.9`$ has practically no variation (cfr. Fig 3). At the higher values of the filling ($`n=0.8`$ and $`n=0.9`$), the number of 2-site states excited at high temperatures is larger than those needed to mime the cluster.

The temperature dependence of the double occupancy $`D`$ at half filling shows a minimum at around $`T1`$ before saturating at the non-interacting value $`n^2/4`$ for very high temperaturesGeorges and Krauth (1993); Schulte and Böhm (1996). This minimum coincides with the kinetic peak in the specific heat. As already discussed for the chemical potential, this peak marks the freedom for the system to occupy further singly occupied states which have, obviously, zero double occupancy and, therefore, lower its total value.

#### V.4.3 The crossing points

In Fig. 13 we report the temperature dependence of the specific heat for $`U5t`$ at half filling (top panel) and $`n=0.75`$ (bottom panel). For these low values of the Coulomb repulsion we identify three crossing points/regions where the specific heat is almost independent from the value of the Coulomb repulsion $`U`$Vollhardt (1997); Chandra et al. (1999); Mancini et al. (1999). To compute the positions and heights of the crossing points we have expanded the specific heat $`C`$ as function of the Coulomb repulsion $`U`$. We have got

$`C(T,U,n)=C_0(T,n)+C_2(T,n)U^2+O\left(U^4\right)`$ (178a)
$`C_0(T,n=1)=2\beta ^2t^2\left[1\mathrm{tanh}^2\left(\beta t\right)\right]`$ (178b)
$$C_2(T,n=1)=\frac{\mathrm{e}^{2\beta t}}{8T^4\left(1+\mathrm{e}^{2\beta t}\right)^6}\left[48t^2\mathrm{e}^{2\beta t}\left(13\mathrm{e}^{2\beta t}+\mathrm{e}^{4\beta t}\right)+2tT\left(1+20\mathrm{e}^{2\beta t}20\mathrm{e}^{6\beta t}\mathrm{e}^{8\beta t}\right)+2T^2\left(1+\mathrm{e}^{2\beta t}\right)^4\right]$$
(179)
The temperatures $`T^{}(n)`$ at which the crossing point can be observed are determined by the equation $`C_2(T,n)=0`$. At half filling and at $`n=0.75`$, this equation gives the following results:

$`T_1^{}(n=1)0.499C_00.563`$ (180a)
$`T_2^{}(n=1)0.997C_00.841`$ (180b)
$`T_3^{}(n=1)2.058C_00.376`$ (180c)
$`T_1^{}(n=0.75)0.951C_00.523`$ (181a)
$`T_2^{}(n=0.75)2.024C_00.342`$ (181b)
Vollhardt found an approximate formula for $`C_0`$ and $`T^{}`$, as function of the dimensionality of the system, for the higher temperature crossing point at half fillingChandra et al. (1999). The values that this formula gives for $`d=1`$, the *dimensionality* of the 2-site system within periodic boundary conditions, are close, as regards $`C_0`$, to the exact ones of Eq. (180c).

It is worth noting that higher is the temperature narrower is the crossing region; ranging from quite wide regions for the two lower temperatures to a very sharp point for the higher temperature. Also the chemical potential (Fig. 3 (top panel) at $`n0.275,0.75,1.25,1.725`$) and the double occupancy show quite clear crossing points in (no dependence on) temperature.

### V.5 Charge, spin and pair correlations

The charge (i.e., $`\chi _c^\alpha `$), spin (i.e., $`\chi _s^\alpha `$) and pair (i.e., $`d`$) correlation functions contain information regarding the spatial distributions of the corresponding quantities. Obviously, we can only consider the quantum fluctuations in the paramagnetic state, which is the only admissible equilibrium state on a finite cluster.

The charge correlation function $`\chi _c^\alpha `$ is shown in Fig. 14 (top) as function of the interaction $`U`$ for different values of the temperature $`T`$ and $`n=1`$. In the non-interacting case ($`U=0`$) and at zero temperature we have

$$\chi _c^\alpha n^2=\{\begin{array}{ccc}\frac{1}{2}n^2\hfill & \text{if}\hfill & n1\hfill \\ \frac{1}{2}(2n)^2\hfill & \text{if}\hfill & n1\hfill \end{array}$$
(182)
as only the singlet states with one electron ($`n1`$) \[one hole ($`n1`$)\] per site contribute. In fact, these are the states which lower the most the internal energy as they maximize the absolute value of the kinetic one. In the strongly interacting limit ($`U\mathrm{}`$) and at zero temperature we have

$$\chi _c^\alpha n^2=\{\begin{array}{ccc}n^2\hfill & \text{if}\hfill & n0.5\hfill \\ (1n)^2\hfill & \text{if}\hfill & 0.5n1.5\hfill \\ (2n)^2\hfill & \text{if}\hfill & n1.5\hfill \end{array}$$
(183)
as only the states with a single electron ($`n0.5`$) \[a single hole ($`n1.5`$)\] or only with one electron per site ($`0.5n1`$) \[one hole per site ($`1n1.5`$)\] contribute. In fact, no double occupancy is allowed and there is no gain in the kinetic energy if the singlet states are used as $`J_U0`$. At intermediate values of the coupling $`U`$ the double occupied states play a relevant role and we found results between the limiting cases reported above: on increasing the coupling $`U`$ the double occupancy diminishes and consequently the states with one electron per site increase their contribution and the value of $`\chi _c^\alpha `$ (see Fig. 14 (top)). On increasing the temperature $`T`$, more and more states become available and the correlation tends to its reducible part, $`n^2`$.

The spin correlation function $`\chi _s^\alpha `$ is reported in Fig. 14 (bottom) for the same parameters chosen for $`\chi _c^\alpha `$. In the non-interacting case ($`U=0`$), we have $`\chi _s^\alpha =\chi _c^\alpha `$: in absence of interaction the charge and spin behave coherently. On increasing the interaction potential $`U`$ the singlet states are favored and the spin correlations get enhanced ($`J_U`$ gets closer and closer to $`4t^2/U`$ which is the scale of energy of the spin excitations, see Fig. 15). Then, further increasing $`U`$ the value of $`J_U`$ tends to zero and the spin correlations get suppressed except at zero temperature \[see Fig. 14 (bottom)\]. A rapid suppression of the spin fluctuations can be caused also by the increment of the temperature.

The pair correlation function $`d`$ is reported in Fig. 16 as function of the temperature $`T`$ at $`n=1`$ and for $`U=0`$, $`2`$, $`4`$, $`8`$, $`12`$ and $`16`$. Obviously, $`d`$ decreases on increasing the Coulomb interaction $`U`$ and is maximum at half-filling where the singlet state is the favorite one. On increasing the temperature $`T`$ (see Fig. 16), more and more states become available and the pair correlation function $`d`$ tends to zero.

### V.6 Charge and spin susceptibilities

The charge and spin susceptibilities $`\chi _\mu (k,\omega )`$ are given by

$$\chi _\mu (k,\omega )=[n_\mu (i)n_\mu (j)]$$
(184)
where no summation is implied by repeated indices. By means of the expression of the causal Green’s function given in Sec. IV.2.1 is immediate to compute the following static susceptibilities through the spectral theorem

$`t\chi _c(k,0)`$ $`=C_{cc}^\alpha {\displaystyle \frac{U}{2t}}d`$ (185)
$`t\chi _s(k,0)`$ $`=C_{cc}^\alpha {\displaystyle \frac{1}{12}}{\displaystyle \frac{U}{t}}\chi _s^\alpha `$ (186)
where $`d`$ and $`\chi _s^\alpha `$ are the first-neighbor pair and spin correlation functions, respectively.

In Fig. 17, the charge $`\chi _c(k,0)`$ and spin $`\chi _s(k,0)`$ susceptibilities are reported as function of the temperature $`T`$ for different values of the Coulomb repulsion $`U`$ and $`n=1`$. As it results clear by Eqs. 185, the two susceptibilities result identical in the non-interacting case $`U=0`$: once more charge and spin behave coherently in absence of interaction. On increasing the Coulomb repulsion, the two susceptibilities behave in opposite manners: the charge susceptibility is strongly suppressed, in particular at half-filling where a charge gap develops in the strongly interacting limit ($`U\mathrm{}`$); the spin susceptibility is greatly enhanced, in particular at half-filling where the singlet state is the favorite one. The behavior of the two susceptibilities as functions of the temperature deserve special attention as it is directly connected with the formation of a gap in the corresponding channel. The spin susceptibility shows the typical paramagnetic behavior: plateau at low temperatures and Curie tail at high temperatures. The charge susceptibility instead shows the presence of a gap in the strongly interacting limit ($`U\mathrm{}`$). The small upturn at low temperatures and medium-high values of the Coulomb repulsion $`U`$ is due to the finite level-spacing of the system. In particular, it arises from the intra-subband gaps characteristic of the scale of energy $`J_U`$.

In the $`t`$-$`J`$ model, we simply have

$`t\chi _c(k,0)`$ $`=C_{11}^\alpha `$ (187)
$`t\chi _s(k,0)`$ $`=C_{11}^\alpha {\displaystyle \frac{1}{3}}{\displaystyle \frac{t}{J}}\chi _s^\alpha `$ (188)
### V.7 Thermal compressibility

The thermal compressibility is defined as

$$\kappa =\frac{1}{n^2}\frac{n}{\mu }$$
(189)
By using some general quantum statistical relations, which can be established between the particle density $`n`$ and the chemical potential $`\mu `$ Feng and Mancini (2001), we can express the thermal compressibility in terms of the density-density correlation function

$$\kappa =\frac{1}{Tn^2}\left[\left(n+2Dn^2\right)+\left(\chi _c^\alpha n^2\right)\right]$$
(190)
In Fig. 18, the thermal compressibility $`\kappa `$ is reported as function of the temperature $`T`$ for different values of the Coulomb repulsion $`U`$ and $`n=1`$. The system is completely incompressible $`\kappa =0`$ when no more particles are allowed to enter the system: the chemical potential, which is a measure of the energy necessary to insert a new particle in the system, diverges. Obviously, the system is extremely eager to accept particles at very low fillings, in order to increase in absolute value the kinetic energy and lower the total one, and absolutely incompressible at $`n=2`$, when all the quantum states are filled. At zero temperature, the Coulomb repulsion makes incompressible also the states with commensurate fractional fillings $`n=0.5`$ and $`n=1.5`$. At half-filling, on increasing the Coulomb repulsion, the compressibility is rapidly suppressed according to the very high price in energy that should be paid to add one particle to the singlet state which is the one favored by the Coulomb repulsion. The most relevant feature for this property is the presence of a well-defined peak when it is plotted versus temperature. A finite, but low in comparison with the actual value of the Coulomb repulsion, temperature permits to overcome the suppression related to the formation of a charge gap. A further increment of the temperature makes available more and more states and the system is driven back to be incompressible.

It is possible to write the compressibility for the Hubbard and $`t`$-$`J`$ models in an identical way

$$\kappa =\frac{2}{Tn^2}\left(\mathrm{\Gamma }_{11c}n^2\right)$$
(191)
where $`\mathrm{\Gamma }_{11c}`$ is the zero frequency function in the charge channel and its ergodic value is exactly $`n^2`$. According to this, the compressibility is a direct measure of the ergodicity of the charge dynamics.

### V.8 Optical conductivity

In the framework of the linear response and by using the Ward-Takahashi identitiesMancini:98a, which relate the current-current propagator to the charge-charge one by exploiting the charge conservation, the optical conductivity is given by

$$\sigma _1(\omega )=D_w\delta (\omega )$$
(192)
with the Drude weight $`D_w`$ given by

$$D_w=4\pi \mathrm{e}^2a^2tC_{cc}^\alpha $$
(193)
For the 2-site system, the incoherent part of optical conductivity is zero as no contribution comes from the imaginary part of retarded current-current propagator.

In Fig. 19, the Drude weight $`D_w`$, normalized by $`2\pi \mathrm{e}^2a^2`$, is reported as function of the filling $`n`$ for different values of the Coulomb repulsion $`U`$ and $`T=1`$. At half-filling, the Coulomb repulsion tends to suppress the Drude weight and drives the system to be insulating. However, the Drude weight vanishes only in the limit $`U`$ infinite. Higher the temperature more states with no contribution to the kinetic energy result available and lower is the Drude weight.

## VI Ergodicity

One of the main issues on which this manuscript wishes to draw attention is the non-ergodicity of the charge and spin dynamics in the two-site Hubbard model. The non-ergodic dynamics in this system is due to the finite number of degenerate states available to the system when dealing with finite temperatures and/or incommensurate fillings and/or no interaction (cfr. Table 1 and Eq. 205). The finite level spacing confines the system to degenerate or non-degenerate states, according to the filling and interaction strength, making the dynamic ergodic or non-ergodic according to the degree of degeneracy of the ground state. The temperature then opens up the possibility for quite complicate mixtures, only possible according to our choice to work in the grand-canonical ensemble in order to get results comparable with larger systems for which the level spacing is thinner. The coupling to an heat and a particle reservoirs have given us the possibility to simulate a continuous tuning of the filling and temperature on a finite system (with a finite degree of freedom) as it is possible only in a bulk system with infinite degrees of freedom.

The zero-frequency constants for charge $`\mathrm{\Gamma }_{11c}`$ and spin $`\mathrm{\Gamma }_{11s}`$ channels have been plotted as functions of the filling $`n`$ for different values of the Coulomb repulsion $`U`$ and $`T=0.1`$ and as functions of the interaction $`U`$ for different values of temperature $`T`$ and $`n=1`$ in Figs. 20 and 21. Their ergodic values are $`n^2`$ and $`0`$, respectively. $`\mathrm{\Gamma }_{11c}`$ assumes its ergodic value only at zero temperature and commensurate fillings; in particular, only at integer commensurate fillings for any interaction strength and also at fractional ones for a finite interaction strength. On the other hand, $`\mathrm{\Gamma }_{11s}`$ assumes its ergodic value only at zero temperature and integer commensurate fillings for any interaction strength. At half-filling, while the charge dynamics is constrained to be ergodic in the strongly interacting limit ($`U\mathrm{}`$) by the formation of a charge gap, the spin dynamics gets more and more non-ergodic as larger the interaction strength is. The states that are left available in this condition have a quite different behavior once probed with respect to the charge or spin response.

## VII Conclusions

We have studied the 2-site Hubbard and $`t`$-$`J`$ models by means of the Green’s function and equations of motion formalism. The main results can be so summarized:

* We have got a complete basis of eigenoperators for the fermionic and bosonic sectors which could be used to get a controlled approximation in the study of the lattice case by means of any approximation that strongly relies on the choice of the basic field.
* We have identified the eigenoperator responsible for the appearance of the exchange scale of energy $`J`$.
* We have illustrated, once moreMancini and Avella (2003), as the local algebra constrains can properly fix the representation and easily give the values of the bosonic correlations that appear as internal parameters in the fermionic sector and of the zero frequency functions that appear in the bosonic sector. This also permits, while studying the fermionic sector, to avoid the opening of the bosonic one (i.e., the charge, the spin and pair channels) and all the heavy calculations required to solve it.
* We have explored the possibility of the 2-site systems to mime the behavior of larger clusters, as regards some of their physical properties, by using an higher temperature and by exploiting the qualitative properties of the level spacing. We have to report that the exact results of the 2-site system managed to reproduce the Mott-Hubbard *MIT* of the bulk, the behavior of some local quantities (i.e., the double occupancy, the local magnetic moment and the kinetic energy) and of some thermodynamic quantities (i.e., the energy, the specific heat and the entropy) of larger clusters. This shows that the necessary energy scales are already present in the 2-site system. According to this, as already said in the first point above, we strongly believe that the operatorial basis that exactly solves this system can also give excellent results if used for the bulk.
* The study of the specific heat has given many valuable information regarding: the scales of energy present in the system, their origin, interaction and possibility of resolution, the range of parameters for which the $`t`$-$`J`$ model faithfully reproduces the low energy/temperature behavior of the Hubbard model, the explanation for the temperature dependence of the local properties, the existence of crossing points.
* We have shown how relevant is the determination of the zero-frequency constants in order to correctly compute the bosonic Green’s functions. In particular, we have shown that, for the 2-site system, they assume values very far from the ergodic ones.

###### Acknowledgements.

We wish to thank A. MoreoDuffy and Moreo (1997) and P. PrelovsekPrelovsek ; Bonca and Prelovsek (2003); Jaklic and Prelovsek (2000) for providing us with the numerical data.

## Appendix A The eigenproblem

The eigenproblem of a given *grand canonical* Hamiltonian $`H`$ is solved once the latter is diagonalized on the Fock space of the system under study, that is

$`H|n`$ $`=E_n|n`$ (194a)
$`N|n`$ $`=N_n|n`$ (194b)
where $`N=_in(i)`$ is the total number operator and $`|n`$ is a complete orthonormal basis.

Once the eigenproblem is solved, we can compute the thermal average of an operator $`\mathrm{\Phi }`$ by means of the following expression

$$\mathrm{\Phi }=\frac{1}{Z}\underset{n}{}n\left|\mathrm{\Phi }\right|n\mathrm{e}^{\beta E_n}$$
(195)
where $`Z=_n\mathrm{e}^{\beta E_n}`$ is the grand canonical partition function and $`\beta `$ is the inverse temperature.

If a finite minimal energy $`E_{\mathrm{min}}`$ exists then

$$\underset{T0}{lim}\mathrm{\Phi }=\frac{1}{M}\underset{n|E_n=E_{\mathrm{min}}}{}n\left|\mathrm{\Phi }\right|n$$
(196)
where $`M`$ is the number of eigenstates $`|n`$ with $`E_n=E_{\mathrm{min}}`$.

On the basis of the knowledge of the set of eigenstates and eigenvalues of $`H`$ (194a) it is possible to derive the expressions for Green’s functions and correlation functions. Let $`\psi (i)`$ be a field operator in the Heisenberg scheme $`\psi (i)=\psi (𝐢,t)=\mathrm{e}^{\mathrm{i}Ht}\psi (𝐢)\mathrm{e}^{\mathrm{i}Ht}`$; we do not specify the nature, fermionic or bosonic, of $`\psi (i)`$ ($`\psi (i)`$ can be, for instance, either $`c(i)`$ or $`n(i)`$). By considering the two-time thermodynamic Green’s functions Bogolubov:59; Zubarev (1960, 1974), let us define the causal function

$$\begin{array}{c}G_{\psi \psi ^{}}^{C(\eta )}(i,j)=\theta \left(t_it_j\right)\psi (i)\psi ^{}(j)\hfill \\ \hfill \eta \theta \left(t_jt_i\right)\psi ^{}(j)\psi (i)\end{array}$$
(197)
the retarded and advanced functions

$$G_{\psi \psi ^{}}^{R,A(\eta )}(i,j)=\pm \theta \left[\pm \left(t_it_j\right)\right][\psi (i),\psi ^{}(j)]_\eta $$
(198)
and the correlation function

$$C_{\psi \psi ^{}}(i,j)=\psi (i)\psi ^{}(j)$$
(199)
Here $`\eta =\pm 1`$; usually, it is convenient to take $`\eta =1`$ ($`\eta =1`$) for a fermionic (bosonic) field $`\psi `$ (i.e., for a composite field constituted of an odd (even) number of original fields) in order to exploit the canonical (anti)commutation relations of the constituting original fields; but, in principle, both choices are possible. Accordingly, we define

$$[A,B]_\eta =\{\begin{array}{cc}\{A,B\}=AB+BA\hfill & \text{for }\eta =1\hfill \\ [A,B]=ABBA\hfill & \text{for }\eta =1\hfill \end{array}$$
(200)
By using the definition of thermal average (195) it is possible to derive the following expressions for the Green’s functions (197, 198) and correlation functions (199) in terms of eigenstates and eigenvalues

$$\begin{array}{c}G_{\psi \psi ^{}}^{C(\eta )}(𝐢,𝐣,\omega )=\mathrm{\Gamma }_{\psi \psi ^{}}(𝐢,𝐣)\left[(1+\eta )𝒫\frac{1}{\omega }\mathrm{i}\pi (1\eta )\delta (\omega )\right]\hfill \\ \hfill +\frac{1}{Z}\underset{\begin{array}{c}n,m\\ E_nE_m\end{array}}{}\left(\frac{A_{\psi \psi ^{}}^{n,m}(𝐢,𝐣)\mathrm{e}^{\beta E_n}}{\omega +E_nE_m+\mathrm{i}\delta }+\frac{\eta A_{\psi \psi ^{}}^{n,m}(𝐢,𝐣)\mathrm{e}^{\beta E_m}}{\omega +E_nE_m\mathrm{i}\delta }\right)\end{array}$$
(201)
$$\begin{array}{c}G_{\psi \psi ^{}}^{R,A(\eta )}(𝐢,𝐣,\omega )=\mathrm{\Gamma }_{\psi \psi ^{}}(𝐢,𝐣)\frac{1+\eta }{\omega \pm \mathrm{i}\delta }\hfill \\ \hfill +\frac{1}{Z}\underset{\begin{array}{c}n,m\\ E_nE_m\end{array}}{}\frac{A_{\psi \psi ^{}}^{n,m}(𝐢,𝐣)\left(\mathrm{e}^{\beta E_n}+\eta \mathrm{e}^{\beta E_m}\right)}{\omega +E_nE_m\pm \mathrm{i}\delta }\end{array}$$
(202)
$$\begin{array}{c}C_{\psi \psi ^{}}(𝐢,𝐣,\omega )=2\pi \mathrm{\Gamma }_{\psi \psi ^{}}(𝐢,𝐣)\delta (\omega )\hfill \\ \hfill +\frac{2\pi }{Z}\underset{\begin{array}{c}n,m\\ E_nE_m\end{array}}{}\mathrm{e}^{\beta E_n}A_{\psi \psi ^{}}^{n,m}(𝐢,𝐣)\delta (\omega +E_nE_m)\end{array}$$
(203)
where

$$A_{\psi \psi ^{}}^{n,m}(𝐢,𝐣)=n|\psi (𝐢)|mm|\psi ^{}(𝐣)|n$$
(204)
and the zero frequency constant $`\mathrm{\Gamma }_{\psi \psi ^{}}(𝐢,𝐣)`$ has the following representation

$$\mathrm{\Gamma }_{\psi \psi ^{}}(𝐢,𝐣)=\frac{1}{Z}\underset{\begin{array}{c}n,m\\ E_n=E_m\end{array}}{}\mathrm{e}^{\beta E_n}A_{\psi \psi ^{}}^{n,m}(𝐢,𝐣)$$
(205)
It is worth noticing that, in this presentation (i.e., in terms of eigenstates and eigenvalues), the Green’s functions and correlation functions are fully determined up to the value of the chemical potential, which is present in the expressions of the eigenvalues. As usual, the computation of the chemical potential requires the inversion of the expression, in terms of eigenstates and eigenvalues, for the number of particle per site (A.3.1). In the main text instead, we presented expressions for the same quantities (Green’s and correlation functions) in terms of eigenenergies of eigenoperators and of correlation functions of these latter. In this case, more parameters appeared ($`\mathrm{\Delta }`$ and $`p`$ in the fermionic sector and $`\chi _s^\alpha `$, $`d`$, $`\gamma `$ and the zero frequency constant $`\mathrm{\Gamma }_{11\mu }`$ in the bosonic sector) that have been computed self-consistently as we usually do for the chemical potential. The two presentations, although equivalent (i.e., they obviously lead to the same results), require two different self-consistent procedures to come to the computation of the physical quantities. The reason of this occurrence resides in the level of knowledge of the representation we have in the two cases. In the first case, we have the full knowledge of the states of the system and we should only fix their occupancy with respect to the average number of particle per site we wish to fix; in the second case, as we have no direct knowledge of the states before fixing the counting we have to reduce the Hilbert space to the correct one (i.e., the one of the system under analysis), that is, we have to impose constraints in order to select only those states that enjoy the correct *symmetry* properties (e.g., we have to discard states with site double occupied by electrons with parallel spins). At any rate, we want to emphasize once more that all the results obtained in this paper by means of the Green’s function formalism exactly coincide, as it should be, with those computable by means of the thermal averages.

The eigenstates of the systems under study are given by linear combinations of the vectors spanning their Fock spaces. These vectors can be displayed as $`|a,b`$ where $`a`$ and $`b`$ denote the occupancy of each site, respectively. In particular, we have $`0`$ for an empty site, $``$ ($``$) for a single occupied site by a spin-up (spin-down) electron and $``$ for a double occupied site. This latter is not allowed in the $`t`$-$`J`$ model (see Tab. 1).

### A.1 The Hubbard model

The eigenstates and eigenenergies of the 2-site Hubbard model are reported in Tab. 1. The coefficients $`\alpha _1`$ and $`\alpha _2`$ are determined by the orthonormality of the eigenstates and have the following expressions:

$`\alpha _1`$ $`={\displaystyle \frac{\left(U+4J_U\right)\sqrt{2}}{2\sqrt{\left(U+4J_U\right)^2+16t^2}}}`$ (206)
$`\alpha _2`$ $`={\displaystyle \frac{2t\sqrt{2}}{\sqrt{\left(U+4J_U\right)^2+16t^2}}}`$ (207)
We also have $`\alpha _1^2+\alpha _2^2=1/2`$.

Expressions (201), (202) and (203) show that the Green’s functions and the correlation functions are completely determined (up to the value of the chemical potential) once the matrices $`A_{\psi \psi ^{}}^{n,m}(𝐢,𝐣)`$ are known. We here give the results for some relevant operators ($`𝐢`$, $`𝐣=\mathrm{𝟏}`$, $`\mathrm{𝟐}`$; $`\mathrm{\Delta }_{a,b}=\delta _{a,n}\delta _{b,m}`$):

$`A_{\xi _{}\xi _{}^{}}^{n,m}(\mathrm{𝟏},\mathrm{𝟏})=B_{\xi _{}\xi _{}^{}}^{n,m}+C_{\xi _{}\xi _{}^{}}^{n,m}`$ (208)
$`A_{\xi _{}\xi _{}^{}}^{n,m}(\mathrm{𝟏},\mathrm{𝟐})=B_{\xi _{}\xi _{}^{}}^{n,m}+C_{\xi _{}\xi _{}^{}}^{n,m}`$ (209)
$`A_{\xi _{}\xi _{}^{}}^{n,m}(\mathrm{𝟏},\mathrm{𝟏})=B_{\xi _{}\xi _{}^{}}^{n,m}+C_{\xi _{}\xi _{}^{}}^{n,m}`$ (210)
$`A_{\xi _{}\xi _{}^{}}^{n,m}(\mathrm{𝟏},\mathrm{𝟐})=B_{\xi _{}\xi _{}^{}}^{n,m}+C_{\xi _{}\xi _{}^{}}^{n,m}`$ (211)
where

$$\begin{array}{c}B_{\xi _{}\xi _{}^{}}^{n,m}=\frac{1}{2}\left(\mathrm{\Delta }_{1,2}+\mathrm{\Delta }_{4,6}+\frac{1}{2}\mathrm{\Delta }_{5,8}+\alpha _1^2\mathrm{\Delta }_{3,9}\right)\hfill \\ \hfill +\frac{1}{2}\left(\alpha _2^2\mathrm{\Delta }_{3,11}+\alpha _2^2\mathrm{\Delta }_{9,14}+\frac{1}{2}\mathrm{\Delta }_{10,12}+\alpha _1^2\mathrm{\Delta }_{11,14}\right)\end{array}$$
(212)
$$\begin{array}{c}C_{\xi _{}\xi _{}^{}}^{n,m}=\frac{1}{2}\left(\mathrm{\Delta }_{1,4}+\mathrm{\Delta }_{2,6}+\frac{1}{2}\mathrm{\Delta }_{3,8}+\alpha _1^2\mathrm{\Delta }_{5,9}\right)\hfill \\ \hfill +\frac{1}{2}\left(\alpha _2^2\mathrm{\Delta }_{5,11}+\alpha _2^2\mathrm{\Delta }_{9,12}+\frac{1}{2}\mathrm{\Delta }_{10,14}+\alpha _1^2\mathrm{\Delta }_{11,12}\right)\end{array}$$
(213)
$$\begin{array}{c}B_{\xi _{}\xi _{}^{}}^{n,m}=\frac{1}{2}\left(\mathrm{\Delta }_{1,3}+\mathrm{\Delta }_{5,7}+\frac{1}{2}\mathrm{\Delta }_{4,8}+\alpha _1^2\mathrm{\Delta }_{2,9}\right)\hfill \\ \hfill +\frac{1}{2}\left(\alpha _2^2\mathrm{\Delta }_{2,11}+\alpha _2^2\mathrm{\Delta }_{9,15}+\frac{1}{2}\mathrm{\Delta }_{10,13}+\alpha _1^2\mathrm{\Delta }_{11,15}\right)\end{array}$$
(214)
$$\begin{array}{c}C_{\xi _{}\xi _{}^{}}^{n,m}=\frac{1}{2}\left(\mathrm{\Delta }_{1,5}+\mathrm{\Delta }_{3,7}+\frac{1}{2}\mathrm{\Delta }_{2,8}+\alpha _1^2\mathrm{\Delta }_{4,9}\right)\hfill \\ \hfill +\frac{1}{2}\left(\alpha _2^2\mathrm{\Delta }_{4,11}+\alpha _2^2\mathrm{\Delta }_{9,13}+\frac{1}{2}\mathrm{\Delta }_{10,15}+\alpha _1^2\mathrm{\Delta }_{11,13}\right)\end{array}$$
(215)
$`A_{\eta _{}\eta _{}^{}}^{n,m}(\mathrm{𝟏},\mathrm{𝟏})=B_{\eta _{}\eta _{}^{}}^{n,m}+C_{\eta _{}\eta _{}^{}}^{n,m}`$ (216)
$`A_{\eta _{}\eta _{}^{}}^{n,m}(\mathrm{𝟏},\mathrm{𝟐})=B_{\eta _{}\eta _{}^{}}^{n,m}+C_{\eta _{}\eta _{}^{}}^{n,m}`$ (217)
$`A_{\eta _{}\eta _{}^{}}^{n,m}(\mathrm{𝟏},\mathrm{𝟏})=B_{\eta _{}\eta _{}^{}}^{n,m}+C_{\eta _{}\eta _{}^{}}^{n,m}`$ (218)
$`A_{\eta _{}\eta _{}^{}}^{n,m}(\mathrm{𝟏},\mathrm{𝟐})=B_{\eta _{}\eta _{}^{}}^{n,m}+C_{\eta _{}\eta _{}^{}}^{n,m}`$ (219)
where

$$\begin{array}{c}B_{\eta _{}\eta _{}^{}}^{n,m}=\frac{1}{2}\left(\mathrm{\Delta }_{7,15}+\mathrm{\Delta }_{15,16}+\frac{1}{2}\mathrm{\Delta }_{5,10}+\alpha _1^2\mathrm{\Delta }_{3,11}\right)\hfill \\ \hfill +\frac{1}{2}\left(\alpha _2^2\mathrm{\Delta }_{3,9}+\alpha _2^2\mathrm{\Delta }_{11,12}+\frac{1}{2}\mathrm{\Delta }_{8,14}+\alpha _1^2\mathrm{\Delta }_{9,12}\right)\end{array}$$
(220)
$$\begin{array}{c}C_{\eta _{}\eta _{}^{}}^{n,m}=\frac{1}{2}\left(\mathrm{\Delta }_{7,13}+\mathrm{\Delta }_{13,16}+\frac{1}{2}\mathrm{\Delta }_{3,10}+\alpha _1^2\mathrm{\Delta }_{5,11}\right)\hfill \\ \hfill +\frac{1}{2}\left(\alpha _2^2\mathrm{\Delta }_{5,9}+\alpha _2^2\mathrm{\Delta }_{11,14}+\frac{1}{2}\mathrm{\Delta }_{8,12}+\alpha _1^2\mathrm{\Delta }_{9,14}\right)\end{array}$$
(221)
$$\begin{array}{c}B_{\eta _{}\eta _{}^{}}^{n,m}=\frac{1}{2}\left(\mathrm{\Delta }_{6,14}+\mathrm{\Delta }_{14,16}+\frac{1}{2}\mathrm{\Delta }_{4,10}+\alpha _1^2\mathrm{\Delta }_{2,11}\right)\hfill \\ \hfill +\frac{1}{2}\left(\alpha _2^2\mathrm{\Delta }_{2,9}+\alpha _2^2\mathrm{\Delta }_{11,13}+\frac{1}{2}\mathrm{\Delta }_{8,15}+\alpha _1^2\mathrm{\Delta }_{9,13}\right)\end{array}$$
(222)
$$\begin{array}{c}C_{\eta _{}\eta _{}^{}}^{n,m}=\frac{1}{2}\left(\mathrm{\Delta }_{6,12}+\mathrm{\Delta }_{12,16}+\frac{1}{2}\mathrm{\Delta }_{2,10}+\alpha _1^2\mathrm{\Delta }_{4,11}\right)\hfill \\ \hfill +\frac{1}{2}\left(\alpha _2^2\mathrm{\Delta }_{4,9}+\alpha _2^2\mathrm{\Delta }_{11,15}+\frac{1}{2}\mathrm{\Delta }_{8,13}+\alpha _1^2\mathrm{\Delta }_{9,15}\right)\end{array}$$
(223)
$`A_{\mathrm{nn}}^{n,m}(\mathrm{𝟏},\mathrm{𝟏})=B_{\mathrm{nn}}^{n,m}+C_{\mathrm{nn}}^{n,m}+D_{\mathrm{nn}}^{n,m}`$ (224)
$`A_{\mathrm{nn}}^{n,m}(\mathrm{𝟏},\mathrm{𝟐})=B_{\mathrm{nn}}^{n,m}C_{\mathrm{nn}}^{n,m}D_{\mathrm{nn}}^{n,m}`$ (225)
where

$$\begin{array}{c}B_{\mathrm{nn}}^{n,m}=\frac{1}{4}\left(\mathrm{\Delta }_{2,2}+\mathrm{\Delta }_{3,3}+\mathrm{\Delta }_{4,4}+\mathrm{\Delta }_{5,5}+4\mathrm{\Delta }_{6,6}\right)\hfill \\ \hfill +\left(\mathrm{\Delta }_{7,7}+\mathrm{\Delta }_{8,8}+\mathrm{\Delta }_{9,9}+\mathrm{\Delta }_{10,10}+\mathrm{\Delta }_{11,11}\right)\\ \hfill +\frac{1}{4}\left(9\mathrm{\Delta }_{12,12}+9\mathrm{\Delta }_{13,13}+9\mathrm{\Delta }_{14,14}+9\mathrm{\Delta }_{15,15}+16\mathrm{\Delta }_{16,16}\right)\end{array}$$
(226)
$$\begin{array}{c}C_{\mathrm{nn}}^{n,m}=\frac{1}{4}\left(\mathrm{\Delta }_{2,4}+\mathrm{\Delta }_{3,5}+\mathrm{\Delta }_{4,2}+\mathrm{\Delta }_{5,3}\right)\hfill \\ \hfill +\frac{1}{4}\left(\mathrm{\Delta }_{12,14}+\mathrm{\Delta }_{13,15}+\mathrm{\Delta }_{14,12}+\mathrm{\Delta }_{15,13}\right)\end{array}$$
(227)
$$\begin{array}{c}D_{\mathrm{nn}}^{n,m}=2\left(\alpha _2^2\mathrm{\Delta }_{9,10}+\alpha _2^2\mathrm{\Delta }_{10,9}+\alpha _1^2\mathrm{\Delta }_{10,11}+\alpha _1^2\mathrm{\Delta }_{11,10}\right)\hfill \end{array}$$
(228)
$`A_{n_3n_3}^{n,m}(\mathrm{𝟏},\mathrm{𝟏})=B_{n_3n_3}^{n,m}+C_{n_3n_3}^{n,m}+D_{n_3n_3}^{n,m}`$ (229)
$`A_{n_3n_3}^{n,m}(\mathrm{𝟏},\mathrm{𝟐})=B_{n_3n_3}^{n,m}C_{n_3n_3}^{n,m}D_{n_3n_3}^{n,m}`$ (230)
where

$$\begin{array}{c}B_{n_3n_3}^{n,m}=\frac{1}{4}\left(\mathrm{\Delta }_{2,2}+\mathrm{\Delta }_{3,3}+\mathrm{\Delta }_{4,4}+\mathrm{\Delta }_{5,5}+4\mathrm{\Delta }_{6,6}\right)\hfill \\ \hfill +\frac{1}{4}\left(4\mathrm{\Delta }_{7,7}+\mathrm{\Delta }_{12,12}+\mathrm{\Delta }_{13,13}+\mathrm{\Delta }_{14,14}+\mathrm{\Delta }_{15,15}\right)\end{array}$$
(231)
$$C_{n_3n_3}^{n,m}=C_{nn}^{n,m}$$
(232)
$$\begin{array}{c}D_{n_3n_3}^{n,m}=2\left(\alpha _2^2\mathrm{\Delta }_{8,11}+\alpha _2^2\mathrm{\Delta }_{11,8}+\alpha _1^2\mathrm{\Delta }_{8,9}+\alpha _1^2\mathrm{\Delta }_{9,8}\right)\hfill \end{array}$$
(233)
### A.2 The $`t`$-$`J`$ model

The eigenstates and eigenenergies of the 2-site $`t`$-$`J`$ model coincides with the first nine of the Hubbard model with $`\alpha _1=1/\sqrt{2}`$, $`\alpha _2=0`$ and $`J_U=J`$. Actually, in the strong coupling regime (i.e., $`Ut`$) we have $`\alpha _11/\sqrt{2}`$, $`\alpha _20`$ and $`J_U4t^2/U`$. In particular, this latter value ($`4t^2/U`$) is the one we get in the strong coupling regime when deriving the $`t`$-$`J`$ model from the Hubbard one. All other states of the Hubbard model can not be realized in the $`t`$-$`J`$ model owing to the exclusion of the double occupancies.

Following the same prescription (only first nine states, $`\alpha _1=1/\sqrt{2}`$, $`\alpha _2=0`$ and $`J_U=J`$), it is possible to get $`A_{\xi \xi ^{}}^{n,m}(𝐢,𝐣)`$, $`A_{\nu \nu }^{n,m}(𝐢,𝐣)`$ and $`A_{\nu _3\nu _3}^{n,m}(𝐢,𝐣)`$ form the corresponding expressions given for the Hubbard model ($`A_{\xi \xi ^{}}^{n,m}(𝐢,𝐣)`$, $`A_{nn}^{n,m}(𝐢,𝐣)`$ and $`A_{n_3n_3}^{n,m}(𝐢,𝐣)`$, respectively).

### A.3 The properties

Given the expressions of $`A_{\psi \psi ^{}}^{n,m}(𝐢,𝐣)`$, we can now compute the relevant physical quantities.

#### A.3.1 The Hubbard model

The partition function is given by

$`Z`$ $`=1+2\mathrm{e}^{\beta E_2}+2\mathrm{e}^{\beta E_4}+3\mathrm{e}^{\beta E_6}+\mathrm{e}^{\beta E_9}+\mathrm{e}^{\beta E_{10}}`$
$`+\mathrm{e}^{\beta E_{11}}+2\mathrm{e}^{\beta E_{12}}+2\mathrm{e}^{\beta E_{14}}+\mathrm{e}^{\beta E_{16}}`$ (234)
The chemical potential $`\mu `$ can be computed inverting the following expression for the particle number per site

$`n`$ $`={\displaystyle \frac{1}{Z}}(\mathrm{e}^{\beta E_2}+\mathrm{e}^{\beta E_4}+3\mathrm{e}^{\beta E_6}+\mathrm{e}^{\beta E_9}+\mathrm{e}^{\beta E_{10}}`$
$`+\mathrm{e}^{\beta E_{11}}+3\mathrm{e}^{\beta E_{12}}+3\mathrm{e}^{\beta E_{14}}+2\mathrm{e}^{\beta E_{16}})`$ (235)
This expression determines the chemical potential as a function of $`n`$, $`T`$ and $`U`$.

The self-consistent parameter $`\mathrm{\Delta }`$ and $`p`$ are given by

$$\mathrm{\Delta }=\frac{1}{2Z}\left(\mathrm{e}^{\beta E_2}\mathrm{e}^{\beta E_4}\mathrm{e}^{\beta E_{13}}+\mathrm{e}^{\beta E_{15}}\right)$$
(236)
$`p`$ $`={\displaystyle \frac{1}{2Z}}(3\mathrm{e}^{\beta E_6}\mathrm{e}^{\beta E_9}+\mathrm{e}^{\beta E_{10}}\mathrm{e}^{\beta E_{11}}+2\mathrm{e}^{\beta E_{12}}`$
$`+2\mathrm{e}^{\beta E_{14}}+2\mathrm{e}^{\beta E_{16}})`$ (237)
The correlation functions $`C_{12}^\alpha `$ and $`C_{cc}^\alpha `$ are given by

$$C_{12}^\alpha =\frac{\alpha _1\alpha _2}{Z}\left(\mathrm{e}^{\beta E_9}\mathrm{e}^{\beta E_{11}}\right)$$
(238)
$`C_{cc}^\alpha `$ $`={\displaystyle \frac{1}{2Z}}[\mathrm{e}^{\beta E_2}\mathrm{e}^{\beta E_4}+4\alpha _1\alpha _2(\mathrm{e}^{\beta E_9}\mathrm{e}^{\beta E_{11}})`$
$`+\mathrm{e}^{\beta E_{12}}\mathrm{e}^{\beta E_{14}}]`$ (239)
The double occupancy per site is given by

D $`={\displaystyle \frac{1}{2Z}}(2\alpha _2^2\mathrm{e}^{\beta E_9}+\mathrm{e}^{\beta E_{10}}+2\alpha _1^2\mathrm{e}^{\beta E_{11}}+2\mathrm{e}^{\beta E_{12}}`$
$`+2\mathrm{e}^{\beta E_{14}}+2\mathrm{e}^{\beta E_{16}})`$ (240)
Spin ($`\chi _s^\alpha `$), charge ($`\chi _c^\alpha `$) and pair ($`d`$) correlation functions are given by

$$\chi _s^\alpha =\frac{3}{Z}\left(\mathrm{e}^{\beta E_6}2\alpha _1^2\mathrm{e}^{\beta E_9}2\alpha _2^2\mathrm{e}^{\beta E_{11}}\right)$$
(241)
$`\chi _c^\alpha `$ $`={\displaystyle \frac{1}{Z}}(3\mathrm{e}^{\beta E_6}+2\alpha _1^2\mathrm{e}^{\beta E_9}+2\alpha _2^2\mathrm{e}^{\beta E_{11}}`$
$`+4\mathrm{e}^{\beta E_{12}}+4\mathrm{e}^{\beta E_{14}}+4\mathrm{e}^{\beta E_{16}})`$ (242)
$$d=\frac{1}{2Z}\left(2\alpha _2^2\mathrm{e}^{\beta E_9}\mathrm{e}^{\beta E_{10}}+2\alpha _1^2\mathrm{e}^{\beta E_{11}}\right)$$
(243)
The zero-frequency constant $`\mathrm{\Gamma }_{110}=\mathrm{\Gamma }_{\mathrm{nn}}(\mathrm{𝟏},\mathrm{𝟐})`$ and $`\mathrm{\Gamma }_{113}=\mathrm{\Gamma }_{n_3n_3}(\mathrm{𝟏},\mathrm{𝟐})`$ are given by

$`\mathrm{\Gamma }_{110}`$ $`={\displaystyle \frac{1}{2Z}}(\mathrm{e}^{\beta E_2}+\mathrm{e}^{\beta E_4}+6\mathrm{e}^{\beta E_6}+2\mathrm{e}^{\beta E_9}+2\mathrm{e}^{\beta E_{10}}`$
$`+2\mathrm{e}^{\beta E_{11}}+9\mathrm{e}^{\beta E_{12}}+9\mathrm{e}^{\beta E_{14}}+8\mathrm{e}^{\beta E_{16}})`$ (244)
$$\mathrm{\Gamma }_{113}=\frac{1}{2Z}\left(\mathrm{e}^{\beta E_2}+\mathrm{e}^{\beta E_4}+4\mathrm{e}^{\beta E_6}+\mathrm{e}^{\beta E_{12}}+\mathrm{e}^{\beta E_{14}}\right)$$
(245)
### A.4 The $`t`$-$`J`$ model

The partition function is given by

$$Z=1+2\mathrm{e}^{\beta E_2}+2\mathrm{e}^{\beta E_4}+3\mathrm{e}^{\beta E_6}+\mathrm{e}^{\beta E_9}$$
(246)
The chemical potential $`\mu `$ can be computed through the following expression for the particle number per site

$$n=\frac{1}{Z}(\mathrm{e}^{\beta E_2}+\mathrm{e}^{\beta E_4}+3\mathrm{e}^{\beta E_6}+\mathrm{e}^{\beta E_9})$$
(247)
which can be inverted and give

$$\mu =T\mathrm{ln}\frac{(2n1)\mathrm{cosh}(2\beta t)+\sqrt{(2n1)^2\mathrm{cosh}^2(2\beta t)+n(1n)(3+\mathrm{e}^{4\beta J})}}{(1n)(3+\mathrm{e}^{4\beta J})}$$
(248)
The self-consistent parameter $`C_{11}^\alpha `$ and $`\chi ^\alpha `$ are given by

$$C_{11}^\alpha =\frac{1}{2Z}\left(\mathrm{e}^{\beta E_2}\mathrm{e}^{\beta E_4}\right)$$
(249)
$$\chi _\mu ^\alpha =\frac{2}{Z}\left(3\mathrm{e}^{\beta E_6}\mathrm{e}^{\beta E_9}\right)$$
(250)
Spin ($`\chi _s^\alpha `$) and charge ($`\chi _c^\alpha `$) correlation functions are given by

$$\chi _s^\alpha =\frac{1}{Z}\left(\mathrm{e}^{\beta E_6}\mathrm{e}^{\beta E_9}\right)$$
(251)
$$\chi _c^\alpha =\frac{1}{Z}\left(3\mathrm{e}^{\beta E_6}+\mathrm{e}^{\beta E_9}\right)$$
(252)
The zero-frequency constant $`\mathrm{\Gamma }_{110}=\mathrm{\Gamma }_{\mathrm{nn}}(\mathrm{𝟏},\mathrm{𝟐})`$ and $`\mathrm{\Gamma }_{113}=\mathrm{\Gamma }_{n_3n_3}(\mathrm{𝟏},\mathrm{𝟐})`$ are given by

$`\mathrm{\Gamma }_{110}`$ $`={\displaystyle \frac{1}{2Z}}\left(\mathrm{e}^{\beta E_2}+\mathrm{e}^{\beta E_4}+6\mathrm{e}^{\beta E_6}+2\mathrm{e}^{\beta E_9}\right)`$ (253)
$$\mathrm{\Gamma }_{113}=\frac{1}{2Z}\left(\mathrm{e}^{\beta E_2}+\mathrm{e}^{\beta E_4}+4\mathrm{e}^{\beta E_6}\right)$$
(254)