**HINT**
If we make the substitution $y = x - b$, we get: \begin{align*} \int_{-\infty}^{+\infty} xe^{-a(x-b)^{2}}\mathrm{d}x & = \int_{-\infty}^{+\infty}(y+b)e^{-ay^{2}}\mathrm{d}y \\\ & = \int_{-\infty}^{+\infty}ye^{-ay^{2}}\mathrm{d}y + \int_{-\infty}^{+\infty}be^{-ay^{2}}\mathrm{d}y\\\ & = \int_{-\infty}^{+\infty}be^{-ay^{2}}\mathrm{d}y \end{align*}
once the function $ye^{-ay^{2}}$ is odd. Thus we have: \begin{align*} \left(\int_{-\infty}^{+\infty}be^{-ay^{2}}\mathrm{d}y\right)^{2} & = \left(\int_{-\infty}^{+\infty}be^{-ay^{2}}\mathrm{d}y\right)\times\left(\int_{-\infty}^{+\infty}be^{-az^{2}}\mathrm{d}z\right)\\\ & = \int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}b^{2}e^{-a(y^{2} + z^{2})}\mathrm{d}y\mathrm{d}z \end{align*}
Hence, if you make the change of variable $y = \rho\cos(\theta)$ and $z = \rho\sin(\theta)$, you are able to obtain the sought result. Hope this helps.