Evaluation of Gaussian integral $\int_{0}^{\infty} \mathrm{e}^{-x^2} dx$ How to prove
 $$\int_{0}^{\infty} \mathrm{e}^{-x^2}\, dx = \frac{\sqrt \pi}{2}$$
 A: Consider the mapping $\eta\! : \mathbb{R}^2\to \mathbb{R}$ given by
$$
\eta((x,y)) = \sqrt{x^2+y^2},\quad (x,y)\in\mathbb{R}^2.
$$
(1) Show that the image-measure $\lambda_2\circ\eta^{-1}$ is the measure on $(\mathbb{R},\mathcal{B}(\mathbb{R}))$ with density
$$
f(z)=2\pi z 1_{(0,\infty)}(z),\quad z\in\mathbb{R},
$$
by using Dynkin's Lemma.
(2) Show that
$$
\int_{\mathbb{R}^2} e^{-x^2-y^2}\,\lambda_2(\mathrm{d}x,\mathrm{d}y)=\pi
$$
by using the formula for integration under measurable transformations.
(3) Use Tonelli's theorem to conclude that
$$
\int_{\mathbb{R}}e^{-x^2}\,\lambda(\mathrm{d}x)=\sqrt{\pi},
$$
and now your result follows.
A: A variation on Ross Millikan's answer. 
We can start again with the observation
$$\left(\int_{-\infty}^{\infty}e^{-x^2}dx\right)\left(\int_{-\infty}^{\infty}e^{-y^2}dy\right)=\left(\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}e^{-(x^2+y^2)}dxdy\right)=V.$$
Now $V$ is simply the volume of the body 
$$-\infty < x,y < \infty,\qquad 0 < z < e^{-(x^2+y^2)},$$
or, equivalently,
$$0 < x^2+y^2 < -\ln z,\qquad  0 < z < 1.$$
This implies that the body is a solid of revolution. Using the
disk integration formula, we have
$$V=\int_{0}^{1}\pi(-\ln z)dz=[\pi(z-z\ln z)]_{0}^1=\pi.$$
A: This is Exercise 7.19 of Apostle's Mathematical Nalysis book (second edition))

Define $f$ and $g$ as:
$$f(x):=\left(\int_0^x e^{-t^2}dt\right)^{2} \quad\text{and}\quad g(x):=\left(\int_{0}^{1}\frac{e^{-x^{2}(t^{2}+1)}}{t^{2}+1}dt\right)$$
Now, $$f'(x)=2e^{-x^{2}}\int_{0}^{x}e^{-t^{2}}dt$$ and
$$g'(x)=\int_0^1 \frac{\partial}{\partial x}\left[\frac{e^{-x^2(t^2+1)}}{t^2+1}\right]dt = -2xe^{-x^{2}}\int_{0}^{1}e^{-x^{2}t^{2}}dt$$
So putting $t=tx$, get  $\displaystyle\int_{0}^{1}e^{-x^{2}t^{2}}dt= \frac{1}{x}\displaystyle\int_{0}^{x}e^{-t^{2}}dt$
Then we get: 
$$g'(x)=-2e^{-x^{2}}\int_{0}^{x}e^{-t^{2}}dt$$
Thus $f'(x)+g'(x)=0$ for all $x$, then $f(x)+g(x)$ is an constant function. Also $$f(0)+g(0)=\displaystyle\int_{0}^{1}\frac{1}{t^{2}+1}dt = \displaystyle\frac{\pi}{4}$$
Then $f(x)+g(x)=\displaystyle\frac{\pi}{4}$ for all $x$.
Now $\lim_{x \to{+}\infty}{g(x)}=0$
So $$\displaystyle\frac{\pi}{4} = \lim_{x \to{+}\infty}{f(x)+g(x)}=\lim_{x \to{+}\infty}{f(x)}= \left(\int_{0}^{\infty}e^{-t^{2}}dt\right)^{2}$$
Thus
$$\int_{0}^{\infty}e^{-t^{2}}dt=\sqrt{\frac{\pi}{4}}= \frac{\sqrt{\pi}}{2}$$
The end.
A: I will add an additional Solution to this problem because it is using a powerful integral that Euler derived in one of his famous books about complex analysis (Euler's derivation).
$$\frac{\Gamma(s)}{n(a^2+b^2)^{s/2}}\cos(\alpha s)=\int_0^{\infty}u^{ns-1}e^{-au^n}\cos(bu^n) du$$, where $\tan(\alpha)=b/a$.
Plug in $a=1$, $b=0$ ($\alpha =0$),$n=2$ and $s=1/2$ to get:
$$\frac{\Gamma(1/2)}{2}=\int_0^{\infty}e^{-u^2}du$$
Evaluate Eulers reflection formula $\Gamma(s)\Gamma(1-s)=\frac{\pi}{\sin(\pi s)}$ for $s=1/2$ to get $\Gamma(1/2)=\sqrt{\pi}$.
With that we conclude: 
$$\int_0^{\infty}e^{-u^2}du=\frac{\sqrt{\pi}}{2}$$


*

*Note that there is a similar formula, containing $\sin$ instead of $\cos$. Both formulas can be used for many special integrals like the Fresnel Integral or the Sinc Integral and much more.

*Note that one could directly derive this result from $\Gamma(1/2)$ by substitution, but the formula is far more powerful.
A: It might be worth mentioning that one also can use spherical coordinates in 3-dimensions analogously to the polar coordinates Ross Millikan used above: If $I$ denotes $\int_{-\infty}^{\infty} e^{-x^2}dx$, then we have
$$I^3 =  \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}e^{-x^2 - y^2 - z^2}\,dx\,dy\,dz$$
Switching to spherical coordinates this becomes
$$\int_0^{2\pi}\int_0^{\pi}\int_0^{\infty}\sin(\phi) \rho^2 e^{-\rho^2}\,d\rho\,d\phi\,d\theta$$
Doing the theta and $\phi$ integrations this becomes 
$$I^3 = 4\pi\int_0^{\infty}\rho^2 e^{-\rho^2}\,d\rho$$
One can then integrate parts in this, differentiating $\rho$ and integrating $\rho e^{-\rho^2}$. This leads us to
$$I^3 = 2\pi \int_0^{\infty} e^{-\rho^2}\,d\rho$$
Note the right-hand side is exactly $2\pi\cdot {I \over 2} = \pi I$. Thus $I^3 = \pi I$ and thus $I = \sqrt{\pi}$ as needed. Obviously polar coordinates are faster. Just sayin'...
A: This is similar to user17762's answer, but uses Plancherel's Theorem instead of Poisson summation.  Define the Fourier transform by
$$\mathcal{F}[f](y) = \hat{f}(y) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty f(x) e^{-i xy}\;dx$$
Now,
$$\mathcal{F}[e^{-\frac{1}{2}x^2}]= \frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty e^{\frac{1}{2}x^2} e^{-i xy}\;dx = \frac{e^{-\frac{1}{2}y^2}}{\sqrt{2\pi}}\int_{-\infty}^\infty e^{\frac{1}{2}(x+iy)^2}\;dx$$
Now, consider the contour integral of $e^{-\frac{1}{2}z^2}$ over the rectangular contour with corners at $\pm R$ and $\pm R + iy$.  This integral must be $0$, since $e^{-\frac{1}{2}z^2}$ is analytic.  Taking the limits as $R \to +\infty$, the contributions from the vertical edges go to $0$, so we find that
$$\int_{-\infty}^\infty e^{-\frac{1}{2}x^2}\;dx = \int_{-\infty}^\infty e^{-\frac{1}{2}(x+iy)^2}\;dx$$
Thus, 
$$\mathcal{F}[e^{-\frac{1}{2}x^2}](y) = \frac{e^{-\frac{1}{2}y^2}}{\sqrt{2\pi}}\int_{-\infty}^\infty e^{-\frac{1}{2}x^2}\;dx$$
By Plancherel's Theorem,
$$\int_{-\infty}^\infty e^{x^2}\;dx = \int_{-\infty}^\infty e^{-y^2} \frac{1}{2\pi}\left(\int_{-\infty}^\infty e^{-\frac{1}{2}x^2}\;dx\right)^2\;dy = \frac{1}{2\pi}\left(\int_{-\infty}^\infty e^{-\frac{1}{2}x^2}\;dx\right)^2\int_{-\infty}^\infty e^{-y^2}\;dy$$
dividing through by $\int_{-\infty}^\infty e^{-x^2}\;dx$, we find
$$1 = \frac{1}{2\pi}\left(\int_{-\infty}^\infty e^{-\frac{1}{2}x^2}\;dx\right)^2$$
so,
$$\sqrt{2\pi} = \int_{-\infty}^\infty e^{-\frac{1}{2}x^2}\;dx$$
Changing variables and observing that $\int_0^\infty e^{-x^2}\;dx = \frac{1}{2}\int_{-\infty}^\infty e^{-x^2}\;dx$, we find that
$$\int_0^\infty e^{-x^2}\;dx = \frac{\sqrt{\pi}}{2}$$
A: An alternative derivation is to show that 
$$\int_{0}^{\infty}xe^{-x^{2}y^{2}}\; \mathrm{d}y=I,$$
where $I$ is your integral:
$$I:=\int_{0}^{\infty}e^{-x^2}\; \mathrm{d}x,$$
and then evaluate $I^2$ by reversing the order of integration. If $x>0$, then
$$\int_{0}^{\infty}xe^{-x^{2}y^{2}}\; \mathrm{d}y=x\int_{0}^{\infty}e^{-{(xy)}^2}\; \mathrm{d}y=x\int_{0}^{\infty}e^{-u^2}\dfrac{\mathrm{d}u}{x}=\int_{0}^{\infty}e^{-u^2}\; \mathrm{d}u=I.$$
Thus
$$\begin{aligned}I^2&=\int_{0}^{\infty}e^{-x^2}\; \mathrm{d}x\int_{0}^{\infty}xe^{-x^{2}y^{2}}\; \mathrm{d}y=\displaystyle\int_{0}^{\infty}\mathrm{d}y\int_{0}^{\infty}xe^{-x^2}e^{-x^{2}y^{2}}\; \mathrm{d}x\\ &=\int_{0}^{\infty}\mathrm{d}y\int_{0}^{\infty}xe^{-x^{2}(1+y^2)}\; \mathrm{d}x=\int_{0}^{\infty}\mathrm{d}y\dfrac{1}{2\left( 1+y^{2}\right) }\left[ -e^{-x^{2}\left( 1+y^{2}\right) }\right] _{x=0}^{\infty }\\ &=\int_{0}^{\infty }\dfrac{1}{2\left( 1+y^{2}\right) }\; \mathrm{d}y=\dfrac{1}{2}\left[ \arctan y\right] _{y=0}^{\infty }=\dfrac{\pi}{4}.\end{aligned}$$
So
$$I=\dfrac{\sqrt{\pi}}{2}.$$
A: Another proof, from G.M. Fichtengoltz, Calculus Course, page 612. 
$$K=\int_{0}^{\infty} e^{-x^2} dx $$
It easy to see (and prove) that,  $\max\{(1+t)e^{-t}\}=1$ at $t=0$, hence for all $t\in\mathbb{R}$:
$$(1+t)e^{-t}<1$$
Substitution of $t=\pm x^2$, leads us to:
$$(1-x^2)e^{x^2}<1 \ \ \ \ \text{and} \ \ \ \ \ (1+x^2)e^{-x^2}<1 $$ 
So,
$${1-x^2} <e^{-x^2}<\frac{1}{1+x^2} \ \ \ \ \ \ (x>0) $$ 
Now, at the left inequality we restrict our $x$ to be in $(0,1)$ (so that, $1-x^2>0$), and in the right inequality let $x>0$. Raising all the inequalities with natural number $n$, we get,
$$\underset{x\in (0,1)}{(1-x^2)^n<e^{-nx^2}} \ \ \ \ \text{and} \ \ \ \ \ \underset{x>0}{e^{-nx^2}<\frac{1}{(1+x^2)^n}}$$
Integrating the first inequality from $0$ to  $1$, and the second inequality from $0$ to $+\infty$ we'll get:
$$\int_0^1({1-x^2})^ndx <\int_0^1 e^{-nx^2} dx<\int_0^{\infty} e^{-nx^2} dx<\int_0^{\infty}\frac{dx}{(1+x^2)^n}$$
But,
$$\int_0^{\infty} e^{-nx^2}dx=\frac{1}{\sqrt{n}}K \ \ \ \ \ \ (\text{substitution} \ \  u=\sqrt{nx}),$$
$$\int_0^1({1-x^2})^ndx=\int_0^{\frac{\pi}{2}}\sin^{2n+1}(v)dv=\frac{(2n)!!}{(2n+1)!!}  \ \ \ (\text{substitution} \ \ x=\cos(v))$$
and, finally,
$$\int_0^{\infty}\frac{dx}{(1+x^2)^n}=\int_0^{\frac{\pi}{2}}\sin^{2n-2}(v)dv=\frac{(2n-3)!!}{(2n-2)!!}\frac{\pi}{2} \ \ \ (\text{substitution} \ \ x=\text{ctg}(v))$$
Hence, our unknown, $K$ is bound:
$$\sqrt{n}\frac{(2n)!!}{(2n+1)!!}<K<\sqrt{n}\frac{(2n-3)!!}{(2n-2)!!}\frac{\pi}{2}$$
or,
$$\frac{n}{2n+1}\frac{((2n)!!)^2}{((2n-1)!!)^2(2n+1)}<K^2<\frac{n}{2n-1}\frac{((2n-3)!!)^2(2n-1)}{((2n-2)!!)^2}(\frac{\pi}{2})^2$$
Now, the final step - Wallis Formula
:
$$\lim_{n\to\infty}\frac{((2n)!!)^2}{((2n-1)!!)^2(2n+1)}=\frac{\pi}{2}$$
Then, when $n$ tends to $\infty$ in our last inequality, we get:
$$K^2=\frac{\pi}{4}$$
and,
$$K=\frac{\sqrt{\pi}}{2}  \ \ \ \ \ \text{as}\ K>0 $$
A: By using Beta and Gamma functions properties we may simply obtain that:
$$\operatorname B\left(\tfrac 12,\tfrac12\right)=\frac{\left[\Gamma(\tfrac{1}{2})\right]^{2}}{\Gamma{(1)}}=\left[\Gamma(\tfrac{1}{2})\right]^{2}$$
$$\operatorname B\left(\tfrac{1}{2},\tfrac{1}{2}\right)=\frac{\pi}{\sin{\frac{\pi}{2}}}=\pi$$
In other words we have that:
$$\Gamma(\tfrac{1}{2})=\sqrt{\pi}\longrightarrow \space\int\limits_0^\infty x^{\frac{-1}{2}} e^{-x} \,\mathrm dx = \sqrt{\pi}$$
By substitution $x=t^2$ we get the final result:
$$2\int\limits_0^\infty e^{-t^2} \,\mathrm dt=\sqrt{\pi} \longrightarrow \int\limits_0^\infty e^{-t^2} \,\mathrm dt = \frac{\sqrt{\pi}}{2}.$$ 
Q.E.D. (Chris)
A: Using Normal density:
$$
\int_0^\infty e^{-x^2}dx=\sqrt \pi \int_0^\infty \frac{1}{\sqrt \pi}e^{-x^2}dx=\sqrt \pi P(X\geq0),
$$
where $X\sim \mathrm{Normal}(0,\frac{1}{2})$. Remember that the normal distribution is symmetric around its mean, so $P(X>0)=P(X>E(X))=\frac{1}{2}$, and the result follows.
A: Here's a proof that only requires elementary, but clever, calculus manipulations.
Let $I_A = \int_0^A e^{-x^2}dx$.
$$\begin{split} I_A^2 &= \int_0^A \int_0^A e^{-(x^2+y^2)}dxdy\\
&= \int_0^A \int_0^A \sum_{n\geq 0} \frac{(-1)^n}{n!}(x^2+y^2)^ndxdy\\
&=\sum_{n\geq 0} \frac{(-1)^n}{n!}\int_0^A \int_0^A\sum_{k=0}^n{n \choose k}x^{2k}y^{2n-2k}dxdy\\
&=\sum_{n\geq 0} \frac{(-1)^n}{n!}\sum_{k=0}^n{n \choose k}\frac{A^{2k+1}}{2k+1}\frac{A^{2n-2k+1}}{2n-2k+1}\\
&= \sum_{n\geq 0} \frac{(-1)^n}{n!}A^{2n+2}\sum_{k=0}^n{n \choose k}\frac{1}{2k+1}\frac{1}{2n-2k+1}
\end{split}$$
Now, note that
$$\begin{split}
\sum_{k=0}^n{n \choose k}\frac{1}{2k+1}\frac{1}{2n-2k+1} &= \frac 1 {2n+2}\sum_{k=0}^n{n \choose k}\left (\frac{1}{2k+1}+\frac{1}{2n-2k+1}\right)\\
& =\frac 1 {n+1}\sum_{k=0}^n{n \choose k}\frac{1}{2k+1}
\end{split}
$$
Thus,
$$\begin{split} I_A^2 &=
\sum_{n\geq 0} \frac{(-1)^n}{(n+1)!}A^{2n+2}\sum_{k=0}^n{n \choose k}\frac{1}{2k+1}\\
&= \sum_{n\geq 0} \frac{(-1)^n}{(n+1)!}A^{2n+2}\sum_{k=0}^n{n \choose k}\int_0^1x^{2k}dx\\
&= \sum_{n\geq 0} \frac{(-1)^n}{(n+1)!}A^{2n+2}\int_0^1(1+x^2)^ndx\\
&= \int_0^1\frac 1 {1+x^2}\sum_{n\geq 0}\frac{(-1)^n}{(n+1)!}A^{2n+2}(1+x^2)^{n+1}dx\\
&= \int_0^1\frac{1-e^{-A^2(1+x^2)}}{1+x^2}dx\\
&= \frac{\pi}4 + \mathcal O\left(e^{-A^2}\right)
\end{split}$$
Taking the limit as $A\rightarrow+\infty$ yields the result.
A: First, we notice that $$n!\ =\ \int_0^\infty e^{-\sqrt[n]x}\ dx\quad\iff\quad\tfrac1n!\ =\ \int_0^\infty e^{-x^n}dx\quad\rightarrow\quad\tfrac12!\ =\ \int_0^\infty e^{-x^2}dx$$ Then we further notice that $$\int_0^1\Big(1-\sqrt[n]x\Big)^m\,dx\ =\ \int_0^1\Big(1-\sqrt[m]x\Big)^n\,dx\ =\ \frac1{C_{m+n}^n}\ =\ \frac1{C_{m+n}^m}\ =\ \frac{m!\,n!}{(m+n)!}$$ From where we deduce that $$\frac\pi4\ =\ \int_0^1\sqrt{1-x^2}\,dx\ =\ \frac{\Big(\frac12!\Big)^2}{\Big(\frac12 + \frac12\Big)!}\ =\ \Big(\tfrac12!\Big)^2$$ Which leads us to conclude that $$\int_0^\infty e^{-x^2}dx\ =\ \tfrac12!\ =\ \sqrt{\pi\over4}\ =\ \frac{\sqrt\pi}2$$ QED
A: The following argument, similar to Bryan Yock's, is a Feynman parameter trick I invented in Integrating $\int^{\infty}_0 e^{-x^2}\,dx$ using Feynman's parametrization trick

Let $$I(b) = \int_0^\infty \frac {e^{-x^2}}{1+(x/b)^2} \mathrm d x = \int_0^\infty \frac{e^{-b^2y^2}}{1+y^2} b\,\mathrm dy$$ so that $I(0)=0$, $I'(0)= \pi/2$ and $I(\infty)$ is the thing we want to evaluate.
Now note that rather than differentiating directly, it's convenient to multiply by some stuff first to save ourselves some trouble. Specifically, note
$$\left(\frac 1 b e^{-b^2}I\right)' = -2b \int_0^\infty e^{-b^2(1+y^2)} \mathrm d y = -2 e^{-b^2} I(\infty)$$
Then usually at this point we would solve the differential equation for all $b$, and use the known information at the origin to infer the information at infinity. Not so easy here because the indefinite integral of $e^{-x^2}$ isn't known. But we don't actually need the solution in between; we only need to relate information at the origin and infinity. Therefore, we can connect these points by simply integrating the equation definitely; applying $\int_0^\infty \mathrm d b$ we obtain
$$-I'(0)= -2 I(\infty)^2 \quad \implies \quad I(\infty) = \frac{\sqrt \pi} 2$$
A: This is an old favorite of mine.
Define $$I=\int_{-\infty}^{+\infty} e^{-x^2} dx$$ 
Then $$I^2=\bigg(\int_{-\infty}^{+\infty} e^{-x^2} dx\bigg)\bigg(\int_{-\infty}^{+\infty} e^{-y^2} dy\bigg)$$  
$$I^2=\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}e^{-(x^2+y^2)} dxdy$$  
Now change to polar coordinates
$$I^2=\int_{0}^{+2 \pi}\int_{0}^{+\infty}e^{-r^2} rdrd\theta$$  
The $\theta$ integral just gives $2\pi$, while the $r$ integral succumbs to the substitution $u=r^2$  
$$I^2=2\pi\int_{0}^{+\infty}e^{-u}du/2=\pi$$ 
So $$I=\sqrt{\pi}$$ and your integral is half this by symmetry
I have always wondered if somebody found it this way, or did it first using complex variables and noticed this would work.
A: No entender mucho el idioma pero aqui dejo mi solucion en una imagen didactica. 
A: Another way is to make use of the Poisson summation formula. I will work with the Fourier transform $$\hat{f}(\xi) = \int_{-\infty}^{\infty} f(x) \exp(-2 \pi i \xi x) dx$$ The Poisson summation formula states that $$\sum_{\xi \in \mathbb{Z}} \hat{f}(\xi) = \sum_{n \in \mathbb{Z}} f(n).$$
Now take $f(x) = \exp(-\pi x^2)$. We then get that
\begin{align}
\hat{f}(\xi) & = \int_{-\infty}^{\infty} \exp(- \pi x^2) \exp(-2 \pi i \xi x) dx\\
& = \int_{-\infty}^{\infty} \exp(-\pi(x+i \xi)^2 - \pi \xi^2) dx\\
& = \exp( - \pi \xi^2)\int_{-\infty}^{\infty} \exp(-\pi(x+i \xi)^2) dx\\
& = \exp( - \pi \xi^2)\int_{-\infty + ic}^{\infty+ic} \exp(-\pi x^2) dx
\end{align}
By integrating from $-\infty+ic$ to $\infty + ic$, I mean integrate along the line Im$(x) = c$ from left to right. Now since the integrand is analytic, we can move this contour to $X$ axis and conclude that
$$\int_{-\infty + ic}^{\infty+ic} \exp(-\pi x^2) dx = \int_{-\infty}^{\infty} \exp(-\pi x^2) dx$$
Hence, we get that $$\hat{f}(\xi) = C \exp( - \pi \xi^2)$$where $C = \displaystyle \int_{-\infty}^{\infty} \exp(-\pi x^2) dx$. Now make use of the Poisson summation formula to get that
$$C \left(\sum_{\xi \in \mathbb{Z}} \exp(-\pi \xi^2) \right) = \left(\sum_{x \in \mathbb{Z}} \exp(-\pi x^2) \right)$$
We can afford to cancel $\displaystyle \left(\sum_{x \in \mathbb{Z}} \exp(-\pi x^2) \right)$ since it converges  and hence we can conclude that $$C = \displaystyle \int_{-\infty}^{\infty} \exp(-\pi x^2) dx = 1$$
Suitable scaling gives you the integral and answer you are looking for.
A: Change variables. Let $z=x^2$. 
We find $\int_{0}^{\infty} e^{-x^2} dx = \frac{1}{2} \Gamma(\frac{1}{2}) = \frac{\sqrt{\pi}}{2}$. 
Addendum: 
Setting $z=1/2$ in Euler's reflection formula,
$\Gamma(1-z)\Gamma(z) = \pi/\sin \pi z$,
we find $\Gamma(1/2) = \sqrt{\pi}$. 
A: Consider the integral: 
$$
\int_0^\infty t^{-1/2}{e^{-t}} dt
$$
We perform a change of variables $u=t^{1/2}$ and $du= \frac{1}{2}t^{-1/2} dt$.
The integral then becomes:
$$
\int_0^\infty t^{-1/2}{e^{-t}} dt=\int_0^\infty 2{e^{-u^2}} du.  
$$ 
Now let us consider the well-known integral:
$$ 
\frac{\pi}{2}=\int_0^\infty \frac{1}{1+x^2} dx
$$ 
We can expand the right hand side into a double integral:
$$
\int_0^\infty \frac{1}{1+x^2} dx= \int_0^\infty \int_0^\infty  e^{-y(1+x^2)}dy dx=\int_0^\infty \int_0^\infty  e^{-y-yx^2}dy dx
$$
Reversing the order of integration: 
$$
\int_0^\infty \int_0^\infty  e^{-y-yx^2}dy dx=\int_0^\infty \int_0^\infty  e^{-y-yx^2}dx dy 
$$
Now, we can perform a change of variables $x^2=\frac{u^2}{y}$ and $2xdx=\frac{2u}{y}du$ or $dx=y^{-1/2}du$ 
$$
\int_0^\infty \int_0^\infty  e^{-y-yx^2}dx dy=\int_0^\infty \int_0^\infty  y^{-1/2} e^{-y-u^2}du dy= \int_0^\infty y^{-1/2} e^{-y} dy\int_0^\infty e^{-u^2}du     
$$
Because of what was established earlier:
$$\int_0^\infty y^{-1/2}{e^{-y}} dy=\int_0^\infty 2{e^{-u^2}} du $$
$$
\frac{\pi}{2}=\int_0^\infty \frac{1}{1+x^2} dx= 2 \left(\int_0^\infty {e^{-u^2}} du\right)^2
$$
Thus,
$$
\frac{\pi}{4}=\left(\int_0^\infty {e^{-u^2}} du\right)^2
$$
The desired result follows upon taking the square root of both sides.
Remarks


*

*The integral considered at the start of the solution is $\Gamma \left(\frac{1}{2} \right).$

*We essentially evaluated $$\int_{0}^{\infty} \frac{1}{1+x^2} \ dx$$ in two different ways; we know the closed form recognizing it is an arctangent integral, but the crux of the proof is to show it is the same as $\frac{\Gamma(\frac{1}{2})^2}{2}.$ This same arctangent integral appears in proofs 2,3, and 4 of http://www.math.uconn.edu/~kconrad/blurbs/analysis/gaussianintegral.pdf. 

*Tonelli's Theorem enables us to reverse the order of integration in the expanded double integral. See https://en.wikipedia.org/wiki/Fubini%27s_theorem.

