Why does combining $\int_{-\infty}^{\infty} e^{-x^2} dx\int_{-\infty}^{\infty} e^{-y^2} dy$ into a $dx\,dy$ integral work?

$$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} e^{-x^2}e^{-y^2} dx dy$$

Because we treat $y$ constant in the first, inner, integral we can pull it out.

$$=\int_{-\infty}^{\infty} e^{-y^2} \int_{-\infty}^{\infty} e^{-x^2} dx dy$$

Now because $\int_{-\infty}^{\infty} e^{-x^2} dx$ is some constant we can pull it out,

$$=\int_{-\infty}^{\infty} e^{-x^2} dx\int_{-\infty}^{\infty} e^{-y^2} dy$$

The result we got is generalizable, given $g(x,y)=f(x)h(y)$ we have,

$$\int_{a}^{b} \int_{c}^{d} g(x,y) dx dy=\int_{a}^{b} h(y) dy \int_{c}^{d} f(x) dx$$

Provided everything we write down exists.

Fubini-Tonelli theorem implies that if $\int_{\mathbb{R}} dx \left(\int_{\mathbb{R}} dy |f(x,y)|\right) < \infty$ then $$\int_{\mathbb{R}^2} f(x,y)dxdy = \int_{\mathbb{R}} dx \left(\int_{\mathbb{R}} dy f(x,y)\right).$$ First, to establish the premise, note that $f(x,y) = e^{-(x^2+y^2)},$ is positive so $|f(x,y)| = f(x,y).$ So we just need to establish the finiteness of the RHS. The inner integral is $$ \int_{\mathbb{R}} dy e^{-(x^2+y^2)} = e^{-x^2}\int_{\mathbb{R}} dy e^{-y^2}. $$ By comparing the integrand with $e^{-|y|}$ which can be integrated to 2 by elementary means, we see that $$\int_{\mathbb{R}} e^{-y^2}dy = C < \infty$$ so that $$\int_{\mathbb{R}} dx \left(\int_{\mathbb{R}} dy e^{-(x^2+y^2)}\right) = C\int_{\mathbb{R}} dx e^{-x^2} = C^2 < \infty.$$

Thus the theorem gives $$ \int_{\mathbb{R}^2} dxdy e^{-(x^2+y^2)} = C^2 $$ where $$ C = \int_{\mathbb{R}} dx e^{-x^2}.$$

The two-dimensional integral can then be transformed to polar coordinates by the change of variables theorem, which (by Fubini again) is a doable iterated integral that comes out to $\pi.$

The integral $\displaystyle \int_{-\infty}^\infty e^{-y^2} \,dy$ does not change as $x$ goes from $-\infty$ to $\infty$, so you have \begin{align} \int_{-\infty}^\infty e^{-x^2} \, dx \cdot \int_{-\infty}^\infty e^{-y^2} \, dy & = \int_{-\infty}^\infty e^{-x^2} \, dx \cdot \text{constant} \\[10pt] & = \int_{-\infty}^\infty \left( e^{-x^2}\cdot\text{constant} \right)\,dx \\[10pt] & = \int_{-\infty}^\infty \left( e^{-x^2} \int_{-\infty}^\infty e^{-y^2}\,dy \right)\,dx \end{align} The factor $e^{-x^2}$ does not change as $y$ goes from $-\infty$ to $\infty$, so you now have \begin{align} \int_{-\infty}^\infty \left( e^{-x^2} \int_{-\infty}^\infty e^{-y^2}\,dy \right) \, dx &= \int_{-\infty}^\infty \left( \text{constant} \cdot \int_{-\infty}^\infty e^{-y^2}\,dy \right) \,dx \\[10pt] & = \int_{-\infty}^\infty \left( \int_{-\infty}^\infty \text{constant} \cdot e^{-y^2}\,dy \right) \,dx \\[10pt] & = \int_{-\infty}^\infty \left( \int_{-\infty}^\infty e^{-x^2} \cdot e^{-y^2}\,dy \right) \,dx \end{align}