Showing that $10!=6!7!$ via Gamma function

It is well-known that $10! = 6! 7! $, I want to prove it via algebraic manipulations of Gamma function, i.e show that: $$ \int_{0}^{+\infty} x^7 e^{-x} \, dx \int_{0}^{+\infty} y^6 e^{-y} \, dy = \int_{0}^{+\infty} z^{10} e^{-z} \, dz \qquad (1) $$ I tried to write the LHS as:

$$ \int_{0}^{+\infty} \int_{0}^{+\infty} x (xy)^6 e^{-(x+y)} \, dx dy \qquad (2)$$

and go for the substitution $u = x + y $, $v=xy$, but that leads to complicated calculations. Is there a reasonable way/substitution to turn the double integral in $(2)$ into the RHS in $(1)$. Thanks in advance for any contribution.

Solution 1:

The goal as I understand it is to use some nifty change of variables to turn this double integral $$ LHS=\int_{0}^{+\infty} \int_{0}^{+\infty} x (xy)^6 e^{-(x+y)} \, dx dy=6!\times 7!, $$ into $$ RHS=\int_0^{+\infty} z^{10} e^{-z} \, dz=10!. $$ This appears to me to be an exercise in manipulating the product of Gamma functions or as it is otherwise known a Beta function. A common substitution made when dealing with such a double integral turns the exponential weighting into a Gaussian $$ x=s^2 , \ y=t^2 , \ \frac{\partial (x,y)}{\partial (s,t)}=4st, \\LHS=\int_{0}^{+\infty} \int_{0}^{+\infty} x (xy)^6 e^{-(x+y)} \, dx dy=4\int_{0}^{+\infty} \int_{0}^{+\infty} s^2 (st)^{13} e^{-(s^2+t^2)} \, ds dt. $$ This particular substitution is valuable because it makes it simpler to use the polar parametrization of the integration region which is desirable because it reduces the Gaussian weight function to single variable dependence $$ s=r\cos{(\theta)} ,\ t=r\sin{(\theta)} , \ \frac{\partial (s,t)}{\partial (r,\theta)}=r, \\LHS=4\int_{0}^{+\infty} e^{-r^2} \, r^{29}dr\int_{0}^{\pi/2} \cos^{15}{(\theta)} \sin^{13}{(\theta)} d\theta. $$ The $\theta$-integral is easily evaluated using the trigonometric substitution we all learn in high school $$ \int_{0}^{\pi/2} \cos^{15}{(\theta)} \sin^{13}{(\theta)} \ d\theta \\=\int_{0}^{\pi/2} \cos^{14}{(\theta)} (1-\cos^{2}{(\theta)})^6\cos{(\theta)}\sin{(\theta)} \ du, \\u=\cos^{2}{(\theta)} , \ \frac{du}{d\theta}=-2\cos{(\theta)}\sin{(\theta)}, \\=\frac{1}{2}\int_{0}^{1} u^7 (1-u)^6 \ du=1/48048. $$ Finally, the radial integral may be shown to be proportional to the RHS by using the first substitution I introduced and integrating by parts $$ r=\sqrt{z} , \ \frac{dr}{dz}=\frac{1}{2}z^{-1/2} , \\\int_{0}^{+\infty} e^{-r^2} \, r^{29}dr=\frac{1}{2}\int_{0}^{+\infty} e^{-z} \, z^{14} dz=\frac{14\times13\times12\times11}{2}\int_{0}^{+\infty} e^{-z} \, z^{10} dz, $$ hence $$ LHS=\frac{4\times14\times13\times12\times11}{2\times48048}\int_{0}^{+\infty} e^{-z} \, z^{10} dz=\int_{0}^{+\infty} e^{-z} \, z^{10} dz=RHS. $$ As required, however it is worth noting that the quickest way to prove this identity is to simply use the definition of factorial and cancel down $$ LHS=6!7!=\frac{6\times5\times4\times3\times2}{8\times9\times10}10!=10!=RHS. $$ Unsurprisingly this is mirrored by using integration by parts on the integral representation of the Gamma function.