Two particle operator Why is the two-particle (fermionic, cause for bosonic operators it is immediately clear that both representations are the same) Hamiltonian given by 
$$ H = \sum_{a,b,c,d} \langle ab|V|cd \rangle a_{a}^{\dagger}a_{b}^{\dagger}a_da_c$$
I would have rather assumed that the order is just the same as in the expectation value
$$ H = \sum_{a,b,c,d} \langle ab|V|cd \rangle a_{a}^{\dagger}a_{b}^{\dagger}a_ca_d.$$
Is there a simple argument why the first version is right and the second one not?
 A: $\renewcommand{\bbraket}[3]{\left\langle #1 | #2 | #3 \right\rangle}$

Is there a simple argument why the first version is right and the
  second one not?

The correct form for a two-body interaction in second quantization is always:
$$
V=\frac{1}{2}\int d^3r \int d^3r' V(|\vec r-\vec r'|)\hat \psi^\dagger(\vec r)\hat \psi^\dagger(\vec r')\hat \psi(\vec r')\hat \psi(\vec r)
$$
Note carefully the order of the $\vec r$s and the $\vec r'$s. This is the correct order. Any other order doesn't reduce to the correct first quantized form, which you can check by acting on a state of fixed particle number. I.e., check that it reduces to:
$$
\frac{1}{2}\sum_{i\neq j=1}^N V(|\vec r_i-\vec r_j|)\;,
$$
where $\vec r_i$ and $\vec r_j$ are position operators in first quantization.
So, now you plug in some orbital basis and corresponding operators like
$$
\hat \psi(\vec r)=\sum_a\phi_a(\vec r)\hat a_a
$$
and you end up with
$$
V=\sum_{a,b,c,d}a_a^\dagger a_b^\dagger a_ca_d\int d^3r d^3r' \phi_a^*(r)\phi^*_b(r')\phi_c(r')\phi_d(r)V(r-r')
$$
and what you mean by
$$
\bbraket{ij}{V}{kl}
$$
is
$$
\int d^3r d^3r' \phi_i^*(r)\phi^*_j(r')\phi_l(r')\phi_k(r)V(r-r')
$$
because you are writing the i and k terms first (in |ij> and |kl>) because they both are integrated over "r" and the j and l term second because they are both integrated over r'.
So, we have
$$
V=\sum_{a,b,c,d}a_a^\dagger a_b^\dagger a_ca_d \bbraket{ab}{V}{dc}
$$
A: Another way of showing this is to consider the definition of the ket 
$$|ab\rangle = a^{\dagger}_ab^{\dagger}_b |0\rangle$$
and the bra
$$\langle ab| = \langle 0 | a_b a_a$$
and to look at the matrix element of $H$ as it is defined.
$$ \langle ab | \hat{H} | cd \rangle = \langle ab | V | cd \rangle
 \times \langle 0 | a_b a_a a^{\dagger}_aa^{\dagger}_ba_da_ca_c^{\dagger}a_d^{\dagger} | 0 \rangle $$
(with three other similar terms that give the same result). Then simplify the expectation value of the creation/annihilation operators using the anti-commutation relation
$$ \{ a_a,a_b^{\dagger}\}=\delta_{ab}$$
to move all the creation operators to the right so that you end up with
$$ \langle 0 | a_da_ca_ba_a a_a^\dagger a_b^\dagger a_c^\dagger a_d^\dagger |0\rangle
= \langle abcd | abcd \rangle =1.$$
To do this, you need an even number of permutations. If you had used your second suggestion with $a^{\dagger}_aa^{\dagger}_ba_ca_d$, you would have needed an odd number of permutations and gotten an overall -1.
The other, more standard, technique is to move the annihilation operators to the right of the creation operators, since this will give zero. We have
\begin{align}
a_b a_a a^{\dagger}_aa^{\dagger}_ba_da_ca_c^{\dagger}a_d^{\dagger} 
&= a_a a^{\dagger}_a a_b  a^{\dagger}_ba_ca_c^{\dagger}a_da_d^{\dagger} \\
&= (1-a^{\dagger}_aa_a) (1-a^{\dagger}_ba_b)(1-a_c^{\dagger}a_c)(1-a_d^{\dagger}a_d)\\
&= 1
\end{align}
since we can see that any term involving the creation/annihilation operators will give zero when acting on the vacuum.
The essence of the matter is that since the bra is the Hermitian conjugate of the ket, it necessarily has the order of the creation/annihilation operators reversed, and so the form of the two-body operator reflects that.
