# Density-matrix functional theory of the Hubbard model: An exact numerical study

## I Introduction

Density functional theory (DFT) has been the subject of remarkable developments since its original formulation by Hohenberg and Kohn (HK) . After formal improvements, extensions, and an uncountable number of applications to a wide variety of physical problems, this theoretical approach has become the most efficient, albeit not infallible, method of determining the electronic properties of matter from first principles . The most important innovation of DFT, which is actually at the origin of its breakthrough, is to replace the wave function by the electronic density $`\rho (\stackrel{}{r})`$ as the fundamental variable of the many-body problem. In practice, density functional (DF) calculations are largely based on the Kohn-Sham (KS) scheme that reduces the many-body $`N`$-particle problem to the solution of a set of self-consistent single-particle equations . Although this transformation is formally exact, the implementations always require approximations, since the KS equations involve functional derivatives of the unknown interaction energy $`W[\rho (\stackrel{}{r})]`$, usually expressed in terms of the exchange and correlation (XC) energy $`E_{\mathrm{XC}}[\rho (\stackrel{}{r})]`$. Therefore, understanding the functional dependence of $`E_{\mathrm{XC}}[\rho (\stackrel{}{r})]`$ and improving its approximations are central to the development of DF methods. The currently most widespread Ansätze for $`E_{\mathrm{XC}}[\rho (\stackrel{}{r})]`$ —the local density approximation (LDA)  with spin polarized  and gradient corrected extensions — were originally derived from exact results for the homogeneous electron gas. It is one of the purposes of this paper to investigate the properties of the interaction-energy functional from an intrinsically inhomogeneous point of view, namely, by considering exactly solvable many-body lattice models.

Despite the remarkable success of the local spin density approximation, present DFT fails systematically in accounting for phenomena where strong electron correlations play a central role, for example, in heavy-fermion materials or high-$`T_c`$ superconductors. These systems are usually described by simplifying the low-energy electron dynamics using parameterized lattice models such as Pariser-Parr-Pople, Hubbard, or Anderson models and related Hamiltonians . Being in principle an exact theory, the limitations of the DF approach have to be ascribed to the approximations used for exchange and correlation and not to the underlying HKS formalism. It would be therefore very interesting to extend the range of applicability of DFT to strongly correlated systems and to characterize the properties $`E_{\mathrm{XC}}`$ in the limit of strong correlations. Studies of the XC functional on simple models should provide useful insights for future extensions to realistic Hamiltonians. Moreover, taking into account the demonstrated power of the DF approach in ab initio calculations, one may also expect that a DFT with an appropriate $`E_{\mathrm{XC}}`$ could become an efficient tool for studying many-body models, a subject of theoretical interest on its own.

Several properties of DFT on lattice models have been already studied in previous works. Gunnarsson and Schönhammer were, to our knowledge, the first to propose a DF approach on a semiconductor model in order to study the band-gap problem. In this case the local site occupancies were treated as the basic variables. Some years later Schindlmayr and Godby  provided a different formulation of DFT on a lattice by considering as basic variables both diagonal elements $`\gamma _{ii}`$ and off-diagonal elements $`\gamma _{ij}`$ of the single-particle density matrix (see also ). Schönhammer et al. then derived a more general framework that unifies the two previous approaches . Using Levy’s constrained search method  they showed that different basic variables and different $`W`$ functionals can be considered depending on the type of model or perturbation under study. Site occupations alone may be used as basic variables, if only the orbital energies are varied (i.e., if all hopping integrals $`t_{ij}`$ are kept constant for $`ij`$). However, off-diagonal elements of the single-particle density matrix must be included explicitly if the functional $`W`$ is intended to be applied to more general situations involving different values of $`t_{ij}`$, for example, the Hubbard model on various lattice structures or for different interaction regimes, i.e., different $`U/t`$.

In this paper we investigate the properties of Levy’s interaction-energy functional $`W`$ as a function of $`\gamma _{ij}`$ by solving the constrained search minimization problem exactly. In Sec. II the basic formalism of density-matrix functional theory (DMFT) on lattice models is recalled and the equations for determining $`W[\gamma _{ij}]`$ are derived. Sec. III presents and discusses exact results for the correlation energy $`E_\mathrm{C}`$ of the Hubbard model, which is given by the difference between $`W`$ and the Hartree-Fock energy $`E_{\mathrm{HF}}`$. These are obtained, either numerically for finite clusters with different lattice structures, or from the Bethe-Ansatz solution for the one-dimensional chain. Finally, Sec. IV summarizes our conclusions and points out some relevant extensions.

## II Theory

In Sec. II A the main results of Levy’s formulation of DMFT are presented in a form that is appropriate for the study of model Hamiltonians such as the Hubbard model. Here, the hopping integrals $`t_{ij}`$ between sites (or orbitals) $`i`$ and $`j`$ play the role given in conventional DFT to the external potential $`V_{ext}(\stackrel{}{r})`$. Consequently, the single-particle density matrix $`\gamma _{ij}`$ replaces the density $`\rho (\stackrel{}{r})`$ as basic variable . In Sec. II B, we derive equations that allow to determine Levy’s interaction-energy functional $`W[\gamma _{ij}]`$ in terms of the ground-state energy of a many-body Hamiltonian with effective hopping integrals $`\lambda _{ij}`$ that depend implicitly on $`\gamma _{ij}`$.

### A DMFT of lattice models

We consider the many-body Hamiltonian

$$H=\underset{ij\sigma }{}t_{ij}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }+\frac{1}{2}\underset{\genfrac{}{}{0pt}{}{ijkl}{\sigma \sigma ^{}}}{}V_{ijkl}\widehat{c}_{i\sigma }^{}\widehat{c}_{k\sigma ^{}}^{}\widehat{c}_{l\sigma ^{}}\widehat{c}_{j\sigma },$$
(1)   
where $`\widehat{c}_{i\sigma }^{}`$ ($`\widehat{c}_{i\sigma }`$) is the usual creation (annihilation) operator for an electron with spin $`\sigma `$ at site (or orbital) $`i`$. $`H`$ can be regarded as the second quantization of Schrödinger’s equation on a basis . However, in the present paper, the hopping integrals $`t_{ij}`$ and the interaction matrix elements $`V_{ijkl}`$ are taken as parameters to be varied independently. The matrix $`t_{ij}`$ defines the lattice (e.g., one dimensional chains, square or triangular two-dimensional lattices) and the range of single-particle interactions (e.g., up to first or second neighbors). From the ab initio perspective $`t_{ij}`$ is given by the external potential and by the choice of the basis . $`V_{ijkl}`$ defines the type of many-body interactions which may be repulsive (Coulomb like) or attractive (in order to simulate electronic pairing) and which are usually approximated as short ranged (e.g., intra-atomic). Eq. (1) is mainly used in this section to derive general results which can then be applied to various specific models by simplifying the interactions. A particularly relevant example, to be considered in some detail in Sec. III, is the single-band Hubbard model with nearest neighbor (NN) hoppings , which can be obtained from Eq. (1) by setting $`t_{ij}=t`$ for $`i`$ and $`j`$ NN’s, $`t_{ij}=0`$ otherwise, and $`V_{ijkl}=U\delta _{ij}\delta _{kl}\delta _{ik}`$ .

In order to apply DMFT to model Hamiltonians of the form (1) we follow Levy’s constrained search procedure  as proposed by Schindlmayr and Godby . The ground-state energy is determined by minimizing the functional

$$E[\gamma _{ij}]=E_K[\gamma _{ij}]+W[\gamma _{ij}]$$
(2)   
with respect to the single-particle density matrix $`\gamma _{ij}`$. $`E[\gamma _{ij}]`$ is physically defined for all density matrices that can be written as

$$\gamma _{ij}=\underset{\sigma }{}\gamma _{ij\sigma }=\underset{\sigma }{}\mathrm{\Psi }|c_{i\sigma }^{}c_{j\sigma }|\mathrm{\Psi }$$
(3)   
for all $`i`$ and $`j`$, where $`|\mathrm{\Psi }`$ is an $`N`$-particle state. In other words, $`\gamma _{ij}`$ must derive from a physical state. It is then said to be pure-state $`N`$-representable . The first term in Eq. (2) is given by

$$E_K=\underset{ij}{}t_{ij}\gamma _{ij}.$$
(4)   
It includes all single-particle contributions and is usually regarded as the kinetic energy associated with the electronic motion in the lattice. Notice that Eq. (4) yields the exact kinetic energy for a given $`\gamma _{ij}`$. There are no corrections on $`E_K`$ to be included in other parts of the functional as in the KS approach. The second term in Eq. (2) is the interaction-energy functional given by 

$$W[\gamma _{ij}]=min\left[\frac{1}{2}\underset{\genfrac{}{}{0pt}{}{nmkl}{\sigma \sigma ^{}}}{}V_{nmkl}\mathrm{\Psi }[\gamma _{ij}]|\widehat{c}_{n\sigma }^{}\widehat{c}_{k\sigma ^{}}^{}\widehat{c}_{l\sigma ^{}}\widehat{c}_{m\sigma }|\mathrm{\Psi }[\gamma _{ij}]\right].$$
(5)   
The minimization in Eq. (5) implies a search over all $`N`$-particles states $`|\mathrm{\Psi }[\gamma _{ij}]`$ that satisfy $`\mathrm{\Psi }[\gamma _{ij}]|_\sigma \widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }|\mathrm{\Psi }[\gamma _{ij}]=\gamma _{ij}`$ for all $`i`$ and $`j`$. Therefore, $`W[\gamma _{ij}]`$ represents the minimum value of the interaction energy compatible with a given density matrix $`\gamma _{ij}`$. $`W`$ is usually expressed in terms of the Hartree-Fock energy

$$E_{\mathrm{HF}}[\gamma _{ij}]=\frac{1}{2}\underset{\genfrac{}{}{0pt}{}{ijkl}{\sigma \sigma ^{}}}{}V_{ijkl}\left(\gamma _{ij\sigma }\gamma _{kl\sigma ^{}}\delta _{\sigma \sigma ^{}}\gamma _{il\sigma }\gamma _{kj\sigma }\right)$$
(6)   
and the correlation energy $`E_\mathrm{C}[\gamma _{ij}]`$ as

$$W[\gamma _{ij}]=E_{\mathrm{HF}}[\gamma _{ij}]+E_\mathrm{C}[\gamma _{ij}].$$
(7)   
$`W`$ and $`E_\mathrm{C}`$ are universal functionals of $`\gamma _{ij}`$ in the sense that they are independent of $`t_{ij}`$, i.e., of the system under study. They depend on the considered interactions or model, as defined by $`V_{ijkl}`$, on the number of electrons $`N_e`$, and on the structure of the many-body Hilbert space, as given by $`N_e`$ and the number of orbitals or sites $`N_a`$. Notice that $`E_\mathrm{C}`$ in Eq. (7) does not include any exchange contributions. Given $`\gamma _{ij}`$ ($`\gamma _{ij\sigma }=\gamma _{ij}/2`$ in nonmagnetic cases) there is no need to approximate the exchange term, which is taken into account exactly by $`E_{\mathrm{HF}}`$ \[Eq. (6)\]. Nevertheless, if useful in practice, it is of course possible to split $`W`$ in the Hartree energy $`E_\mathrm{H}`$ and the exchange and correlation energy $`E_{\mathrm{XC}}`$ is a similar way as in the KS approach.

The variational principle results from the following two relations :

$$E_{gs}\underset{ij}{}t_{ij}\gamma _{ij}+W[\gamma _{ij}]$$
(8)   
for all pure-state $`N`$-representable $`\gamma _{ij}`$ , and

$$E_{gs}=\underset{ij}{}t_{ij}\gamma _{ij}^{gs}+W[\gamma _{ij}^{gs}],$$
(9)   
where $`E_{gs}=\mathrm{\Psi }_{gs}|H|\mathrm{\Psi }_{gs}`$ refers to the ground-state energy and $`\gamma _{ij}^{gs}=\mathrm{\Psi }_{gs}|_\sigma \widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }|\mathrm{\Psi }_{gs}`$ to the ground-state single-particle density matrix.

As already pointed out in previous works , $`W`$ and $`E_\mathrm{C}`$ depend in general on both diagonal elements $`\gamma _{ii}`$ and off-diagonal elements $`\gamma _{ij}`$ of the density-matrix, since the hopping integrals $`t_{ij}`$ are non local in the sites. The situation is similar to the DF approach proposed by Gilbert for the study of non-local potentials $`V_{ext}(\stackrel{}{r},\stackrel{}{r}^{})`$ as those appearing in the theory of pseudo-potentials . A formulation of DFT on a lattice only in terms of $`\gamma _{ii}`$ would be possible if one would restrict oneself to a family of models with constant $`t_{ij}`$ for $`ij`$. However, in this case the functional $`W[\gamma _{ii}]`$ would depend on the actual value of $`t_{ij}`$ for $`ij`$ .

The functional $`W[\gamma _{ij}]`$, valid for all lattice structures and for all types of hybridizations, can be simplified at the expense of universality if the hopping integrals are short ranged. For example, if only NN hoppings are considered, the kinetic energy $`E_K`$ is independent of the density-matrix elements between sites that are not NN’s. Therefore, the constrained search in Eq. (5) may restricted to the $`|\mathrm{\Psi }[\gamma _{ij}]`$ that satisfy $`\mathrm{\Psi }[\gamma _{ij}]|_\sigma \widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }|\mathrm{\Psi }[\gamma _{ij}]=\gamma _{ij}`$ only for $`i=j`$ and for NN $`ij`$. In this way the number of variables in $`W[\gamma _{ij}]`$ is reduced significantly rendering the interpretation of the functional dependence far simpler. While this is a great practical advantage, it also implies that $`W`$ and $`E_\mathrm{C}`$ lose their universal character since the dependence on the NN $`\gamma _{ij}`$ is now different for different lattices. In Sec. III results for one-, two-, and three-dimensional lattices with NN hoppings are compared in order to quantify this effect.

For the applications in Sec. III we shall consider the single-band Hubbard model with NN hoppings, which in the usual notation is given by

$$H=t\underset{i,j\sigma }{}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }+U\underset{i}{}\widehat{n}_i\widehat{n}_i.$$
(10)   
In this case the interaction energy functional reads

$$W[\gamma _{ij}]=min\left[U\underset{l}{}\mathrm{\Psi }[\gamma _{ij}]|\widehat{n}_l\widehat{n}_l|\mathrm{\Psi }[\gamma _{ij}]\right],$$
(11)   
where the minimization is performed with respect to all $`N`$-particle $`|\mathrm{\Psi }[\gamma _{ij}]`$ satisfying $`\mathrm{\Psi }[\gamma _{ij}]|_\sigma \widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }|\mathrm{\Psi }[\gamma _{ij}]=\gamma _{ij}`$ for $`i`$ and $`j`$ NN’s. If the interactions are repulsive ($`U>0`$) $`W[\gamma _{ij}]`$ represents the minimum average number of double occupations which can be obtained for a given degree of electron delocalization, i.e., for a given value of $`\gamma _{ij}`$. For attractive interactions ($`U<0`$) double occupations are favored and $`W[\gamma _{ij}]`$ corresponds to the maximum of $`_l\widehat{n}_l\widehat{n}_l`$ for a given $`\gamma _{ij}`$.

### B Exact XC energy functional

In order to determine $`E_\mathrm{C}[\gamma _{ij}]`$ and $`W[\gamma _{ij}]`$ we look for the extremes of

$`F`$ $`=`$ $`{\displaystyle \underset{\genfrac{}{}{0pt}{}{ijkl}{\sigma \sigma ^{}}}{}}\left[V_{ijkl}\mathrm{\Psi }|\widehat{c}_{i\sigma }^{}\widehat{c}_{k\sigma ^{}}^{}\widehat{c}_{l\sigma }\widehat{c}_{j\sigma }|\mathrm{\Psi }\right]+\epsilon \left(1\mathrm{\Psi }|\mathrm{\Psi }\right)`$ (12)   
$`+`$ $`{\displaystyle \underset{i,j}{}}\lambda _{ij}\left(\mathrm{\Psi }|{\displaystyle \underset{\sigma }{}}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }|\mathrm{\Psi }\gamma _{ij}\right)`$ (13)   
with respect to $`|\mathrm{\Psi }`$. Lagrange multipliers $`\epsilon `$ and $`\lambda _{ij}`$ have been introduced to enforce the normalization of $`|\mathrm{\Psi }`$ and the conditions on the representability of $`\gamma _{ij}`$. Derivation with respect to $`\mathrm{\Psi }|`$, $`\epsilon `$ and $`\lambda _{ij}`$ yields the eigenvalue equations

$$\underset{ij\sigma }{}\lambda _{ij}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }|\mathrm{\Psi }+\underset{\genfrac{}{}{0pt}{}{ijkl}{\sigma \sigma ^{}}}{}V_{ijkl}\widehat{c}_{i\sigma }^{}\widehat{c}_{k\sigma ^{}}^{}\widehat{c}_{l\sigma }\widehat{c}_{j\sigma }|\mathrm{\Psi }=\epsilon |\mathrm{\Psi },$$
(14)   
and the auxiliary conditions $`\mathrm{\Psi }|\mathrm{\Psi }=1`$ and $`\gamma _{ij}=\mathrm{\Psi }|_\sigma \widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }|\mathrm{\Psi }`$. The Lagrange multipliers $`\lambda _{ij}`$ play the role of hopping integrals to be chosen in order that $`|\mathrm{\Psi }`$ yields the given $`\gamma _{ij}`$. The pure-state representability of $`\gamma _{ij}`$ ensures that there is always a solution . In practice, however, one usually varies $`\lambda _{ij}`$ in order to scan the domain of representability of $`\gamma _{ij}`$. For given $`\lambda _{ij}`$, the eigenstate $`|\mathrm{\Psi }_0`$ corresponding to the lowest eigenvalue of Eq. (14) yields the minimum $`W[\gamma _{ij}]`$ for $`\gamma _{ij}=\mathrm{\Psi }_0|_\sigma \widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }|\mathrm{\Psi }_0`$. Any other $`|\mathrm{\Psi }`$ satisfying $`\gamma _{ij}=\mathrm{\Psi }|_\sigma \widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }|\mathrm{\Psi }`$ would have higher $`\epsilon `$ and thus higher $`W`$. The subset of $`\gamma _{ij}`$ which are representable by a ground-state of Eq. (14) is the physically relevant one, since it necessarily includes the absolute minimum $`\gamma _{ij}^{gs}`$ of $`E[\gamma _{ij}]`$. Nevertheless, it should be noted that pure-state representable $`\gamma _{ij}`$ may be considered that can only be represented by excited states or by linear combinations of eigenstates of Eq. (14). In the later case, $`\lambda _{ij}=0ij`$, and $`|\mathrm{\Psi }_0`$ is an eigenstate of the interaction term with lowest eigenvalue. Examples shall be discussed in Sec. III.

For the Hubbard model Eq. (14) reduces to

$$\underset{\genfrac{}{}{0pt}{}{ij}{\sigma }}{}\lambda _{ij}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma }|\mathrm{\Psi }+U\underset{i}{}\widehat{n}_i\widehat{n}_i|\mathrm{\Psi }=\epsilon |\mathrm{\Psi }.$$
(15)   
This eigenvalue problem can be solved numerically for clusters with different lattice structures and periodic boundary conditions. In this case we expand $`|\mathrm{\Psi }[\gamma _{ij}]`$ in a complete set of basis states $`|\mathrm{\Phi }_m`$ which have definite occupation numbers $`\nu _{i\sigma }^m`$ at all orbitals $`i\sigma `$ ($`\widehat{n}_{i\sigma }|\mathrm{\Phi }_m=\nu _{i\sigma }^m|\mathrm{\Phi }_m`$ with $`\nu _{i\sigma }^m=0`$ or $`1`$). The values of $`\nu _{i\sigma }^m`$ satisfy the usual conservation of the number of electrons $`N_e=N_e+N_e`$ and of the $`z`$ component of the total spin $`S_z=(N_eN_e)/2`$, where $`N_{e\sigma }=_i\nu _{i\sigma }^m`$. For not too large clusters, the lowest energy $`|\mathrm{\Psi }_0[\gamma _{ij}]`$ —the ground state of Eq. (15)— can be determined by sparse-matrix diagonalization procedures, for example, by using Lanczos iterative method . $`|\mathrm{\Psi }_0[\gamma _{ij}]`$ is calculated in the subspace of minimal $`S_z`$ since this ensures that there are no a priori restrictions on the total spin $`S`$. In addition, spin-projector operators may be used to study the dependence of $`E_\mathrm{C}(\gamma _{12})`$ on $`S`$.

For a one-dimensional (1D) chain with NN hoppings $`t_{ij}=t`$, translational symmetry implies equal density-matrix elements $`\gamma _{ij}`$ between NN’s. Therefore, one may set $`\lambda _{ij}=\lambda `$ for all NN $`ij`$, and then Eq. (15) has the same form as the 1D Hubbard model for which Lieb and Wu’s exact solution is available . In this case the lowest eigenvalue $`\epsilon `$ is determined following the work by Shiba . The coupled Bethe-Ansatz equations are solved as a function of $`\lambda `$, band-filling $`n=N_e/N_a`$, and for positive and negative $`U`$, by means of a simple iterative procedure.

## III Results and Discussion

In this section we present and discuss exact results for the correlation energy functional $`E_\mathrm{C}[\gamma _{ij}]`$ of the single-band Hubbard Hamiltonian with nearest neighbor hoppings . Given the lattice structure, $`N_a`$ and $`N_e`$, the model is characterized by the dimensionless parameter $`U/t`$ which measures the competition between kinetic and interaction energies \[see Eq. (10)\]. $`U>0`$ corresponds to the usual intra-atomic repulsive Coulomb interaction, while the attractive case ($`U<0`$) simulates intra-atomic pairing of electrons.

### A Repulsive interaction $`U>0`$

In Fig. 1 the correlation energy $`E_\mathrm{C}`$ of the one-dimensional (1D) Hubbard model is shown for half-band filling ($`N_e=N_a`$) as a function of the density-matrix element or bond order $`\gamma _{12}`$ between NN’s. $`\gamma _{ij}=\gamma _{12}`$ for all NN’s $`i`$ and $`j`$. Results are given for rings of finite length $`N_a`$ as well as for the infinite chain. Several general qualitative features may be identified. First of all we observe that on bipartite lattices  $`E_\mathrm{C}(\gamma _{12})=E_\mathrm{C}(\gamma _{12})`$, since the sign of the NN bond order can be changed without affecting the interaction energy $`W(\gamma _{12})`$ by changing the phase of the local orbitals at one of the sublattices ($`c_{i\sigma }c_{i\sigma }`$ for $`iA`$ and $`c_{j\sigma }`$ unchanged for $`jB`$, where $`A`$ and $`B`$ refer to the sublattices). Let us recall that the domain of definition of $`E_\mathrm{C}(\gamma _{12})`$ is limited by the pure-state representability of $`\gamma _{ij}`$. The upper bound $`\gamma _{12}^{0+}`$ and the lower bound $`\gamma _{12}^0`$ for $`\gamma _{12}`$ ($`\gamma _{12}^{0+}=\gamma _{12}^0=\gamma _{12}^0`$ on bipartite lattices) are the extreme values of the bond order between NN’s on a given lattice and for given $`N_a`$ and $`N_e`$ ($`\gamma _{ij}=\gamma _{12}`$ for all NN $`ij`$). They represent the maximum degree of electron delocalization. $`\gamma _{12}^{0+}`$ and $`\gamma _{12}^0`$ correspond to the extremes of the kinetic energy $`E_K`$ \[$`E_K=_{ij}t_{ij}\gamma _{ij}=(zN_a/2)t\gamma _{12}`$, where $`z`$ is the coordination number\] and thus to the ground state of the Hubbard model for $`U=0`$ \[$`\gamma _{12}^{0+}`$ for $`t>0`$ and $`\gamma _{12}^0`$ for $`t<0`$, see Eq. (10)\]. For $`\gamma _{12}=\gamma _{12}^0`$ the underlying electronic state $`|\mathrm{\Psi }_0`$ is usually a single Slater determinant and therefore $`E_\mathrm{C}(\gamma _{12}^0)=0`$. In other words, the correlation energy vanishes as expected in the fully delocalized limit . As $`|\gamma _{12}|`$ decreases $`E_\mathrm{C}`$ decreases ($`E_\mathrm{C}<0`$) since correlations can reduce the Coulomb energy more and more efficiently as the electrons localize. $`E_\mathrm{C}`$ is minimum in the strongly correlated limit $`\gamma _{12}=\gamma _{12}^{\mathrm{}}`$. For half-band filling this corresponds to a fully localized electronic state ($`\gamma _{12}^{\mathrm{}}=0`$). Here, $`E_\mathrm{C}`$ cancels out the Hartree-Fock energy $`E_{\mathrm{HF}}`$ and the Coulomb energy $`W`$ vanishes ($`E_\mathrm{C}^{\mathrm{}}=E_{\mathrm{HF}}`$) . The ground-state values of $`\gamma _{12}^{gs}`$ and $`E_{gs}`$ for a given $`U/t`$ result from the competition between lowering $`E_\mathrm{C}`$ by decreasing $`\gamma _{12}`$ and lowering $`E_K`$ by increasing it ($`t>0`$). The divergence of $`E_\mathrm{C}/\gamma _{12}`$ for $`\gamma _{12}=\gamma _{12}^0`$ is a necessary condition in order that $`\gamma _{12}^{gs}<\gamma _{12}^0`$ for arbitrary small $`U>0`$. On the other side, for small $`\gamma _{12}`$, we observe that $`(E_\mathrm{C}+E_{\mathrm{HF}})\gamma _{12}^2`$. This implies that for $`U/t1`$, $`\gamma _{12}^{gs}t/U`$ and $`E_{gs}t^2/U`$, a well known result in the Heisenberg limit of the Hubbard model ($`N_e=N_a`$) .

A more quantitative analysis of $`E_\mathrm{C}(\gamma _{12})`$ and in particular the comparison of results for different $`N_a`$ is complicated by the size dependence of $`\gamma _{12}^0`$ and $`E_{\mathrm{HF}}`$. It is therefore useful to measure $`E_\mathrm{C}`$ in units of the Hartree-Fock energy and to bring the domains of representability to a common range by considering $`\epsilon _\mathrm{C}=E_\mathrm{C}/E_{\mathrm{HF}}`$ as a function of $`g_{12}=\gamma _{12}/\gamma _{12}^0`$. Fig. 2 shows that $`\epsilon _\mathrm{C}(g_{12})`$ has approximately the same behavior for all considered $`N_a`$. Finite size effects are small except for the very small sizes. The largest deviations from the common trend are found for $`N_a=N_e=4`$. Here we observe a discontinuous drop of $`\epsilon _\mathrm{C}`$ for $`g_{12}=1`$ ($`g_{12}<1`$) which is due to the degeneracy of the single-particle spectrum. In fact in this case two of the four electrons occupy a doubly degenerate state in the uncorrelated limit and the minimum interaction energy $`W(\gamma _{12})`$ does not correspond to a single-Slater-determinant state even for $`\gamma _{12}=\gamma _{12}^0`$ . As $`N_a`$ increases $`\epsilon _\mathrm{C}(g_{12})`$ approaches the infinite-length limit with alternations around the $`N_a=\mathrm{}`$ curve. The strong similarity between $`\epsilon _\mathrm{C}(g_{12})`$ for small $`N_a`$ and for $`N_a=\mathrm{}`$ is a remarkable result. It suggests that good approximations for $`E_\mathrm{C}(\gamma _{12})`$ in extended systems could be derived from finite cluster calculations.

Fig. 3 shows the band-filling dependence of $`E_\mathrm{C}(\gamma _{12})`$ in a 10-site 1D Hubbard ring. Results are given for $`N_eN_a`$, since for $`N_eN_a`$, $`E_\mathrm{C}(\gamma _{ij},N_e)=E_\mathrm{C}(\gamma _{ij},2N_aN_e)`$ as a result of electron-hole symmetry . Although $`E_\mathrm{C}(\gamma _{12})`$ depends strongly on $`N_e`$, several qualitative properties are shared by all band fillings: (i) As in the half-filled band case, the domain of representability of $`\gamma _{12}`$ is bound by the bond orders in the uncorrelated limits. In fact, $`\gamma _{12}^0\gamma _{12}\gamma _{12}^{0+}`$, where $`\gamma _{12}^{0+}`$ ($`\gamma _{12}^0`$) corresponds to the ground state of the $`U=0`$ tight-binding model for $`t>0`$ ($`t<0`$). On bipartite lattices $`\gamma _{12}^{0+}=\gamma _{12}^0=\gamma _{12}^0`$. Notice that $`\gamma _{12}^0`$ increases monotonously with $`N_e`$ as the single-particle band is filled up. This is an important contribution to the band-filling dependence of $`E_\mathrm{C}`$ (see Fig. 3). (ii) In the delocalized limit, $`E_\mathrm{C}(\gamma _{12}^0)=0`$ for all the $`N_e`$ for which $`W(\gamma _{12}^0)`$ derives from a single Slater determinant . Moreover, the divergence of $`E_\mathrm{C}/\gamma _{12}`$ for $`\gamma _{12}=\gamma _{12}^0`$ indicates that $`\gamma _{12}^{gs}<\gamma _{12}^0`$ for arbitrary small $`U>0`$, as expected from perturbation theory. (iii) Starting from $`\gamma _{12}=\gamma _{12}^0`$, $`E_\mathrm{C}(\gamma _{12})`$ decreases with decreasing $`\gamma _{12}`$ reaching its lowest possible value $`E_\mathrm{C}^{\mathrm{}}=E_{\mathrm{HF}}`$ for $`\gamma _{12}=\gamma _{12}^\mathrm{}+`$ ($`N_eN_a`$). The same behavior is of course observed for $`\gamma _{12}<0`$. In particular, $`E_\mathrm{C}=E_{\mathrm{HF}}`$ also for $`\gamma _{12}=\gamma _{12}^{\mathrm{}}`$. As shown in Fig. 3, $`E_\mathrm{C}^{\mathrm{}}`$ decreases rapidly with increasing $`N_e`$, since $`E_{\mathrm{HF}}`$ increases quadratically with the electron density . (iv) On bipartite lattices $`\gamma _{12}^\mathrm{}+=\gamma _{12}^{\mathrm{}}=\gamma _{12}^{\mathrm{}}`$, while on non-bipartite structures one generally has $`|\gamma _{12}^\mathrm{}+||\gamma _{12}^{\mathrm{}}|`$, since the single-particle spectrum is different for positive and negative energies. The decrease of $`E_\mathrm{C}`$ with decreasing $`|\gamma _{12}|`$ shows that the reduction of the Coulomb energy due to correlations is done at the expense of kinetic energy or electron delocalization, as already discussed for $`N_e=N_a`$ (Fig. 1). (v) $`\gamma _{12}^{\mathrm{}}>0`$ for all $`N_e<N_a`$ ($`\gamma _{12}^{\mathrm{}}=0`$ for $`N_e=N_a`$). $`\gamma _{12}^{\mathrm{}}`$ represents the largest NN bond order that can be constructed under the constraint of vanishing Coulomb repulsion energy. A lower bound for $`\gamma _{12}^{\mathrm{}}`$ is given by the bond order $`\gamma _{12}^{FM}`$ in the fully-polarized ferromagnetic state ($`\gamma _{12}^{\mathrm{}}\gamma _{12}^{FM}`$). This is obtained by occupying the lowest single-particle states with all electrons of the same spin ($`N_eN_a`$). Therefore, $`\gamma _{12}^{FM}`$ increases with $`N_e`$ for $`N_eN_a/2`$ and then decreases for $`N_a/2<N_eN_a`$ reaching $`\gamma _{12}^{FM}=0`$ at half-band filling ($`\gamma _{12}^{FM}>0`$ for $`N_e<N_a`$). In this way the non-monotonous dependence of $`\gamma _{12}^{\mathrm{}}`$ on $`N_e`$ can be explained (see Fig. 3). (vi) The correlation energy is constant and equal to $`E_{\mathrm{HF}}`$ for $`\gamma _{12}^{\mathrm{}}\gamma _{12}\gamma _{12}^\mathrm{}+`$. These values of $`\gamma _{12}`$ can never correspond to the ground-state energy of the Hubbard model, since in this range increasing $`\gamma _{12}`$ always lowers the kinetic energy ($`t>0`$) without increasing the Coulomb repulsion ($`\gamma _{12}^{\mathrm{}}\gamma _{12}^{gs}\gamma _{12}^0`$). For $`\gamma _{12}^{\mathrm{}}<\gamma _{12}<\gamma _{12}^\mathrm{}+`$, $`\gamma _{12}`$ cannot be represented by a ground state of Eq. (15). In this range $`\gamma _{12}`$ can be derived from a linear combination of states having minimal Coulomb repulsion .

In order to compare the functional dependences of the correlation energy for different band fillings, it is useful to scale $`E_\mathrm{C}`$ in units of the Hartree-Fock energy and to bring the relevant domains $`\gamma _{12}^{\mathrm{}}\gamma _{12}\gamma _{12}^0`$ of different $`N_e`$ to a common range. In Fig. 4, $`\epsilon _\mathrm{C}=E_\mathrm{C}/E_{\mathrm{HF}}`$ is shown as a function of $`g_{12}=(\gamma _{12}\gamma _{12}^{\mathrm{}})/(\gamma _{12}^0\gamma _{12}^{\mathrm{}})`$. We observe that the results for $`\epsilon _\mathrm{C}(g_{12})`$ are remarkably similar for all band-fillings. The largest deviations from the common trend are found for $`N_e=4`$. As already discussed for $`N_a=N_e=4`$, this anomalous behavior is related to the degeneracy of the single-particle spectrum and to the finite size of system. Fig. 4 shows that for the Hubbard model the largest part of the dependence of $`E_\mathrm{C}(\gamma _{12})`$ on band filling comes from $`E_{\mathrm{HF}}`$, $`\gamma _{12}^0`$ and $`\gamma _{12}^{\mathrm{}}`$. Similar conclusions are derived from the results for the infinite 1D chain presented in Fig. 5. For a given $`g_{12}`$, $`\epsilon _\mathrm{C}(g_{12})`$ depends weakly on $`N_e/N_a`$ if the carrier density is low ($`N_e/N_a0.4`$), and tends to increase as we approach half-band filling \[see Fig. 5(b)\]. For high carrier densities it become comparatively more difficult to minimize the Coulomb energy for a given degree of delocalization $`g_{12}`$. The effect is most pronounced for $`g_{12}0.8`$$`0.9`$, i.e., close to the uncorrelated limit. As we approach the strongly correlated limit ($`g_{12}0.4`$) the dependence of $`\epsilon _\mathrm{C}`$ on $`N_e/N_a`$ is very weak even for $`N_e/N_a1`$. One concludes that $`\epsilon _\mathrm{C}(g_{12})`$ is a useful basis for introducing practical approximations on more complex systems.

The correlation energy $`E_\mathrm{C}`$ is a universal functional of the complete single-particle density matrix $`\gamma _{ij}`$. $`E_\mathrm{C}[\gamma _{ij}]`$ and $`W[\gamma _{ij}]`$ may depend on $`N_a`$ and $`N_e`$ but are independent of $`t_{ij}`$ and in particular of the lattice structure. The functional $`E_\mathrm{C}(\gamma _{12})`$ considered in this paper depends by definition on the type of lattice, since the constraints imposed in the minimization only apply to NN bonds. In order to investigate this problem we have determined $`E_\mathrm{C}(\gamma _{12})`$ for 2D and 3D finite clusters having $`N_a12`$ sites and periodic boundary conditions. In Fig. 6 we compare these results with those of the 1D $`12`$-site periodic ring. As shown in the inset figure, the qualitative behavior is in all cases very similar. The main quantitative differences come from the domain of representability of $`\gamma _{12}`$, i.e., from the values of $`\gamma _{12}^{0+}`$ and $`\gamma _{12}^0`$ ($`\gamma _{12}^0\gamma _{12}\gamma _{12}^{0+}`$). Once scaled as a function of $`\gamma _{12}/\gamma _{12}^0`$, $`E_\mathrm{C}`$ depends rather weakly on the lattice structure. Notice that the Hartree-Fock energy $`E_{\mathrm{HF}}=(U/4)N_a`$ is the same for all structures. However, for the BCC structure we obtain $`W(\gamma _{12}^0)<E_{\mathrm{HF}}`$, i.e., $`E_\mathrm{C}(\gamma _{12}^0)<0`$, due to degeneracies in the single-particle spectrum of the considered finite cluster \[see inset Fig. 6(b)\]. In order to correct for this finite size effect it is here more appropriate to consider $`\epsilon _\mathrm{C}=[E_\mathrm{C}(\gamma _{12})E_\mathrm{C}(\gamma _{12}^0)]/W(\gamma _{12}^0)`$. Still, the differences in $`\epsilon _\mathrm{C}`$ between BCC and FCC structures appear to be more important than between square and triangular 2D lattices. This is probably related to the degeneracies in the spectrum of the BCC cluster, as already observed for rings with $`N_e=4m`$ \[Figs. 2 and 4(a)\].

The largest changes in $`\epsilon _\mathrm{C}`$ for different lattice structures are observed for intermediate degree of delocalization ($`g_{12}0.7`$$`0.9`$, see Fig. 6). Note that there is no monotonic trend as a function of the lattice dimension. For example, for $`g_{12}=0.7`$$`0.9`$, $`\epsilon _\mathrm{C}`$ first increases somewhat as we go from 1D to 2D lattices, but it then decreases coming close to the 1D curve for the 3D FCC lattice \[$`\epsilon _\mathrm{C}(2\mathrm{D})>\epsilon _\mathrm{C}(\mathrm{FCC})>\epsilon _\mathrm{C}(1\mathrm{D})>\epsilon _\mathrm{C}(\mathrm{BCC})`$ for $`0.7g_{12}0.9`$\]. Finally, it is worth noting that in the strongly correlated limit ($`g_{12}0.3`$) the results for $`\epsilon _\mathrm{C}(g_{12})`$ are nearly the same for all considered lattice structures (see Fig. 6). This should be useful in order to develop simple general approximations to $`E_\mathrm{C}(\gamma _{12})`$ in this limit.

### B Attractive interaction $`U<0`$

The attractive Hubbard model describes itinerant electrons with local intra-atomic pairing ($`U<0`$). The electronic correlations are very different from those found in the repulsive case discussed so far. In particular Levy’s interaction energy functional $`W(\gamma _{ij})`$ now correspond to the maximum average number of double occupation for a given $`\gamma _{ij}`$ \[see Eq. (11)\]. Therefore, it is very interesting to investigate the properties of the correlation energy functional $`E_\mathrm{C}(\gamma _{ij})`$ also for $`U<0`$ and to contrast them with the results of the previous section.

In Fig. 7 the correlation energy $`E_\mathrm{C}(\gamma _{12})`$ of the attractive Hubbard model is given at half-band filling for various finite rings ($`N_a12`$) and for the infinite 1D chain ($`N_e=N_a`$). The band-filling dependence of $`E_\mathrm{C}(\gamma _{12})`$ is shown in Fig. 8 for a finite $`12`$-site ring ($`N_eN_a=12`$). As in the repulsive case, $`\gamma _{12}^0\gamma _{12}\gamma _{12}^{0+}`$ since the domain of representability of $`\gamma _{12}`$ is independent of the form or type of the interaction. Moreover, $`E_\mathrm{C}(\gamma _{12})=E_\mathrm{C}(\gamma _{12})`$ due to the electron-hole symmetry of bipartite lattices . Starting from $`\gamma _{12}^{0+}`$ or $`\gamma _{12}^0`$ ($`\gamma _{12}^{0+}=\gamma _{12}^0=\gamma _{12}^0`$ on bipartite lattices), $`E_\mathrm{C}(\gamma _{12})`$ decreases with decreasing $`|\gamma _{12}|`$ reaching the minimum $`E_\mathrm{C}^{\mathrm{}}`$ for $`\gamma _{12}=\gamma _{12}^{\mathrm{}}`$ and for $`\gamma _{12}=\gamma _{12}^{\mathrm{}}`$ ($`\gamma _{12}^\mathrm{}+=\gamma _{12}^{\mathrm{}}=\gamma _{12}^{\mathrm{}}`$ in this case). For $`N_e`$ even, $`W(\gamma _{12}^{\mathrm{}})=N_eU/2`$, and for $`N_e`$ odd, $`W(\gamma _{12}^{\mathrm{}})=(N_e1)U/2`$, which correspond to the maximum number of electron pairs that can be formed. For $`N_e`$ even, the minimum $`E_\mathrm{C}^{\mathrm{}}=U(N_e/2)[1N_e/(2N_a)]`$ is achieved only for a complete electron localization (i.e, $`\gamma _{12}^{\mathrm{}}=0`$). In contrast, for odd $`N_e`$ a finite-size effect is observed. In this case, one of the electrons remains unpaired even in the limit of strong electron correlations and the minimum of $`E_\mathrm{C}`$ is $`E_\mathrm{C}^{\mathrm{}}=U[(N_e1)/2][1(N_e+1)/(2N_a)]`$. Moreover, non-vanishing $`\gamma _{12}^{\mathrm{}}`$ are obtained as a result of the delocalization of the unpaired electron. $`\gamma _{12}^{\mathrm{}}`$ represents the maximum bond order that can be obtained when $`(N_e1)/2`$ electron pairs are formed ($`\gamma _{12}^{\mathrm{}}0`$ for $`N_a\mathrm{}`$, $`N_e`$ odd). Notice that in all cases the ground state $`\gamma _{12}^{gs}`$ is found in the interval $`\gamma _{12}^{\mathrm{}}\gamma _{12}^{gs}\gamma _{12}^0`$.

It is interesting to observe that $`E_\mathrm{C}(\gamma _{12})`$ can be appropriately scaled in a similar way as for $`U>0`$. In Fig. 8(b), $`\epsilon _\mathrm{C}(g_{12})=E_\mathrm{C}/|E_\mathrm{C}^{\mathrm{}}|`$ is shown as a function of the degree of delocalization $`g_{12}=(\gamma _{12}\gamma _{12}^{\mathrm{}})/(\gamma _{12}^0\gamma _{12}^{\mathrm{}})`$. $`\epsilon _\mathrm{C}(g_{12})`$ presents a pseudo-universal behavior in the sense that it depends weakly on $`N_a`$ and $`N_e`$. The main deviations from the common trend are found for $`N_e=N_a=4`$. As already discussed for $`U>0`$, this is a consequence of degeneracies in the single-particle spectrum. In this case, the wave function corresponding to the minimum in Levy’s functional for $`\gamma _{12}\gamma _{12}^0`$ \[Eq. (11)\] cannot be described by a single Slater determinant and $`W(\gamma _{12}\gamma _{12}^0)<E_{\mathrm{HF}}`$.

## IV Conclusion

Density-matrix functional theory has been applied to lattice Hamiltonians taking the Hubbard model as a particularly relevant example. In this framework the basic variable is the single-particle density matrix $`\gamma _{ij}`$ and the key unknown is the correlation energy functional $`E_\mathrm{C}[\gamma _{ij}]`$. The challenge is therefore to determine $`E_\mathrm{C}[\gamma _{ij}]`$ or to provide with useful accurate approximations for it. In this paper we presented a systematic study of the functional dependence of $`E_\mathrm{C}(\gamma _{12})`$ on periodic lattices, where $`\gamma _{12}`$ is the density-matrix element between nearest neighbors ($`\gamma _{ij}=\gamma _{12}`$ for all NN $`ij`$). Based on finite-cluster exact diagonalizations and on the Bethe-Ansatz solution of the 1D chain, we derived rigorous results for $`E_\mathrm{C}(\gamma _{12})`$ of the Hubbard model as a function of the number of sites $`N_a`$, band filling $`N_e/N_a`$ and lattice structure. A basis for applications of density-matrix functional theory to many-body lattice models is thereby provided. The observed pseudo-universal behavior of $`\epsilon _\mathrm{C}(g_{12})=E_\mathrm{C}/E_{\mathrm{HF}}`$ as a function of $`g_{12}=(\gamma _{12}\gamma _{12}^{\mathrm{}})/(\gamma _{12}^0\gamma _{12}^{\mathrm{}})`$ encourages transferring $`\epsilon _\mathrm{C}(g_{12})`$ from finite-size systems to infinite lattices or even to different lattice geometries. In fact, the exact $`E_\mathrm{C}(\gamma _{12})`$ of the Hubbard dimer has been recently used to infer a simple general Ansatz for $`E_\mathrm{C}(\gamma _{12})`$ . With this approximation to $`E_\mathrm{C}(\gamma _{12})`$ the ground-state energies and charge-excitation gaps of 1D and 2D lattices have been determined successfully in the whole range of $`U/t`$. Further investigations, for example, by considering magnetic impurity models or more complex multiband Hamiltonians, are certainly worthwhile.

###### Acknowledgements.

The authors gratefully acknowledge the financial support provided by CONACyT-Mexico (RLS) and by the Alexander von Humboldt Foundation (GMP).