Let $X$ be a Zero-modified Geometric Distribution with parameters $\mu > 0$ and $\pi \in (-1/\mu, 1)$, denoted by $X \sim \text{ZMG}(\pi, \mu)$. $X$ has a p.m.f. given by $$\Pr(X = k)= \left\{ \begin{array}{ll}\pi + (1-\pi)\left(\frac{1}{1+\mu}\right), \; \text{for} \; k=0\\ (1-\pi)\frac{\mu^{k}}{(1+\mu)^{k+1}}, \; \text{for} \; k = 1, 2, \ldots . \end{array}\right.$$
Consider the probability generating function (p.g.f.) of $X$, denoted by $\varphi_{X}(s) := E(s^{X})$. Note that $$\varphi_{X}(s) = \frac{1+\pi\mu(1-s)}{1+\mu(1-s)},$$ for $|s|<\dfrac{1+\mu}{\mu}$.
Now let $\left\{X_i\right\}_{i=1}^{\infty}\stackrel{iid}{\sim} \text{ZMG}(\pi, \mu)$.
We want to find the probability mass function of $S_n = \sum\limits_{i=1}^{n}{X_i}$.
One way we tried was to compute the p.g.f. of $S_n$ and then shows that it was a p.g.f. of some known distribution (we guess it is a Zero-modified Negative Binomial Distribution). But we weren't able to do the calculations.
\begin{equation*} \begin{aligned} \varphi_{S_n}(s) &= E(s^{S_n}) = E\left(s^{ \sum_{i=1}^{n}{X_i} }\right) = E\left( \prod_{i=1}^{n}s^{X_i} \right), \; X_i \; \text{is i.i.d.} \; \forall i\\ &= \prod_{i=1}^{n}E(s^{X}) = (\varphi_{X}(s))^{n} = \left(\dfrac{ 1+\pi\mu(1-s) }{ 1+\mu(1-s) }\right)^{n}. \end{aligned} \end{equation*}
We failed in manipulated the expression above. If we were able to rewrite $\varphi_{S_n}(s)$ in something like this:
$$\left(\dfrac{ 1+\pi\mu(1-s) }{ 1+\mu(1-s) }\right)^{n}=\sum_{k=0}^{\infty}{s^k f(\pi, \mu, n, k}),$$
$f$ could indicate the distribution (the p.m.f.) we want to find.
Any suggestion would be helpful. Thank you!