Question: What is the density at zero for the pdf of the difference of independent gamma random variables when the sum of the shape parameters is less than or equal to one?
Motivation for asking:
Given two independent gamma random variables $X\sim\mathcal{G}(\alpha_{X},\beta_{X})$ and $Y\sim\mathcal{G}(\alpha_{Y},\beta_{Y})$, the PDFs (using the shape-rate parameterization) of $X$ and $Y$ are \begin{align} \tag{1} f_{X}(x) &=\frac{\beta_{X}^{\alpha_{X}}}{\Gamma(\alpha_{X})}x^{\alpha_{X}-1}e^{-\beta_{X}x},\\ \tag{2} f_{Y}(y) &=\frac{\beta_{Y}^{\alpha_{Y}}}{\Gamma(\alpha_{Y})}y^{\alpha_{Y}-1}e^{-\beta_{Y}y}, \end{align} where the parameters $\alpha_{X}$, $\alpha_{Y}$, $\beta_{X}$, and $\beta_{Y}$ are all greater than zero. Now define the random variable $Z=X-Y$. The pdf of $Z$, $f_{Z}$, can be found by the cross-correlation of $f_{X}$ and $f_{Y}$ to give \begin{equation} f_{Z}(z)=\frac{\beta_{X}^{\alpha_{X}}\beta_{Y}^{\alpha_{Y}}}{\Gamma(\alpha_{X})\Gamma(\alpha_{Y})} \begin{cases} e^{\beta_{Y}z}\int_{0}^{\infty}x^{\alpha_{X}-1}(x-z)^{\alpha_{Y}-1}e^{-(\beta_{X}+\beta_{Y})x}\,\mathrm{d}x & z<0\\ e^{-\beta_{X}z}\int_{0}^{\infty}y^{\alpha_{Y}-1}(x+z)^{\alpha_{X}-1}e^{-(\beta_{X}+\beta_{Y})y}\,\mathrm{d}y & z>0. \end{cases}\tag{3} \end{equation}
Now to find $f_{Z}(0)$ one simply substitutes $z=0$ into the above expressions for $f_{Z}$ when $z<0$ and $z>0$. After the substitution, it can be seen that both pieces will give the same answer of \begin{align} \tag{4.a} f_{Z}(0) &=\frac{\beta_{X}^{\alpha_{X}}\beta_{Y}^{\alpha_{Y}}}{\Gamma(\alpha_{X})\Gamma(\alpha_{Y})} \int_{0}^{\infty}x^{\alpha_{X}+\alpha_{Y}-2}e^{-(\beta_{X}+\beta_{Y})x}\,\mathrm{d}x\\ \tag{4.b} &=\frac{\beta_{X}^{\alpha_{X}}\beta_{Y}^{\alpha_{Y}}}{\Gamma(\alpha_{X})\Gamma(\alpha_{Y})(\beta_{X}+\beta_{Y})^{\alpha_{X}+\alpha_{Y}-1}} \int_{0}^{\infty}t^{\alpha_{X}+\alpha_{Y}-2}e^{-t}\,\mathrm{d}t\\ \tag{4.c} &=\frac{\beta_{X}^{\alpha_{X}}\beta_{Y}^{\alpha_{Y}}}{\Gamma(\alpha_{X})\Gamma(\alpha_{Y})(\beta_{X}+\beta_{Y})^{\alpha_{X}+\alpha_{Y}-1}}\Gamma(\alpha_{X}+\alpha_{Y}-1) \end{align}
Here is where my question comes into play... The quantity $\Gamma(\alpha_{X}+\alpha_{Y}-1)$ is not defined for $\alpha_{X}+\alpha_{Y}=1$, and even more disturbing is negative for $0<\alpha_{X}+\alpha_{Y}<1$. How can this be when the quantity $\alpha_{X}+\alpha_{Y}$ can take on these values? How would one find the value of $f_{Z}(0)$?
As it turns out the issue resides at Eq. (4.c). The integral definition for $\Gamma(z)$ is only defined for $\Re z>0$, i.e. \begin{equation} \Gamma(z):=\int_{0}^{\infty}t^{z-1}e^{-t}\,\mathrm{d}t\qquad\text{for}\ \Re z>0. \end{equation}
For $z\leq0$, $\Gamma(z)$ is defined by analytic continuation. Therefore, if $0<\alpha_{1}+\alpha_{2}\leq1$, the integral in Eq. (4.b) is: \begin{equation} \begin{aligned} \int_{0}^{\infty} t^{\alpha_{3}-2}e^{-t}\,\mathrm{d}t &=-\left(\lim_{t\to\infty}\Gamma(\alpha_{3}-1,t)-\lim_{t\to 0^{+}}\Gamma(\alpha_{3}-1,t)\right)\\ &= \begin{cases} \infty &\text{for}\ 0<\alpha_{3}\leq1\,,\\ \Gamma(\alpha_{3}-1) &\text{for}\ \alpha_{3}>1\,. \end{cases} \end{aligned} \end{equation} where $\alpha_{3}=\alpha_{1}+\alpha_{2}$. Using this solution, $f_{Z}(0)$ is
\begin{equation} f_{Z}(0) = \begin{cases} \infty &\text{for}\ 0<\alpha_{3}\leq1\,,\\ \frac{\beta_{1}^{\alpha_{1}}\beta_{2}^{\alpha_{2}}\Gamma(\alpha_{3}-1)}{\beta_{3}^{\alpha_{3}-1}\Gamma(\alpha_{1})\Gamma(\alpha_{2})} &\text{for}\ \alpha_{3}>1\,. \end{cases} \end{equation}