The Odds 3 (or More) Group Elements Commute 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}) $$
 A: This question is studied in the article Isoclinism classes and commutativity degrees of finite groups by P. Lescot. The probability that $n+1$ elements of $G$ pairwise commute is called there the $n$-th commutativity degree $d_n(G)$. A kind of recursive formula for $d_n(G)$ in terms of $d_{n-1}$ of centralizers is proved (Lemma 4.1). Lescot also proves that if $G$ is not abelian then $d_n(G) \leq \frac{3 \cdot 2^n - 1}{2^{2n+1}}$ with equality if and only if $G$ is isoclinic to the quaternion group $Q_8$.
A: There are $k(G)|G|$ commuting pairs $(y,z)$ of elements of $G.$ How many commuting triples $(x,y,z)$ of elements of $G$ are there? If we fix the first component $x$ of the triple, note that $x$ permutes (by conjugation) the commuting pairs of elements of $G,$ and the only such pairs it fixes are the commuting pairs which already have both components in $C_{G}(x).$ Hence $x$ fixes $k(C_{G}(x))|C_{G}(x)|$ 
commuting pairs, and the total number of commuting triples in $G \times G \times G$ is
$\sum_{x \in G} k(C_{G}(x))|C_{G}(x)|$. If $G$ has $k$ conjugacy classes, with representatives $\{ x_{i} : 1 \leq i \leq k \}$, this may be rewritten as 
$|G| \sum_{i=1}^{k} k(C_{G}(x_{i})).$ The probability you require, assuming a uniform distribution, is $\frac{1}{|G|^{2}}\left( \sum_{i=1}^{k} k(C_{G}(x_{i})) \right)$. 
The discussion (and Burnside's counting lemma) makes it clear that this is (number of orbits of $G$ by conjugation on commuting pairs of elements of $G$), divided by $|G|^{2}$.  
Later edit: Since I noticed the "(or more)" in the question title, the pattern is now clear:
let $c_{n}(G)$ denote the number of commuting (ordered) $n$-tuples of elements of $G$. Then we have $c_{n+1}(G) = \sum_{x \in G} c_{n}(C_{G}(x))$.
