What is the distribution of the sum of squared chi-square random variables? What would be the distribution of the following equation:
$$y = a^2 + 2ad + d^2$$
where $a$ and $d$ are independent non-central chi-square random variables with $2 \textbf{M}$ degrees of freedom.
OBS.: The r.v.'s generating both $a$ and $d$ have $\mu = 0$ and $\sigma^2 \neq 1$, let's say $\sigma^2 = c$.
 A: If $a, d\sim\chi^2_{2M}$ are independent, then $X=a+d$ will have $\chi^2_{4M}$ distribution. Since $X$ is non-negative, CDF of $Y=a^2+2ad+d^2=(a+d)^2=X^2$ can be found by noting $$F_Y(y)=P(Y\leq y)=P(X^2\leq y)=P(X\leq \sqrt{y})=F_X(\sqrt{y}).$$ Therefore, $$f_Y(y)=\frac{1}{2\sqrt{y}}f_X(\sqrt{y})=\frac{1}{2^{2M+1}\Gamma(2M)}y^{M-1}e^{-\sqrt{y}/2}.$$
If $a$ and $d$ are correlated then things are much more intricate. See for example N. H. Gordon & P. F. Ramig's Cumulative distribution function of the sum of
correlated chi-squared random variables (1983) for a definition of multivariate chi-squared and distribution of its sum.
If $\mu\neq 2M$ then you are dealing with non-central chi-squared so the above will no longer be valid. This post may provide some insight.
EDIT: Based on the new information it seems $a$ and $d$ are formed by summing up normal r.v. with non-unit variance. Recall if $Z\sim N(0, 1)$ then $\sqrt{c}Z\sim N(0, c)$. Since now $$a=c\sum_{i=1}^{2M}Z_i^2=d,$$ both $a,d$ will have chi-squared distribution scaled by $c$, i.e. $\Gamma(M, 2c)$ distribution. In this case $X=a+d$ will be $\Gamma(2M, 2c)$ distributed. As a result, for $Y=X^2$ we have $$f_Y(y)=\frac{1}{2(2c)^{2M}\Gamma(2M)}y^{M-1}e^{-\sqrt{y}/2c}.$$
A: Since a non-central chi-square is a sum of independent rv's, then the sum of two independent non-central chi-squares $X = a+b$ is also a non-central chi-square with parameters the sum of the corresponding parameters of the two components, $k_x = k_a+k_b$ (degrees of freedom), $\lambda_x = \lambda_a+\lambda_b$ (non-centrality parameter). 
To obtain the distribution function of its square $Y =X^2$ , one can apply the "CDF method" (as in @francis answer), 
$$F_Y(y)=P(Y\leq y)=P(X^2\leq y)=P(X\leq \sqrt{y})=F_X(\sqrt{y})$$
and where 
$$F_X(x)=1 - Q_{k_x/2} \left( \sqrt{\lambda_x}, \sqrt{x} \right)$$
so
$$F_Y(y)=1 - Q_{k_x/2} \left( \sqrt{\lambda_x}, y^{1/4} \right)$$
where $Q$ here is Marcum's Q-function.
The above apply to non-central chi-squares formed as sums of independent squared normals each with unitary variance but different mean.
ADDENDUM RESPONDING TO QUESTION'S EDIT 
If the base rv's are $N(0,c)$, then the square of each is a $Gamma (1/2,2c)$
see https://stats.stackexchange.com/a/122864/28746 . 
So the rv $a \sim Gamma (M, 2c)$ and $b \sim Gamma (M, 2c)$ so also $X = a+b \sim Gamma(2M, 2c)$  (shape-scale parametrization, and see the wikipedia article for the additive properties for Gamma).  
Then one can apply again the CDF method to find the CDF of the square $Y = X^2$
