Consider GUE$(N)$, the set of $N \times N$ random Hermitian matrices where the elements on the diagonal are i.i.d. according to $\mathcal{N}(0,1)$ and both the real and imaginary parts of the upper-triangular elements are i.i.d. according to $\mathcal{N}(0,1/2)$. I am interested in the (suitably normalized) probability distribution on the spacing between the two largest (or smallest) eigenvalues as $N \rightarrow \infty$. Specifically I'd like to know $$F_\Delta(s) = \lim_{N \rightarrow \infty} \text{Prob} \left( N^{1/6} \left( \lambda_N - \lambda_{N-1} \right) \leq s \right) \,,$$ where $\lambda_N$ is the largest eigenvalue of a GUE$(N)$ matrix and $\lambda_{N-1}$ is the next-to-largest.
What I've gotten from the literature so far: The distribution of $\lambda_N$ is well-known to follow the Tracy-Widom law: $$\lim_{N \rightarrow \infty} \text{Prob} \left( N^{2/3} \left( \frac{\lambda_N}{\sqrt{N}} - 2 \right) \leq s \right) = F_2(s)$$ where $F_2$ is the Tracy-Widom distribution. From this formula we can expect that the typical spacing between the two largest eigenvalues in the GUE$(N)$ will be of order $N^{-1/6}$, hence the inverse factor in the first formula above. This post by T. Tao discussing this paper focusses on the spacing between two consecutive eigenvalues in the "bulk" of the spectrum, which is not what we're interested in here. In that post he mentions "In the edge case [...] the distribution is given by the famous Tracy-Widom law." However I have not been able to deduce the result. Also in this paper, at the bottom of p.6, it is stated that the joint density for $(\lambda_N,\lambda_{N-1})$ follows a "multivariate Tracy-Widom" distribution (from which we could deduce the distribution on the difference $\lambda_N-\lambda_{N-1}$), referring to this book. However the only remark about this in that book I found was in Remark 3.7.3 on p.134 where one rather vaguely states that some joint distribution for the suitably normalized pair $(\lambda_N,\lambda_{N-1})$ exists.
