Itzykson-Zuber integral over orthogonal groups

Question

I would like to know is there a closed form expression for the following Itzykson-Zuber integral for the orthogonal case.

$I=\int_{\mathcal{O}(p)} \exp\left(-\frac{1}{2}~\mathrm{tr}(\Sigma^{-1}HLH^{T})\right)~\mathrm{d}H$

where ${\mathcal{O}(p)}$ is a Orthogonal group of $p \times p$ symmetric matrices, $\Sigma \neq aI $ is a covariance matrix of a column vector of a correlated Gaussian matrix $A$ that occurs in the Wishart matrix formation $W=A^{T}A$, $dH$ is normalized Harr measure and $L$ is the diagonal matrix $\mathrm{diag}(l_1,l_2,\cdots,l_p)$

I know that an infinite Zonal polynomial series exists for the integral. But, I am interested to know is there any closed-form or tractable solutions for the integral.

Answer 1

Take a look in the abstracts of the 16 workshop in Bedlew (non coummutativ harmonic..) 6-12 of 7

Answer 2

I present the answer for $p=2$ in here. We write the orthogonal matrix $H$ by using the Cayley' parametrization: \begin{equation} H = \left(1 + A\right) \cdot \left( 1 - A \right)^{-1} \end{equation} where \begin{equation} A = \left( \begin{array}{cc} 0 & a \\ -a & 0 \end{array} \right) \end{equation} is a real skew-symmetric matrix. A straightforward calculation shows that: \begin{equation} H = \frac{1}{1+a^2} \left( \begin{array}{rr} 1-a^2 & 2 a \\ - 2 a & 1- a^2 \end{array} \right) \end{equation} Now we have: \begin{equation} \left(\det \Sigma\right) Tr\left[\Sigma^{-1} H L H^T\right] = \frac{1}{2} Tr[\Sigma](l_1+l_2) + \frac{(l_1-l_2)}{2}\left((\Sigma_{1,1} - \Sigma_{2,2}) \cos(4 y) - 2 \Sigma_{1,2}\sin(4 y)\right) \end{equation} where $a := \tan(y)$. Another straightforward calculation gives the Jacobian: \begin{equation} \frac{\partial(c_{1,1}, c_{1,2}, c_{2,2})}{\partial(l_1,l_2,a)} = \frac{2 (l_1-l_2)}{1+a^2} \end{equation} where \begin{equation} \left( \begin{array}{cc} c_{1,1} & c_{1,2} \\ c_{1,2} & c_{2,2} \end{array} \right) = H L H^T \end{equation} Bringing everything together our integral reads: \begin{eqnarray} I &=& 2 (l_1-l_2) e^{-\frac{1}{4} \frac{Tr\left[\Sigma\right]}{\det \Sigma} (l_1+l_2)} \cdot \int\limits_{-\frac{\pi}{2}}^{\frac{\pi}{2}} e^{\frac{(l_2-l_1)}{4 \det \Sigma} \left((\Sigma_{1,1}-\Sigma_{2,2}) \cos(4 y) - 2 \Sigma_{1,2} \sin(4 y)\right)} d y \\ &=& (l_1-l_2) e^{-\frac{1}{4} \frac{Tr\left[\Sigma\right]}{\det \Sigma} (l_1+l_2)} \cdot I_0\left((l_2-l_1) \frac{\sqrt{(Tr \Sigma)^2 - 4 \det(\Sigma)}}{4 \det \Sigma}\right) \end{eqnarray} where $I_0$ is the modified Bessel function of order zero. As we can see the results depends only on the rotational invariants of the matrix $\Sigma$.

Answer 3

I show a solution for a three by three matrix ($p=3$). The calculations follow the same pattern as in my $p=2$ solution above, but as we will see things get slightly more complicated.

We start, as before, with a Cayley' parametrization of the orthogonal matrix $H$. We have: \begin{equation} H = (1 + A) \cdot (1 - A)^{-1} \end{equation} where \begin{equation} A = \left( \begin{array}{ccc} 0 & a & b \\ -a & 0 & c \\ -b & -c & 0 \end{array} \right) \end{equation} is a skew-symmetric matrix and $a$, $b$ and $c$ are real numbers. Combining the equations above we get: \begin{equation} H = \frac{1}{1+a^2+b^2+c^2} \left( \begin{array}{ccc} 1-a^2-b^2+c^2 & 2(a- b c) & 2(b+a c) \\ 2(-a-b c) & 1-a^2+b^2- c^2 & 2(c - a b) \\ 2(-b + a c) & 2(-c- a b) & 1+a^2-b^2-c^2 \end{array} \right) \end{equation} Now we compute the Jacobian: \begin{equation} \frac{\partial(c_{1,1}, c_{2,2}, c_{3,3}, c_{1,2}, c_{1,3}, c_{2,3})}{\partial(l_1,l_2,l_3,a,b,c)} = - \frac{8(l_1-l_2)(l_1-l_3)(l_2-l_3)}{(1+a^2+b^2+c^2)^2} \end{equation} where \begin{equation} \left( \begin{array}{ccc} c_{1,1} & c_{1,2} & c_{1,3} \\ c_{1,2} & c_{2,2} & c_{2,3} \\ c_{1,3} & c_{2,3} & c_{3,3} \end{array} \right) = H L H^T \end{equation} At this stage we have have reduced the integral over the orthogonal group to an integral over $(a,b,c) \in {\mathcal R}^3$. Since the expressions in the exponential are quadratic forms it is a natural choice to go to polar coordinates. So we set: \begin{eqnarray} a &=& r \cos(\theta) \cos(\phi) \\ b &=& r \cos(\theta) \sin(\phi) \\ c &=& r \sin(\theta) \end{eqnarray} where $r = \tan(\lambda)$. Here $(\lambda, \theta, \phi) \in (-\frac{\pi}{2},\frac{\pi}{2})^2 \times (0, 2 \pi)$.

The Jacobian of this transformation is easy to compute: \begin{equation} \frac{\partial(a,b,c)}{\partial(\lambda,\theta,\phi)} = \frac{\sin(\lambda)^2}{\cos(\lambda)^4} \cos(\theta) \end{equation} Combining the two Jacobians together we get: \begin{equation} \frac{\partial(c_{1,1}, c_{2,2}, c_{3,3}, c_{1,2}, c_{1,3}, c_{2,3})}{\partial(l_1,l_2,l_3,a,b,c)} \frac{\partial(a,b,c)}{\partial(\lambda,\theta,\phi)} = -8 \Delta(\vec{l}) \cdot \sin(\lambda)^2 \cos(\theta) \end{equation} where $\Delta(\lambda) := \prod\limits_{1\le i < j \le 3} (l_i - l_j)$ is the Vandermonde determinant of the matrix $L$. Now we have to write down the the trace in the exponential. Unfortunately this becomes messy. This quantity is a huge sum of terms that contain products of powers of cosines and sines of our angles $\theta$, $\phi$ and $\lambda$. However we can always decompose those non-linear terms into sines and cosines or integer multiples of our angles just by using Euler's formula $\cos(\phi) = (e^{\imath \phi} + e^{-\imath \phi})/2$ and $\sin(\phi) = (e^{\imath \phi} - e^{-\imath \phi})/(2 \imath)$. We have done that using Mathematica and we got the following result: \begin{eqnarray} &&I= 8 \Delta(\vec{l}) \exp(-\frac{1}{2 \det \Sigma} \sum\limits_{j=1}^3 {\mathcal A}_j l_j) \cdot \int\limits_{ (-\frac{\pi}{2},\frac{\pi}{2})^2 \times (0, 2 \pi)} \sin(\lambda)^2 \cos(\theta) \\&& \exp\left(-\frac{1}{2 \det \Sigma} \sum\limits_{i=0}^4 \sum\limits_{j=-4}^4 \sum\limits_{k=-4}^4 \left[ {\mathcal X}_{i,j,k} \cos( i \lambda + j \phi + k \theta) + {\mathcal Y}_{i,j,k} \sin( i \lambda + j \phi + k \theta) \right] \right) d \lambda d \theta d \phi \end{eqnarray} Here: \begin{eqnarray} {\mathcal A}_1 &=& \frac{1}{128} (-28 \Sigma_{1,2}^2 - 28 \Sigma_{1,3}^2 - 72 \Sigma_{2,3}^2 + 28 \Sigma_{1,1} \Sigma_{2,2} + 28 \Sigma_{1,1} \Sigma_{3,3} + 72 \Sigma_{2,2} \Sigma_{3,3}) \\ {\mathcal A}_2 &=& \frac{1}{128} (-41 \Sigma_{1,2}^2 - 59 \Sigma_{1,3}^2 - 28 \Sigma_{2,3}^2 + 41 \Sigma_{1,1} \Sigma_{2,2} + 59 \Sigma_{1,1} \Sigma_{3,3} + 28 \Sigma_{2,2} \Sigma_{3,3}) \\ {\mathcal A}_3 &=& \frac{1}{128} (-59 \Sigma_{1,2}^2 - 41 \Sigma_{1,3}^2 - 28 \Sigma_{2,3}^2 + 59 \Sigma_{1,1} \Sigma_{2,2} + 41 \Sigma_{1,1} \Sigma_{3,3} + 28 \Sigma_{2,2} \Sigma_{3,3}) \end{eqnarray} and \begin{equation} \sum\limits_{j=1}^3 {\mathcal A}_j = \frac{1}{2} \left((Tr[ \Sigma])^2 - Tr[ \Sigma^2]\right) \end{equation} the quantities $ {\mathcal X}_{i,j,k}$ and $ {\mathcal Y}_{i,j,k}$ are all of the form: \begin{eqnarray} {\mathcal X}_{i,j,k} &=& {\mathcal X}^{(1)}_{i,j,k} \cdot (l_1-l_3) + {\mathcal X}^{(2)}_{i,j,k} \cdot (l_2-l_3) \\ {\mathcal Y}_{i,j,k} &=& {\mathcal Y}^{(1)}_{i,j,k} \cdot (l_1-l_3) + {\mathcal Y}^{(2)}_{i,j,k} \cdot (l_2-l_3) \end{eqnarray} where ${\mathcal X}^{(p)}_{i,j,k}$ and ${\mathcal Y}^{(p)}_{i,j,k}$ for $p=1,2$ are some quadratic forms of elements of the matrix $\Sigma$. The three dimensional integral above can be expressed in terms of Generalized Bessel Functions (GBF). Since I only started learning about GBFs I leave the integral as it is for the time being and will provide the final solution later.

Let us define a Generalised modified Bessel function as follows: \begin{eqnarray} &&{\mathcal I}_{n,m,o}^{(d_1,d_2,d_3)} \left( \left\{{\mathcal X}_{i,j,k}\right\},\left\{{\mathcal Y}_{i,j,k}\right\} \right) := \\&& (2\pi)^3 \sum\limits_{\left\{n_{p,q,r}\right\} \in {\mathcal Z}^{\theta}} \prod\limits_{p=0}^{d_1-1} \prod\limits_{q=-d_2}^{d_2} \prod\limits_{r=-d_3}^{d_3} {\mathcal I}_{n_{p,q,o}} \left(\sqrt{(x_{p,q,r})^2 + (y_{p,q,r})^2}\right) e^{\imath n_{p,q,r} \lambda_{p,q,r}} \end{eqnarray} Here $\theta := d_1(2 d_2+1)(2 d_3+1)$ and $\lambda_{p,q,r} := -\arctan(x_{p,q,r}/y_{p,q,r})$ and the sum on the right hand side runs over $\theta$ independent integers $n_{p,q,r}$ subject to following constraints: \begin{equation} \sum\limits_{p=0}^{d_1-1} \sum\limits_{q=-d_2}^{d_2} \sum\limits_{r=-d_3}^{d_3} n_{p,q,r} \left\{ \begin{array}{c} p \\ q \\ r \end{array} \right\} = \left\{ \begin{array}{c} n \\ m \\ o \end{array} \right\} \end{equation} Then the integral in question reads: \begin{eqnarray} I=- \Delta(\vec{l}) e^{-\frac{\sum\limits_{j=1}^3 {\mathcal A}_j l_j}{2 \det \Sigma}} \left( {\mathcal I}_{-2,-1,0} + {\mathcal I}_{-2,1,0} - 2{\mathcal I}_{0,-1,0} - 2{\mathcal I}_{0,1,0} +{\mathcal I}_{2,-1,0}+{\mathcal I}_{2,1,0} \right) \end{eqnarray} Here we have dropped the tensor arguments which are all the same and read: \begin{equation} -\frac{1}{2 \det \Sigma} \left(\left\{{\mathcal X}_{i,j,k}\right\}, \left\{{\mathcal Y}_{i,j,k}\right\}\right) \end{equation} The dimensions reads $\left(d_1,d_2,d_3\right) = \left(5,4,4\right)$.

Answer 4

Without loss of generality we can assume that $\Sigma$ is a diagonal matrix. Indeed since $\Sigma$ is a full rank symmetric matrix it has a spectral decomposition $\Sigma = O^T \cdot {\mathcal D} \cdot O$ where $O\cdot O^T = O^T\cdot O=1$ and ${\mathcal D}:=(s_1,s_2,s_3)$ is diagonal. Then of course $\Sigma^{-1} = O^T \cdot {\mathcal D}^{-1} \cdot O$ and \begin{equation} Tr\left( \Sigma^{-1} \cdot H L H^T \right)= Tr\left( O^T \cdot {\mathcal D}^{-1} \cdot O \cdot H L H^T\right) = Tr\left( {\mathcal D}^{-1} O \cdot H L H^T O^T\right)= Tr\left({\mathcal D}^{-1} \cdot OH \cdot L \cdot (OH)^T\right) \end{equation} where we used the cyclic property of the trace. Now in our integral we can change variables $H \rightarrow OH$ and use the fact that the Jacobian of this mapping is unity.

Now we construct an orthogonal matrix as a product of three rotational matrices each one in three mutually orthogonal planes, $(1,2)$, $(1,3)$ and $(2,3)$ respectively. We have: \begin{eqnarray} &&\left(H_{1,2},H_{1,3},H_{2,3}\right) = \\&& \left( \left(\begin{array}{ccc} \cos(\phi_1) & -\sin(\phi_1) & 0 \\ \sin(\phi_1) & \cos(\phi_1) & 0 \\ 0 & 0 & 1 \end{array} \right), \left(\begin{array}{ccc} \cos(\phi_2) & 0 & -\sin(\phi_2) \\ 0 & 1 & 0 \\ \sin(\phi_2) & 0 & \cos(\phi_2) \end{array} \right), \left(\begin{array}{ccc} 1 & 0 & 0 \\ 0 & \cos(\phi_3) & -\sin(\phi_3) \\ 0 & \sin(\phi_3) & \cos(\phi_3) \end{array} \right) \right) \end{eqnarray} Therefore: \begin{eqnarray} &&H= H_{1,2} H_{1,3} H_{2,3} = \\ && \begin{array}{rrr} \cos(\phi_1) \cos(\phi_2),& -\cos(\phi_3) \sin(\phi_1) - \cos(\phi_1) \sin(\phi_2) \sin(\phi_3),& -\cos(\phi_1) \cos(\phi_3) \sin(\phi_2) + \sin(\phi_1) \sin(\phi_3) \\ \cos(\phi_2) \sin(\phi_1),& \cos(\phi_1) \cos(\phi_3) - \sin(\phi_1) \sin(\phi_2) \sin(\phi_3),& -\cos(\phi_3) \sin(\phi_1) \sin(\phi_2) - \cos(\phi_1) \sin(\phi_3) \\ \sin(\phi_2),&\cos(\phi_2) \sin(\phi_3),&\cos(\phi_2) \cos(\phi_3) \end{array} \end{eqnarray} The Jacobian of the mapping from the matrix elements to the eigenvalues and the angles is quite simple, much simpler than before. It reads: \begin{equation} \frac{\partial(c_{1,1},c_{2,2},c_{3,3},c_{1,2},c_{1,3},c_{2,3})}{\partial(l_1,l_2,l_3,\phi_1,\phi_2,\phi_3)} = (l_1-l_2)(l_1-l_3)(l_2-l_3) \cos(\phi_2) \end{equation} where \begin{equation} \left( \begin{array}{ccc} c_{1,1} & c_{2,2} & c_{1,3} \\ c_{1,2} & c_{2,2} & c_{2,3} \\ c_{1,3} & c_{2,3} & c_{3,3} \end{array} \right) = H L H^T \end{equation} The expression inside the exponential reads: \begin{eqnarray} && Tr\left[\Sigma^{-1} H L H^T\right] =\\&& \left(\sum\limits_{j=1}^3 {\mathcal A}_j l_j\right) +\\ && \sum\limits_{j=1}^3 {\mathcal X}_j \cdot \cos(2 \phi_j) + \sum\limits_{1\le i < j \le 3} {\mathcal X}^+_{i,j} \cdot \cos(2\phi_i+2\phi_j) + \sum\limits_{1\le i < j \le 3} {\mathcal X}^-_{i,j} \cdot \cos(2\phi_i-2\phi_j) +\\ &&\sum_{1\le i<j<k\le 3} {\mathcal X}_{i,j,k} \cos(2(\phi_i+\phi_j+\phi_k))+\\ &&\sum_{1\le i<j<k\le 3} {\mathcal X}^{(-,+,+)}_{i,j,k} \cos(2(-\phi_i+\phi_j+\phi_k))+\\ &&\sum_{1\le i<j<k\le 3} {\mathcal X}^{(+,-,+)}_{i,j,k} \cos(2(\phi_i-\phi_j+\phi_k))+\\ &&\sum_{1\le i<j<k\le 3} {\mathcal X}^{(+,+,-)}_{i,j,k} \cos(2(\phi_i+\phi_j-\phi_k))+\\ &&\sum\limits_{1\le i<j<k\le 3} {\mathcal Y}_{i,j,k} \sin(2 \phi_i+\phi_j+2 \phi_k)+\\ &&\sum\limits_{1\le i<j<k\le 3} {\mathcal Y}^{(-,+,+)}_{i,j,k} \sin(-2 \phi_i+\phi_j+2 \phi_k)+\\ &&\sum\limits_{1\le i<j<k\le 3} {\mathcal Y}^{(+,-,+)}_{i,j,k} \sin(2 \phi_i-\phi_j+2 \phi_k)+\\ &&\sum\limits_{1\le i<j<k\le 3} {\mathcal Y}^{(+,+,-)}_{i,j,k} \sin(2 \phi_i+\phi_j-2 \phi_k) \end{eqnarray} where \begin{eqnarray} {\mathcal A}_1 &=& \frac{1}{4 s_1}+\frac{1}{4 s_2}+\frac{1}{2 s_3} \\ {\mathcal A}_2 &=& \frac{3}{8 s_1}+\frac{3}{8 s_2}+\frac{1}{4 s_3} \\ {\mathcal A}_3 &=& \frac{3}{8 s_1}+\frac{3}{8 s_2}+\frac{1}{4 s_3} \end{eqnarray} and \begin{equation} \sum\limits_{j=1}^3 {\mathcal A}_j = \frac{1}{s_1}+\frac{1}{s_2}+\frac{1}{s_3} = = \frac{1}{2} \frac{\left(Tr[\Sigma]^2 - Tr[\Sigma^2]\right)}{\det(\Sigma)} \end{equation} and \begin{eqnarray} {\mathcal X}_1 &=& \frac{1}{8}((l_2-l_1)+(l_3-l_1)) \frac{s_1-s_2}{s_1 s_2}\\ {\mathcal X}_2 &=& \frac{1}{8}((l_2-l_1)+(l_3-l_1)) \left( \frac{s_2-s_3}{s_2 s_3} + \frac{s_1-s_3}{s_1 s_3}\right) \\ {\mathcal X}_3 &=& \frac{1}{8} (l_3-l_2) \left( \frac{s_2-s_3}{s_2 s_3} + \frac{s_1-s_3}{s_1 s_3}\right)\\ {\mathcal X}^+_{1,2} &=&\frac{1}{16} ((l_2-l_1)+(l_3-l_1)) \frac{s_1-s_2}{s_1 s_2}\\ {\mathcal X}^+_{1,3} &=&\frac{3}{16} (l_2-l_3) \frac{s_1-s_2}{s_1 s_2}\\ {\mathcal X}^+_{2,3} &=& \frac{1}{16} (l_3-l_2) \left( \frac{s_2-s_3}{s_2 s_3} + \frac{s_1-s_3}{s_1 s_3}\right)\\ {\mathcal X}^-_{1,2} &=& {\mathcal X}^+_{1,2}\\ {\mathcal X}^-_{2,3} &=& {\mathcal X}^+_{2,3}\\ {\mathcal X}^-_{1,3} &=& {\mathcal X}^+_{1,3}\\ {\mathcal X}_{1,2,3}&=& \frac{1}{32} (l_3-l_2)\frac{s_1-s_2}{s_1 s_2}\\ {\mathcal X}^{(-,+,+)}_{1,2,3}&=&{\mathcal X}_{1,2,3}\\ {\mathcal X}^{(+,-,+)}_{1,2,3}&=&{\mathcal X}_{1,2,3}\\ {\mathcal X}^{(+,+,-)}_{1,2,3}&=&{\mathcal X}_{1,2,3}\\ {\mathcal Y}_{1,2,3}&=& \frac{1}{8} (l_2-l_3) \frac{s_1-s_2}{s_1 s_2}\\ {\mathcal Y}^{(-,+,+)}_{1,2,3}&=&-{\mathcal Y}_{1,2,3}\\ {\mathcal Y}^{(+,-,+)}_{1,2,3}&=&-{\mathcal Y}_{1,2,3}\\ {\mathcal Y}^{(+,+,-)}_{1,2,3}&=&-{\mathcal Y}_{1,2,3} \end{eqnarray}

Therefore the final result has the same form as before: \begin{eqnarray} &&I = \Delta(\vec{l}) \cdot e^{-\frac{1}{2} \sum\limits_{j=1}^3 {\mathcal A}_j l_j} \cdot \int\limits_{[0,2\pi]^3} \left|\cos(\phi_2)\right| \\&& \exp(-1/2)\left(\right.\\ &&\sum\limits_{j=1}^3 {\mathcal X}_j \cdot \cos(2 \phi_j) + \sum\limits_{1\le i < j \le 3} {\mathcal X}^+_{i,j} \cdot \cos(2\phi_i+2\phi_j) + \sum\limits_{1\le i < j \le 3} {\mathcal X}^-_{i,j} \cdot \cos(2\phi_i-2\phi_j) +\\ &&\sum_{1\le i<j<k\le 3} {\mathcal X}_{i,j,k} \cos(2(\phi_i+\phi_j+\phi_k))+\\ &&\sum_{1\le i<j<k\le 3} {\mathcal X}^{(-,+,+)}_{i,j,k} \cos(2(-\phi_i+\phi_j+\phi_k))+\\ &&\sum_{1\le i<j<k\le 3} {\mathcal X}^{(+,-,+)}_{i,j,k} \cos(2(\phi_i-\phi_j+\phi_k))+\\ &&\sum_{1\le i<j<k\le 3} {\mathcal X}^{(+,+,-)}_{i,j,k} \cos(2(\phi_i+\phi_j-\phi_k))+\\ &&\sum\limits_{1\le i<j<k\le 3} {\mathcal Y}_{i,j,k} \sin(2 \phi_i+\phi_j+2 \phi_k)+\\ &&\sum\limits_{1\le i<j<k\le 3} {\mathcal Y}^{(-,+,+)}_{i,j,k} \sin(-2 \phi_i+\phi_j+2 \phi_k)+\\ &&\sum\limits_{1\le i<j<k\le 3} {\mathcal Y}^{(+,-,+)}_{i,j,k} \sin(2 \phi_i-\phi_j+2 \phi_k)+\\ &&\sum\limits_{1\le i<j<k\le 3} {\mathcal Y}^{(+,+,-)}_{i,j,k} \sin(2 \phi_i+\phi_j-2 \phi_k)\\ &&\left. \right) \prod\limits_{j=1}^3 d\phi_j \quad (ii) \end{eqnarray} This integral defines a Generalized (matrix) Bessel function (see Dattoli, G.; Torre, A.; Lorenzutta, S., Theory of multivariable Bessel functions and elliptic modular functions, Matematiche 53, No. 2, 387-399 (1998). ZBL1156.33303. for example) and therefore can be expressed as follows:

\begin{eqnarray} &&I = \Delta(\vec{l}) \cdot e^{-\frac{1}{2} \sum\limits_{j=1}^3 {\mathcal A}_j l_j} \cdot \underline{(2\pi)^2 \cdot \frac{4 (-1)^{|{\mathcal N}|/2+1}}{{\mathcal N}^2-1} }\cdot \left\{ \right.\\ && \prod\limits_{j=1}^3 \left[\imath ^{n_j} \cdot J_{n_j}\left(\frac{\imath}{2} {\mathcal X}_j\right)\right] \cdot \prod\limits_{\begin{array}{c} 1 \le i < j \le 3 \\ \xi=\pm \end{array} } \left[ \imath^{n^{(\xi)}_{i,j}} \cdot J_{n^{(\xi)}_{i,j}}\left(\frac{\imath}{2} {\mathcal X}^{\xi}_{i,j} \right) \right] \cdot \\ && \prod\limits_{\xi=0}^3 \left[ \imath^{n^{c,\vec{e}_\xi}} \cdot J_{n^{c,\vec{e}_\xi}}\left( \frac{\imath}{2} {\mathcal X}_{1,2,3}^{\vec{e}_\xi}\right) \right] \cdot \prod\limits_{\xi=0}^3 \left[ J_{n^{s,\vec{e}_\xi}} \left( \frac{\imath}{2} {\mathcal Y}_{1,2,3}^{\vec{e}_\xi}\right) \right] \\ \left. \right\} \quad (ii) \end{eqnarray} where the quantity in curly braces in the right hand side is being summed over $ \vec{n}:= \left(n_1,n_2,n_3;n_{1,2}^{+},n_{1,2}^{-},n_{1,3}^{+},n_{1,3}^{-},n_{2,3}^{+},n_{2,3}^{-};(n^{c,\vec{e}_\xi})_{\xi=0}^3; (n^{s,\vec{e}_\xi})_{\xi=0}^3 \right) \in {\mathbb Z}^{17} $ subject to following constraints:

\begin{eqnarray} 2 n_1 &+& 2(\begin{array}{c} n_{1,2}^{+}+n_{1,2}^{-}\\ + n_{1,3}^{+}+n_{1,3}^{-} \end{array}) + 2(\begin{array}{c} n^{c,\vec{e}_0}-n^{c,\vec{e}_1} \\ +n^{c,\vec{e}_2} + n^{c,\vec{e}_3} \end{array}) + 2(\begin{array}{c} n^{s,\vec{e}_0}-n^{s,\vec{e}_1} \\ +n^{s,\vec{e}_2} + n^{s,\vec{e}_3} \end{array}) &=& 0 \\ % 2 n_2 &+& 2(\begin{array}{c} n_{1,2}^{+}-n_{1,2}^{-}\\ + n_{2,3}^{+}+n_{2,3}^{-} \end{array}) + 2(\begin{array}{c} n^{c,\vec{e}_0}+n^{c,\vec{e}_1} \\ -n^{c,\vec{e}_2} + n^{c,\vec{e}_3} \end{array}) + 1(\begin{array}{c} n^{s,\vec{e}_0}+n^{s,\vec{e}_1} \\ -n^{s,\vec{e}_2} + n^{s,\vec{e}_3} \end{array}) &=& {\mathcal N} \in 2{\mathbb N} \\ % 2 n_3 &+& 2(\begin{array}{c} n_{1,3}^{+}-n_{1,3}^{-}\\ + n_{2,3}^{+}-n_{2,3}^{-} \end{array}) + 2(\begin{array}{c} n^{c,\vec{e}_0}+n^{c,\vec{e}_1} \\ +n^{c,\vec{e}_2} - n^{c,\vec{e}_3} \end{array}) + 2(\begin{array}{c} n^{s,\vec{e}_0}+n^{s,\vec{e}_1} \\ +n^{s,\vec{e}_2} - n^{s,\vec{e}_3} \end{array}) &=& 0 \end{eqnarray} In here $ \left( \vec{e}_\xi \right)_{\xi=0}^3 $ are defined as follows:

\begin{eqnarray} \vec{e}_0 &=& (+,+,+) \\ \vec{e}_1 &=& (-,+,+) \\ \vec{e}_2 &=& (+,-,+) \\ \vec{e}_3 &=& (+,+,-) \end{eqnarray}

Note 0: Equation $(ii) $ follows from applying the Jacobi-Anger expansion to each of the seventeen terms inside the exponential in $(i)$ and then using the orthogonality of plane waves along with the following identity:

\begin{equation} \int\limits_0^{2\pi} e^{\imath n \phi} \cdot \left| \cos(\phi) \right| d\phi = \left\{ \begin{array}{lll} \frac{4 (-1)^{|n|/2+1} }{n^2-1} & \quad \mbox{if $n$ is even }\\ 0 & \quad \mbox{otherwise} \end{array} \right. \end{equation}

Note 1: In general the integral over the angles, ie.e equation $(ii)$, depends only on $\left(l_1-l_2, l_1-l_3, l_2-l_3\right)$ on one hand and on $\left(1/s_1-1/s_2, 1/s_1-1/s_3, 1/s_2-1/s_3\right)$ on the other hand.

Note 2: On the face of it equation $(ii)$ does not seem to be useful in the generic case, since it is not quite clear how fast that multiple sum will converge, however in some particular cases that equation might be useful. Let us list some of those cases now.

1. Take $s_1=s_2=s_3= s$. Then $\vec{n} = \left(\underbrace{0,\cdots,0}_{17}\right)$ and ${\mathcal N}= 0$ and $I = \Delta(\vec{l}) \cdot e^{-1/(2 s) \cdot Tr[L]} \cdot (2\pi)^2 \cdot 4 $. This makes sense by comparing with the definition in the body of the question.
2. Take $s_1=s_2 \neq s_3$. Then $\vec{n} = \left(0,n_2,n_3;0,0,0,0,n_{2,3}^+,n_{2,3}^{-};(0)_{\xi=0}^3;(0)_{\xi=0}^3 \right)$. Then the result reads:

\begin{eqnarray} &&I = \Delta(\vec{l}) \cdot e^{-\frac{1}{2} \sum\limits_{j=1}^3 {\mathcal A}_j l_j} \cdot \underline{(2\pi)^2 \cdot \frac{4 (-1)^{|{\mathcal N}|/2+1}}{{\mathcal N}^2-1} }\cdot \left\{ \right.\\ && \prod\limits_{j=2}^3 \left[\imath ^{n_j} \cdot J_{n_j}\left(\frac{\imath}{2} {\mathcal X}_j\right)\right] \cdot \prod\limits_{\begin{array}{c} (i,j)=(2,3) \\ \xi=\pm \end{array} } \left[ \imath^{n^{(\xi)}_{i,j}} \cdot J_{n^{(\xi)}_{i,j}}\left(\frac{\imath}{2} {\mathcal X}^{\xi}_{i,j} \right) \right] \\ \left. \right\} \quad (iia) \end{eqnarray}

where the expression in curly brackets in $(iia)$ is summed over $(n_2,n_3;n_{2,3}^{+},n_{2,3}^{-}) \in {\mathbb Z}^4 $ subject to a constraint $2n_3 + 2(n_{2,3}^+-n_{2,3}^-)=0$ and to ${\mathcal N}:= 2 n_2+2(n_{2,3}^++n_{2,3}^-) $.

3. Take $l_2=l_3 \neq l_1$. Then $\vec{n} = \left(n_1,n_2,0;n_{1,2}^+,n_{1,2}^-,0,0,0,0;(0)_{\xi=0}^3;(0)_{\xi=0}^3 \right)$. Then the result reads:

\begin{eqnarray} &&I = \Delta(\vec{l}) \cdot e^{-\frac{1}{2} \sum\limits_{j=1}^3 {\mathcal A}_j l_j} \cdot \underline{(2\pi)^2 \cdot \frac{4 (-1)^{|{\mathcal N}|/2+1}}{{\mathcal N}^2-1} }\cdot \left\{ \right.\\ && \prod\limits_{j=1}^2 \left[\imath ^{n_j} \cdot J_{n_j}\left(\frac{\imath}{2} {\mathcal X}_j\right)\right] \cdot \prod\limits_{\begin{array}{c} (i,j)=(1,2) \\ \xi=\pm \end{array} } \left[ \imath^{n^{(\xi)}_{i,j}} \cdot J_{n^{(\xi)}_{i,j}}\left(\frac{\imath}{2} {\mathcal X}^{\xi}_{i,j} \right) \right] \\ \left. \right\} \quad (iib) \end{eqnarray}

where the expression in curly brackets in $(iib)$ is summed over $(n_1,n_2;n_{1,2}^{+},n_{1,2}^{-}) \in {\mathbb Z}^4 $ subject to a constraint $2n_1 + 2(n_{1,2}^++n_{1,2}^-)=0$ and to ${\mathcal N}:= 2 n_2+2(n_{1,2}^+-n_{1,2}^-) $.

Answer 5

I an writing an answer for $p=4$. Again, as explained in my answer for $p=3$, we can assume without loss of generality that the matrix $\Sigma$ is diagonal having eigenvalues $\left(s_i\right)_{i=1}^p$.

We construct the orthogonal matrix as product of $\binom{4}{2}$ rotational matrices in planes $(i,j)$ for $1\le i < j \le 4$. \begin{eqnarray} &&H=\\&&\left( \begin{array}{cccc} \cos (\phi_1) & -\sin (\phi_1) & 0 & 0 \\ \sin (\phi_1) & \cos (\phi_1) & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right) \left( \begin{array}{cccc} \cos (\phi_2) & 0 & -\sin (\phi_2) & 0 \\ 0 & 1 & 0 & 0 \\ \sin (\phi_2) & 0 & \cos (\phi_2) & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right)\\ && \left( \begin{array}{cccc} \cos (\phi_3 & 0 & 0 & -\sin (\phi_3) \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ \sin (\phi_3) & 0 & 0 & \cos (\phi_3) \\ \end{array} \right) \left( \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & \cos (\phi_4) & -\sin (\phi_4) & 0 \\ 0 & \sin (\phi_4) & \cos (\phi_4) & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right)\\ && \left( \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & \cos (\phi_5) & 0 & -\sin (\phi_5) \\ 0 & 0 & 1 & 0 \\ 0 & \sin (\phi_5) & 0 & \cos (\phi_5) \\ \end{array} \right) \left( \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & \cos (\phi_6) & -\sin (\phi_6) \\ 0 & 0 & \sin (\phi_6) & \cos (\phi_6) \\ \end{array} \right) \end{eqnarray} where each angle $\phi_j \in [0,2\pi]$ for $j=1,\cdots,6$. The Jacobian equals: \begin{equation} \frac{\partial(c_{1,1},c_{2,2},c_{3,3},c_{4,4},c_{1,2},c_{1,3}c_{1,4},c_{2,3},c_{2,4},c_{3,4})}{\partial(l_1,l_2,l_3,l_4,\phi_1,\phi_2,\phi_3,\phi_4,\phi_5,\phi_6)} = \cos(\phi_2) \cos(\phi_3)^2 \cos(\phi_5)\cdot \prod\limits_{1 \le i < j \le 4} (l_i-l_j) \end{equation} where \begin{equation} \left( c_{i,j} \right)_{i,j=1}^4 = H L H^T \end{equation} Note, that we have guessed the result and verified it by numerical calculations. The expression inside the exponential reads: \begin{eqnarray} &&Tr\left[\Sigma^{-1} H L H^T\right] =\\ &&\sum\limits_{j=1}^4 {\mathcal A}_j l_j + \sum\limits_{i_1=0}^2 \prod\limits_{\xi=2}^6 \sum\limits_{i_\xi=-2}^2 \left( {\mathcal X}_{\vec{i}} \cos(\vec{i} \cdot \vec{\phi}) + {\mathcal Y}_{\vec{i}} \sin(\vec{i} \cdot \vec{\phi}) \right) \end{eqnarray} where $\vec{\phi} = \left(\phi_j\right)_{j=1}^6$ and $\vec{i} := \left(i_\xi\right)_{\xi=1}^6$. The coefficients in the first term on the right hand side read: \begin{eqnarray} {\mathcal A}_1 &=& \frac{1}{8 s_1} + \frac{1}{8 s_2} + \frac{1}{4 s_3} + \frac{1}{2 s_4} \\ {\mathcal A}_2 &=& \frac{1}{4 s_1} + \frac{1}{4 s_2} + \frac{1}{4 s_3} + \frac{1}{4 s_4} \\ {\mathcal A}_3 &=& \frac{5}{16 s_1} + \frac{5}{16 s_2} + \frac{1}{4 s_3} + \frac{1}{8 s_4} \\ {\mathcal A}_4 &=& \frac{5}{16 s_1} + \frac{5}{16 s_2} + \frac{1}{4 s_3} + \frac{1}{8 s_4} \end{eqnarray} and we have: \begin{equation} \sum\limits_{j=1}^4 {\mathcal A}_j = \frac{1}{s_1} + \frac{1}{s_2} + \frac{1}{s_3} + \frac{1}{s_4} = \frac{1}{3!} \frac{\left(Tr\left[\Sigma\right]^3 - 3 Tr\left[\Sigma\right]Tr\left[\Sigma^2\right] + 2 Tr\left[\Sigma^3\right]\right)}{\det(\Sigma) } \end{equation} as it should be.

Therefore the integral reads: \begin{eqnarray} &&I = \left|\prod\limits_{1 \le i < j \le 4} (l_i-l_j)\right| e^{-\frac{1}{2} \sum\limits_{j=1}^4 {\mathcal A}_j l_j } \cdot \\ &&\int\limits_{[0,2\pi]^6} \left|\cos(\phi_2) \cos(\phi_3)^2 \cos(\phi_5)\right| \exp\left((-\frac{1}{2})\sum\limits_{i_1=0}^2 \prod\limits_{\xi=2}^6 \sum\limits_{i_\xi=-2}^2 \left( {\mathcal X}_{\vec{i}} \cos(\vec{i} \cdot \vec{\phi}) + {\mathcal Y}_{\vec{i}} \sin(\vec{i} \cdot \vec{\phi}) \right)\right) \prod\limits_{j=1}^6 d\phi_j \end{eqnarray} where the coefficients ${\mathcal X}_{\vec{i}}$ and ${\mathcal Y}_{\vec{i}}$ along with the number #$(\vec{i})$ of corresponding index values are given in the tables below:

\begin{array}{cc} {\mathcal X}_{\vec{i}} & \#(\vec{i}) \\ -\frac{\left(l_3-l_4\right) \left(s_1 \left(2 s_2-s_4\right)-s_2 s_4\right)}{16 s_1 s_2 s_4} & 1 \\ -\frac{\left(2 l_2-l_3-l_4\right) \left(s_1 \left(2 s_2-s_4\right)-s_2 s_4\right)}{16 s_1 s_2 s_4} & 1 \\ -\frac{\left(2 l_2-l_3-l_4\right) \left(s_1 \left(2 s_2-s_3\right)-s_2 s_3\right)}{32 s_1 s_2 s_3} & 1 \\ \frac{\left(4 l_1-2 l_2-l_3-l_4\right) \left(s_2 s_3 s_4+s_1 \left(s_3 s_4+s_2 \left(2 s_4-4 s_3\right)\right)\right)}{32 s_1 s_2 s_3 s_4} & 1 \\ -\frac{\left(2 l_1-l_3-l_4\right) \left(s_1 \left(2 s_2-s_3\right)-s_2 s_3\right)}{16 s_1 s_2 s_3} & 1 \\ -\frac{\left(2 l_1-l_3-l_4\right) \left(s_1-s_2\right)}{16 s_1 s_2} & 1 \\ -\frac{\left(l_3-l_4\right) \left(s_1 \left(2 s_2-s_4\right)-s_2 s_4\right)}{32 s_1 s_2 s_4} & 2 \\ \frac{3 \left(l_3-l_4\right) \left(s_1 \left(2 s_2-s_3\right)-s_2 s_3\right)}{64 s_1 s_2 s_3} & 2 \\ \frac{\left(l_3-l_4\right) \left(s_2 s_3 s_4+s_1 \left(s_3 s_4+s_2 \left(2 s_4-4 s_3\right)\right)\right)}{64 s_1 s_2 s_3 s_4} & 2 \\ \frac{\left(2 l_2-l_3-l_4\right) \left(s_2 s_3 s_4+s_1 \left(s_3 s_4+s_2 \left(2 s_4-4 s_3\right)\right)\right)}{64 s_1 s_2 s_3 s_4} & 2 \\ \frac{\left(l_3-l_4\right) \left(s_1 \left(2 s_2-s_3\right)-s_2 s_3\right)}{32 s_1 s_2 s_3} & 2 \\ \frac{\left(2 l_2-l_3-l_4\right) \left(s_1 \left(2 s_2-s_3\right)-s_2 s_3\right)}{32 s_1 s_2 s_3} & 2 \\ \frac{\left(2 l_2-l_3-l_4\right) \left(s_1-s_2\right)}{32 s_1 s_2} & 2 \\ -\frac{\left(4 l_1-2 l_2-l_3-l_4\right) \left(s_1 \left(2 s_2-s_3\right)-s_2 s_3\right)}{64 s_1 s_2 s_3} & 2 \\ \frac{3 \left(2 l_2-l_3-l_4\right) \left(s_1-s_2\right)}{64 s_1 s_2} & 2 \\ -\frac{\left(4 l_1-2 l_2-l_3-l_4\right) \left(s_1-s_2\right)}{64 s_1 s_2} & 2 \\ -\frac{\left(2 l_1-l_3-l_4\right) \left(s_1-s_2\right)}{32 s_1 s_2} & 2 \\ \frac{\left(2 l_2-l_3-l_4\right) \left(s_1 \left(2 s_2-s_3\right)-s_2 s_3\right)}{64 s_1 s_2 s_3} & 4 \\ \frac{\left(l_3-l_4\right) \left(s_2 s_3 s_4+s_1 \left(s_3 s_4+s_2 \left(2 s_4-4 s_3\right)\right)\right)}{128 s_1 s_2 s_3 s_4} & 4 \\ \frac{\left(l_3-l_4\right) \left(s_1 \left(2 s_2-s_3\right)-s_2 s_3\right)}{64 s_1 s_2 s_3} & 4 \\ \frac{3 \left(l_3-l_4\right) \left(s_1 \left(2 s_2-s_3\right)-s_2 s_3\right)}{128 s_1 s_2 s_3} & 4 \\ -\frac{9 \left(l_3-l_4\right) \left(s_1-s_2\right)}{128 s_1 s_2} & 4 \\ -\frac{\left(l_3-l_4\right) \left(s_1-s_2\right)}{128 s_1 s_2} & 4 \\ \frac{3 \left(2 l_2-l_3-l_4\right) \left(s_1-s_2\right)}{128 s_1 s_2} & 4 \\ \frac{\left(2 l_2-l_3-l_4\right) \left(s_1-s_2\right)}{64 s_1 s_2} & 4 \\ -\frac{\left(4 l_1-2 l_2-l_3-l_4\right) \left(s_1-s_2\right)}{128 s_1 s_2} & 4 \\ \frac{\left(l_3-l_4\right) \left(s_1 \left(2 s_2-s_3\right)-s_2 s_3\right)}{128 s_1 s_2 s_3} & 8 \\ -\frac{\left(l_3-l_4\right) \left(s_1-s_2\right)}{32 s_1 s_2} & 8 \\ \frac{\left(2 l_2-l_3-l_4\right) \left(s_1-s_2\right)}{128 s_1 s_2} & 8 \\ \frac{\left(l_3-l_4\right) \left(s_1-s_2\right)}{128 s_1 s_2} & 8 \\ -\frac{\left(2 l_2-l_3-l_4\right) \left(s_1 \left(2 s_2-s_3\right)-s_2 s_3\right)}{128 s_1 s_2 s_3} & 8 \\ -\frac{\left(2 l_2-l_3-l_4\right) \left(s_1 \left(2 s_2-s_3\right)-s_2 s_3\right)}{64 s_1 s_2 s_3} & 8 \\ \frac{\left(l_3-l_4\right) \left(s_1-s_2\right)}{32 s_1 s_2} & 10 \\ \frac{\left(l_3-l_4\right) \left(s_1-s_2\right)}{256 s_1 s_2} & 16 \\ -\frac{\left(l_3-l_4\right) \left(s_1-s_2\right)}{64 s_1 s_2} & 16 \\ -\frac{\left(l_3-l_4\right) \left(s_1 \left(2 s_2-s_3\right)-s_2 s_3\right)}{256 s_1 s_2 s_3} & 16 \\ \frac{3 \left(l_3-l_4\right) \left(s_1-s_2\right)}{256 s_1 s_2} & 16 \\ -\frac{\left(2 l_2-l_3-l_4\right) \left(s_1-s_2\right)}{256 s_1 s_2} & 16 \\ -\frac{\left(l_3-l_4\right) \left(s_1 \left(2 s_2-s_3\right)-s_2 s_3\right)}{128 s_1 s_2 s_3} & 16 \\ -\frac{\left(2 l_2-l_3-l_4\right) \left(s_1-s_2\right)}{128 s_1 s_2} & 16 \\ \frac{\left(l_3-l_4\right) \left(s_1-s_2\right)}{64 s_1 s_2} & 24 \\ -\frac{\left(l_3-l_4\right) \left(s_1-s_2\right)}{512 s_1 s_2} & 32 \\ -\frac{\left(l_3-l_4\right) \left(s_1-s_2\right)}{256 s_1 s_2} & 32 \\ \end{array} and \begin{array}{cc} {\mathcal Y}_{\vec{i}} & \#(\vec{i}) \\ -\frac{\left(l_3-l_4\right) \left(s_1 \left(2 s_2-s_3\right)-s_2 s_3\right)}{32 s_1 s_2 s_3} & 2 \\ \frac{\left(l_3-l_4\right) \left(s_1 \left(2 s_2-s_3\right)-s_2 s_3\right)}{32 s_1 s_2 s_3} & 2 \\ -\frac{\left(2 l_2-l_3-l_4\right) \left(s_1-s_2\right)}{32 s_1 s_2} & 2 \\ \frac{\left(2 l_2-l_3-l_4\right) \left(s_1-s_2\right)}{32 s_1 s_2} & 2 \\ \frac{3 \left(l_3-l_4\right) \left(s_1-s_2\right)}{64 s_1 s_2} & 8 \\ -\frac{3 \left(l_3-l_4\right) \left(s_1-s_2\right)}{64 s_1 s_2} & 8 \\ \frac{\left(l_3-l_4\right) \left(s_1 \left(2 s_2-s_3\right)-s_2 s_3\right)}{64 s_1 s_2 s_3} & 12 \\ -\frac{\left(l_3-l_4\right) \left(s_1 \left(2 s_2-s_3\right)-s_2 s_3\right)}{64 s_1 s_2 s_3} & 12 \\ -\frac{\left(2 l_2-l_3-l_4\right) \left(s_1-s_2\right)}{64 s_1 s_2} & 12 \\ \frac{\left(2 l_2-l_3-l_4\right) \left(s_1-s_2\right)}{64 s_1 s_2} & 12 \\ \frac{\left(l_3-l_4\right) \left(s_1-s_2\right)}{128 s_1 s_2} & 48 \\ -\frac{\left(l_3-l_4\right) \left(s_1-s_2\right)}{128 s_1 s_2} & 48 \\ \end{array} Note 1: The integral above is again,per definition, equal to some generalized Bessel function.

Note 2: We refrain from writing out all the indices $\vec{i}$ explicitly since there is no room for this in here and, in addition, it is still hard to see any pattern in them.

Note 3: Interestingly enough the coefficients are functions of differences of eigenvalues only rather than the eigenvalues themselves. In other words the coefficients, and what follows the whole integral, depend on $\left(l_i-l_j\right)_{1\le i < j \le 4}$ and on $\left(1/s_i - 1/s_j\right)_{1\le i < j \le 4}$ only.