Apparent problem in using Wick's theorem to calculate matrix elements of two-body operators In the second quantized notation, a two-body operator $\hat{O}$ can be written as
$$\hat{O} = \sum\limits_{x_1,x_2,x_3,x_4} O_{x_1,x_2,x_3,x_4} a^\dagger_{x_1}a^\dagger_{x_2}a_{x_4}a_{x_3} ,$$ where $$ O_{x_1,x_2,x_3,x_4} = \langle x_1x_2 \rvert\hat{O}\rvert x_3x_4\rangle .$$
However, if we write $\langle ij \rvert\hat{O}\rvert kl\rangle$ as
$$\langle \Omega \rvert a_{i}a_{j} \hat{O}a^\dagger_{k}a^\dagger_{l} \rvert \Omega\rangle = \sum\limits_{x_1,x_2,x_3,x_4} O_{x_1,x_2,x_3,x_4} \langle \Omega \rvert a_{i}a_{j}a^\dagger_{x_1}a^\dagger_{x_2} a_{x_4}a_{x_3}  a^\dagger_{k}a^\dagger_{l} \rvert \Omega\rangle $$
The vacuum expectation can be calculated using Wick's theorem to be
$$\langle \Omega \rvert a_{i}a_{j}a^\dagger_{x_1}a^\dagger_{x_2} a_{x_4}a_{x_3}  a^\dagger_{k}a^\dagger_{l} \rvert \Omega\rangle = \delta_{i x_1}\delta_{j x_2}\delta_{k x_4}\delta_{l x_3}-\delta_{i x_1}\delta_{j x_2}\delta_{k x_3}\delta_{l x_4}-\delta_{i x_2}\delta_{j x_1}\delta_{k x_4}\delta_{l x_3}+\delta_{i x_2}\delta_{j x_1}\delta_{k x_3}\delta_{l x_4}$$
Then $$\langle ij \rvert\hat{O}\rvert kl\rangle = O_{ijlk}- O_{ijkl}-O_{jilk}+O_{jikl} ,$$ not just $O_{ijkl}$!
Where am I committing a mistake that results in this apparent contradiction?
 A: In fact by using the symmetry properties of the two-particle states you can prove that all the terms in the right are exactly equal, which gives a paradoxical result: $\langle ij |O|kl \rangle =4 \langle ij |O|kl \rangle $. I think there might be an additional factor of one fourth in the original definition if the expectation value is to be taken seriously. In some texts (Coleman, Mahan, Girvin & Yang) they define the coefficient in front of the four operators as $O_{ijkl} =\int d1 d2 \overline{\phi}_i (1) \overline{\phi}_j (2)\phi_l (2) \phi_k (1)O(1,2)$, it can just be abuse of notation to write it explicitly as the expectation value for the properly symmetrized two-particle states.
A: The answer by @juanmendez is correct - these terms are equal due to the symmetries of the Coulomb interaction.
Your calculation is correct, and it is an important feature of the second quantization that you obtain results like that. This is the same result that you would obtain without the second quantization, using correctly symmetrized wave functions:
$$\phi_{ij}(x_1,x_2) =\frac{1}{\sqrt{2}}\left[\phi_i(x_1)\phi_j(x_2)-\phi_i(x_2)\phi_j(x_1)\right].$$
