A $k$-periodic matrix $M$ is an $\mathbb{R}^{n \times n}$ matrix such that $M^{k+1} = M$
Is it possible to characterize this set of $k$-periodic matrices? My goal is to exploit this characterization to randomly generate an arbitrary number of these that are all different from one another. Permutation matrices are periodic but otherwise don't fit the bill - there is a finite set of them ($n!$) and they are not necessarily $k$-peroidic - I am looking for an infinite set.
For $k=1$, every matrix is $k$-periodic, and so you can simply pick some arbitrary $n\times n$ matrices. We'll henceforth assume $k>1$.
Here is a characterization over $\mathbb C$. A matrix is $k$-periodic if and only if its Jordan normal form is, i.e. if and only if every Jordan block in its Jordan normal form is.
For blocks with nonzero eigenvalue, this eigenvalue must satisfy $\lambda^{k-1}=1$, and thus must be a $(k-1)$st root of unity. For a diagonal block, this is also sufficient; for non-diagonal blocks $J$ with nonzero eigenvalue, $J^{k-1}$ will have nonzero entries above the diagonal, and thus $J^k$ will not equal $J$ (since $J$ is invertible). So, the only possible blocks with nonzero eigenvalue are the blocks $[\lambda]$ where $\lambda^{k-1}=1$.
For a block $J$ of size $m$ with eigenvalue $0$, the matrix $J^k$ has entries of $1$ on the diagonal $k$ entries above the main diagonal, and zeros everywhere else (or is $0$ for $k\geq m$); this does not equal $J^k$ if $m>1$ and $k>1$. So, the only possible Jordan block with $k>1$ and eigenvalue $0$ is the block $[0]$.
In other words, a $k$-periodic matrix must be diagonalizable, with eigenvalues $\lambda$ that satisfy $\lambda^k=\lambda$. This means that one can generate $k$-periodic $n\times n$ matrices $M$ by first choosing an arbitrary invertible $n\times n$ matrix $P$, then choosing an arbitrary $n\times n$ diagonal matrix $D$ with diagonal entries satisfying $\lambda^k=\lambda$, and finally set $M=PDP^{-1}$.
If you're looking only for real matrices, you can do this by more carefully constructing $P$ to be a matrix defined by eigenvectors, so that the eigenvector $v$ corresponding to a complex entry $\lambda$ of $D$ is the conjugate of the eigenvector $v'$ corresponding to an entry $\overline\lambda$.