Number of $3\times 3$ anticommuting matrices over finite fields $\mathbb{F}_p$ is (or is not?) polynomial in $p$? There are rare algebraic varieties such that the number of points over finite fields $\mathbb{F}_p$ is given by a polynomial in $p$. One notable series of examples is the commuting variety: $[A,B]=0$ of $n\times n$ matrices $A,B$ over finite field. The computation was obtained by Feit and Fine in 1960 and many generalizations have been obtained recently. But it seems results in natural generality are not yet achieved (see below).
Question 1: Consider pairs of $3\times 3$ anticommuting $AB+BA = 0 $  over finite fields $\mathbb{F}_p$, is true that their number is polynomial in $p$, for $p>2$ ? May be one needs to exclude some other primes, not only  $p=2$ or $p=2$ is the only exception ?
Question 2: If the number of points is indeed given by a polynomial ($p>2$),  then it is given by the polynomial found by Roland Bacher and Peter Taylor in MO 404760:
$$2p^{10}+7p^9-3p^8-6p^7-4p^6+3p^5+4p^4-2p^3$$
(My direct calculation yields 221157 and 31511625 matrices for $p=3,5$ respectively,
and colleagues found that it is the only polynomial of degree 10 which satisfies these conditions and has minimal possible coefficients. Heuristic to search for polynomials with the smallest possible coefficients works quite fine in my experience for such questions.)
Question 3: Bonus question. It might be count is polynomial for any $n$ and $n$x$n$ anticommuting $[A,B]=0$, and there is nice generating function over $n$ for such polynomials - similar to Feit,Fine result for commuting matrices ? (Well, it might be better to leave it for separate question).
======================
PS
For 2x2 case the polynomial seems to be given by $+p^5+3p^4-2p^3-2p^2+p$ for $p>2$, checked till $p=19$.
Anticommuting variety has been studied recently e.g.: Anti-commuting varieties.
======================
Context
We might expect similar results not only for pairs, but triples, n-tuples of commuting/anticommuting matrices: MO271752, but commuting/anticommuting is just an example, it should be true for much wider class of algebras - "categorical exponential formula" might be the right context for such questions MO272045, MO275524, however presently known forms seems does not cover even Feit,Fine case. See also very nice results and connections with Hasse-Weil zeta function in  Yifeng Huang 2021. The general question about varieties which are polynomial count seems also not so simple as disproof of Kontsevich conjecture indicates.
Such varieties thought to be defined over the mysterious "field with one element".
 A: Here is an explicit formula. Let $C_n(q)$ be the number of pairs of anticommuting matrices over $F_q$ of size $n$. Let
$$
\phi(q,T)=\sum_{n=0}^\infty \frac{T^n}{\prod_{i=1}^n(1-q^{-i})}.
$$
Then
$$
\sum_{n=0}^\infty \frac{C_n(q) T^n}{q^{n^2}\prod_{i=1}^n(1-q^{-i})} = \prod_{n=1}^\infty \frac{\phi(q,T^n)}{(1-q T^{2n})(1-T^{2n-1})}.
$$
The calculation is a bit tedious. Some of the ideas:
We have $|GL_n(F_q)|=q^{n^2}\prod_{i=1}^n(1-q^{-i})$, so when we divide $C_n(q)$ by this number we are counting pairs $A,B$ up to conjugation with weight $\frac{1}{|\text{automorphisms(A,B)}|}$.
The generating function splits as a product of the corresponding functions for $3$ cases:

*

*$A$ invertible, $B$ invertible. Here the computation is reduced to the number of conjugacy classes $C$ in $GL_n$ satisfying $-C=C$. The answer is
$$
\prod_{n=1}^\infty \frac{1-T^{2n}}{1-q T^{2n}}.
$$


*$A$ invertible, $B$ nilpotent. It is clear that $-B$ is conjugate to $B$, so the corresponding generating function counts nilpotent conjugacy classes
$$
\prod_{n=1}^\infty \frac{1}{1-T^{n}}.
$$


*$A$ nilpotent, $B$ arbitrary. Here we compute the dimension of the kernel of the operator $B\to BA+AB$, obtain $\sum_{i} \lambda_i^2$ where $\lambda$ is the conjugate partition to the partition giving the sizes of Jordan blocks of $A$. The result is
$$
\prod_{n=1}^\infty \phi(q,T^n).
$$
Finally, we multiply all the generating functions together. and obtain the claimed result.
By the way, using the fact that the number of nilpotent matrices is $q^{n^2-n}$ one can check that the counting function for pairs nilpotent, nilpotent is given by
$$
F_3(q,q^{-1} T),
$$
where $F_3$ is the function for the case 3. So we obtain
$$
F_3(q,T)=F_3(q,q^{-1} T) \prod_{n=1}^\infty \frac{1}{1-T^{n}},
$$
from which it is possible to deduce our formula for $F_3$ using
$$
\phi(q,T)=\phi(q,q^{-1} T) \frac{1}{1-T}.
$$
But the formula for $F_3$ can be also obtained by a direct computation.
A: I think one can run a similar argument as YCor's to show the number of anticommuting pairs of $n\times n$ matrices over $\mathbb F_q$ is polynomial in $q$ for $q$ odd, for any $n$. This argument is not as detailed as I would like, but I believe they can all be filled in.
To do this, note that $AB + BA =0 $ is a system of linear equations in $B$. So we're looking for $\sum_a q^{ \dim \{ B \mid AB+BA= 0\}}$.
Now the operator $B \mapsto AB + BA$ can be written as the operator from a vector space $V \otimes V^{\vee}$ to itself given by $A \otimes I + I \otimes A^T$. Expressed this way, it's clear that its kernel is invariant under conjugating $A$ and $A^T$ separately. Noting that $A^T$ is conjugate to $A$, we may assume that $A$ and $A^T$ are both in Jordan normal form, and it's the same Jordan normal form. (Here we use that the dimension of the kernel may be calculated over a larger field.)
The kernel now may be derived from a standard computation of the kernel of $A_1 \otimes I + I \otimes A_2$, which may be done one Jordan block at a time. It is the sum over ordered pairs of a Jordan block of $A$ of size $n_1$ and eigenvalue $\lambda_1$ and a Jordan block fo $A$ of size $n_2$ and eigenvalue $\lambda_2$ of $$\begin{cases} n_1 +n_2 -1 & \lambda_2 =-\lambda_1 \\ 0 & \textrm{otherwise}\end{cases}$$
This depends only on the following features of $A$: The number of distinct prime factors of its characteristic polynomial, the degree of each prime factor, for each prime factor, the number of associated Jordan canonical form / rational canonical form blocks of each size, which prime factors are related to which other prime factors by the relation $f(x) = \pm f(-x)$, and which prime factors satisfy $f(x) = f(-x)$. Indeed, given such data, we know how many eigenvalues it has, for each eigenvalue, the sizes of its Jordan blocks, and for each eigenvalue, which other eigenvalue, if any, is its negative, and this can be used to calculate $\sum_{\lambda_1,\lambda_2, n_1,n_2}\begin{cases} n_1 +n_2 -1 & \lambda_2 =-\lambda_1 \\ 0 & \textrm{otherwise}\end{cases}$.
So it suffices to show the number of matrices consistent with one fixed list of the number of prime factors, the number of Jordan blocks of each size for each prime factor, etc. is a polynomial in $q$. This data determines for similar reasons the size of the centralizer of $A$, which is a polynomial in $q$ dividing the order of $GL_n(\mathbb F_q)$, and thus the size of the conjugacy class, so it suffices to show the number of conjugacy classes consistent with this data is polynomial in $q$.
(Why does the size of the centralizer divide the size of $GL_n$? Because the centralizer looks like an extension of a unipotent group by a product of things of the form $GL_a(\mathbb F_{q^b})$ where the sum of $ab$ is $n$. The order of $GL_a(\mathbb F_{q^b})$ divides the order of $GL_{ab}(\mathbb F_{q})$ by matching up the factors, and then the product of these divides $GL_n$ by the usual flag variety calculation. The unipotent part gives a power of $q$, but if the power of $q$ doesn't match up as abstract polynomials, it also won't match up for integers, because the other factors never accidentally give us powers of $q$, so that's ok also.)
Choosing the prime factors one at a time, we see that it suffices that the number of primes of degree $d$ over $\mathbb F_q$ is polynomial in $q$, as is the number of primes satisfying $f(x) =\pm f(-x)$ and $f(x) \neq \pm f(-x)$, because then the number of choices for each prime is given by a polynomial minus a constant representing forbidden values. The polynomiality of the primes is standard, so it suffices to show polynomiality of primes satisfying $f(x) = \pm f(-x)$.
In degree $1$ there is one such prime, in every other odd degree there are none, and in even degree $2n$, every such prime arises as $g(x^2)$ for $g$ prime of degree $n$. However, not all $g(x^2)$ are irreducible - they can also be $h(x) h(-x)$ for $h$ of degree $n$ irreducible not equal to $\pm h(-x)$. However, because the number of such $h$ is polynomial by induction on $n$, we only need to subtract a polynomial term from the number of irreducible $g$ of degree $n$ to get the answer, giving a polynomial, as desired.
A: If you were counting matrices $A$, $B$ in $GL_n$ with $AB=-BA$, this could be attacked using the methods from
Hausel, Tamás; Rodriguez-Villegas, Fernando, Mixed Hodge polynomials of character varieties. With an appendix by Nicholas M. Katz. (arXiv version), Invent. Math. 174, No. 3, 555-624 (2008). ZBL1213.14020 .
First of all, taking determinants of both sides, notice that $n$ must be even. We can rewrite the relation as $ABA^{-1} B^{-1} = - \mathrm{Id}_n$. Hausel and Rodriguez-Villegas spend most of the paper studying the count of solutions to $ABA^{-1}B^{-1} = \zeta \mathrm{Id}_n$, for $\zeta$ a primitive $n$-th root of unity but, as they note in Section 3.8, their methods can also be adapted to non-primitive $n$-th roots of unity, such as $-1$.
