Let $W=[X_1,Y_1,\cdots,X_N,Y_N]^T$, and $Z=[X_1+iY_1,\cdots,X_N+iY_N]^T$, where $W\sim~N(0,\Sigma)$, and $W\sim~CN(0,M)$. According to Wiki, the pdf of $2N$ multivariate of zero mean Gaussian distribution is: $$f_W(w)=\frac{1}{\sqrt{(2\pi)^{2N} |\Sigma|}} \text{exp}\left(\frac{-w^H \Sigma^{-1} w}{2} \right),\tag{1}$$ and the pdf of circular symmetric complex Gaussian distribution is: $$f_Z(z)=\frac{1}{(\pi)^{N}|M|}\text{exp}\left(-z^H M^{-1} z \right).\tag{2}$$
I understand the earlier since everything is real. But how to prove (2) from (1)? And how do I plot (2), i.e., do I have to consider the complex domain since Z is complex?
Let's introduce in the game the complex conjugate variable \begin{align} z_{n}^{\star} = x_{n} - jy_{n} \end{align} in addition to \begin{align} z_{n} = x_{n} + jy_{n} \end{align}
We can now write the following identity: \begin{align} \begin{pmatrix} z_{n} \\ z_{n}^{\star} \end{pmatrix}=J \begin{pmatrix} x_{n} \\ y_{n} \end{pmatrix} \end{align} where it is easy to see that \begin{align} J= \begin{pmatrix} 1&j \\ 1&-j \end{pmatrix} \end{align}
We can now define the vector $v = [z_{1}z_{1}^{\star},\dots,z_{N}z_{N}^{\star}]^{T}$
Similar to what we did above we can now extend the identity and write $v=\tilde{J}w$, where \begin{align} \tilde{J}= \begin{pmatrix} J & 0 & \dots \\ \vdots & \ddots & \\ 0 & & J \end{pmatrix} \end{align} is a diagonal block matrix with $J$ on the main diagonal.
If we now write the covariance of $v$, $T=\mathbf{E}[vv^{H}]$. Hence $T = \tilde{J}\Sigma\tilde{J}^{H}$
The quadratic form $w^{T}\Sigma^{-1}w = v^{H}T^{-1}v$, considering $\tilde{J}^{-1} = 1/2\tilde{J}^{H}$, and $w=\tilde{J}^{-1}v$
We further need $\det(\Sigma) = \det(1/2\tilde{J}^{H}T1/2\tilde{J}) = (1/2)^{2N}\det(T)$
All of the above gives us the pdf of $v$: \begin{align} f(v) = \frac{1}{\pi^{N}\det(T)^{1/2}}\exp(-1/2v^{H}T^{-1}v) \end{align}
Now let's get to our special case. It should be clear by now that the expression above should be a generalization of the initial problem.
Notice that with an appropriate choice of a permutation matrix $\Pi$ one can write $u = \Pi v = [z,z^{\star}]^{T}$
Thus $v^{H}T^{-1}v = u^{H}(\Pi T \Pi^{T})^{-1}u = u^{H}(U)^{-1}u$ where because of complex circularity: \begin{align} U= \begin{pmatrix} M&0 \\ 0&M^{\star} \end{pmatrix} \end{align}
which thus leads to $v^{H}T^{-1}v = 2z^{H}M^{-1}z$. With $\det(T) = \det(\Pi)\det(U)\det(\Pi) = \det(U) = \det(M)\det(M^{\star}) = \det(M)^2$
Substituting all the element back in $f(v)$ we get: \begin{align} f(z) = \frac{1}{\pi^{N}\det(M)}\exp(-z^{H}M^{-1}z) \end{align} for a circular symmetric random vector.
For plotting one can just observe that with the derivation above, the 1-D complex distribution is identical to the joint distribution of the i.i.d real and imaginary parts (that only depends thus on the single parameter which is the variance of the unidimensional real variable).