Mismatch in the subscripts of two-body operators in occupation number formalism I'm reading Prof. Mattuck's "A Guide to Feynman Diagram" and came across this rather simple deviation of the "two-body" operator in the occupation number formalism. How the author labelled the subscripts of creation/destruction operators makes me puzzled a bit. The part that I found the most puzzling is listed here:

it can be shown that the "two-body" operator
  \begin{equation}\mathcal{O}=\frac{1}{2} \sum_{\substack{i, j=1\\(i \neq j)}}^{N} \mathcal{O}\left(\mathbf{r}_{l}, \mathbf{p}_{i}, \mathbf{r}_{j}, \mathbf{p}_{j}\right)\end{equation}
  For instance the interaction potential
  \begin{equation}V\left(\mathbf{r}_{1}, \dots, \mathbf{r}_{N}\right)=\frac{1}{2} \sum_{\substack{i, j=1\\ (i \neq j)}} V\left(\mathbf{r}_{i}-\mathbf{r}_{j}\right)\end{equation}
  becomes
  \begin{equation}\mathcal{O}^{\mathrm{occ}}=\frac{1}{2} \sum_{k l m n} \mathcal{O}_{klmn} c_{l}^{\dagger} c_{k}^{\dagger} c_{m} c_{n}\quad (1)\end{equation}
  where
  \begin{equation}\mathcal{O}_{klmn}=\int d^{3} \mathbf{r} \int d^{3} \mathbf{r}^{\prime} \phi_{k}^{*}(\mathbf{r}) \phi_{i}^{*}\left(\mathbf{r}^{\prime}\right) \mathcal{O}\left(\mathbf{r}, \mathbf{r}^{\prime} ; \mathbf{p}, \mathbf{p}^{\prime}\right) \phi_{m}(\mathbf{r}) \phi_{n}\left(\mathbf{r}^{\prime}\right)\end{equation}

It is obvious that the transition amplitude ($\mathcal{O}_{klmn}$) in (1) is related to the Bhabha scattering diagram by labeling the momenta using the rule of "left out-right out-left in-right in":

Now my question is, why don't we write the operator as $\sum_{k l m n} \mathcal{O}_{klmn} c_{k}^{\dagger} c_{l}^{\dagger} c_{m} c_{n}$ instead of $\sum_{k l m n} \mathcal{O}_{klmn} c_{l}^{\dagger} c_{k}^{\dagger} c_{m} c_{n}$? According to the "anti-commutation" rule we know that the latter one can be obtained from the former one by adding a (-1) factor, but what would be missing if we use the former one at the first place?
 A: The important point here is that the order of indices in the matrix element is not the same as in the operator product.
Indeed, if we have fields represented by
$$\hat{\psi}(x) = \sum_n c_n\phi_n(x),$$
Then the Coulomb term is written as
$$\hat{V} = \frac{1}{2}\int dx dx'\hat{\psi}^\dagger(x)\hat{\psi}^\dagger(x')v(x-x')\hat{\psi}(x')\hat{\psi}(x) =  
\frac{1}{2}\sum_{k, l, m,n}\langle k, l|v|m, n\rangle c_k^\dagger c_l^\dagger c_n c_m,$$ that is 
$$\hat{\psi}^\dagger(x)\hat{\psi}^\dagger(x')\hat{\psi}(x')\hat{\psi}(x) =
\hat{\psi}^\dagger(x)\hat{\psi}(x)\hat{\psi}^\dagger(x')\hat{\psi}(x') = \hat{n}(x)\hat{n}(x'),$$
so that the product of the charge densities has the sign corresponding to the repulsive interaction, whereas the order of indices in the matrix element is the same for its bra and ket vectors, as it should be:
$$ \langle k, l|v|m, n\rangle = \frac{1}{2}\int dx dx'\phi_k^*(x)\phi_l^*(x')v(x-x')\phi_m(x')\phi_n(x).$$
A: Thanks @Vadim for making a good start. I think this is just a matter of making things streamlined. Using the Dirac notation, we have, for a state in the occupation number formalism:
$$
\langle n_1,n_2,\ldots,n_i\ldots|=\overline{| n_1,n_2,\ldots,n_i\rangle}
$$
where the overline means the Hermitian adjoint. However, for a product of operators like the ones in @Vadim's answer, we have:
$$
(\hat{\phi}_k\hat{\phi}_l)^{\dagger}=\hat{\phi}_l^{\dagger}\hat{\phi}_k^{\dagger}
$$
Thus, the second equation in @Vadim's answer becomes, in Dirac's notation:
$$
\begin{aligned}
\hat{V}&=\frac{1}{2}\sum_{klmn}\langle k,l|V|m,n\rangle=\frac{1}{2}\sum_{klmn}\int\int d\mathbf{r}d\mathbf{r}^{\prime}(\hat{\phi}_k\hat{\phi}_l)^{\dagger}V(\mathbf{r}-\mathbf{r}^{\prime})\hat{\phi}_m\hat{\phi}_n\\
&=\frac{1}{2}\sum_{klmn}\int\int d\mathbf{r}d\mathbf{r}^{\prime}\hat{\phi}_l^{\dagger}\hat{\phi}_k^{\dagger}V(\mathbf{r}-\mathbf{r}^{\prime})\hat{\phi}_m\hat{\phi}_n=\frac{1}{2}\sum_{klmn}V_{klmn}c^{\dagger}_lc^{\dagger}_kc_mc_n
\end{aligned}
$$
where
$$
V_{klmn}=\int\int d\mathbf{r}d\mathbf{r}^{\prime}\phi^{\dagger}_k\phi^{\dagger}_lV(\mathbf{r}-\mathbf{r}^{\prime})\phi_m\phi_n
$$
Thus we don't have mismatches in Dirac notation to keep things intuitive, but we do need to take care of the effects of " ${}^{\dagger}$ " when we translate the notation into integrals.
