It's well known in random matrix theory that locally the eigenvalues of a random matrix from the Gaussian unitary ensemble tend to a sine-kernel determinantal point process. Likewise, locally the eigenangles of a random unitary matrix (the 'circular unitary ensemble') also tend to a sine-kernel determinanantal point process. (An accessible introductory account of these facts, with definitions, is given in section 1.3 of this survey of Bourgade and Keating.)
One usually proves these facts by explicitly computing the joint intensity of the eigenvalues of an $N\times N$ GUE matrix -- extracting in this way the square of a Vandermonde determinant term -- then using orthogonal polynomials and Gaudin's lemma to find the correlation functions of GUE, and finally using something like Plancherel-Rotach asymptotics to evaluate the correlation function in the limit. A similar process holds for unitary matrices, with the first step just being the Weyl integration formula. This all is outlined in Bourgade and Keating's note.
My question is whether there is a less computational approach to see that the eigenvalues of a random $GUE(N)$ matrix and the eigenangles of a random $U(N)$ matrix tend locally to the same distribution (bypassing that this distribution is the sine-kernel process). To make the question more precise, is there for instance a simple demonstration of the fact that for nice test functions $\eta$, and $E\in(-2,2)$,
$$ \lim_{N\rightarrow\infty}\mathbb{E}_{U(N)} \sum_{i\neq j} \eta(N\theta_i,N\theta_j) = \lim_{N\rightarrow\infty} \mathbb{E}_{GUE(N)} \sum_{i\neq j}\eta\big( N\rho_{sc}(E)(\lambda_i-E), N\rho_{sc}(E)(\lambda_j-E) \big) $$
without actually computing what either limit is? (Or indeed, even doing so much computation as to extract a Vandermonde squared term.)
Here $\theta_i\in[-1/2,1/2)$ are the eigenangles of a random unitary matrix, and $\lambda_i$ are the eigenvalues of a random GUE matrix, normalized as in Bourgade and Keating's survey (this is the normalization such that nearly all eigenvalues lie in the interval $(-2,2)$). For the purposes of this question, define the weight $\rho_{sc}$ by the relation $$ \lim_{N\rightarrow\infty}\mathbb{E}_{GUE(N)} \frac{1}{N}\sum_i f( \lambda_i ) = \int f(x) \rho_{sc}(x)\,dx. $$
I'm aware the question is very much a naive one, but I have never gotten an answer that was entirely satisfying to me. I know enough representation theory to recite here the fact that the space of positive definite Hermitian matrices is in correspondence with $GL(N,\mathbb{C})/U(N)$, which is dual to $U(N)$, with a vague sense acquired from multiple sources that this is the relevant fact. I don't understand the theory of symmetric spaces well enough however to see why this implies the formula highlighted above. In general the less algebraic sophistication required to understand an answer the easier it will be for me (and perhaps others) to assimilate.
