Generating function of $SO(N)$ random matrix 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.
 A: For the $n=3$ case, it is best to use quaternionic methods to express $SO(3)$ as a quotient of $S^3$.  Explicitly, for $a\in S^3$ one can check that the matrix 
$$ A = \left[\begin{matrix}
    a_1^2-a_2^2-a_3^2+a_4^2 & 2a_1a_2-2a_3a_4 & 2a_1a_3+2a_2a_4 \\
    2a_1a_2+2a_3a_4 & -a_1^2+a_2^2-a_3^2+a_4^2 & -2a_1a_4+2a_2a_3 \\
    2a_1a_3-2a_2a_4 & 2a_1a_4+2a_2a_3 & -a_1^2-a_2^2+a_3^2+a_4^2
   \end{matrix}\right]
$$
lies in $SO(3)$.  (This is based on an identification $\mathbb{H}\simeq\mathbb{R^4}$ with $1\in\mathbb{H}$ corresponding to $(0,0,0,1)\in\mathbb{R}^4$.)  This construction comes from a group homomorphism, so the Haar measure on $SO(3)$ corresponds to the natural round measure on $S^3$, suitably rescaled.  Thus, we just want to integrate the function $\exp(f)$ over $S^3$, where
$$ f(a) = (a_1^2-a_2^2-a_3^2+a_4^2)\lambda_1 + 
          (-a_1^2+a_2^2-a_3^2+a_4^2)\lambda_2 +
          (-a_1^2-a_2^2+a_3^2+a_4^2)\lambda_3. 
$$
Using $\sum_ia_i^2=1$, we get $f(a)=\sum_i\lambda_i-2g(a)$, where
$$ g(a) = (a_2^2+a_3^2)\lambda_1 + (a_1^2+a_3^2)\lambda_2 + (a_1^2+a_2^2)\lambda_3.$$
This is at least visibly invariant under permutations of $\lambda$.  Next, because $a_4$ does not appear in $g(a)$, it is natural to use a stereographic parametrisation $\sigma\colon\mathbb{R}^3\to S^3$, as follows:
$$ \sigma(x_1,x_2,x_3) = \frac{(2x_1,2x_2,2x_3,\|x\|^2-1)}{\|x\|^2+1} $$
The Jacobian determinant of $\sigma$ turns out to be $8/(\|x\|^2+1)^3$.  This leads to an expression for $Z_3$ as an integral over $\mathbb{R}^3$ that retains all the natural symmetries, but it does not seem easy to evaluate either numerically or symbolically.
A: Singular value decomposition: The real $N\times N$ matrix $J$ has singular values $\sigma_1,\sigma_2\ldots\sigma_N\geq 0$ in the decomposition $J=U\,{\rm diag}\,(\sigma_1,\sigma_2\ldots\sigma_N)V$ with $U,V\in{\rm O}(N)$. We need to restrict $U,V$ to ${\rm SO}(N)$, which we can do by allowing for one of the singular values to become negative:
$$J=U\,{\rm diag}\,(\lambda_1,\lambda_2\ldots\lambda_N)V,\;\;U,V\in{\rm SO}(N),$$
$$\lambda_n=\sigma_n\geq 0,\;\;{\rm for}\;\;n=1,2,\ldots N-1,$$
$$\lambda_N=\begin{cases}
\sigma_N\;{\rm if}\;{\rm det}\,M\geq 0,\\
-\sigma_N\;{\rm if}\;{\rm det}\,M<0.
\end{cases}$$
Note that this need to consider the $\lambda_n$'s implies that the integral $Z_N$ is not only a function of the invariants ${\rm Tr}((J^TJ)^n)$.   
Without loss of generality we may now take a diagonal $J={\rm diag}\,(\lambda_1,\lambda_2,\ldots \lambda_N)$ and then we need the integral
$$Z_N=\int_{{\rm O}(N)} dM \exp\left(\sum_{n=1}^N \lambda_n M_{nn}\right).$$
Large-$N$ limit: I do not think there is a closed form expression for arbitrary $N$, but for $N\gg 1$ the diagonal matrix elements of $M$ are independent Gaussians (mean zero, variance $1/N$), resulting in
$$Z_N\rightarrow \exp\left(\frac{1}{2N}\sum_{n=1}^N \lambda_n^2\right),\;\;N\gg 1.$$
This result for $Z_N$ in the large-$N$ limit depends only on the singular values (because $\lambda_n$ enters only quadratically), but for small-$N$ this is no longer true.    
Small-$N$ results:
For small $N$ we can use the parametrization of the orthogonal group given in Appendix B of arXiv:1405.3115. 
• $N=2$
$$Z_2=\frac{1}{\pi}\int_0^{\pi} \exp\bigl[(\lambda_1+\lambda_2)\cos\theta\bigr]\,d\theta=I_0(\lambda_1+\lambda_2)$$
Note that in terms of the singular values considered in the OP, this would read $Z_2=I_0(\sigma_1\pm\sigma_2)$, where the $\pm$ sign is the sign of ${\rm det}\,J$.
• $N=3$
$$Z_3=\frac{1}{8\pi^2}\int_0^{2\pi} d\alpha\int_0^{2\pi} d\alpha'\int_0^\pi \sin\theta\, d\theta$$
$$\qquad\qquad\times\exp\left[\cos \alpha \cos \alpha' (\lambda_1+\lambda_2 \cos \theta)-\sin \alpha \sin \alpha' (\lambda_1 \cos \theta+\lambda_2)+\lambda_3\cos \theta\right]$$
The two integrals over $\alpha,\alpha'$ can be carried out (by rewriting them as integrals over $\alpha\pm\alpha'$), with the result
$$Z_3=\frac{1}{2}\int_{-1}^1   I_0\left[\tfrac{1}{2}(\lambda_1-\lambda_2) (1-x)\right] I_0\left[\tfrac{1}{2} (\lambda_1+\lambda_2)(1+x)\right]\,e^{\lambda_3 x}\,dx.$$
As shown by Paul Enta, this integral can be written as the inverse Laplace transform
${\cal L}^{-1}[F(s)](t)$ for $t=1$ of the function 
$$F_3(s)=\frac{1}{\sqrt{(s-\lambda_1-\lambda_2-\lambda_3) (s+\lambda_1+\lambda_2-\lambda_3) (s+\lambda_1-\lambda_2+\lambda_3) (s-\lambda_1+\lambda_2+\lambda_3)}},$$
to demonstrate that it is invariant under permutation of the $\lambda_n$'s.
Notice also that $F_3(s)$ is invariant under a sign change of two of the $\lambda_n$'s, but not under a sign change of one single $\lambda_n$. This is consistent with the definition of the $\lambda_n$'s in terms of the singular values, given above. 
A closed-form expression of either the Bessel-function integral or the inverse Laplace transform seems not forthcoming, except for some special cases (see the comments by Adam).
Average over the full orthogonal group: If we integrate over the full group ${\rm O}(N)$, instead of only over ${\rm SO}(N)$, we should combine the results for $\pm\lambda_N$,
$$Z_{{\rm O}(N)}(\lambda_1,\lambda_2,\ldots \lambda_{N-1},\lambda_N)=
\tfrac{1}{2}Z_{{\rm SO}\,(N)}(\lambda_1,\lambda_2,\ldots \lambda_{N-1},\lambda_N)+
\tfrac{1}{2}Z_{{\rm SO}(N)}(\lambda_1,\lambda_2,\ldots \lambda_{N-1},-\lambda_N).$$
There is now no need to distinguish the $\lambda_n$'s from the singular values $\sigma_n$. In particular, for $N=2$ one has
$$Z_{{\rm O}(2)}=\tfrac{1}{2}I_0(\sigma_1+\sigma_2)+\tfrac{1}{2}I_0(\sigma_1-\sigma_2).$$
