Disclaimer: This is a cross-post from math.stackexchange. Given that there is little activity on the subject (random-matrice) on the aformentioned site, and given that many interesting discussion on this exact subject have taken place here, I feel that this might be a suitable place for this question, althought it is not exactly 'research level'.
Many texts (e.g. 1-2-3) on random matrices start with some variation of the identity:
$$\rho_1(\lambda) = \frac{1}{\pi} \text{Im}\{\langle\text{Tr}(\mathbf{X}-\lambda\mathbf{I})^{-1}\rangle\}$$
where the brakets denote the average over the ensemble of eigvenvalues and 'Tr' is the trace of the resolvant matrix $(\mathbf{X}-\lambda\mathbf{I})^{-1}$, where $\lambda\in\mathbb{C}\backslash\mathbb{R}$. The left side is the 1 point correlation function, i.e. the density, which is also defined as
$$\rho_1(\lambda)= \big\langle\sum_i \delta(\lambda-\lambda_i)\big\rangle $$
(same ensemble). I would like to understand this identity better, so I re-did the complete derivation, as I could not find it anywhere, but I'm off by some $2N$ factor. And I cannot seem to identify a mistake. So here goes nothing. (disclaimer II: I'm a physicist, so some arguments might not be as sound as they could be, forgive me)
We first express the matrices $\mathbf{X},\mathbf{I}_{N\times N}$ and $R_{\mathbf{X}}^{-1}(\lambda):=\mathbf{X}-\lambda\mathbf{I}_{N\times N}$ in term of their spectra \begin{align} &x_{mn} = \sum_{j=1}^N\lambda_j v_m^{(j)} v_{n}^{(j)} & &\Longleftrightarrow& &\mathbf{X} = \sum_{j=1}^N\lambda_j \mathbf{v}^T_j\mathbf{v}_j &\\ &I_{mn} = \sum_{j=1}^Nv_m^{(j)} v_{n}^{(j)} & &\Longleftrightarrow& &\mathbf{I} = \sum_{j=1}^N\mathbf{v}^T_j\mathbf{v}_j&\\ &R_{\mathbf{X}}^{-1}(\lambda)_{mn} = \sum_{j=1}^Nv_m^{(j)} v_{n}^{(j)} (\lambda_j -\lambda)& &\Longleftrightarrow& &R_{\mathbf{X}}^{-1}(\lambda)= \sum_{j=1}^N \mathbf{v}^T_j\mathbf{v}_j(\lambda_j-\lambda) \end{align} Notationally, I use $\mathbf{v}_j^T$ to express the $j^{th}$ eigenvetor as a column matrice, and $v_m^{(j)}$ to denote the m$^{th}$ element of the $j^{th}$ eigenvector. Moreoever, we assume that the norm of the basis is chosen to be equal to one such that all the above identities hold.
By inspection, it is evident that the spectral representation of the resolvent must be \begin{align} R_{\mathbf{X}}(\lambda)_{mn}= \sum_{j=1}^N\frac{v_m^{(j)} v_{n}^{(j)}}{(\lambda_j-\lambda)} \quad\Longleftrightarrow\quad R_{\mathbf{X}}(\lambda)= \sum_{j=1}^N\frac{\mathbf{v}^T_j\mathbf{v}_j}{ (\lambda_j-\lambda)}, \end{align} which is a straightforward thing to verify \begin{align*} R_{\mathbf{X}}^{-1}(\lambda) &R_{\mathbf{X}}(\lambda) \\ &=\left(\sum_{j=1}^N \mathbf{v}^T_j\mathbf{v}_j(\lambda_j-\lambda)\right)\left(\sum_{k=1}^N\frac{\mathbf{v}^T_k\mathbf{v}_k}{ (\lambda_k-\lambda)}\right)\\ &=\sum_{j,k} \frac{(\lambda_j-\lambda)}{(\lambda_k-\lambda)} \bigl(\mathbf{v}_j^T\mathbf{v}_j\bigr)\mathbf{v}_k^T\mathbf{v}_k\\ &=\sum_{j,k} \frac{(\lambda_j-\lambda)}{(\lambda_k-\lambda)} \bigl( \delta_{jk}\bigr)\mathbf{v}_j^T\mathbf{v}_k\\ &=\sum_{j} \mathbf{v}_j^T\mathbf{v}_j \equiv\mathbf{I}.\\ \end{align*} The trace of the resolvent is directly obtained from the above expression for its element $(m,n)$ \begin{align} \text{Tr}R_{\mathbf{X}}(\lambda+i\varepsilon) \equiv \sum_{k=1}^N R_{\mathbf{X}}(\lambda+i\varepsilon)_{kk}=\sum_{j,k=1}^N \frac{v^{(j)}_kv^{(j)}_k}{\bigl(\lambda_j-(\lambda+i\varepsilon)\bigr)}. \end{align} where the constraint that the argument of the resolvent must be complex and not on the real line is now explicit ($\in \mathbb{C}\backslash\mathbb{R}$, where $\lambda$ is the real part, and $i\varepsilon$ the imaginary part).
The average of this trace for some properly normalized and arbitrary pdf $\rho(\lambda_1,\lambda_2,...,\lambda_N)$ is thus \begin{align} \langle \text{Tr}(R_{\mathbf{X}}(\lambda+i\varepsilon))\rangle &= \int_{-\infty}^\infty d\lambda_1 ...\int_{-\infty}^\infty d\lambda_N\ \rho(\lambda_1,...,\lambda_N)\sum_{j,k=1}^N\frac{v^{(j)}_k v^{(j)}_k}{\bigl(\lambda_j-(\lambda+i\varepsilon)\bigr)}\notag\\ &=\sum_{j,k=1}^Nv^{(j)}_kv^{(j)}_k\int_{-\infty}^\infty d\lambda_1... \int_{-\infty}^\infty d\lambda_j\ ...\int_{-\infty}^\infty d\lambda_N\ \frac{\rho(\lambda_1,...,\lambda_j,...,\lambda_N)}{\bigl(\lambda_j-(\lambda+i\varepsilon)\bigr)} \end{align}
The integrals over $d\lambda_\ell$ for $\ell\neq j$ are trivial since $\{\lambda_\ell\}$ only appear in the probability density. Introducting the marginal probability $\rho_{(1)}(x)$, we thus have \begin{align} \langle \text{Tr}(R_{\mathbf{X}}(\lambda+i\varepsilon))\rangle &=\sum_{j,k=1}^Nv^{(j)}_kv^{(j)}_k \int_{-\infty}^\infty d\lambda_j\ \frac{\rho_{(1)}(\lambda_j)}{\bigl(\lambda_j-(\lambda+i\varepsilon)\bigr)}, \end{align} where we supposed that the variables $\lambda_1,...,\lambda_N$ are indistinguishable (i.e. same marginal probability density $\forall j$).
We now consider the complex path integral over $C$ -- a closed semi-circle that goes to infinity in the upper half of the complex plane. The marginal pdf is normalizable by definition and there must therefore exist a region at $\infty$ where the integral over the open semi-circle $C_R$ is null, leaving us with an integral over the real line only. Mathematically: \begin{align} &\oint_C d z \frac{\rho_{(1)}(z)}{\bigl(z-(\lambda+i\varepsilon)\bigr)} = \\& \int_{-\infty}^\infty d\lambda_j\ \frac{\rho_{(1)}(\lambda_j)}{\bigl(\lambda_j-(\lambda+i\varepsilon)\bigr)} + \int_{C_R} d z \frac{\rho_{(1)}(z)}{\bigl(z-(\lambda+i\varepsilon)\bigr)} =\int_{-\infty}^\infty d\lambda_j\ \frac{\rho_{(1)}(\lambda_j)}{\bigl(\lambda_j-(\lambda+i\varepsilon)\bigr)} \end{align} But, this integral is also proportionnal to the sum of the residue in $C$. Here, the pdf cannot feature poles--yet again, by definition--which leaves us with a single pole of order 1 at $z=\lambda+i\varepsilon$, with corresponding residue \begin{align} \text{Res}\left(\frac{\rho_{(1)}(z)}{z-\lambda-i\varepsilon},\lambda+i\varepsilon\right)&=\lim_{z\to\lambda+i\varepsilon } \bigl(z-(\lambda+i\varepsilon)\bigr)\frac{\rho_{(1)}(z)}{z-(\lambda+i\varepsilon)}\notag\\ &=\rho_{(1)}(\lambda+i\varepsilon) \end{align} i.e. \begin{align} \oint_C d z \frac{\rho_{(1)}(z)}{\bigl(z-(\lambda+i\varepsilon)\bigr)} = 2\pi i\text{Res}\left(\frac{\rho_{(1)}(z)}{z-\lambda-i\varepsilon},\lambda+i\varepsilon\right) \end{align}
Putting all these results together yields \begin{align} \langle \text{Tr}(R_{\mathbf{X}}(\lambda+i\varepsilon))\rangle =2\pi i \rho_{(1)}(\lambda+i\varepsilon) \sum_{j,k=1}^Nv^{(j)}_kv^{(j)}_k = 2N\pi i \rho_{(1)}(\lambda+i\varepsilon) \end{align} Or in the limit $\varepsilon\to0$ \begin{align} \lim_{\varepsilon\to0}\langle \text{Tr}(R_{\mathbf{X}}(\lambda+i\varepsilon))\rangle = 2N\pi i \rho_{(1)}(\lambda) \end{align}
Taking the imaginary part finally leads to
\begin{align} \rho_{(1)}(\lambda) = \frac{1}{2N\pi} \lim_{\varepsilon\to0}\ \mathrm{Im} \bigg\langle \text{Tr} \big(R_{\mathbf{X}}(\lambda+i\varepsilon)\big)\bigg\rangle \end{align} which is the desired result divided by an extra 2N factor!
