# 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=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.

• If you expand $e^{{\rm Tr}(J^TM)}=\sum_n \frac{1}{n!}[{\rm Tr}(J^TM)]^n$, you can use the results discussed in this question – Marcel Jun 7 '18 at 11:29
• @Marcel : Thank you for the reference. In addition of being quite hard to decipher for a physicist like me, I was hopping to get an explicit solution (for $N=3$, say). This solution exists at least for $N=2$: $Z_2=I_0\left(\sqrt{{\rm Tr}(J J^T)+2\det J}\right)$, with $I_\nu$ the modified Bessel function of the first kind. It is quite unclear to me how to recover this results from the answers to the question you linked to. – Adam Jun 7 '18 at 12:04

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).$$

• There is at least a closed form for $N=2$, see the edit part of my question, and my hope is that it can be compute for some specific $N$ (say $N=3$). Thanks for the large $N$ limit, though. – Adam Jun 7 '18 at 14:07
• Sure, it is not hard to find an integral form of $Z_3$, but I'm looking for a more explicit form. One of the main issue with this representation is that the symmetry under permutation of $\lambda_i$ is not explicit (this would be needed for further calculation), see also math.stackexchange.com/questions/2808873/… – Adam Jun 7 '18 at 14:34
• Your integral for $Z_3$ seems to give $Z_3(0.3,0.4,0.5)\neq Z_3(0.4,0.3,0.5)$, but it seems clear that the true $Z_3$ should be symmetric in $\lambda$. – Neil Strickland Jun 7 '18 at 16:30
• @NeilStrickland --- I corrected a typo (the integrals over $\alpha$ and $\alpha'$ should run from $0$ to $2\pi$, not from $0$ to $\pi$), and now it is symmetric; thanks for catching this. – Carlo Beenakker Jun 7 '18 at 20:47
• Using this representation, one can also compute the case $\lambda_3=0$ and $\lambda_1=\lambda_2=\lambda$, which reads $I_0(2\lambda)+\frac\pi2 L_1(2\lambda) I_0(2\lambda)-\frac\pi2 L_0(2\lambda) I_1(2\lambda)$ with $L_\nu(x)$ the Struve L function. – Adam Jun 11 '18 at 14:49

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.

• Yes, written in this way, the integral can be seen to be invariant under permutation. But unfortunately, it does not help to write $Z$ as a function explicitely invariant (i.e. without using a change of variable in the integral). I was hoping that maybe a clever parametrization of SO(3) might do the job (see also this question : math.stackexchange.com/questions/2808873/…), but I haven't got any luck about that yet. – Adam Jun 8 '18 at 8:31