I happened to come across the following term while doing an excercise on perturbation theory
\begin{equation} H^2=J^2\sum_{<i,j>}\sum_{<k,l>}\sigma_i\ \sigma_j\ \sigma_k\ \sigma_l \end{equation}
Where $H$ is the Ising Hamiltonian, $J$ is the coupling constant, $\sigma_I=\pm1$ and $<i,j>$ means sum over neighbouring sites.
When computing the mean value $\langle H^2 \rangle$
\begin{equation} \langle H^2\rangle=J^2\frac{1}{Z}\sum_{\{\sigma\}}e^{-\beta H}\sum_{<i,j>}\sum_{<k,l>}\sigma_i\ \sigma_j\ \sigma_k\ \sigma_l \end{equation}
I started to wonder if every term in this sum is non-zero. for example, terms where $i=k, j=l$ give a contribution of $2J^2N$ because $\sigma_i \sigma_j \sigma_k \sigma_l=1$ and one is left out with a sum over pairs of sites. Using the fact that $\sum_{\sigma=-1}^1\sigma^2=2$ and $\sum_{\sigma=-1}^1\sigma=0$ is there any argument to determine if some combinations of $i,j,k,l$ don't contribute to the mean value?
Edit: I would be really interested in knowing, for example, how much this term is
\begin{equation} \langle H^2\rangle_{i=j}= J^2\frac{1}{Z}\sum_{\{\sigma\}}e^{-\beta H}\sum_{<i,j>}\sum_{<k,l>}\sigma_i\ \sigma_j\ \sigma_k\ \sigma_l\ \delta_{ik} \end{equation}
