During reading proof of Stirling's formula $n! \sim \sqrt{2n\pi}\left(\frac ne\right)^n$, I started finding the reason of $\pi$'s presence in the expression; which led to Gaussian Integral, i.e. $\int^{\infty}\_{-\infty}{e^{-x^2}}dx = \sqrt{\pi}$.

The main suspect according to me in its derivation is:

$$ \left(\int\_{-\infty}^{\infty} e^{-x^{2}} \, dx\right)^{2}
= \int\_{-\infty}^{\infty} e^{-x^{2}} \, dx \int\_{-\infty}^{\infty} e^{-y^{2}} \, dy
= \int\_{-\infty}^{\infty} \int\_{-\infty}^{\infty} e^{-(x^2+y^2)} \, dxdy. $$

This has that term $( x^2 + y^2 )$ which seems to be taking a sort of circle and consequently $\pi$ in the formula. 
I am not finding any bigger picture of how exactly this arrival is taking place here.

**Question**:

How can I visualize the Gaussian integral to get the intuition for $\pi$'s presence in it? Also,

$$ \left(\int\_{-\infty}^{\infty} e^{-x^{2}} \, dx\right)^{2} = \int\_{-\infty}^{\infty} e^{-x^{2}} \, dx \int\_{-\infty}^{\infty} e^{-y^{2}} \, dy $$ 

is quite confusing for me from perspective of $x$ and $y$ axis' emergence out of area's square. (I know its correct mathematically, but this is something that has boggled the noggin of mine.)

Thanks.