prove Gaussian integral using polar coordinates 
The proof method is to equate expression$\mathrm{\iint_{-\infty}^\infty\,e^{-(x^2+y^2)}}$ (Cartesian)with $\mathrm{\int_0^{2\pi}\int_0^{\infty}e^{-r^2}drd\theta}$(polar)
however, the answer goes into great length to prove that the integral is bounded.

to evaluate $\mathrm{\lim_{a \to \infty}\int_{-a}^ae^{-x^2}dx\int_{-a}^ae^{-y^2}dy}$ we are using the fact that these integrals are bounded...

why does the integral needs to be bounded? what would happen if it is not bounded? 
 A: In my view, the posted proof (from a solution manual?) leaves something to be desired in terms of mathematical care. In lieu of a point-by-point critique, here's an independent conceptual overview about why care is needed with this type of question.

When you write an improper double integral
$$
\iint_{\mathbf{R}^{2}} f(x, y)\, dx\, dy,
\tag{1}
$$
you're implicitly making a big assumption, that if $D_{a}$ is an "increasing" family of closed, bounded sets that "converges to the entire plane",[1] then
$$
\lim_{a \to \infty} \iint_{D_{a}} f(x, y)\, dx\, dy
$$
exists and has the same value independently of the family $D_{a}$.
Here, $f(x, y) = e^{-(x^{2} + y^{2})}$, a positive, "rapidly-decaying" function, and there are two particular families of sets: the closed disks $D_{a} = \{(x, y) : x^{2} + y^{2} \leq a^{2}\}$, and the closed squares $D_{a}' = [-a, a] \times [-a, a]$. The integrals over $D_{a}$ can be calculated easily in polar coordinates. The integrals over $D_{a}'$ are conveniently interpreted as squares of truncated one-variable Gaussian integrals. To make a mathematical deduction about the Gaussian, it's important to prove the starred equality in
\begin{align*}
2\pi = \lim_{a \to \infty} 2\pi\int_{0}^{a} re^{-r^{2}}\, dr
  &= \lim_{a \to \infty} \iint_{D_{a}} e^{-(x^{2} + y^{2})}\, dx\, dy \\
  &\stackrel{*}{=} \lim_{a \to \infty} \iint_{D_{a}'} e^{-(x^{2} + y^{2})}\, dx\, dy \\
  &= \lim_{a \to \infty} \left(\int_{-a}^{a} e^{-x^{2}}\, dx\right)^{2}
   = \left(\int_{-\infty}^{\infty} e^{-x^{2}}\, dx\right)^{2}.
\end{align*}
There are many ways to do this, boiling down to "absolute convergence" of the improper integral. For example, on $D_{a}'' := D_{a}' \setminus D_{a}$, the portion of the square outside the disk, the integrand is bounded by $e^{-a^{2}}$, while the area of $D_{a}''$ is $(4 - \pi) a^{2} \leq a^{2}$; consequently, for each positive real number $a$ we have
\begin{align*}
0 &\leq \iint_{D_{a}''} e^{-(x^{2} + y^{2})}\, dx\, dy \\
  &= \iint_{D_{a}'} e^{-(x^{2} + y^{2})}\, dx\, dy - \iint_{D_{a}} e^{-(x^{2} + y^{2})}\, dx\, dy \\
  &\leq a^{2} e^{-a^{2}},
\end{align*}
and the upper bound converges to $0$ as $a \to \infty$, proving the two families of integral have the same limit.

If you try to run this argument with a suitably-chosen double integral that converges non-absolutely, and with suitable sets $D_{a}$ and $D_{a}'$, it may turn out that the difference
$$
\iint_{D_{a}'} f(x, y)\, dx\, dy - \iint_{D_{a}} f(x, y)\, dx\, dy
$$
does not converge at all, much less converge to $0$. For example, letting $k$ denote a positive integer, if
$$
f(x, y) = \begin{cases}
   \frac{1}{2k-1} & \text{if $k - 1 \leq x < k$ and $-1 \leq y < 0$,} \\
  -\frac{1}{2k} & \text{if $k - 1 \leq y < k$ and $-1 \leq x < 0$,} \\
  0 & \text{otherwise},
\end{cases}
$$
and if
$$
D_{a} = [-a, a] \times [-a, a],\qquad
D_{a}' = [-a, a] \times [-2a, 2a],
$$
then for positive integer $N$,
$$
\iint_{D_{N}} f(x, y)\, dx\, dy = \sum_{k=1}^{2N} \frac{(-1)^{k}}{k},\quad
\iint_{D_{N}'} f(x, y)\, dx\, dy
  = \sum_{k=1}^{N} \frac{1}{2k - 1} - \sum_{k=1}^{2N} \frac{1}{2k}.
$$
These are rearranged partials sums of the alternating harmonic series, well-known not to have the same limit. It follows that
$$
\iint_{\mathbf{R}^{2}} f(x, y)\, dx\, dy
$$
is not well-defined.



*

*Precisely, a family of closed, bounded sets $D_{a}$ depending on a positive real number $a$ converges to the entire plane if "for every closed, bounded set $K$, there exists a real number $M$ such that $K \subset D_{a}$ for all $a \geq M$".

