What is $\pi$ doing here in Gaussian integral? 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.
 A: One way to visualize this is to note that the integral
$$
I^2=\int_{-\infty}^\infty\int_\infty^\infty e^{-(x^2+y^2)}\,\mathrm{d}x\,\mathrm{d}y
$$
is the volume of the region between $z=e^{-(x^2+y^2)}=e^{-r^2}$ and the $xy$-plane.  So taking slices parallel to the $xy$-plane, we have
$$
I^2=\int_0^1 A(z)\,\mathrm{d}z
$$
where $A(z)$ is the area of the slice at height $z$ intersecting the region.  The slice is a disc of radius $r$ where $r^2=-\log z$.  So
$$
I^2=\int_0^1 -\pi\log z\,\mathrm{d}z
$$
and there is your $\pi$.  We must also check the integral doesn't end up cancelling this $\pi$ out:
$$
\frac{I^2}{\pi}=-\int_0^1\log z\,\mathrm{d}z=-\Big[z\log z-z\Big]_{z=0}^{z=1}=1+\underbrace{\lim_{z\to 0+} z\log z}_{=0}.
$$
A: I understand your frustration; the calculation tells us $\pi$ appears in the answer, but it seems to come from nowhere. What does exponentiation have to do with circles/rotation? Weirder still, why is it square-rooted?
Surprisingly, the best place to get an intuition for this isn't geometry; it's statistics.
One way to think of it is this. Since $r^2=x^2+y^2$, $\exp -r^2$ is separable in $x,\,y$ (meaning it's a function of $x$, times a function of $y$). Similarly, the Jacobian $r$ in $dxdy=rdrd\theta$ is separable (indeed, the "function of $\theta$" we'd use is constant, but that still counts). And where do separable functions famously come up? The "joint" distribution of independent random variables.
Now $2r\exp -r^2$ is a pdf on $[0,\,\infty)$, and $\frac{1}{2\pi}$ is a pdf on $[0,\,2\pi)$. Multiplying two expressions for $1$ to get another expression for $1$, we have a joint distribution for polar variables $r,\,\theta$ whereby they're also independent:
$$1=\int_{r\ge 0,\,0\le\theta\le2\pi}f(r,\,\theta)drd\theta,\,f:=\frac{1}{\pi}r\exp-r^2.$$But deliciously, this distribution also makes the Cartesian coordinates $x\leftrightarrow y$ independent! (In fact, you can show that, to within scaling, this is the only way for a distribution in the plane to satisfy both independence conditions.) Obviously we can rewrite the above integral as $$1=\int_{x,\,y\in\Bbb R}\frac{f(r,\,\theta)}{r}dxdy=\int_{x,\,y\in\Bbb R}\frac{1}{\pi}\exp\left(-x^2\right)\exp\left(-y^2\right)dxdy.$$The Cartesian formalism has the added beautiful consequence that the distributions of $x,\,y$ are identical. Indeed, $X$ has pdf $$\frac{1}{\sqrt{\pi}}\exp-x^2,$$which is equivalent to your original result. So in short, the reason $\sqrt{\pi}$ comes up is because of the very special way you can simultaneously make $x,\,y$ independent and $r,\,\theta$ independent.
A: First of all, this is not an area:
$$ \int _{-\infty}^\infty e^{-x^{2}}\,dx . $$
It is a definite integral that evaluates to a number. You might associate the number with an area, but only under specific interpretations of what the integral is for. (For example, if someone asks for the area of the region below the graph of the function $e^{-x^{2}}$ and above the $x$-axis, then this integral is how you compute that area.)
In this exercise the fact that there are such interpretations is mostly spurious and misleading.
Just to be clear, you really should not get hung up on having to imagine an area every time you see an integral. There are many applications 
(for example, electromagnetism)
that are chock full of integrals, most of which are only related to any kind of "area" in the most abstract way.
So when you see
$$ \int_{-\infty}^\infty e^{-x^{2}}\,dx  \int_{-\infty}^\infty e^{-y^{2}}\,dy, $$
it's just the product of two numbers, nothing more.
Certainly not the product of two areas; that makes no sense (at least in this context).
The "two" numbers in fact are just one number:
$$ \int_{-\infty}^\infty e^{-x^{2}}\,dx = \int_{-\infty}^\infty e^{-y^{2}}\,dy.$$
Whether we write $x$ or $y$ or $\theta$ inside the integral really doesn't matter (yet).
It's when we write
$$ \int_{-\infty}^\infty \int_{-\infty}^\infty e^{-(x^2+y^2)}\,dx\,dy $$
that it finally makes sense to regard the $x$ and $y$ as $x$ and $y$ coordinates on a Cartesian plane, because that interpretation helps in visualizing the transformation to polar coordinates (if that's what comes next) or in reinterpreting the volume under $e^{-(x^2+y^2)}$ as a set of concentric shells or stacked disks.
So the progression is number, number times itself (i.e. number squared),
number times itself written slightly differently, 
rearrange the order of integration so we have a double integral instead of two single ones, and then
reinterpret the double integral as a volume between a curved surface and the Cartesian plane.
The last step should be the first time we invoke any geometric intuition.
The step 
$$ \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)}\,dx\,dy $$
is admittedly a bit much to swallow all at once if you're not used to this sort of thing. Step by step, it can be done like this:
$$ \int_{-\infty}^\infty e^{-x^2}\,dx = K. $$
This requires that we know (or assume) the integral exists.
But if it exists, it's just a number, so we can call it $K.$
$$ \int_{-\infty}^\infty e^{-x^2}\,dx  \int_{-\infty}^\infty e^{-y^2}\,dy 
= K \int_{-\infty}^\infty e^{-y^2}\,dy .$$
Simple substitution using the previous equation.
$$ K \int_{-\infty}^\infty e^{-y^2}\,dy = \int_{-\infty}^\infty Ke^{-y^2}\,dy.$$
Since $K$ is simply a fixed number (though we haven't computed it yet), we get the same result multiplying by it before or after integrating.
$$\int_{-\infty}^\infty Ke^{-y^2}\,dy =
\int_{-\infty}^\infty \left(\int_{-\infty}^\infty e^{-x^2}\,dx\right)e^{-y^2}\,dy.
$$
All we did here was reverse the substitution we did earlier.
Now we look at what is inside the $dy$ integral:
$$ \left(\int_{-\infty}^\infty e^{-x^2}\,dx\right)e^{-y^2}.$$
In this expression $e^{-y^2}$ is simply an (unknown) number that we are multiplying the integeral $\int_{-\infty}^\infty e^{-x^2}\,dx$ by, so again we can get the same result multiplying inside the integral as outside:
$$ e^{-y^2}\int_{-\infty}^\infty e^{-x^2}\,dx = 
\int_{-\infty}^\infty e^{-y^2}e^{-x^2}\,dx.$$
Finally, $e^{-y^2}e^{-x^2} = e^{-(x^2+y^2)}.$
Putting it all together we get
$$
\int_{-\infty}^\infty \left(\int_{-\infty}^\infty e^{-x^2}\,dx\right)e^{-y^2}\,dy
= \int_{-\infty}^\infty \left(\int_{-\infty}^\infty e^{-(x^2+y^2)}\,dx\right)dy.
$$
And then Fubini's Theorem lets us treat
$$
\int_{-\infty}^\infty \left(\int_{-\infty}^\infty e^{-(x^2+y^2)}\,dx\right)dy
= \int_{-\infty}^\infty \int_{-\infty}^\infty e^{-(x^2+y^2)}\,dx\,dy
$$
as an integral over the combined $(x,y)$ coordinates of a Cartesian plane
instead of just one integral inside another,
which we need in order to find the circles in the shape of this integral.

Conceptually, the direction in which I would approach the calculation of the constant in the Gaussian distribution is exactly the opposite of the direction taken in the previous part of this answer.
That is, I would start by setting up a joint distribution of two iid Gaussian variables as a function over an $x,y$ Cartesian plane,
using an as-yet-unknown constant factor to make this a probability distribution.
That is, I would start with a distribution that is already two-dimensional.
Then I would show that the integral of that distribution can be written as the product of two integrals, which can be rewritten as the square of one integral.
