For a Hubbard hamiltonian, $$\hat{H} = \sum_{i,j,\sigma} t_{i j} \hat{c}^\dagger_{i \sigma} \hat{c}_{j \sigma} + \sum_i U_i \hat{n}_{i \uparrow} \hat{n}_{i \downarrow}, $$ the mean-field solution can be obtained by using: $ \hat{n}_{i \uparrow} \hat{n}_{i \downarrow} \approx n_{i \uparrow} \hat{n}_{i \downarrow} + n_{i \downarrow} \hat{n}_{i \uparrow} - n_{i \uparrow} n_{i \downarrow} $ and solving self-consistenly the resulting one-electron hamiltonian. In general, for a generic many-body term $\hat{c}^\dagger_{i \sigma}\hat{c}^\dagger_{j \sigma^\prime}\hat{c}_{k \sigma^\prime}\hat{c}_{l \sigma},$ the approximation is: $$\hat{c}^\dagger_{i \sigma}\hat{c}^\dagger_{j \sigma^\prime}\hat{c}_{k \sigma^\prime}\hat{c}_{l \sigma} \approx \langle \hat{c}^\dagger_{i \sigma} \hat{c}_{l \sigma} \rangle \hat{c}^\dagger_{j \sigma^\prime}\hat{c}_{k \sigma^\prime} + \langle \hat{c}^\dagger_{j \sigma^\prime}\hat{c}_{k \sigma^\prime} \rangle \hat{c}^\dagger_{i \sigma} \hat{c}_{l \sigma} - \langle \hat{c}^\dagger_{j \sigma^\prime}\hat{c}_{k \sigma^\prime} \rangle \langle \hat{c}^\dagger_{i \sigma} \hat{c}_{l \sigma} \rangle \\ -\langle \hat{c}^\dagger_{i \sigma} \hat{c}_{k \sigma^\prime} \rangle \hat{c}^\dagger_{j \sigma^\prime}\hat{c}_{l \sigma} - \langle \hat{c}^\dagger_{j \sigma^\prime} \hat{c}_{l \sigma}\rangle \hat{c}^\dagger_{i \sigma}\hat{c}_{k \sigma^\prime} + \langle \hat{c}^\dagger_{j \sigma^\prime} \hat{c}_{l \sigma} \rangle \langle \hat{c}^\dagger_{i \sigma}\hat{c}_{k \sigma^\prime} \rangle $$

(we assume the particle number is conserved, so that terms like $\langle \hat{c}^\dagger \hat{c}^\dagger \rangle $ are not important). Now, $S_z$ also being a good quantum number, if $\sigma \neq \sigma^\prime$ we will also ignore the terms on the second line, as done for the approximation above for Hubbard model. The Hubbard model on a lattice can be solved easily on mean-field thanks to the fact that on the self-consistency loop we only need the charges, and never non-diagonal charges (other elements of the density matrix) like $\langle \hat{c}^\dagger_{j \sigma^\prime} \hat{c}_{l \sigma} \rangle$.

However think now of a slightly more complicated model on a lattice. It's like the hamiltonian $\hat{H}$ above but now we add many-body interactions between neighboring sites; something like: $$ \dfrac{1}{2}\sum_{i,j, \sigma,\sigma^\prime} J_{ij} \hat{n}_{i \sigma} \hat{n}_{j \sigma^\prime}.$$ Now, for the terms with $\sigma \neq \sigma^\prime$ we have the simple approximation $\hat{n}_{i \uparrow} \hat{n}_{j \downarrow} \approx n_{i \uparrow} \hat{n}_{j \downarrow} + n_{j \downarrow} \hat{n}_{i \uparrow} - n_{i \uparrow} n_{j \downarrow}$. However, for terms like $\hat{n}_{i \sigma} \hat{n}_{j \sigma}$, when doing the mean-field approximation we will unavoidable find terms like those of the second line for the generic equation above; that is, the term: $$-\langle \hat{c}^\dagger_{i \sigma} \hat{c}_{j \sigma} \rangle \hat{c}^\dagger_{j \sigma}\hat{c}_{i \sigma} - \langle \hat{c}^\dagger_{j \sigma} \hat{c}_{i \sigma}\rangle \hat{c}^\dagger_{i \sigma}\hat{c}_{j \sigma} + \langle \hat{c}^\dagger_{i \sigma} \hat{c}_{j \sigma} \rangle \langle \hat{c}^\dagger_{j \sigma}\hat{c}_{i \sigma} \rangle,$$ and here is the problem. We need to include somehow in the self-consistency the non-diagonal charges $n_{ij\sigma}:=\langle \hat{c}^\dagger_{i \sigma} \hat{c}_{j \sigma} \rangle. $

In order to do self-consistency on the charges we should express the $n_{ij \sigma}$ in terms of the $n_{i \sigma},$ but I don't think there is any analytical way of doing this even when all expected values are taken over Slater deterinants. Am I right? I know that I can express sums of non-diagonal charges, $\sum_j n_{i j \sigma}$ comfortably as sums of electron-hole interactions, $\sum_{i \neq j} n_{i \sigma} \left( 1 - n_{i \sigma} \right)$, which is convenient to write the Hartree-Fock energy of lattice models when the $J_{ij}$ parameters are constant. But in general, an isolated $n_{ij \sigma}$ seems untreatable. Two possible solutions come to mind:

1. Finding approximations for $n_{ij \sigma}$ in terms of $n_{k \sigma}$. How can this be done? Can the identity with $\sum_{i \neq j} n_{i \sigma} \left( 1 - n_{i \sigma} \right)$ (that is exact over Slater determinants) be used?

2. Doing self-consistency in terms of $n_{ij \sigma}$ instead of just the charges. I have tried this last option for things like PPP parametrization on hydrocarbons and I've encountered cases on which, even using a mixing algorithm on all the density matrix, the convergence is quite problematic. Has this be done before and how? (meaning, Lattice Hartree-Fock where self-consistency is done in terms of the full density matrix)

For a Hubbard hamiltonian, $$\hat{H} = \sum_{i,j,\sigma} t_{i j} \hat{c}^\dagger_{i \sigma} \hat{c}_{j \sigma} + \sum_i U_i \hat{n}_{i \uparrow} \hat{n}_{i \downarrow}, $$ the mean-field solution can be obtained by using: $ \hat{n}_{i \uparrow} \hat{n}_{i \downarrow} \approx n_{i \uparrow} \hat{n}_{i \downarrow} + n_{i \downarrow} \hat{n}_{i \uparrow} - n_{i \uparrow} n_{i \downarrow} $ and solving self-consistenly the resulting one-electron hamiltonian. In general, for a generic many-body term $\hat{c}^\dagger_{i \sigma}\hat{c}^\dagger_{j \sigma^\prime}\hat{c}_{k \sigma^\prime}\hat{c}_{l \sigma},$ the approximation is: $$\hat{c}^\dagger_{i \sigma}\hat{c}^\dagger_{j \sigma^\prime}\hat{c}_{k \sigma^\prime}\hat{c}_{l \sigma} \approx \langle \hat{c}^\dagger_{i \sigma} \hat{c}_{l \sigma} \rangle \hat{c}^\dagger_{j \sigma^\prime}\hat{c}_{k \sigma^\prime} + \langle \hat{c}^\dagger_{j \sigma^\prime}\hat{c}_{k \sigma^\prime} \rangle \hat{c}^\dagger_{i \sigma} \hat{c}_{l \sigma} - \langle \hat{c}^\dagger_{j \sigma^\prime}\hat{c}_{k \sigma^\prime} \rangle \langle \hat{c}^\dagger_{i \sigma} \hat{c}_{l \sigma} \rangle \\ -\langle \hat{c}^\dagger_{i \sigma} \hat{c}_{k \sigma^\prime} \rangle \hat{c}^\dagger_{j \sigma^\prime}\hat{c}_{l \sigma} - \langle \hat{c}^\dagger_{j \sigma^\prime} \hat{c}_{l \sigma}\rangle \hat{c}^\dagger_{i \sigma}\hat{c}_{k \sigma^\prime} + \langle \hat{c}^\dagger_{j \sigma^\prime} \hat{c}_{l \sigma} \rangle \langle \hat{c}^\dagger_{i \sigma}\hat{c}_{k \sigma^\prime} \rangle $$

(we assume the particle number is conserved, so that terms like $\langle \hat{c}^\dagger \hat{c}^\dagger \rangle $ are not important). Now, $S_z$ also being a good quantum number, if $\sigma \neq \sigma^\prime$ we will also ignore the terms on the second line, as done for the approximation above for Hubbard model. The Hubbard model on a lattice can be solved easily on mean-field thanks to the fact that on the self-consistency loop we only need the charges, and never non-diagonal charges (other elements of the density matrix) like $\langle \hat{c}^\dagger_{j \sigma^\prime} \hat{c}_{l \sigma} \rangle$.

However think now of a slightly more complicated model on a lattice. It's like the hamiltonian $\hat{H}$ above but now we add many-body interactions between neighboring sites; something like: $$ \dfrac{1}{2}\sum_{i,j, \sigma,\sigma^\prime} J_{ij} \hat{n}_{i \sigma} \hat{n}_{j \sigma^\prime}.$$ Now, for the terms with $\sigma \neq \sigma^\prime$ we have the simple approximation $\hat{n}_{i \uparrow} \hat{n}_{j \downarrow} \approx n_{i \uparrow} \hat{n}_{j \downarrow} + n_{j \downarrow} \hat{n}_{i \uparrow} - n_{i \uparrow} n_{j \downarrow}$. However, for terms like $\hat{n}_{i \sigma} \hat{n}_{j \sigma}$, when doing the mean-field approximation we will unavoidable find terms like those of the second line for the generic equation above; that is, the term: $$-\langle \hat{c}^\dagger_{i \sigma} \hat{c}_{j \sigma} \rangle \hat{c}^\dagger_{j \sigma}\hat{c}_{i \sigma} - \langle \hat{c}^\dagger_{j \sigma} \hat{c}_{i \sigma}\rangle \hat{c}^\dagger_{i \sigma}\hat{c}_{j \sigma} + \langle \hat{c}^\dagger_{i \sigma} \hat{c}_{j \sigma} \rangle \langle \hat{c}^\dagger_{j \sigma}\hat{c}_{i \sigma} \rangle,$$ and here is the problem. We need to include somehow in the self-consistency the non-diagonal charges $n_{ij\sigma}:=\langle \hat{c}^\dagger_{i \sigma} \hat{c}_{j \sigma} \rangle. $

In order to do self-consistency on the charges we should express the $n_{ij \sigma}$ in terms of the $n_{i \sigma},$ but I don't think there is any analytical way of doing this even when all expected values are taken over Slater deterinants. Am I right? I know that I can express sums of non-diagonal charges, $\sum_j n_{i j \sigma}$ comfortably as sums of electron-hole interactions, $\sum_{i \neq j} n_{i \sigma} \left( 1 - n_{i \sigma} \right)$, which is convenient to write the Hartree-Fock energy of lattice models when the $J_{ij}$ parameters are constant. But in general, an isolated $n_{ij \sigma}$ seems untreatable. Two possible solutions come to mind:

1. Finding approximations for $n_{ij \sigma}$ in terms of $n_{k \sigma}$. How can this be done? Can the identity with $\sum_{i \neq j} n_{i \sigma} \left( 1 - n_{i \sigma} \right)$ (that is exact over Slater determinants) be used?

2. Doing self-consistency in terms of $n_{ij \sigma}$ instead of just the charges. I have tried this last option for things like PPP parametrization on hydrocarbons and I've encountered cases on which, even using a mixing algorithm on all the density matrix, the convergence is quite problematic. Has this be done before and how? (meaning, Lattice Hartree-Fock where self-consistency is done in terms of the full density matrix)