Let $q\in \mathbb{C}$ be a primitive $n$th root of unity, for some $n>1$. Consider the $n^2\times n^2$-matrix $$M=\left( q^{ki+lj}\right)_{(k,l),(i,j)}$$ indexed by all pairs $(k,l), (i,j) \in \{0,\cdots, n-1\}^2$ using lexicographic order.
This means for $n=2$ and $q=-1$: $$M = \left(\begin{array}{rrrr} 1&1&1&1\\1 & 1 &-1&-1\\ 1&-1&1&-1\\1&-1&-1&1\end{array}\right)$$
And for $n=3$: $$M = \left(\begin{array}{rrrrrrrrr} 1&1&1&1&1&1&1&1&1\\ 1 & q & q^2 & 1 & q & q^2 & 1 & q & q^2\\ 1 & q^2 & q & 1 & q^2 & q & 1 & q^2 & q\\ 1 & 1 & 1 & q & q & q & q^2 & q^2 & q^2\\ 1 & q & q^2 & q & q^2 & 1 & q^2 & 1 & q\\ 1 & q^2 & q & q & 1 & q^2 & q^2 & q & 1\\ 1 & 1 & 1 & q^2 & q^2 & q^2 & q & q & q\\ 1 & q & q^2 & q^2 & 1 & q & q & q^2 & 1\\ 1 & q^2 & q & q^2 & q & 1 & q & 1 & q^2 \end{array}\right)$$
My question is: Are these matrices $M$ always invertible? I guess yes, but I was not able of showing it. However I calculated their determinants for the cases $n<17$ and found $\mathrm{det}(M) = \lambda n^{n^2}$, where $\lambda\in\{-1,1\}$. So for $n=2$ it is $16$ and for $n=3$ it is $-19683 = -3^9$. Actually it turned out that the determinant is negative iff $n=4k-3$ for some $k$.
Any help or suggestion are very welcome.
Let us have a look at your $9\times 9$ example. It has the form $$\left(\begin{matrix}V & V & V\\V & qV & q^2V\\V & q^2V & qV\end{matrix}\right),$$where $V$ is the $3\times 3$ (rearranged) Vandermonde matrix in the upper left. This is actually the Kronecker product of $V$ with itself. I think, this is a general pattern, so that the $n^2\times n^2$ matrix is the Konecker product of $V = (q^{ij})_{i,j=0}^{n-1}$ with itself. The determinant is $\det^{2n}(V)$. And since $V$ is of (rearranged) Vandermonde form, your matrix is always invertible.
EDIT: Note that you might have to rearrange your rows to get that form (see your $4\times 4$ example). This gives you the minus sign for some dimensions.
EDIT: Look here: https://www.udacity.com/course/viewer#!/c-ud061/l-3495828730/m-1178758803. It shows that $V(q)V(q^{-1}) = nI$. Taking into account that $q^{-1} = \overline q$, we have $V(q^{-1}) = V(\overline q) = \overline{V(q)}$. Therefore, $$ n^n = \det(nI) = \det(V)\det(\overline V) = \det V \overline{\det V} = |\det V|^2. $$ Thus, $|\det V| = n^{n/2}$. So, the abolute value of the determinant of your matrix is $\det^{2n}(V) = n^{n^2}$, which agrees with your result.