How many pairs of nilpotent, commuting matrices are there in $M_n(\mathbb{F}_q)$?

I can give a formula which shows that the answer is polynomial in $q$. Let $\lambda$ be a partition, with path $\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_a$. More specifically, let $$\begin{array}{lclclclc} \lambda_1 &=& \lambda_2 &=& \cdots &=& \lambda_{b_1} &\geq \\ \lambda_{b_1+1} &=& \lambda_{b_1+2} &=& \cdots &=& \lambda_{b_1+b_2} &\geq \\ \cdots \\ \lambda_{b_1+\cdots+b_{r-1}+1} &=& \lambda_{b_1+\cdots+b_{r-1}+1+2} &=& \cdots &=& \lambda_{b_1+\cdots+b_{r-1}+b_r} & \\ \end{array}$$ so $b_1+b_2+\cdots + b_r = a$. The number of pairs is $$\sum_{|\lambda|=n} q^{\binom{n}{2} + \sum_{j=1}^r \binom{b_j}{2}} \frac{\prod_{k=1}^{n} (q^k-1)}{\prod_{j=1}^r \prod_{i=1}^{b_j} (q^i-1)} = q^{\binom{n}{2}} \sum_{|\lambda|=n} q^{\sum_{j=1}^r \binom{b_j}{2}} \left[ \begin{matrix} a \\ b_1, b_2, \ldots, b_r \end{matrix} \right]_q \prod_{k=a+1}^n (q^k-1)$$ where the $\left[ \begin{matrix} a \\ b_1, b_2, \ldots, b_r \end{matrix} \right]_q$ is the $q$-deformation of the multinomial coefficient.

For example, here is the computation for $n=3$: $$\begin{array}{|l|l|l|c@{}c@{}c|} \hline \lambda & a & (b_1, \ldots, b_r) & q^{\binom{b_j}{2}} & \left[ \begin{matrix} a \\ b_1, b_2, \ldots, b_r \end{matrix} \right]_q & \prod_{k=a+1}^n (q^k-1) \\ \hline (3) & 1 & (1) & & & (q^2-1)(q^3-1) \\ (2,1) & 2 & (1,1) & & (q+1) & (q^3-1) \\ (1,1,1) & 3 & (3) & q^3 & & \\ \hline \end{array}$$ and, sure enough $$q^3 \left( (q^2-1)(q^3-1) \ + \ (q+1) (q^3-1) \ + \ q^3 \right) = q^8+q^7+q^6-q^5-q^4$$

Proof As you would expect, the partitions $\lambda$ are encoding Jordan normal forms. For a partition $\lambda$ of $n$, let $J(\lambda)$ be a $n \times n$ nilpotent matrix with Jordan blocks of size $\lambda$. Let $E(\lambda)$ be the set of $n \times n$ matrices which commute with $J(\lambda)$, so $E(\lambda)$ is a ring. Let $A(\lambda)$ denote the unit group in $E(\lambda)$ and let $N(\lambda)$ denote the nilpotents of $E(\lambda)$. Let $GL_n$ be the group of $n \times n$ invertible matrices over $\mathbb{F}_q$.

The number of $n \times n$ nilpotent matrices with Jordan type $\lambda$ is $\frac{|GL_n|}{|A(\lambda)|}$, since $GL_n$ acts transitively with stabilizer $A(\lambda)$. For any matrix in this orbit, the number of nilpotent matrices which commute with it is $|N(\lambda)|$. So we want to compute $$|GL_n| \sum_{|\lambda|=n} \frac{|N(\lambda)|}{|A(\lambda)|} = |GL_n| \sum_{|\lambda|=n} \frac{|N(\lambda)|/|E(\lambda)|}{|A(\lambda)|/|E(\lambda)|}$$ I like to think of the numerator and denominator of the fraction as the probability that a matrix in $E(\lambda)$ will be nilpotent or invertible (respectively).

It is now useful to use the ring structure on $E(\lambda)$. Let $M(\lambda)$ be the $\mathbb{F}_q[t]$ module which is $\mathbb{F}_q^n$ as a vector space, and where $t$ acts by $J(\lambda)$. The next paragraph is very similar to my answer here, where the role of $M(\lambda)$ is played by a finite abelian group $\bigoplus \mathbb{Z}/p^{\lambda_i}$.

The ring $E(\lambda) = \mathrm{End}_{\mathbb{F}_q[t]} M(\lambda)$. A matrix $X$ in $E(\lambda)$ is nilpotent (respectively invertible) if and only if $X$ acts nilpotently (respectively invertibly) on $M(\lambda) / t M(\lambda)$. As a vector space, $M(\lambda)/t M(\lambda)$ has dimension $a$, and $E(\lambda)$ acts on it by block upper-triangular matrices with block structure $(b_1, b_2, \ldots b_r)$. A block upper triangular matrix is nilpotent (resp. invertible) if and only if each of the blocks is nilpotent (resp. invertible), and each of the blocks is independent.

The probability that a $b \times b$ matrix is invertible if $\prod_{i=1}^b (1-q^{-i})$. From your earlier question, the probability that a $b \times b$ matrix is nilpotent is $q^{-b}$. Now plug in the formula for $|GL_n|$ and rearrange to get the above.

The answer to these kinds of questions is usually nicely formulated in terms of generating functions. So you are trying to get a formula for the formal sum $$\Lambda^{n,n}(z,q)=\sum_{d \geq 0} \frac{X^{n,n}_d(q)}{|GL_d(\mathbb{F}_q)|}z^d$$ where $X^{n,n}_d(q)$ stands for the number of pairs of commuting nilpotent matrices in $gl_d(\mathbb{F}_q)$. Here the exponents n,n in the notation is to recall that we are considering nilpotent matrices. One may consider also the problem of computing the number of pairs of commuting matrices in which the first is nilpotent, or the number of pairs of commuting matrices (in with no nilpotency condition). This leads to two variants

$$\Lambda^n(z,q)=\sum_{d \geq 0} \frac{X^n_d(q)}{|GL_d(\mathbb{F}_q)|}z^d$$

$$\Lambda(z,q)=\sum_{d \geq 0} \frac{X_d(q)}{|GL_d(\mathbb{F}_q)|}z^d$$

with obvious notations.

In this last case (pairs of arbitrary commuting matrices, the answer was given (in these terms) by Feit and Fine in the paper :

W. Feit, N. J. Fine, Pairs of commuting matrices over a finite field, Duke Math. J 27, 1960 91--94.

In terms of plethystic exponentials, the answers are as follows :

$$\Lambda(z,q)=Exp\left(\sum_{d \geq 1} \frac{q^2z}{(q-1)(1-z)}\right)$$ $$\Lambda^n(z,q)=Exp\left(\sum_{d \geq 1} \frac{qz}{(q-1)(1-z)}\right)$$ $$\Lambda^{n,n}(z,q)=Exp\left(\sum_{d \geq 1} \frac{z}{(q-1)(1-z)}\right)$$

The plethystic exponentials are defined as follows : for any formal power series $f(q,z)$ in q and z with rational coefficients, set $$Exp(f(q,z))=exp\left(\sum_{l \geq 1} \frac{1}{l} f(q^l,z^l)\right).$$

The above formulas may be proved using a direct argument as explained above by David, (together with some clever combinatorics). There are also so more fancy ways of proving these, which can be generalized to the setting of an arbitrary quiver (the present case corresponds to the case of the Jordan quiver); in that case the answer is given in terms of the so-called `Kac polynomials' of the quiver. When no nilpotency condition is required, this is done by Mozgovoy in

S. Mozgovoy, Motivic Donaldson-Thomas invariants and McKay correspondence, arXiv:1107.6044 (2011)

When one imposes a 'half' nilpotency condition, this is done in my paper

O. Schiffmann, On the number of points of the Lusztig nilpotent variety over a finite field, arXiv:1212.3772

(this paper is in the process of being revised and improved -- there are some typos inside).