The number of commuting m-tuples is divisible by order of group: Improvements? The number of commuting pairs of elements in finite group G is equal to the product $k(G)*|G|$ (see MO271757 ) where $k(G)$ is the number of conjugacy classes. Thus it is is divisible by $|G|$ (the number of elements of $G$). 
That divisibility also follows from a theorem by L. Solomon, stated below.
Question 0: It seems that the number of commuting $m$-tuples $c_m(G)$ is also divisible by $|G|$ for any $m$; is this correct? It appears to follow from a result cited in  Klyachko, Mkrtchyan (details below).
Question 1: Does the ratio $c_m(G)/|G|$ have some group theoretic interpretation for $m>2$?  (When $m=2$, this is the number of conjugacy classes).
If the group is abelian, then obviously $c_m(G) = |G|^m$, so it is divisible by a very high power of $|G|$. 
Question 2: Is there any improvement possible for this type of divisibility by $|G|$ for  nilpotent or $p$-groups ? 
Remark: Any improvement cannot contradict the analogues of the 5/8
bound for general $m$: $c_{m+1}(G) \leq \frac{3 \cdot 2^m - 1}{2^{2m+1}} |G|^{m+1}$ by  Lescot (see MO108392). 

Reminder 
Let me state theorems cited in Klyachko & Mkrtchyan 2012 (found in MO98639), and apply it to our situation.


*

*Solomon theorem [1969]. In any group, the number of solutions to a system of coefficient-free equations is divisible
by the order of this group if the number of equations is less than the number of unknowns.

*Application Consider the equation  $xy=yx$ in a group---one equation, two unknowns---thus the number of solutions should be divisible by the order of the group. Hence the number of commuting pairs is divisible by the order of group.
Note that we cannot apply that theorem for the number of commuting triples, $m$-tuples, since the number of equations exceeds the number of unknowns.
So one needs a refinement, and that seemed to be known:


*

*Gordon,Rodriguez-Villegas theorem arXiv:1105.6066. In any group, the number of solutions of a system of coefficient free
equations is divisible by the order of this group if the rank of the matrix composed of exponent sums of $i$th
unknown in $j$th equation is less than the number of unknowns.


(It is presented like this by Klyachko & Mkrtchyan 2012, it is not immediately clear (to me) how to extract this formulation from the original paper). 


*

*Application Consider equations defining commuting $m$-tuples: $x_ix_jx_i^{-1}x_j^{-1} = 1$. The sums the of exponents is ZERO! Thus the rank of the matrix is zero and hence the theorem ensures that the number of solutions is divisible by order of $G$.


I hope that this correct, and that an expert can confirm it.

Motivation
I hope (see MO271752) that there should be a nice generating function,
 $$ \sum_{n\geq 0} \frac{|\mathrm {
commuting~} m\mathrm{-tuples~  in~ GL}(n,F_q)|}{|GL(n,F_q)|}  x^n = ??? $$
This is similar to the known,
  $$ \sum_{n\geq 0} \frac{|\mathrm {
commuting~} \mathrm{pairs~  in~ GL}(n,F_q)|}{|GL(n,F_q)|} x^n = \prod_{j\geq 1}
    \frac{1-x^j}{1-qx^j} $$
which would imply  divisibility results at least for 
$GL(n,F_q)$ (and some other groups too), so it is 
nice to have  support for such a belief. 
About other groups, it would be  quite interesting for me to know especially
about the group $UT(n,q)$ (unitriangular matrices over a finite field)---
by what power of $q$ the numbers $c_m(G)$
are divisible.
 A: As Alexander requested, I am writing some details of my comment.

Gordon–Rodriguez-Villegas Theorem (the original form, see J.Algebra 2012 or arXiv 2011).
  If $\Gamma$ is a finitely generated group and the index of its commutator subgroup is infinite, then the number of homomorphisms $\Gamma\to G$ is divisible by $|G|$ for any group $G$.

You may or may not assume that $G$ is finite. In the latter case, the divisibility is in the sense of cardinal arithmetics (an infinite cardinal is divisible by all less or equal nonzero cardinals).

Gordon–Rodriguez-Villegas Theorem (an ``equational form'', see [KM2012]).
  Suppose that 
  $$
\{w_1(x_1,\dots,x_n)=1,\ w_2(x_1,\dots,x_n)=1,\dots\}
$$
  is a system of coefficient-free equations over a group $G$, i.e. $w_i$ are elements of the free group $F(x_1,\dots,x_n)$ (the number of equations may be infinite, but the number of unknowns is finite). 
  Then the number of solutions to this system is divisible by $|G|$ whenever the rank of the matrix $A$ composed of the exponent sums of $i$th unknown in $j$th equation is less than the number of unknowns.

Clearely, this is an equivalent statement, because


*

*the number of homomorphisms
$$
\Gamma=\langle x_1,\dots,x_n\ |\ w_1=1,\ w_2=1,\dots\rangle\to G
$$
is equal to the number of solutions to the system of equations
$\{w_1(x_1,\dots,x_n)=1,\ w_2(x_1,\dots,x_n)=1,\dots\}$
in $G$;

*the index of the commutator subgroup of $\Gamma$ is infinite if and only
if $\mbox{rank}\, A<n$.  
Thus, Alexander is right — the the matrix $A$ in his case is zero.
The divisibility part in Qiaochu Yuan's Claim is an immediate corollary of 
the Gordon–Rodriguez-Villegas Theorem (as the commutator subgroup of 
$\mathbb Z\times (\mbox{something})$ is of infinite index).
As for the ratio
$$
\frac{|\{\Gamma\to G\}|}{|G|}=\frac{\mbox{the number of solutions to the system above}}{|G|},
$$
it has a natural interpretation in the general case, but perhaps it was not stated explicidly anywhere. Anyway, everybody can obtain such an interpretation from a short and quite elementary proof of a generalisation of the GRV theorem that can be found in [this self-advertisement] (which contains also some other surprising divisibility corollaries).
A: Another strategy is to use induction on $m,$ there being nothing to do if $m=1.$ Note that $(x_{1},x_{2},\ldots,x_{m})$ is a commuting $m$-tuple if and only if $(x_{2},\ldots,x_{m})$ is a commuting $(m-1)$-tuple of $C_{G}(x_{1}).$ By induction then, for every choice of $x_{1},$ the number of completions of $x_{1}$ to a commuting $m$-tuple with first component $x_{1}$ is divisible by $|C_{G}(x_{1})|.$ It follows easily that given a conjugacy class $C$ of $G,$ the number of commuting $m$-tuples whose first component is an element of $C$ is divisible by $|G|.$ Since the conjugacy classes partition $G,$ the number of commuting $m$-tuples is indeed divisible by $|G|.$
Notice also that if $p$ is a prime which does not divide $|G|$, then the number of mutually commuting $p$-tuples of elements of $G$ is congruent to $|G|$ (mod $p$), since the cyclic permutation of order $p$ acts on such commuting $p$-tuples with $|G|$ fixed points (the ordered $p$-tuples with all components equal).
A: The answer to questions 0 and 1 is yes. Here is a generalization.

Claim: Let $\pi$ be a finitely generated group and $G$ be a finite group. Then
$$\frac{|\text{Hom}(\pi \times \mathbb{Z}, G)|}{|G|}$$
is equal to the number of conjugacy classes of homomorphisms $\pi \to G$.

We get an answer to questions 0 and 1 by setting $\pi = \mathbb{Z}^n$: $\frac{c_{n+1}(G)}{|G|}$ is equal to the number of conjugacy classes of homomorphisms $\mathbb{Z}^n \to G$.
A proof can be given along the lines of this blog post. The idea is to consider the mapping groupoid $X = [B \pi, BG]$ of homomorphisms $\pi \to G$, then to take its "loop groupoid"
$$LX = [B \mathbb{Z}, X] \cong [B(\pi \times \mathbb{Z}), BG].$$
It's an easy computation to verify that the groupoid cardinality of a loop groupoid $LX$ is equal to the number of connected components (isomorphism classes) of $X$. On the other hand, $LX$ is the quotient groupoid obtained from the conjugation action of $G$ on $\text{Hom}(\pi \times \mathbb{Z}, G)$, and hence its groupoid cardinality is also the above expression.

Edit: It is maybe worth pointing out the following connection to TQFT, which makes the above proof somewhat more natural than just using Burnside's lemma. In short, in $n$-dimensions the (untwisted) Dijkgraaf-Witten TQFT $Z_G$ with finite gauge group $G$ assigns to

*

*a closed $n$-manifold $X$ the number $\frac{\text{Hom}(\pi_1(X), G)|}{|G|}$ (the "volume" of the "moduli space" of $G$-bundles on $X$), and to

*a closed $(n-1)$-manifold $X$ the vector space of functions on the set of isomorphism classes of $G$-bundles on $X$.

It's a general feature of any TQFT that if $X$ is $(n-1)$-dimensional, then
$$Z(X \times S^1) = \dim Z(X).$$
Applying this to Dijkgraaf-Witten theory recovers a special case of the above result for $\pi = \pi_1(X)$. In particular the original desired result about $\pi = \mathbb{Z}^n$ can be seen as coming from Dijkgraaf-Witten theory as applied to tori.
A: As @AlexanderChervov requested, I am writing the details of the elementary argument for Qiaochu's answer as a community wiki.  First note that $\hom(\pi\times \mathbb Z, G)$ is in bijection with the set $X$ of pairs $(\varphi,g)$ where $\varphi\in \hom(\pi,G)$ and $g\in G$ with $g\varphi g^{-1}=\varphi$.  Viewing $\pi\times \mathbb Z$ as an internal direct product, the bijection takes $\psi\colon \pi\times \mathbb Z\to G$ to $(\psi|_\pi,\psi(1))$ and the inverse takes $(\varphi,g)$ to $\psi(x,k) = \varphi(x)g^k$.  
The group $G$ acts on $\hom(\pi, G)$ by conjugation and the orbits are the conjugacy classes of homomorphisms.  We can count $|X|$ as $$|\hom(\pi\times \mathbb Z,G)|=|X|=\sum_{g\in G}|\{\varphi\in \hom(\pi,G)\mid g\varphi g^{-1}=\varphi\}=\sum_{g\in G}|\mathrm{Fix}(g)|$$ and so by the Cauchy-Frobenius-Burnside orbit counting lemma, we get that $$\frac{|\hom(\pi\times \mathbb Z,G)|}{|G|}$$ is the number of conjugacy classes of homomorphisms $\pi\to G$.
