# Can this integral be computed in closed-form?

I am interested in computing a normalizing constant (of a Gaussian density in dimension $3$). Such normalizing constants often do not have a closed form. In dimension $2$, this normalizing constant can be computed in closed form.

The integral I would like to compute is :

$$\int_{\mathbb{R}^{3}} e^{-(r_{1}^{2} + r_{2}^{2} + r_{3}^{3})/2\sigma^{2}} \sinh\Big( \frac{\vert r_{1}-r_{2} \vert}{2} \Big)\sinh\Big( \frac{\vert r_{1}-r_{3} \vert}{2}\Big)\sinh\Big( \frac{\vert r_{2}-r_{3} \vert}{2} \Big)dr_{1}dr_{2}dr_{3}$$

Using the spherical coordinates $(r_{1},r_{2},r_{3}) = (r\sin\varphi\cos\theta,r\sin\varphi\sin\theta,r\cos\varphi)$, I find myself trying to compute the following integral :

$$\int_{0}^{+\infty}\int_{0}^{2\pi}\int_{0}^{\pi} e^{-r^{2}/2\sigma^{2}} \sinh\Big( \frac{r\sin\varphi\vert \cos\theta - \sin\theta\vert}{2}\Big)\sinh\Big( \frac{r\vert\sin\varphi\cos\theta-\cos\phi\vert}{2}\Big) \times \sinh\Big(\frac{r\vert\sin\varphi\sin\theta-\cos\varphi\vert}{2}\Big)r^{2}\sin\varphi \, drd\theta d\varphi$$

I'm not an expert but, looking at this integral, I doubt there exist a closed expression.

Since you have a full symmetry by exchanging $r_1$, $r_2$ and $r_3$, the integral can be computed as $$3!\int_{-\infty}^{+\infty}\mathrm dx \int_{-\infty}^x\mathrm dy\int_{-\infty}^y\mathrm dz \mathrm e^{-(x^2+y^2+z^2)/2\sigma^2}\sinh\left(\frac{x-y}2\right)\sinh\left(\frac{x-z}2\right)\sinh\left(\frac{y-z}2\right).$$ Then writing the $\sinh$ in the exponential form, you end up with a sum of product integrals that are all of Gaussian type.