Given a prime power $q$, I would like to enumerate (preferably up to isomorphism*) all the permutation polynomials $f(x)$ on $K = GF(q^3)$ satisfying the following conditions:

- $f(ax) = af(x)$ for all $a \in GF(q)$,
- $Tr_{K/F}(x f(x)) = 0$. ($F := GF(q)$)

Such polynomials correspond to perfect matchings in the incidence graph of the projective plane $PG(2,q)$ as follows: Look at $GF(q^3)$ as a $3$-dimensional vector space over $GF(q)$ equipped with the symmetric bilinear form $(x,y) \mapsto Tr_{K/F}(xy)$. Let the points of the projective plane be $1$-dimensional subspaces. And let the lines be given by the $2$-dimensional null spaces of the maps $x \mapsto Tr(ax)$, i.e. $a^\perp$, for a given $a \in GF(q^3)$. Therefore, a point $p$ is incident with a line $L$ if $Tr(x_p x_L) = 0$ where $x_p$ and $x_L$ are corresponding members of $GF(q^3)$. Now, every permutation polynomial on $GF(q^3)$ defined by the two conditions above maps $1$-dimensional subspaces to $1$-dimensional subspaces giving rise to a bijection between points and lines of $PG(2,q)$, i.e., a perfect matching in the incidence graph of $PG(2,q)$. We can see that this is an onto map where fiber of each perfect matching is a set of $(q-1)^{1+q+q^2}$ permutation polynomials.

By [1] a $k$-regular bipartite graph of size $2n$ has at least $$\left( \frac{(k-1)^{(k-1)}}{k^{(k-2)}} \right)^n$$ perfect matchings. Therefore there is a lower bound on the number of perfect matchings given by $\frac{q^{qn}}{(q+1)^{n(q-1)}}$ where $n = 1 + q + q^2$.

I also know that for for $q = 3$ there are exactly $3852$ perfect matchings giving rise to $31555584 = 3852 \times 2^{13}$ such permutations. This is from computer computations of the perfect matchings. For $q = 4$ there are $18534400$ of them.

To the finite field experts: Is there a way of theoretically classifying such polynomials? Can we give some bounds on the total number? Any suggestions on what approach might work?

*Look at the action of $\Gamma L(3,q)$ on $GF(q^3)$ viewed as a vector space over $GF(q)$. For a permutation polynomial $f$ and a group element $\sigma$ define $f^\sigma (x) = \hat{\sigma} (f(\sigma(x))$ where $\hat{\sigma}$ is the adjoint of $\sigma$ w.r.t. the symmetric bilinear form $(x, y) \mapsto Tr_{K/F}(xy)$.