# Hartree-Fock decoupling of Hubbard model

Hartree-Fock approximation requires wavefunctions be as separable as possible. I know the basic idea of Hartree-Fock but having some trouble in formalism of second quantization. I am trying to replicate this paper in which they are taking hubbard model in 2D and there is a wavefunction given in equ.3 (also given below) that is being used to find variational energy of Hubbard model.
Hubbard model (in fourier space) is given as: $$H=H_{1}+H_{2}=\sum_k \epsilon(k)c_k^\dagger c_k -\frac{V}{N}\sum_{k,k',q}\epsilon(q)c_{k+q}^\dagger c_kc_{k'-q}^\dagger c_{k'}$$ where $\epsilon(k)=-\cos(k_x)-\cos(k_y)$ and variational wavefunction is given as: $$|\psi\rangle =\underset{k\epsilon RBZ}{\Pi}\frac{\exp{[\tilde{\alpha}\epsilon(k)]}c_k^\dagger+\exp{[-\tilde{\alpha}\epsilon(k)]}c_{k-Q}^\dagger}{\sqrt{2\cosh{[2\alpha\epsilon(k)]}}}|0\rangle$$ where $\tilde{\alpha}=\alpha+i\eta$ is complex variational parameter and $Q=(\pi,\pi)$ . Varitional groudstate energy is $E_0=\underset{\tilde{\alpha}}{min}\langle \psi|H|\psi\rangle =\langle H_{1}\rangle +\langle H_{2}\rangle$.
$\langle H_{1}\rangle$ can be solved simply but I am not understanding how to solve $\langle H_{2}\rangle$. According to this article they are using Hartree-Fock decoupling. System they are studying is spinless fermionic at half-filling which means we can have anomalous ($\langle c_k^\dagger c_{k-Q}\rangle \neq0$) expectation values.

How can I decouple $c_{k+q}^\dagger c_kc_{k'-q}^\dagger c_{k'}$ operators to get the results which are given in equ.6 and equ.7 in this article?

My attempt: (Edited)

$$\langle H_2\rangle =\langle c_{k+q}^\dagger c_{k'-q}^\dagger c_{k'} c_k\rangle \approx \langle c_{k+q}^\dagger c_k\rangle \langle c_{k'-q}^\dagger c_{k'}\rangle -\langle c_{k+q}^\dagger c_{k'}\rangle \langle c_{k'-q}^\dagger c_k\rangle$$ First term is Hartree term and second is Fock term. Let's first take Hartree term at $k'=k+q$ and take $q=0$: $$\langle Hartree Term\rangle _{q=0}\approx \epsilon(0)\sum_{k,k'}[\langle c_{k}^{\dagger} c_k\rangle \langle c_{k'}^\dagger c_{k'}\rangle ]$$ $$=\epsilon(0)[\sum_k \langle \psi|c_k^\dagger c_k|\psi\rangle \sum_{k'} \langle \psi|c_{k'}^\dagger c_{k'}|\psi\rangle ]$$ As $c_k^\dagger c_k=n_k=$number operator and $\sum_k n_k=\sum_{k'}n_{k'}=N/2$ for half filling system where $N=$ total number of sites. So above equation can be written as: $$\langle Hartree Term\rangle _{q=0}=\epsilon(0)N^2/4$$ First term is exactly similar to first term of equ.6 in the article I am replicating. Now let's repeat all setup but with $q=Q$: $$\langle HartreeTerm\rangle _{q=Q}=-\epsilon(0)[\sum_k\frac{\cos(2\eta\epsilon(k))}{2\cosh(2\alpha\epsilon(k))}]^2$$ This is also similar to one given in that article. So, far so good. Now time for Fock terms. Let's take $k'=k+q$ but this time $q\ne0$ so: $$\langle FockTerm\rangle _{q=k'-k}=\sum_{k,k'}\epsilon(k'-k)\langle c_{k'}^\dagger c_{k'}\rangle \langle c_k^\dagger c_k\rangle$$ I am not sure how to solve this one because $k$ and $k'$ are entangled with each other in $\epsilon(k'-k)$ factor. NEED HELP ._. I am quite sure that this term is the one which will give equ.7b and equ.7c can be obtained from Fock term by taking $k'=k+q+Q$

Someone please help me or give me some references from which I can understand this kind of calculations.