By the spectral decomposition of $M$, the distribution of $y^T My$ is the same as that of $\sum_i\sigma^2\mu_i Z_i^2$, where the $\mu_i$'s are the eigenvalues of $M$ and the $Z_i$'s are iid standard normal random variables (r.v.'s). By rescaling, without loss of generality $\sigma=1$.
In the paper by Székely & Bakirov, exact lower and upper bounds on the cumulative distribution function of quadratic forms $Q:=y^T My$ are given. In particular, it is shown that
$$\inf_{n\ge1,M\ge0,EQ=1}P(Q\le x)=
\inf_{n\ge1}P(\chi^2_n/n\le x)=P(\chi^2_{n(x)}/n(x)\le x),
$$
where $x>0$, $n(x)$ depends explicitly and only on $x$, and $\chi^2_n$ is a r.v. with the chi-squared distribution with $n$ degrees of freedom.
Among other results, it is also shown in that paper that
$$\sup_{M\ge0,EQ=1}P(Q\le x)=P(\chi^2_n/n\le x)I\{x\ge y(n)\}
+s(x)I\{1\le x<y(n)\}+P(\chi^2_1\le x)I\{x<1\},
$$
where $I$ is the indicator, $x>0$, $y(n)$ depends explicitly and only on $n$, and
$$s(x):=\sup_{0\le\mu\le1/2}P(\mu Z_1^2+(1-\mu)Z_2^2\le x).
$$