I'm confused on the thinking behind one of the steps for integrating a Gaussian. Namely, the following:
$$I = \int_{-\infty}^\infty e^{-x^2}dx$$ $$I^2 = \left(\int_{-\infty}^\infty e^{-x^2}dx \right) \cdot \left(\int_{-\infty}^\infty e^{-y^2}dy \right)$$
From here, I know how it's technically done. However, I don't see how this connection is fair.
To me, it should be:
$$I^2 = \left(\int_{-\infty}^\infty e^{-x^2}dx \right) \cdot \left(\int_{-\infty}^\infty e^{-x^2}dx \right)$$
I've heard it argued that it's due to the fact that $x$ and $y$ here are arbitrary variables, so it's fine in this case. But, we use this to turn $x^2 + y^2$ to $r^2$, which implies that $x^2$ and $y^2$ are fundamentally different. Namely, orthogonal.
So I don't see how this is fair to do. We require the implication that $x^2 = y^2$ to compute the integral this way but that simply can't be as we contradict that argument by acknowledging orthogonality when switching to polar.
$I$ is a number. You obtain $I$ by integrating the Gaussian over the real line. That is, $$ I := \int_{\mathbb{R}} \mathrm{e}^{-x^2}\,\mathrm{d}x. $$ Notice that after the integration, the variable $x$ is gone. It plays no role, and is just a name for a quantity which varies over $\mathbb{R}$ in the integration. Before going any farther, perhaps it would be useful to consider the following exercise:
What you should notice is that the name of the variable with respect to which we are integrating doesn't matter. In this sense, both $x$ and $y$ are "dummy" variables in the integration. They are not elements of $\mathbb{R}^2$ and are not orthogonal to each other. They are simply labels.
So, we know how $I$ is defined, but the name of the variable doesn't matter. Hence we can write $$ I = \int_{\mathbb{R}} \mathrm{e}^{-y^2}\,\mathrm{d}y $$ just as easily as the original definition. Or even $$ I =\int_{\mathbb{R}} \mathrm{e}^{-s^2}\,\mathrm{d}s =\int_{\mathbb{R}} \mathrm{e}^{-t^2}\,\mathrm{d}t.\tag{1} $$ The actual names are completely irrelevant. Taking this last version, we can write \begin{align*} I^2 &= \left( \int_{\mathbb{R}} \mathrm{e}^{-s^2}\,\mathrm{d}s \right)^2 \\ &= \left( \int_{\mathbb{R}} \mathrm{e}^{-s^2}\,\mathrm{d}s \right)\left( \int_{\mathbb{R}} \mathrm{e}^{-s^2}\,\mathrm{d}s \right) \\ &= \left( \int_{\mathbb{R}} \mathrm{e}^{-s^2}\,\mathrm{d}s \right)\left( \int_{\mathbb{R}} \mathrm{e}^{-t^2}\,\mathrm{d}t \right) && \tag{using (1)} \\ &= \int_{\mathrm{R}} \mathrm{e}^{-s^2} \left( \int_{\mathbb{R}} \mathrm{e}^{-t^2}\,\mathrm{d}t \right) \, \mathrm{d}s \tag{2} \\ &= \int_{\mathrm{R}} \left( \int_{\mathbb{R}} \mathrm{e}^{-s^2} \mathrm{e}^{-t^2}\,\mathrm{d}t\right) \, \mathrm{d}s \tag{3} \\ &= \iint_{\mathbb{R}^2} \mathrm{e}^{-(s^2+t^2)}\,\mathrm{d}(x,y). \end{align*} At (2), we are using the fact that $$ \left( \int \mathrm{e}^{-s^2} \,\mathrm{d}s\right) k = \int_{\mathbb{R}} \mathrm{e}^{-s^2}k \,\mathrm{d}s $$ for any constant $k$ which does not depend on $s$, and the fact that $$ \int_{\mathbb{R}} \mathrm{e}^{-t^2} \,\mathrm{d}t $$ is a constant (since it doesn't depend on $s$). At (3), we are playing the same game, but this time using the fact that the value of $\mathrm{e}^{-s^2}$ does not depend on $t$. Once we get the integral at (3), we can apply the Fubini-Tonelli in order to interpret the iterated integral as an integral over a region, then make a change of variables and work with polar (rather than Cartesian) coordinates
To give a simpler example, observe that (via similar arguments) $$ 1 = \int_{0}^{1} 2s \,\mathrm{d}s. $$ Now, observe that \begin{align} 1 = 1^2 &= \left( \int_{0}^{1} 2s \,\mathrm{d}s \right)^2 \\ &= \left( \int_{0}^{1} 2s \,\mathrm{d}s \right) \left( \int_{0}^{1} 2t \,\mathrm{d}t \right) \\ &= \int_{0}^{1} 2s \left( \int_{0}^{1} 2t \,\mathrm{d}t \right) \,\mathrm{ds} \\ &= \int_{0}^{1} \left( \int_{0}^{1} (2s)(2t)\,\mathrm{d}t\right)\, \mathrm{d}s \\ &= 4 \int_{0}^{1} \left( \int_{0}^{1} st \,\mathrm{d}t\right)\, \mathrm{d}s \tag{4} \\ &= 4 \int_{0}^{1} \left. \frac{1}{2}st^2 \right|_{t=0}^{1} \,\mathrm{ds} \\ &= 2 \int_{0}^{1} s\, \mathrm{d}s \\ &= 2 \cdot \left.\frac{1}{2}s^2 \right|_{s=0}^{1} \\ &= 1, \end{align} which is the expected result. Here, we can evaluate the iterated integral obtained at (4) directly (i.e. no change of variables is required), but otherwise this example follows the same principles as the argument for integrating the Gaussian.