Solution (source code)

= Solution

At a Hartree-Fock stationary point the reference-single terms vanish. Expanding the normalized expectation through second order gives
$$
\mathcal E^{(2)}
=\sum_{ia,jb}\lambda_i^a\lambda_j^b
\left[
\langle\Phi_i^a|(\hat H-E_0)|\Phi_j^b\rangle
+\langle\Phi_0|\hat H|\Phi_{ij}^{ab}\rangle
\right].
$$
The reference-double term arises twice in the numerator and cancels the factor $1/2$ in the <wavefunction>. The normalization subtracts $E_0\sum_{ia}(\lambda_i^a)^2$.

Using <Slater-Condon rules> in canonical orbitals,
$$
\langle\Phi_i^a|(\hat H-E_0)|\Phi_j^b\rangle
=(\gamma_a-\gamma_i)\delta_{ij}\delta_{ab}
+\langle aj|ib\rangle-\langle aj|bi\rangle,
$$
and $\langle\Phi_0|\hat H|\Phi_{ij}^{ab}\rangle
=\langle ab|ij\rangle-\langle ab|ji\rangle$. Combining them gives the requested $X_{ij}^{ab}$. The orbital energies are
$$
\gamma_p=h_{pp}+
\sum_{k\ {\rm occ}}\left(\langle pk|pk\rangle-\langle pk|kp\rangle\right),
$$
for either occupied or virtual $p$. Diagonalizing the occupied and virtual Fock blocks is the suitable orbital choice that makes the delta form valid for coincident as well as distinct indices.