Inconsistency of single-particle approach with non-interacting many-body systems Suppose that we have a two-particle non-interacting system and we are to solve the following Hamiltonian for this system:
$$H=\epsilon_1 c_1^{\dagger}c_1 +\epsilon_2 c_2^{\dagger}c_2 +vc_1^{\dagger}c_2+v^* c_2^{\dagger}c_1$$
If we solve this Hamiltonian by two different approach, single-particle and many particle, we obtain different results as follows:
1- In single particle approach we have a $2\times 2$ matrix in the bases $c_1^{\dagger}|0>$ and $c_2^{\dagger}|0>$ where $|0>$ is the vacuum state. Solution of this matrix gives the two new states,$\psi_1$ and $\psi_2$, with energies:
$$E_1=\frac{1}{2} [\epsilon_1+\epsilon_2 +\sqrt{(\epsilon_1+\epsilon_2)^2-4(\epsilon_1 \epsilon_2 -|v|^2)}]$$ and
$$E_1=\frac{1}{2} [\epsilon_1+\epsilon_2 -\sqrt{(\epsilon_1+\epsilon_2)^2-4(\epsilon_1 \epsilon_2 -|v|^2)}]$$
2- If we solve this problem in many-body approach, we have the Hamiltonian in the four bases; $|0>$, $c_1^{\dagger}|0>$, $c_2^{\dagger}|0>$, $c_1^{\dagger}c_2^{\dagger}|0>$ and we have three blocks of the Hamiltonian, one for zero particle, another for one particle in each state and the other for two particles in both states.
It is clear that the two-particle block gives one of the solution of the Hamiltonian (with $c_1^{\dagger}c_2^{\dagger}|0>$ as its eigenstate) when we have two particles. However, if we solve this two-particle problem with single-particle approach in such a way that we put one particle in the state $\psi_1$ and the other in the state $\psi_2$, the solution is different from $c_1^{\dagger}c_2^{\dagger}|0>$.
Any help is appreciated.
 A: the states $c^{\dagger}_1 c^{\dagger}_2 |0\rangle$ and $\psi^{\dagger}_1 \psi^{\dagger}_2 |0\rangle$ should be identical, up to a global phase. You can see that both of them are eigenstates of the Hamiltonian with the same energy $\epsilon_1 + \epsilon_2$.
You can see that if you act directly on the state with the Hamiltonian you'll get
$$H c^{\dagger}_1 c^{\dagger}_2 |0\rangle = (\epsilon_1 + \epsilon_2)c^{\dagger}_1 c^{\dagger}_2 |0\rangle$$
(notice that the $\nu$-terms give zero as $c^{\dagger}_2 c_1 c^{\dagger}_1 c^{\dagger}_2$ has two $c^{\dagger}_2$ with no $c_2$ between them, so it will vanish, and similarly $c^{\dagger}_1 c_2 c^{\dagger}_1 c^{\dagger}_2$).
For the composition of the two single particle eigenstates you found, the energy of the doubly occupied state is the sum of the energies (as they are non-interacting) which again gives $\epsilon_1 + \epsilon_2$.
If you'll write the single states explicitly and carry the multiplication of them (taking care to normalize both of them) you'll see that $\psi^{\dagger}_1 \psi^{\dagger}_2 \propto c^{\dagger}_1 c^{\dagger}_2$
You can see that even without finding the expansion coefficients explicitly. Just by the notion of orthogonality we know that
$$\psi^{\dagger}_1 = e^{i\varphi}\cos(\theta) c^{\dagger}_1 + \sin(\theta) c^{\dagger}_2 \\
\psi^{\dagger}_2 = \cos(\theta) c^{\dagger}_2 - e^{i\varphi}\sin(\theta) c^{\dagger}_1 $$
which is the only way to get $\left\{ \psi_i, \psi^{\dagger}_j \right\} = \delta_{i,j}$, up to arbitrary phase assigned to the states (here I just chose the relative phase to always assigned to $c^{\dagger}_1$, but any other choice is valid). So taking this form
$$ \psi^{\dagger}_1 \psi^{\dagger}_2 = \left[e^{i\varphi}\cos(\theta) c^{\dagger}_1 + \sin(\theta) c^{\dagger}_2\right]\left[\cos(\theta) c^{\dagger}_2 - e^{i\varphi}\sin(\theta) c^{\dagger}_1\right] = e^{i\varphi}\left[\cos^{2}(\theta) + \sin^{2}(\theta)\right]c^{\dagger}_1c^{\dagger}_2 = e^{i\varphi}c^{\dagger}_1 c^{\dagger}_2$$
which is what we wanted to show. The global phase of course doesn't matter at all. Note that here I used $c^{\dagger}_i c^{\dagger}_i = 0$ and also $c^{\dagger}_1 c^{\dagger}_2 = -c^{\dagger}_2 c^{\dagger}_1$.
