We square an integral, but why change a variable? 
If we square an integral, we also change the integration variable in one of the integrals. But why is this actually correct?

For example, say I have the following:
Solve $\int_{-\infty}^\infty e^{-x^2} dx$. Let $I=\int_{-\infty}^\infty e^{-x^2} dx$, so
\begin{align}
I^2 &=\bigg(\int_{-\infty}^\infty e^{-x^2} dx\bigg)^2\\
&= \bigg(\int_{-\infty}^\infty e^{-x^2} dx\bigg) \times 
\bigg(
\underbrace{
\int_{-\infty}^{\infty} e^{-y^2}dy
}_{\text{Why}?}
\bigg) 
\\
&=\int_{-\infty}^\infty\bigg(\int_{-\infty}^\infty e^{-(x^2+y^2)} dx\bigg)dy 
\end{align}
But why is the following wrong:
\begin{align}
I^2 &=\bigg(\int_{-\infty}^\infty e^{-x^2} dx\bigg)^2\\
&= \bigg(\int_{-\infty}^\infty e^{-x^2} dx\bigg) \times \bigg(\int_{-\infty}^{\infty} e^{-x^2}dx\bigg) \\
&=\int_{-\infty}^\infty\bigg(\int_{-\infty}^\infty e^{-x^2-x^2} dx\bigg)dx \\
&=\int_{-\infty}^\infty\bigg(\int_{-\infty}^\infty e^{-2x^2} dx\bigg)dx \qquad ?
\end{align}
 A: The short answer is that in your second to last identity you are neglecting cross terms in your multiplication.
A simple example may be used to demonstrate the error, and we will use a strict summation instead of an integration for clarity.  Consider the sum
\begin{equation}
  \sum_{x=1}^3 x = 1 + 2 +3 = 6.
\end{equation}
Now, squaring the summation yields
\begin{align}
  \left(\sum_{x=1}^3x\right)^2 &= (1+2+3)^2 = (1+2+3)\times(1+2+3)
\end{align}
In order to properly calculate this quantity (long-hand) requires one to completely distribute including cross terms:
\begin{align}
  &\ {\color{white}+}\ 1\times 1 + 1\times 2 + 1\times 3 \notag \\
  (1+2+3)\times(1+2+3) =& +2\times 1 + 2\times 2 + 2\times 3 = 36. \notag \\
  & + 3\times 1 + 3\times 2 + 3\times 3
\end{align}
However, if the cross terms were neglected one would obtain
\begin{equation}
  (1+2+3)\times(1+2+3) \ne (1^2+2^2+3^2) = 1\times 1 + 2\times 2 + 3\times 3 = 14.
\end{equation}
Written another way, we conclude that
\begin{equation}
  \left(\sum_{x=1}^3x\right)^2 = \left(\sum_{x=1}^3x\right)\left(\sum_{x=1}^3x\right) \ne \sum_{x=1}^3x^2.
\end{equation}
However, as noted in the comments, the $x$ in the summation (or the in the definite integral) is just a dummy variable, and may be replaced with another symbol, such as $y$:
\begin{equation}
  \left(\sum_{x=1}^3x\right)^2 = \left(\sum_{x=1}^3x\right)\left(\sum_{y=1}^3y\right).
\end{equation}
The reason for changing the dummy variable is that we may now say that the summation over $x$ and $y$ are independent, and so we may rearrange the summations:
\begin{equation}
  \left(\sum_{x=1}^3x\right)^2 = \left(\sum_{x=1}^3x\right)\left(\sum_{y=1}^3y\right) = \sum_{x=1}^3\sum_{y=1}^3x\times y.
\end{equation}
Again, integration is just summation, and so the same fact holds true for integrals.  Furthermore, it doesn't matter what the function inside the summation/integral is.  Thus, we may write
\begin{equation}
  \left(\int_a^b f(x)\ dx\right)^2 = \left(\int_a^b f(x)\ dx\right)\left(\int_a^b f(y)\ dy\right) = \int_a^b\int_a^b f(x)\times f(y)\ dx\ dy.
\end{equation}
Using $f(x)=e^{-x^2}$ results in your example problem.
A: Because changing the integrand's variable doesn't change the value of the integral:$$\int_{a}^{b}f(x)dx=\int_{a}^{b}f(y)dy$$the same argument holds for summation as following$$\sum_{n=a}^{b}k_n=\sum_{m=a}^{b}k_m$$also we know that$$\left(\sum_{n=a}^{b}k_n\right)^2=\left(k_a+...+k_b\right)^2=k_a^2+...+k_b^2+2k_ak_{a+1}+...+2k_{b-1}k_b$$and$$\sum_{n=a}^{b}k_n\sum_{m=a}^{b}k_m=(k_a+...+k_b)(k_a+...+k_b)$$which by expanding the terms and rearranging them leads to$$\left(\sum_{n=a}^{b}k_n\right)^2=\sum_{n=a}^{b}k_n\sum_{m=a}^{b}k_m$$since integral is intrinsicly a summation, this can be generalized to integral operator either.
A: I think the best way of visualizing what is happening intuitively is thinking in terms of programming. 
If you make a program and you have two functions (in the programming sense) which do not interact, you can call "$x$" the variable in both functions. They are isolated.
However, if a function is supposed to handle two different variables independently in the same context, you cannot name them the same thing.
Formally and mathematically (and also changing notation from $\int f(x)dx$ to $\int f$ so that the point may become clearer), we have the following chain of equalities.
\begin{align*}
\underline{\left(\int_{\mathbb{R}}  x \mapsto e^{-x^2}\right)} \cdot \underline{\left(\int_{\mathbb{R}}  x \mapsto e^{-x^2} \right)}&\stackrel{Linearity}=\left(\int_{\mathbb{R}} \underline{\left(\int_{\mathbb{R}}  x \mapsto e^{-x^2}\right)}\left(x \mapsto e^{-x^2} \right)\right) \\
&\underset{\text{fctn mult.}}{\overset{\text{Def. of}}{=}} \int_{\mathbb{R}}  x \mapsto \left(e^{-x^2} \underline{\left(\int_{\mathbb{R}}  y \mapsto e^{-y^2}\right)}\right) \\
&\underset{\text{}}{\overset{\text{Linearity}}{=}} \int_{\mathbb{R}} x \mapsto \left(\int_{\mathbb{R}} y \mapsto e^{-x^2-y^2}\right) \\
&\stackrel{Fubini}{=}\int_{\mathbb{R}^2}(x,y) \mapsto e^{-x^2-y^2}.
\end{align*}
Using the same name ($x$) for the variables in both integrations is "fine" up to the third equality (I changed it in the second equality for the sake of clarity, otherwise the third would be a little mystic). Using $x$ instead of $y$ in the second line would be extremely bad taste, but not explicitly wrong. Using $x$ instead of $y$ in the right side of the first equality is edgy, but not so much.
Indeed, each underline integral is a closed box with respect to the rest: they are fixed numbers, which don't interact at all with the rest of the environment other than via the fact that they are numbers. It is in the third equality that we are multiplying a constant (relative to the inner function which we are integrating) which depends on $x$. It does not depend on the function inside the integral at all, but it is invading its space and interacting with it. To be extremely clear, call the function $x \mapsto e^{-x^2}$ by $f$. The second line is then
$$\int_{\mathbb{R}}  x \mapsto \left(e^{-x^2} \left(\int_{\mathbb{R}}  f\right)\right) $$
$$=\int_{\mathbb{R}} \left(  x \mapsto  \int_{\mathbb{R}} (e^{-x^2}\cdot f)\right),$$
where we multiplied a constant to the inner integral.
What you are proposing is that the equality should read
$$\int_{\mathbb{R}}  x \mapsto \left(e^{-x^2} \left(\int_{\mathbb{R}}  f\right)\right) $$
$$=\int_{\mathbb{R}} \left(  x \mapsto  \int_{\mathbb{R}} g\right),$$
where $g(t)=e^{-t^2}f(t)$. This doesn't make sense. Essentially, what you are implying is that
$$f \cdot\int g=\int f\cdot g,$$
which is not true.
