Two-level correlation function of eigenvalues for large random matrices One can define the density of eigenvalues of a $N\times N$ Hermitian random matrix $H$ as:
\begin{equation}
\rho(\lambda)=\left \langle\frac{1}{N} \operatorname{Tr} \delta(\lambda-H)\right\rangle
\end{equation}
Where $\langle \dots \rangle$ denotes the average over the distribution of $H$.
In the large $N$ limit, it is famously known that $\rho(\lambda)$ will approach the Wigner semi circular law (given some conditions on the moments of the distribution of $H$). This can be shown with various methods, one that I favor is using the resolvent and computing its schur complements.
One can define the two-level correlation function:
\begin{equation}
\rho^{(2)}(\lambda, \mu)=\left \langle\frac{1}{N} \operatorname{Tr} \delta(\lambda-H) \frac{1}{N} \operatorname{Tr} \delta(\mu-H)\right\rangle
\end{equation}
The following paper 1 provides a method to compute this quantity, and provides an exhaustive list of the existing methods to compute eq $(2)$. However the paper is now $25$ years old: is there any known results that extends their method? Can we compute this quantity using the schur complement and the resolvent? What about non-hermitian matrices?
Edit:
My initial query was whether we could could compute equation $(2)$ using the resolvents and the Schur complements.
The paper mentions that contrary to $\rho(\lambda)$ there is no universality for $\rho^{(2)}(\lambda, \mu)$: it would depend on the choice of distribution. However in a certain regime, with large $N$ and small $\lambda- \mu$ then universal properties can be derived. Has this been we shown this using the Schur complements of the resolvents?
1 : Brézin, E., & Hikami, S. (1996). Correlations of nearby levels induced by a random potential. Nuclear Physics B, 479(3), 697-706. link: https://arxiv.org/pdf/cond-mat/9605046.pdf
 A: Let me address the issue raised by the OP of the universality of the two-point correlation function.
The universality of $\rho^{(2)}(\lambda,\mu)$ does exist if one considers the correlations locally, on the scale of the mean eigenvalue spacing. This is relevant for many applications, because the correlations decay quickly with increasing $|\lambda-\mu|$.
If the eigenvalue distribution has the Gibbs form
$$P(\lambda_1,\lambda_2,\ldots \lambda_N)\propto e^{-\beta W},$$
$$W=\sum_{i<j}u(\lambda_i-\lambda_j)+\sum_i V(\lambda_i),$$
then the two-point correlation function
$$K(\lambda,\mu)=N^2\rho^{(2)}(\lambda,\mu)-N^2\rho(\lambda)\rho(\mu)$$
is given in the large-$N$ local limit by the functional inverse $u^{\text{inv}}$ of $u$,
$$K(\lambda,\mu)=\frac{1}{\beta}u^{\text{inv}}(\lambda,\mu).$$
See section 1.D of arXiv:cond-mat/9612179 .
Notice that this is independent of the "potential" $V$, only the eigenvalue interaction $u$ enters. That is the sense in which the two-point correlation function is "universal" on the local scale for a broad class of RMT ensembles with a logarithmic eigenvalue repulsion, $u(\lambda,\mu)=-\log|\lambda-\mu|$, see for example arXiv:cond-mat/9310010.
