General form of two particle operators is $$\hat{B} = \sum\limits_{ij,kl} B_{ij,kl} \hat{a}^{\dagger}_i\hat{a}^{\dagger}_j\hat{a}_k\hat{a}_l$$
and I am asked to compute the expected value of $\hat{B^2}$. Some algebra gives us:
$$\hat{B^2} = \sum\limits_{ij,kl, ab, cd} B_{ij,kl}B_{ab, cd} \hat{a}^{\dagger}_i\hat{a}^{\dagger}_j\hat{a}_k\hat{a}_l\hat{a}^{\dagger}_a\hat{a}^{\dagger}_b\hat{a}_c\hat{a}_d$$
When computing $\langle\{ n_k \} |\hat{B}^2 | \{n_k\} \rangle$ one gets nonzero terms only if the number of particles in every state doesn't change. Thus for every annihilation operator $\hat{a}$ we must have an creation operator $\hat{a}^\dagger$ acting on the same index.
The problem boils down to pairing the set $\{i,j,a,b\}$ with $\{k,l,c,d\}$ and splitting the sum over $N^8$ to subsets of that space, rearrange the operators to get the number operators $\hat{n}_i := \hat{a}_i^\dagger \hat{a}_i$ at which point such "subset" of $B^2$ acting on $|\{n_k\}\rangle$ will give us products of corresponding occupation numbers.
Yet, there are $4!$ pairs to be made, and we must also counts things like $i =j = a = k = l = c \neq b =d$ making it too hard to "brute force" by hand.
Since we are effectively looking at permutations there must exist a smart way of going on about this! Thus my question is: how does one efficiently solve this? To make matters worse we have two cases, one for bosons where $B_{ij,kl} = B_{ji,kl} = B_{ji,lk} =B_{ij,lk}$ and the operators commute, second for fermions with $B_{ij,kl} = -B_{ji,kl} = B_{ji,lk} =-B_{ij,lk}$ and the operators anticommute.