Compute integral about cosh How to simplify following form?
$$(2\pi)^{-d/2} \det{(\Sigma)}^{-1/2} e^{-\mu^T \Sigma^{-1} \mu}  \int_{x \in R^d} e^{-\frac{1}{2}  x^T \Sigma^{-1}x } \cosh(\mu^T \Sigma^{-1} x) \ln(\cosh(\mu^T \Sigma^{-1} x)) dx $$
where $x \in R^d$, $\Sigma \in R^{d \times d}$ and $\Sigma > 0$ is a symmetric positive-definite matrix.
Above form is a multivariate version of Eq(4) in the paper http://mdpi.org/entropy/papers/e10030200.pdf
The original form is:
$$A = \frac{1}{\sqrt{2\pi} \sigma} e^{-\mu^2 / (2 \sigma ^2)} \int_{-\infty}^{\infty} e^{-x^2/(2\sigma^2) } \cosh(\mu x/\sigma^2) \ln(\cosh(\mu x/\sigma^2)) dx$$
let $y=\mu x/\sigma^2$, we have 
$$A =  \frac{1}{\sqrt{2\pi} \sigma} e^{-\mu^2 / (2 \sigma ^2)} (\sigma^2/\mu) \int_{-\infty}^{\infty} e^{-\sigma^2 y^2 / (2\mu^2)} \cosh(y)  \ln(\cosh(y)) dy \\
   = \frac{2}{\sqrt{2\pi} \sigma} e^{-\mu^2 / (2 \sigma)} (\sigma/\mu) \int_{0}^{\infty} e^{-\sigma^2 y^2 / (2\mu^2)} \cosh(y)  \ln(\cosh(y)) dy$$
Further let $\alpha = \sigma / \mu $, we get
$$A = \frac{2}{\sqrt{2\pi} \alpha} e^{-\alpha^2 / 2 }  \int_{0}^{\infty} e^{-y^2 / (2\alpha^2)} \cosh(y)  \ln(\cosh(y)) dy$$
I wanna ask whether we can get a similar simplification for the multivariate form.
 A: Let's write the integral as:
$$
I = \frac{
\exp\left\{-\frac{ \boldsymbol\mu\cdot\mathbf\Sigma^{-1}\cdot\boldsymbol\mu}{2}\right\}
 }{\sqrt{(2\pi)^n\textrm{det}\;(\mathbf\Sigma)}} \int_{\mathbb{R}^n}\exp\left\{-\frac{\boldsymbol{x}\cdot\mathbf\Sigma^{-1}\cdot\boldsymbol x}{2}\right\}\cosh\left(\boldsymbol\mu\cdot\mathbf\Sigma^{-1}\cdot\boldsymbol x\right)\;\ln\cosh\left(\boldsymbol\mu\cdot\mathbf\Sigma^{-1}\cdot\boldsymbol x\right)\;\textrm{d}\boldsymbol{x}
$$
Since $\mathbf{\Sigma}$ is a symmetric positive definite matrix, you can decompose it as:
$$ \mathbf{\Sigma} = \mathbf{Q}\cdot\mathbf{\Lambda}\cdot\mathbf{Q}^\textrm{T}$$
With $\textbf{Q}$ an orthonormal matrix ($\mathbf{Q}\cdot\mathbf{Q}^\textrm{T}=\textbf{I}$) and $\mathbf\Lambda$ diagonal with only strict positive diagonal elements.
Now define $\boldsymbol y = \mathbf{Q}\cdot\boldsymbol x$ and you can merge the integrals back together:


*

*$\boldsymbol y = \mathbf{Q}\cdot\boldsymbol x$

*$\mathrm{d}\boldsymbol y = \mathrm{d}\boldsymbol x$

*$\boldsymbol{\xi} = \mathbf{Q}\cdot\boldsymbol{\mu}$
leading to:
$$
I = \frac{
\exp\left\{-\frac{ \boldsymbol\xi\cdot\mathbf\Lambda^{-1}\cdot\boldsymbol\xi}{2}\right\}
 }{\sqrt{(2\pi)^n\textrm{det}\;(\mathbf\Lambda)}} \int_{\mathbb{R}^n}\exp\left\{-\frac{\boldsymbol{y}\cdot\mathbf\Lambda^{-1}\cdot\boldsymbol y}{2}\right\}\cosh\left(\boldsymbol\xi\cdot\mathbf{\Lambda}^{-1}\cdot\boldsymbol y\right)\;\ln\cosh\left(\boldsymbol\xi\cdot\mathbf{\Lambda}^{-1}\cdot\boldsymbol y\right)\;\textrm{d}\boldsymbol{y}
$$
Since $\mathbf\Lambda$ is diagonal with strict positive elements, we can define $\boldsymbol z=\sqrt{\mathbf\Lambda^{-1}}\boldsymbol y$ and $\boldsymbol\zeta=\sqrt{\mathbf\Lambda^{-1}}\boldsymbol\xi$ which gives:
$$
I = \frac{
\exp\left\{-\frac{ \boldsymbol\zeta\cdot\boldsymbol\zeta}{2}\right\}
 }{\sqrt{(2\pi)^n}} \int_{\mathbb{R}^n}\exp\left\{-\frac{\boldsymbol z\cdot\boldsymbol z}{2}\right\}\cosh\left(\boldsymbol\zeta\cdot\boldsymbol z\right)\;\ln\cosh\left(\boldsymbol\zeta\cdot\boldsymbol z\right)\;\textrm{d}\boldsymbol{z}
$$
Now we transform to $n$-dimensional spherical coordinates using the $d$-dimensional spherical volume element. Note that we take the axis such that $\boldsymbol\zeta\cdot\boldsymbol z = \zeta r \cos(\varphi_1)$:
$$
I = \frac{
\exp\left\{-\frac{ \zeta^2}{2}\right\}
 }{\sqrt{(2\pi)^n}} \idotsint\exp\left\{-\frac{r^2}{2}\right\}\cosh\left(\zeta r \cos\varphi_1\right)\;\ln\cosh\left(\zeta r \cos\varphi_1\right)r^{n-1}\sin^{n-2}(\varphi_1)\sin^{n-3}(\varphi_2)\cdots \sin(\varphi_{n-2})\,
dr\,d\varphi_1 \, d\varphi_2\cdots d\varphi_{n-1}
$$
The final steps are now to use:
$$\int_0^\pi \sin^n\varphi\;\textrm{d}\varphi = \sqrt{\pi}\frac{\Gamma\left(\frac{n+1}{2}\right)}{\Gamma\left(\frac{n}{2}+1\right)}$$
and realize that:
$$\idotsint \sin^{n-3}(\varphi_2)\cdots \sin(\varphi_{n-2})\,
d\varphi_2\cdots d\varphi_{n-1} = \frac{2\pi^{(n-1)/2}}{\Gamma\left(\frac{n-1}{2}\right)}$$
note: don't forget that $\varphi_{n-1}$ is integrated over $[0,2\pi]$ while all other angles only from $[0,\pi]$
which reduces the integral to the final form:
$$
\bbox[5px,border:2px solid #00A000]{
\begin{align}
I = \sqrt{\frac{2^{n+1}}{2\pi}}
\frac{\exp\left\{-\frac{ \zeta^2}{2}\right\}
 }{\Gamma\left(\frac{n-1}{2}\right)} \int_0^\infty \int_0^\pi&\exp\left\{-\frac{r^2}{2}\right\}\cosh\left(\zeta r \cos\varphi_1\right)\;\ln\cosh\left(\zeta r \cos\varphi_1\right)\\
&\times r^{n-1}\sin^{n-2}(\varphi_1)\,
dr\,d\varphi_1
\end{align}
}
$$
