Sum of generalized gamma variables

Question

This question is in a way duplicating this question, but I'm not happy with the final answer and feel like I need a deeper dive.

Assume that $X_j \sim \Gamma(\alpha_j,\beta_j)$, where $\alpha, \beta$ denotes shape and rate.

I'm interested in deriving the form of $GG(\gamma, \alpha, \beta)$ of the following sum, where $\gamma$ is the form/power parameter:

$$S = \sum_{j = 1}^{d} X_j^2$$

From the previous answer I understand, that $Y_j=X_j^2$ has a Generalized Gamma distribution with $\gamma = 1/2$. It follows, that $\sum_{j=1}^{d} Y_j^{1/2}$ has gamma distribution, therefore $S$ is also Generalized Gamma variable.

However, how would the final parameters of the distribution look like?

Assuming $\beta_j = \beta, \forall j$, it is known that $Z = \sum_{j=1}^{d} Y_j^{1/2} \sim \Gamma(\sum \alpha_j, \beta)$. Knowing $Z$, how should we proceed in deriving the form of $S$? The previous answer claims that it is straightforward from there as if it's $GG(\gamma, \sum \alpha_j, \beta)$, which doesn't feel true since $S^{1/2} \neq Z$.

Answer 1

For the sum of independent Gamma random variables with different shape and rate/scale parameters, the answer has been given in

Moschopoulos, P. G. (1985). The distribution of the sum of independent gamma random variables. Annals of the Institute of Statistical Mathematics, 37(3), 541-544.

The density can be expressed as an infinite sum involving some complicated expressions of the underlying parameters.

Perhaps this can be a base to work the sum of squared variables.