From $$\Gamma(x)\Gamma(y)=\int_0^\infty e^{-t}t^{x-1} \left( \int_0^\infty e^{-t} s^{y-1} ds \right) dt,$$ use a change of variable $s=ut$ to show $$\Gamma(x)\Gamma(y)=\Gamma(x+y)\beta(x,y).$$
Let $s=ut$ so $ds = udt + tdu$ and then
\begin{align*} \Gamma(x)\Gamma(y) &=\int_0^\infty e^{-t}t^{x-1} \left( \int_0^\infty e^{-s} s^{y-1} ds \right) dt \\ &= \int_0^\infty e^{-t} t^{x-1} \left( \int_0^{\infty} e^{-ut} (ut)^{y-1}(udt+tdu)\right)dt \\ &= \int_0^\infty e^{-t} t^{x-1} \left( \int_0^{\infty} (e^{-ut} u^y t^{y-1})dt+ \int_0^{\infty} (e^{-ut}u^{y-1}t^y)du\right)dt \\ \end{align*} My next thought was to use integration by parts, but that didn't pan out. Any suggestions?
Note that Gamma Function for $x>0$ we define: $$\Gamma(x):=\int_0^\infty e^{-t}t^{x-1}dt.$$
And for the Beta Function for $x>0$, $y>0$, we define $$\beta(x,y):=\int_0^1 t^{x-1}(1-t)^{y-1}dt.$$
