Riemann-Hilbert approach to Selberg integral 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?
 A: 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.
