# 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.