What is the error in my derivation of the sum of squares of two standard normal iid variables? I am trying to derive the probability distribution of $Z=Z_1^2 + Z_2^2$. I am doing this in the following way,
$$ \begin{align}p(Z=x) = &\int_0^x p(Z_1^2 = p)\times p(Z_2^2 = x-p) dp\\ \Leftrightarrow ~ &\int_0^x p(Z_1 = \pm\sqrt{p})\times p(Z_2 = \pm\sqrt{x-p})dp \\ \Leftrightarrow &\int_0^x [2p(Z_1 = \pm\sqrt{p})]\times [2p(Z_2 = \pm\sqrt{x-p})]dp \\ \Leftrightarrow 4&\int_0^x (\frac{1}{\sqrt{2\pi}} e^{-0.5p}) \times (\frac{1}{\sqrt{2\pi}} e^{-0.5(x-p)})dp\\ \Leftrightarrow \frac{4}{2\pi}&\int_0^x e^{-0.5p} \times  e^{-0.5(x-p)}dp \\ \Leftrightarrow \frac{4}{2\pi}&\int_0^x e^{-0.5x}dp \\ \Leftrightarrow \frac{4}{2\pi} &[p\times e^{-0.5x}]_0^x  \\ \Leftrightarrow \frac{4}{2\pi}&\int_0^x e^{-0.5x}dp \\ \Leftrightarrow \frac{2}{\pi} &x\times e^{-0.5x}\end{align}$$
However, this should be the same as the gamma distribution, $\operatorname{Gamma}(1,\frac{1}{2})$, which has pdf
$$ \frac{1}{2}e^{-0.5t} $$
With the value $x=0.3 = t$, we get two different answers, so the pdfs cannot be the same. What am I doing wrong? The limits should be fine as any values of $p$ outside of the interval $[0,x]$ would give a negative inside the square roots.

You are missing a term $\frac{1}{\sqrt{x^2-p^2}}$ and you should be ending up with using an integral like this
$$\int_0^x \frac{1}{\sqrt{x^2-p^2}} e^{-x} dp = \frac{\pi}{2} e^{-x}$$
instead of
$$\int_0^x e^{-x} dp = x e^{-x}$$
This error happens in the beginning when you change the variables from $Z_1^2$ and  $Z_2^2$ to $Z_1$ and $Z_2$, your term $dp$ got changed here. The extra term is how the area changes (you integrate along a 1 dimensional line but this is a shortcut for integrating the area for an infinitesimal change of $dZ$ along this line).
It is more foolproof when you compute $p[Z≤x]$ instead of $p[Z=x]$.
