Some time ago I asked about the odds 2 group elements commute. I wonder about the odds that 3 group elements commute. Is there a "closed" formula for the sum

$$ \frac{1}{|G|^3} \sum_{g,h,k} \delta([g,h]=1)\delta([h,k]=1)\delta([k,g]=1)$$

One approach, as mentioned in Kefeng Liu, might be to use the "Heat Kernel" for finite groups.

$$ H(t,x,y) = \frac{1}{|G|} \sum_{\lambda \in \mathrm{Irr}(G)} \mathrm{dim}(\lambda) \chi_\lambda(xy^{-1}) e^{-t f(\lambda)}$$

If I'm not mistaken $f(\lambda)$ is the quadratic Casimir, but not sure. Really, for $t=0$ it reduces to the group theory identity:

$$ \delta(xy^{-1})= \frac{1}{|G|} \sum_{\lambda \in \mathrm{Irr}(G)} \chi_\lambda(x) \overline{\chi_\lambda(y)} = \frac{1}{|G|} \sum_{\lambda \in \mathrm{Irr}(G)} \chi_\lambda(1) \chi_\lambda(xy^{-1}) $$