Marcel's solution provides a  good approach for understanding the marginal statistics.  Here are a few supplementary comments that might give a little insight into the joint distribution.

Let $D$ be the discrete Fourier transform matrix, i.e. the $j,k$-th entry is:
$$D_{j,k}=e^{-2\pi i jk/N}/\sqrt{N}$$
Consider the discrete Fourier transform of the first row of $M$, i.e.
$$(G_1,...,G_N)=(1/\sqrt{N})(X_1,...,X_N)D$$
Let $G$ be a diagonal matrix with diagonal entries $G_1,...,G_N$.
Then $D$ diagonalizes $M$ (see, e.g., [this description][1]), that is:
$$M = D G D^{-1}$$

If we took the $X_i$ to be complex-valued Gaussian variables, then we would be essentially done at this point: since $D$ is unitary, and the Gaussian is spherically symmetric, then the $G_i$ would be i.i.d. (complex) Gaussian random variables with mean 0 and variance $1/n$.  (That is, if we sampled the $G_i$ as i.i.d complex $\mathcal{N}(0,1/n)$, then $DGD^{-1}$ has the same distribution as samples from the original circulant matrix.)  It follows that
$$MM^*=(DGD^{-1})(DG^*D^{-1})=DGG^{-1}D^*$$
It then follows that the eigenvalues of $MM^{*}$ (or, if you prefer, the squared singular values of $M$) are distributed like $n$ i.i.d draws from a rescaled $\chi^2$ distribution with 2 degrees of freedom (which happens to simplify to an exponential distribution), where we rescale by dividing by $n$.  From an eigenvalue perspective, this is a complete characterization of the "orthonormality" of $M$.

However, your $X_i$ are probably real-valued.  This causes a small book-keeping headache, but doesn't really change much.  In brief: since the $X_i$ are real-valued, the $G_i$ will be symmetric.  We can convert our $n\times n$ complex-valued matrices into corresponding $2n \times 2n$ real-valued matrices (say, $D'$ and $G'$). Then observe that there are only $n$ of the $2n$ diagonal elements in $G$ are free, and only $n$ of the $2n$ inputs to $D'$ are non-zero (namely, the ones corresponding to the real values of the inputs in $D$).  Noting that the resulting decimated matrices are still unitary, we then reduce to the previous case, except with $\chi^2/n$ distributions with only 1 degree of freedom.


  [1]: https://www.cs.unm.edu/~williams/cs530/dft.pdf