Sum of two gamma/Erlang random variables $\Gamma(m,\lambda)$ and $\Gamma(n, \mu)$ with integer numbers $m \neq n, \lambda \neq \mu$

A closed form expression is provided in the following paper.

SV Amari, RB Misra, Closed-form expressions for distribution of sum of exponential random variables, IEEE Transactions on Reliability, 46 (4), 519-522.

Update:

We summarize the development in follows:

Step 1: We simplify the case to be $l(\Gamma(m, 1)+k\Gamma(n,1))$ by choosing appropriate $k,l$ for a scale transformation.

Step 2. We want to calculate $$ \Gamma(m,1)+k\Gamma(n,1)=\sum^{m}_{i=1}X_{i}+k\sum^{n}_{i=1}Y_{j} $$

Step 3. For $m=n$ case, we only need to calculate $X+kY$, where $X,Y\sim \Gamma(1,1)$. In the case of $k=1$, we let $$ Z=X+Y,W=\frac{X}{X+Y}, X=ZW, Y=Z-ZW $$ The Jacobian is $Z$. Therefore we have $$ f_{Z,W}(z,w)=ze^{-zw}e^{zw-z}dzdw=ze^{-z}dzdw $$ and $$ Z\sim \Gamma(2,1) $$ as desired. In the general case we have $$ Z=X+kY, W=\frac{X}{X+kY}, X=ZW, Y=\frac{1}{k}(Z-ZW) $$ The Jacobian is $\frac{1}{k}Z$. We thus have $$ f_{Z,W}(z,w)=\frac{1}{k}ze^{-zw}e^{\frac{1}{k}(zw-z)}=\frac{1}{k}ze^{-\frac{k-1}{k}zw-\frac{1}{k}z} $$ and I do not have a good way to factorize it.

A reason this technique might not work in general is the moment generating function does not change when we use the scale transformation, and for different $\beta$ the moment generating function is different. Thus the problem may be better to be attacked numerically.