For the potential term in the Hamiltonian, I understand that we go through the same process as with the kinetic energy term. On the RHS of the TDSE, we get something like $\frac{1}{2}\sum_{i}\sum_{j\neq i}\sum_{k}\sum_{l\neq k}\langle ij|V|kl\rangle C(\dots E_{i-1}E_kE_{i+1}\dots E_{j-1}E_lE_{j+1}\dots;t)$.
Reordering the argument list of $C(\dots)$ on the RHS incurs a phase factor of $(-1)^P$, where $P$ depends on the order of $i, j, k, l$.
Some examples: if their order of ascent is ...
$kilj$: $P = n_{k+1}+\dots+n_{i-1}+n_{l+1}+\dots+n_{j-1}$
$klij$: $P = n_{k+1}+\dots+n_{i-1}+n_{l+1}+\dots+n_{j-1}$
$likj$: $P=n_{i+1}+\dots+n_{k-1}+n_{l+1}+\dots+n_{j-1}+1$ [the extra $+1$ comes from the fact that $i$ and $j$ are positioned with opposite 'polarity' to $k$ and $l$, therefore the $E_k,E_l$ have to be swapped past each other to get into the correct positions]
As above, by my reckoning, $P(klij)=P(kilj)$, however if we calculate $(-1)^P$ by operating on $|\lbrace n_m\rbrace\rangle$ like $a_i^{\dagger} a_j^{\dagger} a_l a_k|\lbrace n_m\rbrace\rangle$, I get $(-1)^{S_k+S_{l}-1+S_{j}-2+S_{i}-2}|\dots n_k-1\dots n_l-1 \dots n_j+1\dots n_i+1\rangle \\ = (-1)^{S_k+S_{l}+S_{j}+S_{i}-1}|\dots n_k-1\dots n_l-1\dots n_j+1\dots n_i+1\dots\rangle$
if the order is $klij$ but
$(-1)^{S_k+S_{l}-1+S_{j}-2+S_{i}-1}|\dots n_k-1\dots n_l-1\dots n_j+1\dots n_i+1\dots\rangle\\= (-1)^{S_k+S_{l}+S_{j}+S_{i}}|\dots n_k-1\dots n_l-1\dots n_j+1\dots n_i+1\dots\rangle$
if the order is $kilj$. The order of $i$ and $l$ seems to matter when I calculate the phase factor by successive application of the $a, a^{\dagger}$ but not when I calculate it by reordering the arguments of the $C(\dots)$ coefficients. Where am I going wrong?