Let $n:=M$. Since you say "the denominator follows a Chi-square distribution with $2M$ degrees of freedom", it appears that you assume the $x_i$'s to be iid standard complex Gaussian; anyway, if you need $a>1$, then you can do a simple rescaling.
Thus, let us assume indeed that the $x_i$'s are iid standard complex Gaussian. It is then clear that the distribution of $z$ in $\mathbb C=\mathbb R^2$ is rotation invariant. So, it is enough to find the distribution of $|z|^2$, which is the same as that of \begin{equation} R:=\frac X{(X+Y)^2}, \end{equation} where $X$ and $Y$ are independent random variables such that $X\sim\text{Gamma}(1,2)$ and $Y\sim\text{Gamma}(n-1,2)$.
Considering now the transformation from $(X,Y)$ to $(R,S)$ with $S:=X$ (say) and using the corresponding formula -- see e.g. formula (41) in Transformations, we obtain the joint density of $(R,S)$, say $f_{R,S}$. Then the density of $R$ is given by the formula \begin{equation} f_R(r)=\int_0^\infty f_{R,S}(r,s)\,ds. \end{equation} Using in this integral the substitution $v=\sqrt{rs}$, we find \begin{equation} f_R(r)=\frac1{(n-2)!2^n r^{n+1}}\,\int_0^1 v^n(1-v)^{n-2}\exp\Big\{-\frac v{2r}\Big\}\,dv \tag{1} \end{equation} for $r>0$, with $f_R(r)=0$ for $r\le0$.
Using the binomial expansion for $(1-v)^{n-2}$ and then integrating by parts, we can express the latter integral as a finite sum (with the number of summands depending on $n$) of elementary functions. For instance, for $n=3$ and $r>0$ \begin{equation} f_R(r)=\exp\Big\{-\frac1{2r}\Big\}\left(\frac{1}{2 r^2}+96 r+\frac{6}{r}+36\right)-96 r+12; \end{equation} here is the graph of this function:
However, by the dominated convergence theorem, the integral in (1) converges to $c_n:=\int_0^1 v^n(1-v)^{n-2}\,dv=n!(n-2)!/(2n-1)!$ as $r\to\infty$, and so, $f_R(r)\sim \frac{c_n}{(n-2)!2^n r^{n+1}}$ as $r\to\infty$. So, the tail of the distribution of $|z|^2$ is power-like, which shows that "the resulting distribution [of $z$]" is not Gaussian.