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

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 2

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 3

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 4

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

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.

Answer 6

Let us change the notation $p \rightarrow d$. Here we provide an answer for arbitrary $d\ge 3$ in case when the spectrum of the matrix $\Sigma$ is degenerate such that only two eigenvalues are distinct, i.e. $s_1=s_2=\cdots=s_{d-2}$ and $s_{d-1} \neq s_{1}$ and $s_{d}\neq s_1$.

The matrix of rotations reads: \begin{eqnarray} &&H= \prod\limits_{i=1}^d \prod\limits_{j=i+1}^d\\ &&\left( \delta_{\xi,\eta} + (\cos(\phi_{\nu_{i,j}})-1) \delta_{(\xi,\eta),(i,i)} - \sin(\phi_{\nu_{i,j}}) \delta_{(\xi,\eta),(i,j)} + \sin(\phi_{\nu_{i,j}}) \delta_{(\xi,\eta),(j,i)} + (\cos(\phi_{\nu_{i,j}})-1) \delta_{(\xi,\eta),(j,j)} \right)_{\xi=1,\eta=1}^{d,d} \end{eqnarray} where $\nu_{i,j}:= \binom{d}{2} - \binom{d-i+1}{2} + j-i$. where each angle $\phi_j \in [0, 2\pi]$ for $j=1,\cdots,\binom{d}{2}$.

Now, the Jacobian reads: \begin{eqnarray} &&\frac{\partial(\overbrace{c_{1,1},c_{2,2},\cdots,c_{d,d}}^{d},\overbrace{c_{1,2},\cdots,c_{1,d}}^{d-1},\overbrace{c_{2,3},\cdots,c_{2,d}}^{d-2},\cdots,\overbrace{c_{d-2,d-1},c_{d-2,d}}^2,\overbrace{c_{d-1,d}}^1)}{\partial(\underbrace{l_1,\cdots,l_d}_d,\underbrace{\phi_1,\phi_2,\cdots,\phi_{\binom{d}{2}}}_{\binom{d}{2}})}=\\ &&\prod\limits_{j=0}^{d-3} \prod\limits_{\xi=1}^{d-2-j} [\cos(\phi_{j d-\binom{j+1}{2}+1+\xi})]^\xi \cdot \prod\limits_{1 \le i < j \le d} (l_i - l_j) \end{eqnarray} where \begin{equation} \left( c_{i,j} \right)_{i,j=1}^{d} = H L H^T \end{equation}

Again, we just guessed this result by generating the Jacobi matrix in Mathematica and computing the determinant for several values of $d$. Then we checked the result by randomly sampling some integer $d$.

It turns out you can actually find closed form expressions for the matrix elements of $H$. The last $d$th row has the simplest form and then the complexity increases as one moves bottom-up along the rows. Let us define $f_{p,\xi} := \xi d - p - \binom{\xi+1}{2}$ for $p=0,1,2$ and $\xi=1,\cdots,d-1-p$. We have: \begin{eqnarray} H_{d,j} &=& \left(\prod\limits_{\xi=1}^{j-1} \cos \phi_{f_{0,\xi}}\right) \cdot \left\{\begin{array}{rr} \sin\phi_{f_{1,j}} & \mbox{if $j < d$}\\ 1 & \mbox{othw} \end{array}\right\} \\ H_{d-1,j} &=& (-1)^{\delta_{j,d}} \left(\prod\limits_{\xi=1}^{j-1-\delta_{j,d-1}} \cos\phi_{f_{2,\xi}}\right) \cdot \left\{\begin{array}{rr} \sin\phi_{f_{2,j}} & \mbox{if $j \le d-3$}\\ 1 & \mbox{othw} \end{array}\right\} \cdot \left\{\begin{array}{rr} \cos\phi_{f_{1,j}} & \mbox{if $j \le d-2$}\\ \sin\phi_{f_{1,j-1}} & \mbox{othw} \end{array}\right\}\cdot \cos\phi_{f_{0,j}}-\\ && \sum\limits_{l=1}^{j-1-\delta_{j,d}} \left(\prod\limits_{\xi=1}^{l-1} \cos\phi_{f_{1,\xi}}\right) \sin\phi_{f_{1,l}}\cdot \sin\phi_{f_{0,l}} \left(\prod\limits_{\xi=l+1}^{j-1} \cos\phi_{f_{0,\xi}}\right) \left\{\begin{array}{rr} \sin\phi_{f_{0,j}} & \mbox{if $j \le d-1$}\\ 1 & \mbox{othw} \end{array}\right\} \end{eqnarray} Interestingly enough it turns out that the bottom row $H_d$ depends on $\left\{ \phi_{f_{0,\xi}}\right\}_{\xi=1}^{d-1}$ only whereas the row just above it $H_{d-1}$ depends on $\left\{ \phi_{f_{0,\xi}}\right\}_{p=0,\xi=1}^{1,d-1-p}$ only. In other words out of the collection of $\binom{d}{2}$ angles that fully determine the orthogonal matrix the bottom row depends only on a subset of $d-1$ angles whereas the row just above it depends on a subset of $2d-3$ angles.

Now we have enough information to write down the expression in the exponential. We have: \begin{eqnarray} &&Tr\left(\Sigma^{-1} H L H^T \right) = Tr \left( L H^T \cdot \Sigma^{-1} H\right) = \\ &&\frac{1}{s_1} \left(\sum\limits_{j=1}^d l_j\right) + \left(\frac{1}{s_d} - \frac{1}{s_1} \right) \left( \sum\limits_{j=1}^d (H_{d,j})^2 l_j \right) + \left(\frac{1}{s_{d-1}} - \frac{1}{s_1} \right) \left( \sum\limits_{j=1}^d (H_{d-1,j})^2 l_j \right) \end{eqnarray}

Out of the three terms on the right hand side above it is only the last two that depend on the angles. We will handle those terms by adding and subtracting the average over a hypercube. We have:

\begin{eqnarray} \left< (H_{d,j})^2\right>_{O(d)} &=& \frac{1}{2^j} 1_{j\le d-1} + \frac{1}{2^{d-1}} 1_{j=d} \\ \left< (H_{d-1,j})^2\right>_{O(d)} &=& \frac{j}{2^{j+1}} 1_{j \le d-2} + \frac{d}{2^d} 1_{d-1 \le j \le d} \end{eqnarray}

Therefore the integral reads: \begin{eqnarray} &&\left< \exp\left( -\beta Tr\left(\Sigma^{-1} H L H^T \right) \right)\right>_{O(d)} = \left(\prod\limits_{1 \le i < j \le d} \left| l_i - l_j\right|\right) \cdot e^{-\beta \sum\limits_{j=1}^d {\mathcal A}_j l_j} \cdot \\ &&\int\limits_{[0,2\pi]^{\binom{d}{2}}} e^{-\beta (\frac{1}{s_d}-\frac{1}{s_1}) \sum\limits_{j=1}^d((H_{d,j})^2 - \left< (H_{d,j})^2\right>) l_j} \cdot e^{-\beta (\frac{1}{s_{d-1}}-\frac{1}{s_1}) \sum\limits_{j=1}^d((H_{d-1,j})^2 - \left< (H_{d-1,j})^2\right>) l_j} \cdot\\ && \prod\limits_{\xi=0}^{d-2} \prod\limits_{j=1}^{d-1-\xi} [\cos \phi_{f_{\xi,j}}]^{d-1-j-\xi} \frac{d \phi_{f_{\xi,j}}}{2 \pi} \end{eqnarray}

where \begin{eqnarray} {\mathcal A}_j &=&\frac{1}{s_1} \left( 1- \frac{j+2}{2^{j+1}} 1_{j \le d-2} - \frac{d+2}{2^d} 1_{d-1\le j \le d}\right)+ \frac{1}{s_{d-1}} \left( \frac{j}{2^{j+1}} 1_{j\le d-2} + \frac{d}{2^d} 1_{d-1\le j \le d}\right) + \frac{1}{s_d} \left( \frac{1}{2^j} 1_{j \le d-2} + \frac{1}{2^{d-1}} 1_{d-1\le j \le d}\right) \end{eqnarray} for $j=1,\cdots,d$.

Note 1: The integral runs essentially over a $(2d-3)$-dimensional hypercube since the integrand depends only on a subset $\left\{ \phi_{f_{\xi,j}} \right\}_{\xi=0,j=1}^{1,d-1-\xi}$ of all angles. The integration over the remaining angles produces some prefactor only.

Note 2: The integral over the hypercube in question depends only on differences of respective eigenvalues, i.e. on $\left(l_j - l_d\right)_{j=1}^{d-1}$ and on $(1/s_d-1/s_1)$ and $(1/s_{d-1}-1/s_1)$, rather than on the eigenvalues themselves. Indeed this follows from the fact that the sum of the coefficients at the $l$-variables in the exponential is equal to zero as seen below: \begin{equation} \sum\limits_{j=1}^d \left( (H_{d-\xi,j})^2 - \left< (H_{d-\xi,j})^2 \right>_{O(d)} \right) = 1 - 1 = 0 \end{equation} for $\xi=0,1$.