Problem using spin-restricted form of the second-quantized nonrelativistic Hamiltonian
I have a problem that confuses me a lot. The two-electron part of the electronic nonrelativistic Hamiltonian can be written \begin{equation} \frac{1}{2}\sum_{pqrs} (pq|rs) [a^\dagger_{p\alpha}a^\dagger_{r\alpha}a_{s\alpha}a_{q\alpha} + a^\dagger_{p\alpha}a^\dagger_{r\beta}a_{s\beta}a_{q\alpha} + a^\dagger_{p\beta}a^\dagger_{r\alpha}a_{s\alpha}a_{q\beta} + a^\dagger_{p\beta}a^\dagger_{r\beta}a_{s\beta}a_{q\beta}], \end{equation} where the two-electron integral is defined via \begin{equation} (pq|rs) = \iint p^(\mathbf r_1) q(\mathbf r_1) \frac{1}{r_{12}} r^(\mathbf r_2) s(\mathbf r_2) d\mathbf r_1 d\mathbf r_2. \end{equation} I will assume that the single particle wavefunctions (orbitals) are a complete set of eigenfunctions of an effective one-electron Hamiltonian (for example the Fock operator of Hartree-Fock theory) for the hydrogen molecule, written as the product of an orbital part and an eigenfunction of $S_z$ ($\alpha$ or $\beta$).
Then the ground state of the hydrogen molecule can approximately be written as \begin{equation} |\Psi_0\rangle = a^\dagger_{1\alpha}a^\dagger_{1\beta}|\mathrm{vac}\rangle, \end{equation} where 1 is the index of the lowest-energy orbital. So far for the introduction. 656
My problem arises when I consider an alternative form of the two-electron part of the Hamiltonian, \begin{equation} a^\dagger_{p\alpha}a^\dagger_{r\alpha}a_{s\alpha}a_{q\alpha} + a^\dagger_{p\alpha}a^\dagger_{r\beta}a_{s\beta}a_{q\alpha} + a^\dagger_{p\beta}a^\dagger_{r\alpha}a_{s\alpha}a_{q\beta} + a^\dagger_{p\beta}a^\dagger_{r\beta}a_{s\beta}a_{q\beta} = E_{pq}E_{rs} - \delta_{qr}E_{ps}, \end{equation} where the so-called singlet excitation operators are defined via \begin{equation} E_{pq} = a^\dagger_{p\alpha}a_{q\alpha} + a^\dagger_{p\beta}a_{q\beta}. \end{equation} When I calculate the expectation value of the two-electron part in the molecule's ground state with the first formula, I get the result $(11|11)$, which I know is correct. For the form using singlet excitation operators, however, I obtain \begin{equation} \langle \Psi_0 | E_{pq} E_{rs} | \Psi_0\rangle = 2\delta_{qr} \delta_{s1} \delta_{p1} \end{equation} and \begin{equation} \langle \Psi_0 | E_{ps} | \Psi_0\rangle = 2 \delta_{s1} \delta_{p1}. \end{equation} When combining these two parts, one sees that they exactly cancel, meaning that the total expectation value would vanish. Can anyone see where the error is? I have obviously become blind for it, I have calculated it several times forwards and backwards and always get to the same results. 402