Question on integrating a Gaussian function 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.
 A: $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:

Exercise:  Compute
  $$\int_{0}^{1} x^2\,\mathrm{d}x
\qquad\text{and}\qquad
\int_{0}^{1} y^2\, \mathrm{d}y.
$$
  What do you notice about these two quantities?

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.
A: The fact that $\int_{-\infty}^\infty e^{-x^2} \, \mathrm{d}x = \int_{-\infty}^\infty e^{-y^2} \, \mathrm{d}y$ implies absolutely nothing about any relationship, or lack thereof, between $x$ and $y$. They are, as you've heard, dummy variables; they exist as a means to turn a function $x \mapsto e^{-x^2}$ into a definite integral, in this case, the number $\sqrt{\pi}$. There's no hidden assumptions here that $x$ and $y$ are related, or that if plotted, they must be plotted parallel, or anything like that. The definition of integration simply does not rely on any such notions; it's just a way to turn certain functions over certain sets into numbers.
I'm wondering if your confusion might stem from substitution rule, and having to change the boundaries. For example,
$$\int_0^{\frac{\pi}{2}} 2\sin(2x) \, \mathrm{d}x = \int_0^{\pi} \sin(2x) \, \mathrm{d}(2x) = 2.$$
I'm going out on a limb here, but I'm guessing you've seen examples like the above before, seeing how, when changing from $x$ to $2x$, the limits on the integral had to be changed? Surely if $x$ ranges from $-\infty$ to $\infty$, the corresponding bounds for $y$ must be $0$ to $0$, since they are orthogonal? (I apologise if this isn't on the mark, but it can be really hard to figure out exactly what's confusing a student!)
Well, in response to that particular reasoning, I'd point out that, when performing a substitution, you're changing the variable relative to the integrand. By absorbing the $2$ into the $\mathrm{d}x$, suddenly we're integrating $\sin$ instead of $x \mapsto 2\sin(2x)$. These are different functions we are integrating, not merely the same function, expressing it with a different dummy variable. I would agree, for example, that
$$\int_0^{\frac{\pi}{2}} 2\sin(2x) \, \mathrm{d}x = \int_0^{\frac{\pi}{2}} 2\sin(x) \, \mathrm{d}\left(\frac{x}{2}\right),$$
which is just a change in dummy variable. Note that the bounds haven't changed (yet). But, if you were to modify the differential, by pulling the factor of $\frac{1}{2}$ out, then that's when you change the bounds.
