The easiest way to solve this would be compute characteristic function $\mathbb{E}(\mathrm{e}^{i t Y}) = \mathbb{E}(\mathrm{e}^{i t (X_1^2 + X_2^2)}) = \mathbb{E}(\mathrm{e}^{i t X_2^2}) \mathbb{E}(\mathrm{e}^{i t X_1^2})$.
For normal variate:
$ \mathbb{E}(\mathrm{e}^{i t x^2}) = \int_{-\infty}^{\infty} \frac{\mathrm{d}x}{\sqrt{\pi} \sigma} \mathrm{e}^{i t x^2} \mathrm{e}^{- \frac{x^2}{\sigma^2}} = \frac{1}{\sigma} \left( \frac{1}{\sigma^2} - i t\right)^{-\frac{1}{2}} = \frac{1}{\sqrt{1- i t \sigma^2}} $
Then you would compare this to the characteristic function of $\Gamma$ distribution to determite that it is the characteristic function of of $\Gamma(\alpha = \frac{1}{2}, \beta = \sigma^2)$.
Added: I should add, after seeing Michael's response that combining we get $\mathbb{E}(\mathrm{e}^{i t Y}) = (1-i t \sigma^2)^{-1}$, and this is gamma distribution with $\alpha=1$, also known as exponential distribution.
Alternatively you should use the fact that
$X_1 =^d \frac{\sigma}{\sqrt{2}} Z_1$, where
$Z_1$ is standard normal variate. Then
$Y = X_1^2+X_2^2 = \frac{\sigma^2}{2} (Z_1^2+Z_2^2)$. Notice that
$Z_1^2 + Z_2^2$ follows
$\chi^2_2$, i.e chi-squared with two degrees of freedom, hence
$Y$ is a rescaled
$\chi^2_2$ variate.