quantum mechanics - Fetter & Walecka's derivation of second quantised potential term in many-particle TDSE

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?

Answer

Worked out where I've been going wrong. Referring to the content of my question, in the $klij$ case, i.e. $k

$\begin{split} P(klij) &= n_{k+1} + \dots + n_{i-1} + \dots + n_{l+1} + \dots + n_{j-1} - 1 \\ &= S_i - S_k - n_k + S_j - S_l - n_l - 1 \end{split}$

and so the phase factor, $(-1)^{P(klij)} = (-1)^{S_i + S_k - n_k + S_j + S_l - n_l - 1}$, which, with the $\delta_{n_k0}\delta_{n_l0}$ on the RHS, gives $(-1)^{S_i+S_k+S_j+S_l-1}$. This is the same phase factor you get when you operate on the state like $a^{\dagger}_ia^{\dagger}_j a_la_k|\lbrace n_m\rbrace\rangle$ as I mentioned in the original question.