A standard technique of mine when evaluating iterated integrals, such as double integrals, is to invoke a substitution that relates the integration variables of both integrals. To demonstrate and clarify what I mean, I shall present three examples for double integrals:
Example 1: Gaussian integral, $\Gamma\left(\frac{1}{2}\right)=\sqrt{\pi}$
$$I=\Gamma \left(\frac{1}{2}\right)^2 = \int_{0}^{\infty} \int_{0}^{\infty} \frac{1}{\sqrt{x t}}e^{-(t+x)}\,dx\,dt$$
We now enforce the substituion $t=sx \implies dt=x\,ds$:
$$\implies I = \int_{0}^{\infty}\int_{0}^{\infty}\frac{1}{\sqrt{s}}e^{-x (s+1)}\,dx\,ds=\int_{0}^{\infty}\frac{1}{\sqrt{s}(s+1)} \,ds=\pi$$
$$\implies \Gamma\left(\frac{1}{2}\right)=\int_{-\infty}^{\infty}e^{-x^2}\,dx=\sqrt{\pi}$$
This is a much more elementary approach than the common Jacobian polar coordinates method of evaluating the Gaussian integral and is how I originally discovered the solution for myself a few years ago prior to learning about the Jacobian.
Example 2: In a previous answer of mine, I determined the following result:
$$J=\frac{1}{2\pi}\int_{-\infty}^{\infty} \exp\left(\frac{x^2}{2}-\frac{y^2}{2}\right)K_0\left(\sqrt{(x-y)^2+\left(\frac{y^2}{2}-\frac{x^2}{2}\right)^2}\right)\,dy = \frac{e}{2} \sqrt{\pi}\operatorname{erfc}(1)$$ and is independent of $x$. Here $K_0$ is the modified Bessel function of the second kind of order $0$.
A crucial part of my method of evaluation was to invoke a particular substitution on the integral:
$$J=\frac{1}{2\pi}\int_0^\infty\int_0^\infty \frac{\cosh\left(xv\left(\frac{v^2}{4t}+1\right)\right)}{t}\exp\left(-t-\frac{1+2t+x^2}{4t}v^2-\frac{v^4}{16t}\right)\,dv\,dt$$
We now enforce the substitution $t = s v^2 \implies dt = v^2 \,ds$:
$$\begin{align}J&=\frac{1}{2\pi} \int_{0}^{\infty} \int_{0}^{\infty}\frac{\cosh\left(x v \left(\frac{1}{4s}+1\right)\right)}{s} \exp\left(-sv^2-\frac{1+2sv^2+x^2}{4s}-\frac{v^2}{16s}\right) \, dv \, ds\\&=\frac{1}{2\pi}\int_{0}^{\infty}\int_{0}^{\infty}\frac{\cosh\left(xv\left(\frac{1}{4s}+1\right)\right)}{s}\exp\left(-\frac{1+x^2}{4s}\right)\exp\left(-v^2 \frac{(1+4s)^2}{16s}\right)\,dv\,ds\end{align}$$
This allows the result to fall out readily, since:
$$\int_{0}^{\infty}\cosh(a v)\exp\left(-v^2 b^2\right) \, dv=\frac{\sqrt{\pi}}{2b}\exp\left(\frac{a^2}{4b^2}\right)$$
$$\implies J = \frac{1}{2\pi} \int_{0}^{\infty} \frac{\exp \left(-\frac{1+x^2}{4s}\right)}{s} \cdot \frac{\sqrt{\pi}}{2\left(\frac{1+4s}{4\sqrt{s}}\right)}\exp \left(\frac{4s x^2 \left(\frac{1}{4s}+1\right)^2}{(1+4s)^2}\right)\, ds$$
This results in the integral being independent as we wanted since the $x^2$ terms cancel.
$$\implies J = \frac{1}{\sqrt{\pi}} \int_{0}^{\infty} \frac{\exp\left(-\frac{1}{4s}\right)}{\sqrt{s} (1+4s)} \, ds\stackrel{s\,\mapsto\frac{1}{s}}{=}\frac{1}{\sqrt{\pi}} \int_{0}^{\infty} \frac{\exp\left(-\frac{s}{4}\right)}{\sqrt{s} (4+s)} \, ds=\frac{e}{2} \sqrt{\pi} \, \operatorname{erfc} (1)$$
Example 3: Derivation of the relationship between the Beta function and Gamma function. $\Gamma(x)\Gamma(y)=\Gamma(x+y)B(x,y)$:
$$\Gamma(x)\Gamma(y)=\int_{0}^{\infty}\int_{0}^{\infty} e^{-u-v}u^{x-1}v^{y-1}\,du\,dy$$
We now enforce the substitution $u=zt$ and $v=z(1-t)$:
$$\begin{align}\Gamma(x)\Gamma(y)&=\int_{0}^{\infty}\int_{0}^{1}e^{-z}(zt)^{x-1}(z(1-t))^{y-1} z\,dt\,dz\\ &=\int_{0}^{\infty}e^{-z}z^{x+y-1}\,dz \cdot \int_{0}^{1}t^{x-1}(1-t)^{y-1}\\&=\Gamma(x+y)B(x,y)\end{align}$$
This is interestingly the method mentioned in the derivation of this property on the Wikipedia page of the Beta function. In this question about this derivation, the comments and answer invoke the Jacobian, suggesting that such substitutions can be made rigorous using the Jacobian, however, I would appreciate seeing exactly how that works as I only know about the Jacobian based on my own self-teaching.
So to summarise is it always valid to perform substitutions on iterated integrals, such as double or triple integrals in this way, whereby one or more integration variables are written in terms of an integration variable of a different integral in the iteration. I have never come across any issue in using this technique in the past, but have never found rigorous proof in any literature that states that such substitutions are valid; it surprises me that such a simple solution to the Gaussian integral is not often used.
Best Answer
$\newcommand{\d}{\mathrm{d}}$The general substitution theorem is of this form:
$$\int_{\varphi(\Omega)}f\,d\mu=\int_{\Omega}(f\circ\varphi)\cdot|\det\varphi'|\,\d\mu$$
With hypotheses (for some $n$): $\Omega\subseteq\Bbb R^n$ is open, $\varphi$ is a $C^1$ diffeomorphism between $\Omega$ and its image in $\Bbb R^n$, and $f$ is Lebesgue measurable. $\mu$ denotes $n$-dimensional Lebesgue measure (so, just "normal" integration). A good reference for this is peek-a-boo's great answer.
There is actually not all that much to say here. You pass an iterated integral into a product integral via Fubini’s theorem - since the measure space under consideration is complete and $\sigma$-finite, the only consideration is that $f(x,y)$ should be integrable in $x$ for almost every $y$. In almost every “normal” case you deal with this shall be true. From the product integral, you can use the substitution formula. Any of these substitutions which are intertwining two or more of the variables of integration are fine, so long as the substitution function satisfies the above (or some other equivalent) hypotheses. I'll show this for your examples, but for the sake of getting the rigorous argument you do not actually need to read any further.
For simplicity consider two dimensions. You're integrating functions $f(x,t)$ over some such open set in $\Bbb R^n$, in both examples $(0,\infty)\times(0,\infty)$.
Example 1:
Example 2: Exactly the same story
Example 3: ... exactly the same story essentially. You only need to ensure a inverse exists and that the inverse is also continuously differentiable (the inverse function theorem will usually be helpful for this, if the derivative is never null, i.e. the determinant is never $0$). We are integrating $f(u,v)=e^{-(u+v)}u^{x-1}v^{y-1}$. Define $\varphi(z,t)=(zt,z(1-t))$. $\varphi((0,\infty)\times(0,1))=(0,\infty)\times(0,\infty)$ is not too hard to see, and some simple algebra shows it is indeed invertible, the inverse is also $C^1$ (or use the inverse function theorem). You have: $$\int_{\varphi((0,\infty)\times(0,1))}f(u,v)\,\d u\,\d v=\int_{(0,\infty)\times(0,1)}f(\varphi(z,t))\cdot|\det\varphi'|\,\d t\d z$$We get: $$\varphi'=\begin{pmatrix}t&z\\1-t&-z\end{pmatrix},\,|\det\varphi'|=|-zt-z+zt|=z$$So your integral, entirely formally, goes to: $$\int_0^\infty\int_0^1f(zt,z(1-t))\cdot z\,\d t\,\d z$$
And you proceed.