Riemann-Hilbert approach to Selberg integral

Question

I am interested in matrix integrals, and I have seen many mentions to a certain Riemann-Hilbert approach that indicate that this is a very powerful tool to can be used in this area, when coupled with the theory of orthogonal polynomials. I would like to understand it, but the sources I found were hard to follow.

It would be helpful if I could see it in action in a simple case, so I suggest what follows. Consider the matrix integral $$ f(a,N)=\int_{0}^\infty e^{-a\sum_{m=1}^\infty \frac{1}{m}{\rm Tr}(X^m)}|\Delta(X)|^2dX,$$ where $a>0$, $\Delta(X)$ is the Vandermonde and $X$ is diagonal of dimension $N$. I would like to know if this kind of integral is amenable to the R-H approach (with a non-polynomial potential).

The series in the exponent is only finite if $0< X< 1$, in which case it gives $-{\rm Tr}\log(1-X)$ so using that $e^{-\infty}=0$ (I am not sure how rigorous this can be made): $$f(a,N)=\int_{0}^1 \det(1-X)^a|\Delta(X)|^2dX,$$ which is a particular case of the Selberg integral and there is an explicit solution to it, which is $$\prod_{j=0}^{N-1}\frac{\Gamma(a+j+1)j!(j+1)!}{\Gamma(a+N+j+1)}.$$

My question is: how can the R-H approach be applied to the first integral in order to produce the Selberg result?

Answer 1

Let me formulate the problem in a slightly more general way: We seek to evaluate the large-$N$ limit of the matrix integral $$\int e^{-\beta\,{\rm Tr}\,V(X)}|\Delta(X)|^\beta dX\equiv e^{-\beta N^2 F},$$ integrated over $N\times N$ Hermitian matrices $X$. In the OP the index $\beta=2$ and $V(X)=(a/2)\sum_m m^{-1}X^m$, which creates convergence problems, here I assume $V(\lambda)$ is a well defined function of the eigenvalues $\lambda$ of $X$. My goal here is to reduce the calculation of the matrix integral in the large-$N$ limit to the solution of a nonlinear integral equation which can be solved by the Riemann-Hilbert technique.

_{Note for the OP: Selberg integrals are exact finite-$N$ results. The Riemann-Hilbert technique evaluates these integrals in the large-$N$ limit, by the method of stationary phase (or steepest descent). There is therefore no direct connection between the two techniques.}

Because of the invariance under unitary transformations the integral over the matrix elements can be reduced to an integral over the eigenvalues, $$e^{-\beta N^2 F}=\int_{-\infty}^\infty d\lambda_1\cdots \int_{-\infty}^\infty d\lambda_N \,e^{-\beta U(\lambda_1,\ldots\lambda_N)},$$ $$U(\lambda_1,\ldots\lambda_N)=\sum_{i=1}^N V(\lambda_i)-\sum_{i<j}\ln|\lambda_i-\lambda_j|,$$ where I have used that $\Delta(X)=\prod_{i<j}|\lambda_i-\lambda_j|$.

In the large-$N$ limit the integral can be evaluated by the method of stationary phase: $$F=\lim_{N\rightarrow\infty} N^{-2} U(\lambda^\ast_1,\ldots\lambda^\ast_N),$$ where the $\lambda^*_i$'s satisfy the stationarity equations $$\frac{\partial U}{\partial\lambda_i}=0\Rightarrow V'(\lambda_i)=\sum_{j\neq i}\frac{1}{\lambda_i-\lambda_j},\;\;i=1,2,\ldots N.$$ Intuitively, this equation expresses a condition of mechanical equilibrium of $N$ particles on a line, coordinates $\lambda_i$, moving in a potential $V(\lambda)$ and repelling each other pairwise with a logarithmic interaction.

The next step is replace the sum over the discrete index $i$ by an integral over the continuous variable $x$. For that purpose we define a function $\lambda(x)$ by $\lambda_i=N^{1/2}\lambda(i/N)$ and we define the rescaled potential $$V_\infty(\lambda)=\lim_{N\rightarrow\infty}N^{-1}V(N^{1/2}\lambda),$$ assuming that this limit exists. Upon replacement of $\sum_i f(\lambda_i)\mapsto N\int dx\, f(N^{1/2}\lambda(x))$ we thus arrive at $$F=\int_{-\infty}^\infty V_\infty(\lambda(x))\,dx-\tfrac{1}{2}\int_{-\infty}^\infty dx\int_{-\infty}^\infty dy\,[\ln |\lambda(x)-\lambda(y)|+\tfrac{1}{2}\ln N].$$

The unknown function $\lambda(x)$ is determined by the stationarity equation $$\frac{d}{d\lambda}V_\infty(\lambda(x))=\int_{-\infty}^{\infty}dy\frac{1}{\lambda(x)-\lambda(y)}.$$ Defining the density function $\rho(\lambda)=(d\lambda/dx)^{-1}$, with support $(a,b)$, the stationarity equation can be rewritten as $$V'_\infty(\lambda)=\int_{a}^{b}d\mu\frac{\rho(\mu)}{\lambda-\mu},\;\;\lambda\in(a,b).$$ The Cauchy principal value of the integral is intended.

This integral equation can now be solved by the Rieman-Hilbert technique. For a worked out example for $V(\lambda)=\tfrac{1}{2}\lambda^2+(g/N)\lambda^4$ see the paper Planar diagrams by Brézin, Itzykson, Parisi, and Zuber.

---

For the quadratic case, $V(\lambda)=\tfrac{1}{2}\lambda^2=V_\infty(\lambda)$, the solution of the integral equation is $$\rho(\lambda)=\frac{1}{\pi}\sqrt{2-\lambda^2},\;\;\lambda\in(-\sqrt 2,\sqrt 2),$$ which gives $$F=\int_{-\sqrt 2}^{\sqrt 2}\tfrac{1}{2}\lambda^2\rho(\lambda)\,d\lambda-\tfrac{1}{2}\int_{-\sqrt{2}}^{\sqrt{2}}d\lambda\int_{-\sqrt 2}^{\sqrt 2} d\mu\,\rho(\lambda)\rho(\mu)[\ln|\lambda-\mu|+\tfrac{1}{2}\ln N]$$ $$\qquad=\frac{3+\ln 4}{8}-\frac{1}{4}\ln N.$$ This can then be compared with the corresponding Selberg integral, which reads $$S_N=\int e^{-\tfrac{1}{2}\beta\,{\rm Tr}\,X^2}|\Delta(X)|^\beta dX=(2\pi)^{N/2}\beta^{-N/2-\beta N(N-1)/4}\prod_{j=1}^N\frac{\Gamma(1+\beta j/2)}{\Gamma(1+\beta/2)}.$$ So I should check that $$\lim_{N\rightarrow\infty}N^{-2}\ln S_N=-\frac{3+\ln 4 }{8}\beta+\frac{1}{4}\beta\ln N.\;\;(1)$$ I find $$\lim_{N\rightarrow\infty}N^{-2}\ln S_N=-\frac{\beta}{4}\ln \beta+\lim_{N\rightarrow\infty}N^{-2}\sum_{j=1}^N\ln\Gamma(1+\beta j/2)$$ $$=-\frac{\beta}{4}\ln \beta+\lim_{N\rightarrow\infty}N^{-2}\int_1^N \tfrac{1}{2} j \bigl( \beta \ln \beta j- \beta-\beta \ln 2\bigr)\,dj$$ $$=-\frac{3+\ln 4}{8}\beta+\frac{1}{4}\beta\ln N.\;\;(2)$$ It agrees.