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}$$