How does one compute the convolution of two hyperbolic distributions?

Let $X$ and $Y$ be hyperbolically distributed, and let $f_{X} (x)$ and $f_{Y} (x)$ be their associated density functions. Then $$f_{X} (x) = \frac{\gamma}{2 \alpha \delta K_{1} ( \delta \gamma) } e^{- \alpha \sqrt{\delta^{2}+(x-\mu)^{2}} + \beta (x - \mu)} ,$$ where $K_{1} ( \cdot )$ is a Bessel function of the second kind (link), $\alpha, \beta, \delta \in \mathbb{R}$ and $\gamma = \sqrt{\alpha^{2}-\beta^{2}}$. Let's denote this by saying $X \sim H( \alpha, \beta, \delta, \mu)$. You can find more information on the hyperbolic distribution over here.

I am wondering how to compute the convolution $(f_{X} * f_{Y})(x)$ . Whether it is done directly or perhaps by means of Laplace/Fourier transformations and the convolution theorem. Do you have a source addressing this problem? Or do you know how to compute this convolution yourself?

If it is not possible to compute the convolution analytically, what would be a good way to approximate it?