The Gamma function is $\Gamma(\alpha)=\int_0^\infty x^{\alpha-1} e^{-x}\,dx$
Therefore
$\Gamma\left( \frac{1}{2}\right) =\int_0^\infty \frac{1}{\sqrt{x}} e^{-x}\,dx$
Thus, after the change of variable $t=\sqrt{x}$, this turns into the Euler integral $t=\sqrt{x} \implies x = t^2$ $\frac{dt}{dx} = \frac{1}{2\sqrt{x}} = \frac{1}{2t} \iff dx = 2t\,dt$
We have $\int_0^\infty \frac{1}{t} e^{-t^2}\,2t\,dt = 2\int_0^\infty e^{-t^2}\,dt $
And following holds: $\Gamma\left( \frac{1}{2}\right) = 2\int_0^\infty e^{-t^2}\,dt = \int_{-\infty}^\infty e^{-t^2}\,dt$
Which $\int_{-\infty}^\infty e^{-t^2}\,dt$ is the Gaussian integral. It follows:
$\left[ \Gamma\left( \frac{1}{2}\right)\right]^2 = \bigg(\int_{-\infty}^{+\infty} e^{-x^2} dx\bigg)\bigg(\int_{-\infty}^{+\infty} e^{-y^2} dy\bigg)$ $\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}e^{-(x^2+y^2)} \,dx\,dy$
Now change to polar coordinates
$\int_0^{+2 \pi}\int_0^{+\infty}e^{-r^2} r\,dr\,d\theta$
The $\theta$ integral just gives $2\pi$, while the $r$ integral succumbs to the substitution $u=r^2$
$\left[ \Gamma\left( \frac{1}{2}\right)\right]^2 = 2\pi\int_{0}^{+\infty}e^{-u} \,du/2=\pi$ $\Gamma\left( \frac{1}{2}\right)=\sqrt{\pi}$