I am interested in the generating function of $SO(N)$ random matrix, that is, I want to compute $$ Z_N[J]=\int dM e^{{\rm Tr} (J^T M)}, $$ where $dM$ is the $SO(N)$ Haar measure, and $J$ is an arbitrary $N\times N$ matrix. From this generating function, I can generate all correlations $\langle M_{ij}M_{kl}\cdots\rangle$ by taking derivatives with respect to the elements of $J$.

Due to the invariance of the measure, one sees that $Z[J]=Z[U^TJV]$ with $U,V\in SO(N)$, and thus $Z$ only depends on the singular values of $J$. (Stated otherwise, $Z$ only depends on the $N$ invariants ${\rm Tr}((J^TJ)^n)$, $n=1,...,N$).

Finally, at least for $N=2$ and $N=3$, one can show that $Z$ is also invariant under permutations of the singular values of $J$ (maybe it can be generalized for all $N$ ?).

It is not too hard to compute explicitly $Z_2[J]$, which is given in terms of a Bessel function of the sum of the two singular values of $J$.

Is there a way (or has it been done in the literature) to compute $Z_N$ for any $N$ ? I would already be happy with $Z_3$, which I cannot manage to compute explicitly.

EDIT : Here is an attempt, which kind of works for $N=2$, but for which I am stuck for $N=3$. If we define a Laplacian $\Delta=\sum_{ij}\frac{\partial^2}{\partial J_{ij}^2}$ (with $J_{ij}$ the elements of $J$), one shows easily that $$ \Delta Z[J]=N Z[J]. \tag{1} $$ If we call $\lambda_i$ the singular values of $J$ (with $\lambda_1>\lambda_2>\ldots$), using the fact that $Z[J]=Z[\lambda_1,\lambda_2,\ldots]$, one shows (at least for $N=2$ and $N=3$, but it might be generalizable to $N\geq4$) that $$ \Delta Z=\frac{1}{D}\sum_{i}\frac{\partial}{\partial \lambda_i}\left(D\frac{\partial}{\partial \lambda_i}Z\right), $$ where $D=\prod_{i< j}(\lambda_i^2-\lambda_j^2)$ is related to the Jacobian to go from $J_{ij}$ to $\lambda_i$. This equation looks nice enough, so my hope is that a solution exists, I am not quite sure how to find it for $N=3$.

In the case $N=2$, we can compute $Z_2$ exactly via its definition, and it reads $Z_2[\lambda_1,\lambda_2]=I_0(\lambda_1+\lambda_2)$, with $I_\nu$ the modified Bessel function of the first kind. One checks that this is indeed a solution of Eq. (1).

Unfortunately, even in that case, it is not clear to me how to find this solution starting from Eq. (1) only. Given that $D=\lambda_1^2-\lambda_2^2$, it is tempting to defined $u=(\lambda_1+\lambda_2)/2$ and $v=(\lambda_1-\lambda_2)/2$. Then Eq. (1) is solved by separation of variables and we find a family of solution $Z_{2,\mu}$ (where I have already use the fact that $Z_N[0]=1$) : $$ Z_{2,\mu}[u,v]=I_0(\sqrt{\mu}u)I_0(\sqrt{4-\mu}v). $$ Clearly the solution to my problem corresponds to $\mu=4$, but it is not clear to me what is the rigorous argument to pick this value of $\mu$ (since $Z_N>0$ $\forall J$, we must have $\mu\geq 4$ as $I_0$ can be negative for imaginary variables; but how to select $\mu=4$ as the only viable solution ?).

One way to solve this issue of $\mu$ is to use the fact that $Z_2[J=\lambda Id_2]$ can be computed explicitly: $Z_2[J=\lambda Id_2]=I_0(2\lambda)$, which unambiguously selects $\mu=4$.

Since we can always compute the expansion of $Z_N[\lambda Id_N]$ explicitly at least for small $\lambda$, this kind of argument might be enough to fix the constant also for $N>2$. For $N=3$, a few special cases can be computed explicitly, which might help too.

Any insight for the solution of Eq. (1) would be greatly appreciated.