Geometric sum of shape-two gamma variables ($S\stackrel{d}{=}\mathrm{Exp}(a)+\mathrm{Exp}(b)$)

Let the summands have gamma distribution with shape two and rate $\beta$, and let the independent positive-support geometric distribution count have parameter $p$. Put $a=\beta (1-\sqrt{1-p})$ and $b=\beta (1+\sqrt{1-p})$. The aggregate moment-generating function factors as $ab/((a-t)(b-t))$, so it is the sum of two independent exponential distributions of rates $a,b$. Its probability density function is $ab(e^{-as}-e^{-bs})/(b-a)$ for $s>0$, nonnegative and integrating to one. This contrasts with the geometric sum of exponential variables, which is itself exponential.