Let n=2 and C be a symmetric two by two matrix and let $X = \left(X_{1,j}; X_{2,j}\right)_{j=1}^T$ be a two by T real matrix matrix. We consider a following generating function.
\begin{equation} U(t) := \frac{1}{n} \int\limits_{{\mathbb R}^{2 T}} Tr\left[ e^{\frac{\imath t}{T} X X^T} \right] \cdot e^{-\frac{1}{2} Tr\left[ C^{-1} X X^T \right]} \frac{\prod_{j=1}^T dX_{1,j} dX_{2,j}}{\sqrt{(2\pi)^{2T}} \sqrt{\det(C)^T}} \end{equation}
By expanding the first term in the integrand in a Taylor series in t we and integrating term by term we easily establish that: \begin{eqnarray} U(0) &=& 1 \\ \left.\frac{d U(t)}{d \imath t} \right|_{t=0} &=& Tr[C] \\ \left.\frac{d^2 U(t)}{d (\imath t)^2} \right|_{t=0} &=& Tr[C^2](1+\frac{1}{T}) + Tr[C]^2 \frac{2}{T} \\ \end{eqnarray}
Now, let us try to compute the generating function in a different way. We perform a singular value decomposition of the matrix X X^T. We have: \begin{equation} X X^T = U {\mathcal D} U^{-1} \end{equation} where $U$ is a unitary matrix and ${\mathcal D}$ is a diagonal matrix with real entries $s_1,s_2$. We have: \begin{equation} U(t) = \frac{1}{((2 \pi)^2 c_1 c_2)^{T/2}n}\int\limits_{{\mathbb R}_+^2} d s_1 d s_2 (s_2-s_1)^2 (s_1 s_2)^{\frac{T-3}{2}} \left(e^{\frac{\imath t s_1}{T}}+e^{\frac{\imath t s_2}{T}}\right) \cdot \underbrace{\int\limits_{U(2)} dU e^{-\frac{1}{2}Tr\left(C^{-1} U {\mathcal D} U^{-1}\right)}}_{J} \end{equation} Here I used the formula for the joint-probability density of eigenvalues of an ensemble of random matrices (see http://en.wikipedia.org/wiki/Wishart_distribution for example). The last integral on the right-hand side is an integration over the group of two dimensional unitary matrices. We can do that integral owing to the Harish-Chandra-Itzykson-Zuber formula. We have: \begin{equation} J = {\mathcal A}_T \frac{(e^{-\frac{s_1}{2 c_1}-\frac{s_2}{2 c_2}} - e^{-\frac{s_1}{2 c_2}-\frac{s_2}{2 c_1}})}{(\frac{1}{c_1}-\frac{1}{c_2})(s_1-s_2)} \end{equation} where the proportionality constant ${\mathcal A}_T$ depends on $T$ only. Inserting the above into the equation for $U(t)$ we get: \begin{equation} U(t) = \frac{{\mathcal A}_T}{((2 \pi)^2 c_1 c_2)^{T/2}(\frac{1}{c_1}-\frac{1}{c_2})n}\int\limits_{{\mathbb R}_+^2} d s_1 d s_2 (s_2-s_1) (s_1 s_2)^{\frac{T-3}{2}} \left(e^{\frac{\imath t s_1}{T}}+e^{\frac{\imath t s_2}{T}}\right) \cdot (e^{-\frac{s_1}{2 c_1}-\frac{s_2}{2 c_2}} - e^{-\frac{s_1}{2 c_2}-\frac{s_2}{2 c_1}}) \end{equation} Now we set ${\mathcal A}_T=1$ and we use Mathematica to compute the value of our generating function and its derivatives at zero. We have: \begin{eqnarray} U(0) &=& {\mathcal B} \\ \frac{1}{\imath} U^{'}(0) &=& {\mathcal B} \frac{T+1}{2 T} Tr[C] \\ \frac{1}{\imath^2} U^{''}(0) &=& {\mathcal B} \frac{T+1}{2 T} \left[(1+\frac{1}{T})Tr[C^2] + \frac{4}{T} (Tr[C])^2 \right] \end{eqnarray}
where ${\mathcal B} = 4 \sqrt{c_1 c_2} \pi^{-T} \Gamma[(T-1)/2] \Gamma[(T+1)/2]$. These results are clearly different from those stated on the very top. What is the reason for this? I suspect that there is something wrong with the joint-denstity of eigenvalues in the integral over $s_1$ and $s_2$ but I cannot figure out what is wrong. Can anybody help me?