$\newcommand{\ep}{\epsilon}
\newcommand{\ga}{\gamma}
\newcommand{\la}{\lambda}
\newcommand{\La}{\Lambda}
\newcommand{\Si}{\Sigma}
\newcommand{\R}{\mathbb{R}}
\newcommand{\E}{\operatorname{\mathsf E}}
\newcommand{\PP}{\operatorname{\mathsf P}}$
Let
\begin{equation}
R:=\frac{w^T M^3 w}{w^T M^2 w}.
\end{equation}
The distribution of none of the random variables (r.v.'s) $w^T M^3 w$, $w^T M^2 w$, $R$ is Gaussian -- because each of them is positive.
However, $R$ is the ratio of two quadratic forms in iid standard normal r.v.'s $Z_1,\dots,Z_d$. Indeed, by the spherical symmetry of the distribution of the standard normal random vector $Z:=(Z_1,\dots,Z_d)$, the random vector $w$ equals $Z/\|Z\|$ in distribution, where $\|Z\|$ is the Euclidean norm of $Z$. So, we may write $R$ as the ratio of two quadratic forms in $Z_1,\dots,Z_d$:
\begin{equation}
R=\frac{Q_3}{Q_2},\quad Q_j:=Z^T M^j Z.
\end{equation}
Consider now the spectral decomposition $M=O^T\La O$ of the positive definite matrix $M$, where $O$ is an orthogonal matrix and $\La$ is the diagonal matrix with the eigenvalues $\la_1\ge\cdots\ge\la_d>0$ of $M$ on the diagonal. Then
\begin{equation}
Q_j=Z^T O^T \La^j O Z=V^T \La^j V=\sum_1^d\la_i^j V_i^2,
\end{equation}
where $V=(V_1,\dots,V_d):=OZ$, so that $V$ is also a standard normal random vector, and hence $V_1,\dots,V_d$ are iid standard normal r.v.'s.
Note also that $Q_3=\sum_1^d\la_i^3 V_i^2\le\la_1\sum_1^d\la_i^2 V_i^2=\la_1Q_2$ and, similarly, $Q_3\ge\la_dQ_2$, so that we have the specific bounds on the values of $R$, namely $\la_d\le R\le\la_1$.
In particular, if $\la_1=\la_d=:\la$, then $R=\la$, a constant. To avoid this triviality, let us assume that from now on we have $\la_d<\la_1$.
Then, using transformations of distributions of random vectors (see e.g. formula (10)), one can find integral expressions for the density (pdf) $f_R$ of $R$.
For $d=2$, this pdf can be found explicitly. Indeed, then we can write
\begin{equation}
R=\frac{\la_1^3 V_1^2+\la_2^3 V_2^2}{\la_1^2 V_1^2+\la_2^2 V_2^2}
=\frac{\la_1^3 F+\la_2^3}{\la_1^2 F+\la_2^2}=:\rho(F),\quad
F:=\frac{V_1^2}{V_2^2}.
\end{equation}
Let $u:=\rho^{-1}$, so that
\begin{equation}
u(r)=\frac{\la_2^2 (r-\la_2)}{\la_1^2 \left(\la_1-r\right)}
\end{equation}
for $r\in(\la_2,\la_1)$. Note also that $F$ has the $F_{1,1}$ distribution, with pdf $f_F$ given by
\begin{equation}
f_F(x)=\frac1{\pi(1+x)\sqrt x}
\end{equation}
for $x>0$.
So, by formula (5), the pdf $f_R$ of $R$ is given by
\begin{equation}
f_R(r)=f_V(u(r))|u'(r)|=
\frac{\la_1 \la_2}{\pi \sqrt{\left(\la_1-r\right) \left(r-\la_2\right)}
\,[\la_1^2+(\la_1+\la_2) \left(\la_2-r\right)]}
\end{equation}
for $r\in(\la_2,\la_1)$.
Here is the graph of $f_R$ for $\la_1=2$ and $\la_2=1$; the shape is certainly not Gaussian.
