Suppose that $X$ and $Y$ are both gamma distributed with shapes $a,b$ and unit scales / unit rates. To fix ideas, X has Laplace transform given by:
$$L_X(t) = \mathbb{E}(e^{-tX}) = (1 + t)^{-a}$$
How can we compute the Laplace transform of the independent product $Z=XY$, that is the integral:
$$L_{XY}(t) = \mathbb{E}(e^{-tXY}) = \mathbb{E}(L_Y(tX)) = \frac{1}{\Gamma(a)}\int\limits_{0}^{+\infty} (1+tx)^{-b} x^{a-1} e^{-x} \partial x \text{ ?}$$