\chapter{Introduction}\label{s:intro}

\section{Background}\label{s:background}

The gaps between consecutive eigenvalues are among the central objects of random matrix theory. In the bulk of the spectrum, the limiting distribution of an individual gap is known for a wide range of models \cite{ErdosYau2015,Zhang2025}. The smallest and largest gaps within a fixed interval in the bulk are understood for far fewer models. Unlike an individual gap, they are not determined by the local eigenvalue statistics on any fixed scale. The motivation for their study is discussed in the introduction of \cite{BenArousBourgade2013}, which includes a comparison with numerical data for the Riemann zeta function.

We first describe what is known about the largest gaps. For the Gaussian unitary ensemble (GUE) and the circular unitary ensemble (CUE) \cite{Dyson1962,Mehta2004}, Ben Arous and Bourgade found the first-order asymptotics of the largest gaps in the bulk, which are of order $\sqrt{\log N}/N$, and left their fluctuations open \cite{BenArousBourgade2013}. Feng and Wei identified these fluctuations \cite{FengWei2025}. They proved that the largest gaps of the CUE and the GUE, after centering and rescaling, converge to a Poisson point process with an exponential intensity, so that the largest gap has a Gumbel limit. Charlier recently extended this result to unitary-invariant ensembles for a broad class of regular analytic potentials \cite{Charlier2026}. Bourgade proved universality of the smallest gaps in both symmetry classes, and of the largest gaps in the complex Hermitian class, for generalized Wigner matrices satisfying a smoothness assumption \cite{Bourgade2022}. His proof combines new quantitative relaxation estimates for Dyson Brownian motion with reverse heat flow. Landon, Lopatto, and Marcinek proved comparison theorems for the largest gaps of real symmetric and complex Hermitian generalized Wigner matrices without a smoothness assumption \cite{LandonLopattoMarcinek2020}. In the complex Hermitian case, both works use the GUE result of Feng and Wei \cite{FengWei2025} to identify the limiting distribution. Both approaches also apply to the real symmetric largest-gap problem, but identifying its limiting distribution required the corresponding GOE result, which was not known.

The smallest gaps were studied earlier. Vinson's dissertation \cite{Vinson2001} was the first rigorous work on extremal gaps. It determined the limiting distribution of the smallest gap of the CUE and the GUE. Ben Arous and Bourgade found the joint limit of the smallest gaps of these two ensembles \cite{BenArousBourgade2013}, building on Soshnikov's work on determinantal point processes \cite{Soshnikov2005}. The limiting distribution of the smallest gaps is now also known for the Gaussian orthogonal ensemble (GOE) \cite{FengTianWei2019}, for the circular $\beta$-ensemble \cite{Forrester2010} with $\beta$ a positive integer \cite{FengWei2021}, for the Gaussian symplectic ensemble \cite{FengLiYao2024}, and for complex Wishart ensembles and one-cut unitary-invariant ensembles with real-analytic confining potentials \cite{ShiJiang2012}. Bourgade's universality theorem \cite{Bourgade2022} covers the smallest gaps of real symmetric and complex Hermitian generalized Wigner matrices with smooth entries, and Figalli and Guionnet proved limiting laws for the smallest gaps and first-order asymptotics for the largest gaps in perturbative several-matrix models with $\beta = 2$ \cite{FigalliGuionnet2016}. For Wigner matrices and related models, lower bounds for the smallest gaps were proved in \cite{NguyenTaoVu2017}, and for sparse random matrices in \cite{LopattoLuh2021}. Extremal gaps have also been studied for non-Hermitian matrices \cite{LopattoMeeker2026,LopattoOtto2025,Moreillon2026} and for the zeros of random functions \cite{FengGotzeYao2024,FengYao2026,FengMuirhead2026}.

The results above that identify limiting laws or first-order asymptotics for the largest eigenvalue gaps rest on a determinantal structure \cite{Soshnikov2000}, either of the ensemble itself or of the Gaussian ensemble with which it is compared. The GOE does not have this structure, and the limiting distribution of its largest gaps has remained open; the problem was raised in \cite{LandonLopattoMarcinek2020,Bourgade2022,LopattoMeeker2026}. For the circular orthogonal and symplectic ensembles (COE and CSE) \cite{Dyson1962,Mehta2004}, Feng and Wei noted that even the order of the largest gaps is unknown \cite{FengWei2025}. Our main result, \zcref{t:main}, determines the limiting distribution of the largest gaps in the bulk of the spectrum of the GOE.

The analysis of the largest gaps rests on precise asymptotics for the probability that a long interval contains no eigenvalue. For the sine process, which is the limit of the GUE eigenvalues after rescaling in the bulk \cite{AndersonGuionnetZeitouni2010}, this probability is a Fredholm determinant of the sine kernel, and its asymptotics were predicted by des Cloizeaux and Mehta \cite{DesCloizeauxMehta1973} and Dyson \cite{Dyson1976}. The leading term was proved by Widom \cite{Widom1994}, and the constant term, known as the Widom--Dyson constant, by Krasovsky \cite{Krasovsky2004}, Ehrhardt \cite{Ehrhardt2006}, and Deift, Its, Krasovsky, and Zhou \cite{DeiftItsKrasovskyZhou2007}. For $\beta = 1$ and $\beta = 4$, the corresponding asymptotics, including the constant terms, were proved by Ehrhardt \cite{Ehrhardt2007}. For general $\beta > 0$, the bulk scaling limit of the Gaussian $\beta$-ensemble is the $\operatorname{Sine}_\beta$ process \cite{ValkoVirag2009}, and the probability that it has no point in a long interval was determined up to a constant factor by Valk\'o and Vir\'ag \cite{ValkoVirag2010}.

The eigenvalues of the GOE form a Pfaffian point process rather than a determinantal one: the $k$-point correlation function is the Pfaffian of a $2k \times 2k$ matrix built from a $2 \times 2$ matrix kernel, and one entry of this kernel contains the sign function $\sgn(x - y)/2$ \cite{Mehta2004}. For even matrix dimensions, Tracy and Widom \cite{TracyWidom1996} expressed the square of the probability that the GOE has no eigenvalue in a union of intervals as the corresponding probability for the GUE, multiplied by the determinant of a matrix of fixed size; see also \cite{TracyWidom1998,Widom1999,AdlerForresterNagaoVanMoerbeke2000,ForresterMays2009,Sinclair2009,Forrester2010} for the correlation kernels of orthogonal and symplectic ensembles, and \cite{DeiftGioev2007a,DeiftGioev2007b,DeiftGioev2009} for their use in proving universality. Our proof builds on the identity of Tracy and Widom. However, we derive it differently, for even and odd matrix dimensions, and we analyze it on intervals whose length in units of the local mean spacing grows with $N$; see \zcref{s:overview}.

\section{Main results}\label{s:mainresults}

Let $H_N$ be an $N \times N$ matrix from the GOE, normalized so that its distribution has density proportional to $\exp(- N \Tr H^2 / 4)$ with respect to Lebesgue measure on real symmetric matrices. Equivalently, the entries $H_{ij}$ with $i \le j$ are independent centered Gaussian random variables, with variance $2/N$ if $i = j$ and $1/N$ if $i < j$. We let $\lambda_1 < \lambda_2 < \dots < \lambda_N$ denote the eigenvalues of $H_N$, and we let $Z_N$ denote the normalizing constant in their joint density, which is
\begin{equation}\label{e:goedensity}
\frac{1}{Z_N} \prod_{1 \le i < j \le N} |\lambda_j - \lambda_i| \prod_{i = 1}^N e^{- N \lambda_i^2 / 4}
\end{equation}
on $\{ \lambda_1 < \dots < \lambda_N \}$; see \cite{AndersonGuionnetZeitouni2010}. We define
\begin{equation}\label{e:sigmadef}
\sigma(y) = \sqrt{4 - y^2} , \qquad |y| \le 2 .
\end{equation}
With this normalization, the empirical eigenvalue distribution converges to the semicircle law. This is the probability measure on $[-2, 2]$ with density $\sigma(y)/(2\pi)$. The mean spacing between consecutive eigenvalues near $y \in (-2, 2)$ is approximately $2\pi/(N \sigma(y))$.

We define the constant $\kappa_1$ and the function $F$ by
\begin{equation}\label{e:kappa1def}
\kappa_1 = 2^{13/24} \exp \Big( \frac32 \zeta'(-1) \Big) , \qquad F(s) = \frac{s^2}{64} + \frac s8 , \quad s \ge 0 ,
\end{equation}
where $\zeta$ denotes the Riemann zeta function. Both come from the asymptotics of the probability that the $\operatorname{Sine}_1$ process has no point in a long interval, which we recall after \zcref{c:ordered}. For $N \ge 2$, we let $b_N > 0$ denote the solution of $F(b_N) = \log N$, and we define $a_N = F'(b_N)$. Solving the quadratic equation $F(b_N) = \log N$, we find
\begin{equation}\label{e:abdef}
b_N = 8 \sqrt{\log N + \frac14} - 4 , \qquad a_N = \frac{b_N + 4}{32} , \qquad a_N b_N = 2 \log N - \frac{b_N}{8} .
\end{equation}

We fix an interval $I = [y_-, y_+]$ with $-2 < y_- < y_+ < 2$, and we consider the gaps $\lambda_{i + 1} - \lambda_i$ with both endpoints in $I$. We define
\begin{equation}\label{e:sigmaIdef}
\sigma_I = \min_{y \in I} \sigma(y) , \qquad \partial_* I = \{ \xi \in \{ y_-, y_+ \} : \sigma(\xi) = \sigma_I \} .
\end{equation}
Observe that $\partial_* I$ is the set of endpoints of $I$ at which the semicircle density is smallest. It consists of the endpoint with the larger absolute value if $y_- \neq - y_+$, and of the two endpoints if $y_- = - y_+$. In particular, $\xi \neq 0$ for all $\xi \in \partial_* I$. We define
\begin{equation}\label{e:LambdaIdef}
\Lambda_I = \kappa_1 \sum_{\xi \in \partial_* I} \frac{\sigma_I^3}{|\xi|} .
\end{equation}
We let $\mathcal G_N(I)$ denote the set of indices $i \in \{1, \dots, N - 1\}$ such that $\lambda_i, \lambda_{i + 1} \in I$, and for $i \in \mathcal G_N(I)$ we define the \emph{rescaled gap}
\begin{equation}\label{e:taudef}
\tau_{i, N} = a_N \big( N \sigma_I (\lambda_{i + 1} - \lambda_i) - b_N \big) + \frac98 \log b_N .
\end{equation}
Finally, we let $\mathcal X_N$ denote the point process on $\R$ whose points are the numbers $\tau_{i, N}$, $i \in \mathcal G_N(I)$. For a Borel set $A \subset \R$, we let $\mathcal X_N(A)$ denote the number of points of $\mathcal X_N$ in $A$, which is the number of indices $i \in \mathcal G_N(I)$ such that $\tau_{i, N} \in A$. Our main result is the following theorem.

\begin{theorem}\label{t:main}
Let $m \ge 1$, and let $A_1, \dots, A_m \subset \R$ be disjoint intervals, each of the form $(a, b]$ with $a < b$ or $(a, \infty)$. Then, as $N \to \infty$, the vector $(\mathcal X_N(A_1), \dots, \mathcal X_N(A_m))$ converges in distribution to a vector of independent Poisson random variables with means $\Lambda_I \int_{A_l} e^{-x} \dd x$, $1 \le l \le m$. In particular, $\mathcal X_N$ converges in distribution to the Poisson point process on $\R$ with intensity $\Lambda_I e^{-x} \dd x$, with respect to the vague topology.
\end{theorem}

We let $\tau^\downarrow_{1, N} \ge \tau^\downarrow_{2, N} \ge \cdots$ denote the values $\tau_{i, N}$, $i \in \mathcal G_N(I)$, arranged in decreasing order, and we set $\tau^\downarrow_{k, N} = - \infty$ if $\mathcal G_N(I)$ has fewer than $k$ elements. We also let $m_{1, N}$ denote the largest gap $\lambda_{i + 1} - \lambda_i$ with $i \in \mathcal G_N(I)$ (and we set $m_{1, N} = 0$ if $\mathcal G_N(I)$ is empty). For random variables $Y_N$ and positive numbers $\theta_N$, we write $Y_N = O_\Prob(\theta_N)$ if the sequence $(Y_N / \theta_N)_N$ is tight.

\begin{corollary}\label{c:ordered}
For every $k \ge 1$ and $x \in \R$,
\begin{equation}\label{e:ordered}
\lim_{N \to \infty} \Prob\big( \tau^\downarrow_{k, N} \le x \big) = \exp\big( - \Lambda_I e^{-x} \big) \sum_{j = 0}^{k - 1} \frac{(\Lambda_I e^{-x})^j}{j!} .
\end{equation}
In particular, $\tau^\downarrow_{1, N} - \log \Lambda_I$ converges in distribution to the standard Gumbel distribution, which has the distribution function $x \mapsto \exp(-e^{-x})$:
\begin{equation}\label{e:gumbel}
\lim_{N \to \infty} \Prob\big( \tau^\downarrow_{1, N} - \log \Lambda_I \le x \big) = \exp\big( -e^{-x} \big) , \qquad x \in \R .
\end{equation}
Moreover,
\begin{equation}\label{e:firstorder}
N \sigma_I m_{1, N} = 8 \sqrt{\log N} - 4 - \frac94 \frac{\log \log N}{\sqrt{\log N}} + O_\Prob\big( (\log N)^{-1/2} \big) .
\end{equation}
\end{corollary}

We now explain where the constant $\kappa_1$ and the function $F$ come from. For a Borel set $A \subset \R$, we define the \emph{hole probability}
\begin{equation}\label{e:E1def}
E_{1, N}(A) = \Prob\big( \lambda_i \notin A \text{ for all } i \in \{1, \dots, N\} \big) .
\end{equation}
The limit of the hole probability of a short interval is expressed through the sine kernel. We define the \emph{sine kernel} $S$ by
\begin{equation}\label{e:sineintro}
S(u) = \frac{\sin u}{\pi u} , \quad u \neq 0 , \qquad S(0) = \frac1\pi .
\end{equation}
It is the kernel of the sine process, rescaled to have mean spacing $\pi$. For $t > 0$, we let $S_{+, t}$ and $S_{-, t}$ denote the integral operators on $L^2(0, t)$ with kernels $S(u - v) + S(u + v)$ and $S(u - v) - S(u + v)$, respectively, and we define the Fredholm determinants
\begin{equation}\label{e:Dpmintro}
D_+(t) = \det(\Id - S_{+, t}) , \qquad D_-(t) = \det(\Id - S_{-, t}) , \qquad q_1(s) = D_+(s/4) ,
\end{equation}
for $t, s > 0$. We will show in \zcref{c:onehole} that $q_1(s)$ is the limit of $E_{1, N}((y, y + s/(N \sigma(y))))$ as $N \to \infty$, for all $y \in (-2, 2)$ and $s > 0$. The interval $(y, y + s/(N \sigma(y)))$ has length $s/(2\pi)$ in units of the mean spacing near $y$. So $q_1(s)$ is the probability that the $\operatorname{Sine}_1$ process, normalized to have mean spacing $2\pi$, has no point in a given interval of length $s$; this expression for the gap probability is classical \cite{Mehta2004,TracyWidom1996}. The asymptotics of $D_\pm$ were determined by Ehrhardt \cite[Theorem~4.8]{Ehrhardt2007}, confirming predictions of Dyson \cite{Dyson1976}. Ehrhardt proved that
\begin{equation}\label{e:ehrhardtintro}
\log D_\pm(t) = - \frac{t^2}{4} \mp \frac t2 - \frac{\log t}{8} + \Big( \frac{1}{24} \pm \frac14 \Big) \log 2 + \frac32 \zeta'(-1) + o(1) , \qquad t \to \infty .
\end{equation}
Setting $t = s/4$ in \eqref{e:ehrhardtintro} (see \zcref{l:q1}), we have
\begin{equation}\label{e:q1asymp}
q_1(s) = \kappa_1 s^{-1/8} e^{-F(s)} \big( 1 + o(1) \big) , \qquad \frac{D_+(s/4)}{D_-(s/4)} = \sqrt 2 \, e^{-s/4} \big( 1 + o(1) \big) , \qquad s \to \infty ,
\end{equation}
with $\kappa_1$ and $F$ as in \eqref{e:kappa1def}.

\begin{remark}\label{r:constants}
We comment on the constants in \zcref{t:main} and \zcref{c:ordered}.
\begin{enumerate}[label=(\alph*)]
\item By \eqref{e:q1asymp}, the function $F$ is the exponent in the probability that the $\operatorname{Sine}_1$ process has no point in a long interval. Its quadratic term, its linear term, and the power $s^{-1/8}$ in \eqref{e:q1asymp} agree with the asymptotics of Valk\'o and Vir\'ag for the $\operatorname{Sine}_\beta$ process at $\beta = 1$ \cite{ValkoVirag2010}, which are stated for intervals of length $s$ in units in which the mean spacing is $2\pi$. The constant $\kappa_1$ is the constant of Ehrhardt \cite{Ehrhardt2007}, after the change of variables $t = s/4$.
\item The constant $-4$ in \eqref{e:firstorder} comes from the linear term $s/8$ in $F$. The coefficient $9/4$ of the third term comes from the power $s^{-1/8}$ in \eqref{e:q1asymp} and from the fact that the largest gaps are located within distance of order $(\log N)^{-1}$ of the points of $\partial_* I$; see \zcref{r:heuristic}.
\item We consider the GUE with eigenvalue density proportional to $\prod_{i < j} |\lambda_j - \lambda_i|^2 \prod_i e^{-N \lambda_i^2/2}$. Its empirical eigenvalue distribution also converges to the semicircle law on $[-2, 2]$. We let $m^{\mathrm{GUE}}_{1, N}$ denote the largest gap $\lambda_{i + 1} - \lambda_i$ of the GUE with $\lambda_i, \lambda_{i + 1} \in I$. It follows from \cite[Theorem~2]{FengWei2025} that
\begin{equation}
N \sigma_I m^{\mathrm{GUE}}_{1, N} = \sqrt{32 \log N} - \frac{5}{2 \sqrt 2} \frac{\log \log N}{\sqrt{\log N}} + O_\Prob\big( (\log N)^{-1/2} \big) .
\end{equation}
The leading term for the GOE is larger by a factor $\sqrt 2$. It agrees with the prediction $\sqrt{(64/\beta) \log N}$ at $\beta = 1$, which Feng and Wei formulated for the circular $\beta$-ensemble \cite{FengWei2025}.
\end{enumerate}
\end{remark}

The following remark records a consequence of \zcref{c:ordered} for real symmetric generalized Wigner matrices, which was anticipated in \cite{LandonLopattoMarcinek2020} and \cite{Bourgade2022}. We prove it in \zcref{s:factorial}.

\begin{remark}\label{r:universality}
Let $W_N$ be an $N \times N$ real symmetric generalized Wigner matrix satisfying \cite[Definition~1.1]{LandonLopattoMarcinek2020}; in particular, the entries of $\sqrt N W_N$ have finite moments of all orders. For $k \ge 1$, let $\widehat T_{k, N}$ denote the $k$-th largest of the gaps $\mu_{i + 1} - \mu_i$ between consecutive eigenvalues $\mu_1 \le \dots \le \mu_N$ of $W_N$ such that $\mu_i \in I$ (this is the convention of \cite{LandonLopattoMarcinek2020}), and let $\widehat T_{k, N} = 0$ if there are fewer than $k$ such gaps. We define
\begin{equation}
\widehat \tau_{k, N} = a_N \big( N \sigma_I \widehat T_{k, N} - b_N \big) + \frac98 \log b_N .
\end{equation}
Then \eqref{e:ordered} holds with $\widehat \tau_{k, N}$ in place of $\tau^\downarrow_{k, N}$. This follows from \zcref{c:ordered} and the comparison theorem \cite[Corollary~1.6]{LandonLopattoMarcinek2020}.
\end{remark}

\section{Notation}\label{s:notation}

We use the following notation throughout the paper. For a statement $P$, we let $\1_P$ denote $1$ if $P$ holds and $0$ otherwise, and for a set $A \subset \R$, we let $\1_A$ denote its indicator function. For a Borel set $A \subset \R$, we let $|A|$ denote its Lebesgue measure, and for $c > 0$, we write $c A = \{ c a : a \in A \}$. We write $\| \cdot \|_2$ for the norm of $L^2(A)$, and $\| T \|$ for the operator norm of a bounded operator $T$ on $L^2(A)$. We identify $L^2(A)$ with the space of functions in $L^2(\R)$ that vanish outside $A$. For $t \in \R$, we write $(t)_+ = \max(t, 0)$. For $t, y \in \R$, we let $\int_t^y$ denote the oriented integral, so that $\int_t^y = - \int_y^t$ if $y < t$. In \zcref{s:overview}, $C$ and $c$ denote positive constants that do not depend on $N$, and whose values may change from line to line. They may depend on the fixed parameters $\delta$, $\gamma$, $r_0$, $k$, and $x_1, \dots, x_k$ introduced there, and the implicit constants in $O(\cdot)$ have the same property. Limits are taken as $N \to \infty$.

\section{Overview of the proof}\label{s:overview}

We now describe the main steps of the proof. This overview can be read independently of the later sections; all objects in it are defined here, in \zcref{s:mainresults}, or in \zcref{s:notation}. The proof has two parts. The first part, which follows the strategy of Feng and Wei for the GUE \cite{FengWei2025}, reduces \zcref{t:main} to three estimates for hole probabilities. The second part proves these estimates in four steps, and Steps~1--3 are new compared with the proof for the GUE. Where the GUE proof uses that the hole probability is a Fredholm determinant, Step~1 uses an exact identity, valid for even and odd $N$, that writes the square of the GOE hole probability as the GUE hole probability times a determinant of fixed size. Step~2 replaces the sign function in this determinant, which does not decay, by functions that are concentrated near the endpoints of the intervals. Step~3 replaces the elementary bound on the inverse of the operator in the Fredholm determinant, which is a power of $N$ in our setting, by a bound that is exponential in the expected number of eigenvalues in the intervals. Step~4 uses these results and the asymptotics of Hermite functions to compare the GOE hole probability with its limit.

\paragraph{Reduction to hole probabilities.}
For the Poisson point process with intensity $\Lambda_I e^{-x} \dd x$, the sum of $\prod_{j = 1}^k (\tau_j - x_j)_+$ over all ordered $k$-tuples $(\tau_1, \dots, \tau_k)$ of distinct points has expectation
\begin{equation}
\prod_{j = 1}^k \int_\R (t - x_j)_+ \Lambda_I e^{-t} \dd t = \prod_{j = 1}^k \Lambda_I e^{-x_j} .
\end{equation}
By a criterion of the type used in \cite[Lemma~1]{FengWei2025}, which we state in \zcref{l:poisson}, it therefore suffices to show that for every $k \ge 1$ and all $x_1, \dots, x_k \in \R$,
\begin{equation}\label{e:overviewmoments}
\lim_{N \to \infty} \E \sum_{\substack{i_1, \dots, i_k \in \mathcal G_N(I) \\ \text{distinct}}} \prod_{j = 1}^k \big( \tau_{i_j, N} - x_j \big)_+ = \prod_{j = 1}^k \Lambda_I e^{-x_j} .
\end{equation}
The positive parts $(\tau - x)_+$ are used instead of the indicators $\1_{\tau > x}$ because the left side of \eqref{e:overviewmoments} can be expressed through hole probabilities, as follows. For $x \in \R$, we define
\begin{equation}\label{e:overviewell}
\ell_N(x) = \frac{1}{N \sigma_I} \Big( b_N + \frac{x - \frac98 \log b_N}{a_N} \Big) , \qquad \chi_N = a_N N \sigma_I .
\end{equation}
By \eqref{e:taudef}, we have $\tau_{i, N} - x = \chi_N (\lambda_{i + 1} - \lambda_i - \ell_N(x))$ for all $i \in \mathcal G_N(I)$. Further, by Taylor's theorem and the definitions of $a_N$ and $b_N$, we have
\begin{equation}\label{e:overviewFell}
F\big( N \sigma_I \ell_N(x) \big) = \log N + x - \frac98 \log b_N + o(1) .
\end{equation}
We also call the interval $(\lambda_i, \lambda_{i + 1})$ a gap. For a gap $(\lambda_i, \lambda_{i + 1})$ and $\ell > 0$, the set of $y \in \R$ such that $[y, y + \ell] \subset (\lambda_i, \lambda_{i + 1})$ is an interval of length $(\lambda_{i + 1} - \lambda_i - \ell)_+$. For $y_1, \dots, y_k \in \R$, we let $\mathcal W_N(y_1, \dots, y_k)$ denote the event that the intervals $[y_j, y_j + \ell_N(x_j)]$, $1 \le j \le k$, are contained in distinct gaps $(\lambda_i, \lambda_{i + 1})$ with $i \in \mathcal G_N(I)$. Integrating over the left endpoints $y_1, \dots, y_k$, we show in \zcref{l:volume} that, for $N$ large enough that $\ell_N(x_j) > 0$ for all $j$, the expectation in \eqref{e:overviewmoments} equals
\begin{equation}\label{e:overviewvolume}
\chi_N^k \int_{y_-}^{y_+ - \ell_N(x_1)} \cdots \int_{y_-}^{y_+ - \ell_N(x_k)} \Prob\big( \mathcal W_N(y_1, \dots, y_k) \big) \dd y_k \cdots \dd y_1 .
\end{equation}

If two of the intervals $[y_j, y_j + \ell_N(x_j)]$ intersect, then they cannot lie in distinct gaps, so $\Prob(\mathcal W_N(y_1, \dots, y_k)) = 0$. Otherwise, the probability $\Prob(\mathcal W_N(y_1, \dots, y_k))$ is bounded by the hole probability of the union of the intervals $(y_j, y_j + \ell_N(x_j))$. Conversely, suppose that these intervals are separated from each other, and from the endpoints of $I$, by distances at least $3 N^{-1/3}$. Then $\Prob(\mathcal W_N(y_1, \dots, y_k))$ is at least this hole probability minus the probability that some interval of length $N^{-1/3}$ in $I$ contains no eigenvalue. The latter probability is $O(N^{-D})$ for every fixed $D > 0$. Indeed, every interval of length $N^{-1/3}$ in $I$ contains one of $O(N^{1/3})$ fixed intervals of the form $(y, y + s/(N \sigma(y)))$ with $s = \gamma_0 \sqrt{\log N}$, where $\gamma_0 > 0$ is a constant with $\gamma_0^2 / 64 > D + 1/3$. By \eqref{e:overviewonehole} below, applied with $\gamma = \gamma_0$, the hole probability of each of these intervals is at most $2 q_1(\gamma_0 \sqrt{\log N}) \le N^{-\gamma_0^2/64}$, and a union bound gives the claim.

It remains to estimate the hole probabilities of unions of at most $k$ intervals in $I$ with total length at most $C \sqrt{\log N}/N$. We need three estimates. The first is an asymptotic formula for one interval, which we prove in \zcref{c:onehole}. For all $\delta \in (0, 1)$ and $\gamma > 0$, there exists $C > 0$ such that
\begin{equation}\label{e:overviewonehole}
E_{1, N} \Big( \Big( y, y + \frac{s}{N \sigma(y)} \Big) \Big) = q_1(s) \big( 1 + O(e^{C s} N^{-1/2}) \big)
\end{equation}
uniformly in $y \in [-2 + \delta, 2 - \delta]$ and $0 < s \le \gamma \sqrt{\log N}$. The second and third estimates concern open intervals $A_1, \dots, A_r \subset [-2 + \delta, 2 - \delta]$ with disjoint closures and total length $\sum_j |A_j| \le \gamma \sqrt{\log N}/N$, where the number $r$ of intervals is at most a fixed integer $r_0$. We prove in \zcref{c:several} that
\begin{equation}\label{e:overviewquasi}
E_{1, N}\Big( \bigcup_{j = 1}^r A_j \Big) \le e^{C \sqrt{\log N}} \prod_{j = 1}^r E_{1, N}(A_j) ,
\end{equation}
and that if, in addition, $\dist(A_i, A_j) \ge N^{-1/2}$ for all $i \neq j$, then
\begin{equation}\label{e:overviewfactor}
E_{1, N}\Big( \bigcup_{j = 1}^r A_j \Big) = \big( 1 + O(e^{C \sqrt{\log N}} N^{-1/4}) \big) \prod_{j = 1}^r E_{1, N}(A_j) .
\end{equation}

We now explain how these estimates give \eqref{e:overviewmoments}. We fix $\delta > 0$ such that $I \subset [-2 + \delta, 2 - \delta]$, and we take $r_0 = k$. Since $N \ell_N(x) = (8 + o(1)) \sqrt{\log N} / \sigma_I$ for every fixed $x$, we may also fix $\gamma > 0$ such that $N \sum_j \ell_N(x_j) \le \gamma \sqrt{\log N}$ and $2 N \ell_N(x_j) \le \gamma \sqrt{\log N}$ for all $j$ and all large $N$. We first consider the points $(y_1, \dots, y_k)$ in the domain of integration in \eqref{e:overviewvolume} for which the intervals $[y_j, y_j + \ell_N(x_j)]$ are separated from each other, and from the endpoints of $I$, by distances at least $3 N^{-1/3}$. For these points, we replace the probability $\Prob(\mathcal W_N(y_1, \dots, y_k))$ by the hole probability of $\bigcup_j (y_j, y_j + \ell_N(x_j))$, with an error $O(N^{-D})$, where we choose $D > k + 1$ so that the error remains negligible after multiplication by $\chi_N^k \le (N \log N)^k$. Next, we use \eqref{e:overviewfactor} with $A_j = (y_j, y_j + \ell_N(x_j))$ and $r = k$ to replace this hole probability by the product of the hole probabilities $E_{1, N}((y_j, y_j + \ell_N(x_j)))$. For the remaining points at which the intervals $[y_j, y_j + \ell_N(x_j)]$ are disjoint, we use \eqref{e:overviewquasi} to bound $\Prob(\mathcal W_N(y_1, \dots, y_k))$ by $e^{C \sqrt{\log N}}$ times the same product. By \eqref{e:overviewonehole}, \eqref{e:q1asymp}, and \eqref{e:overviewFell}, and since $q_1$ is nonincreasing and $\sigma(y) \ge \sigma_I$ on $I$, each factor $\chi_N E_{1, N}((y_j, y_j + \ell_N(x_j)))$ is $O(\log N)$. Further, the set of these remaining points has $k$-dimensional Lebesgue measure $O(N^{-1/3})$. Their contribution is therefore $O(e^{C \sqrt{\log N}} N^{-1/3} (\log N)^k) = o(1)$, and \eqref{e:overviewmoments} reduces to the statement
\begin{equation}\label{e:overviewfirst}
\lim_{N \to \infty} \chi_N \int_{y_-}^{y_+ - \ell_N(x)} E_{1, N}\big( (y, y + \ell_N(x)) \big) \dd y = \Lambda_I e^{-x} , \qquad x \in \R ,
\end{equation}
which we prove in \zcref{s:firstmoment}. To compute this integral, we apply \eqref{e:overviewonehole} with $s = N \sigma(y) \ell_N(x)$. Note that
\begin{equation}
\begin{gathered}
N \sigma(y) \ell_N(x) - N \sigma_I \ell_N(x) = N \sigma_I \ell_N(x) \Big( \frac{\sigma(y)}{\sigma_I} - 1 \Big) , \\
F'(b_N) b_N = a_N b_N = 2 \log N \big( 1 + o(1) \big) ,
\end{gathered}
\end{equation}
by \eqref{e:abdef}. Set $s_N = N \sigma_I \ell_N(x)$ and $\Psi(y) = \sigma(y) / \sigma_I$. Using \eqref{e:overviewonehole} and \eqref{e:q1asymp}, we find that the integrand is
\begin{equation}
\chi_N \, q_1(s_N) \, \Psi(y)^{-1/8} \exp\Big( - \frac{s_N^2}{64} \big( \Psi(y)^2 - 1 \big) - \frac{s_N}{8} \big( \Psi(y) - 1 \big) \Big) \big( 1 + o(1) \big) ,
\end{equation}
uniformly in $y \in [y_-, y_+ - \ell_N(x)]$. Since $s_N^2/64 \sim \log N$, this shows that the integral is concentrated within distance of order $(\log N)^{-1}$ of the points of $\partial_* I$, where $\sigma(y) = \sigma_I$. Near $\xi \in \partial_* I$, we have
\begin{equation}
\frac{\sigma(y)}{\sigma_I} - 1 = \frac{|\xi|}{\sigma_I^2} \, |y - \xi| + O\big( |y - \xi|^2 \big) ,
\end{equation}
since $\sigma'(\xi) = - \xi / \sigma_I$. A Laplace-type computation then shows that the contribution of $\xi$ to \eqref{e:overviewfirst} is
\begin{equation}
\chi_N \, q_1\big( N \sigma_I \ell_N(x) \big) \, \frac{\sigma_I^2}{2 |\xi| \log N} \big( 1 + o(1) \big) = \frac{a_N b_N}{2 \log N} \, \kappa_1 \frac{\sigma_I^3}{|\xi|} e^{-x} \big( 1 + o(1) \big) ,
\end{equation}
where we used the asymptotics
\begin{equation}\label{e:overviewq1}
q_1\big( N \sigma_I \ell_N(x) \big) = \kappa_1 e^{-x} b_N N^{-1} \big( 1 + o(1) \big) ,
\end{equation}
which follow from \eqref{e:q1asymp} and \eqref{e:overviewFell}. In this computation, the factor $b_N^{9/8}$ coming from $e^{-F(s)}$ is compensated by the factor $s^{-1/8} = b_N^{-1/8} (1 + o(1))$ in \eqref{e:q1asymp} and by the width $(\log N)^{-1}$ of the region near $\xi$, since $\chi_N / (N \log N) = 2 \sigma_I b_N^{-1} (1 + o(1))$. Since $a_N b_N / (2 \log N) \to 1$, summing over $\xi \in \partial_* I$ gives the constant $\Lambda_I$.

The estimates \eqref{e:overviewonehole}, \eqref{e:overviewquasi}, and \eqref{e:overviewfactor} are the main work of this paper. Two features make \eqref{e:overviewonehole} delicate. First, at $s = N \sigma_I \ell_N(x)$, the probability $q_1(s)$ is of order $N^{-1} \sqrt{\log N}$ by \eqref{e:overviewq1}. So we need a relative error bound for a polynomially small probability. Second, the interval contains of order $\sqrt{\log N}$ eigenvalues on average, so the convergence of local statistics on fixed scales does not suffice. For the GUE, the analogous estimates are proved in \cite{FengWei2025} using the determinantal structure. The asymptotics of the Fredholm determinant can then be derived from those of Toeplitz determinants \cite{DeiftItsKrasovskyZhou2007}, and the analogue of \eqref{e:overviewquasi} holds without the factor $e^{C \sqrt{\log N}}$ by the negative association of determinantal point processes \cite[Lemma~3.8]{BenArousBourgade2013}. Neither tool is available for the GOE. The GOE hole probability is the square root of the Fredholm determinant of an operator with a $2 \times 2$ matrix kernel \cite{TracyWidom1998,ForresterMays2009,Sinclair2009}, and one entry of this kernel contains the sign function $\sgn(x - y)/2$, which does not decay as $|x - y| \to \infty$.

\paragraph{Hermite coordinates.}
For the analysis, it is convenient to rescale the eigenvalues by the factor $\sqrt{N/2}$. By \eqref{e:goedensity}, the vector $(\sqrt{N/2} \, \lambda_1, \dots, \sqrt{N/2} \, \lambda_N)$ has a density proportional to $\prod_{i < j} |x_j - x_i| \prod_i e^{-x_i^2/2}$ on $\{x_1 < \dots < x_N\}$, and its empirical distribution concentrates on $[-\sqrt{2N}, \sqrt{2N}]$. We call $x = \sqrt{N/2} \, \lambda$ the \emph{Hermite coordinate}. In the rest of the overview, the variables $x$, $y$, $z$, and $t$ denote points in Hermite coordinates, and $J$ denotes a subset of $\R$ in Hermite coordinates. We let $\Pgoe_N(J)$ denote the probability that $\sqrt{N/2} \, \lambda_i \notin J$ for all $i$, so that $\Pgoe_N(\sqrt{N/2} \, A) = E_{1, N}(A)$ for every Borel set $A \subset \R$. For $\delta \in (0, 1/2]$, we define the \emph{bulk}
\begin{equation}\label{e:overviewbulk}
\bulk_N(\delta) = \big[ -(1 - \delta) \sqrt{2N}, (1 - \delta) \sqrt{2N} \big] ,
\end{equation}
and we note that $\sqrt{N/2} \, [-2 + \delta, 2 - \delta] = \bulk_N(\delta/2)$. We say that an estimate holds in the bulk if, for every fixed $\delta \in (0, 1/2]$, it holds uniformly on $\bulk_N(\delta)$, with constants depending on $\delta$. In Hermite coordinates, the mean spacing between consecutive eigenvalues near a point $x$ in the bulk is approximately $\pi / \sqrt{2N - x^2}$.

The Hermite functions $h_0, h_1, \dots$ are obtained by applying the Gram--Schmidt procedure in $L^2(\R)$ to the functions $x \mapsto x^n e^{-x^2/2}$, $n \ge 0$, and they form an orthonormal basis of $L^2(\R)$. We define
\begin{equation}\label{e:overviewKN}
K_N(x, y) = \sum_{n = 0}^{N - 1} h_n(x) h_n(y) , \qquad (K_N f)(y) = \int_\R K_N(y, x) f(x) \dd x ,
\end{equation}
for $x, y \in \R$ and for $f \in L^2(\R)$ or $f$ bounded. The restriction of $K_N$ to $L^2(\R)$ is the orthogonal projection onto the span of $h_0, \dots, h_{N - 1}$. For a measurable set $J \subset \R$ of finite measure, we let $K_J$ denote the integral operator on $L^2(J)$ with kernel $K_N(x, y)$, $x, y \in J$, and we define $\Dgue_N(J) = \det(\Id - K_J)$. This is the probability that the GUE of \zcref{r:constants}(c), with its eigenvalues multiplied by $\sqrt{N/2}$, has no eigenvalue in $J$ \cite{AndersonGuionnetZeitouni2010}. The eigenvalues of $K_J$ lie in $[0, 1)$, so $\Id - K_J$ is invertible.

\paragraph{Step 1: An exact formula and the endpoint matrix.}
Let $J = \bigcup_{j = 1}^r (a_j, b_j)$ with $a_1 < b_1 < a_2 < \dots < a_r < b_r$. For $t \in \R$, we define $\eps_t(x) = \sgn(x - t)/2$ and
\begin{equation}\label{e:overviewg}
g_t(y) = (K_N \eps_t)(y) - \eps_t(y) , \qquad y \in \R .
\end{equation}
The function $g_t$ measures how far the projection $K_N$ is from preserving the step function $\eps_t$. For $1 \le i \le r$, we define $G_i = (\Id - K_J)^{-1} (\1_J g_{a_i}) \in L^2(J)$. The identity $G_i = \1_J g_{a_i} + K_N(\1_J G_i)$ holds almost everywhere on $J$, and we use its right side to define $G_i$ at every point of $J$. For $t$ in the closure of $J$, we write $G_i(t^+)$ and $G_i(t^-)$ for the one-sided limits of $G_i$ at $t$, whenever they exist. We define the $r \times r$ \emph{endpoint matrix} $\Amat_N(J)$ by
\begin{equation}\label{e:overviewA}
\Amat_N(J)_{ij} = \1_{i = j} - \big( G_i(b_j^-) - G_i(a_j^+) \big) , \qquad 1 \le i, j \le r .
\end{equation}
We measure the size of $J$ by its \emph{microscopic length} $L(J) = \sqrt{2N} |J|$. For $J = \sqrt{N/2} \, A$, we have $L(J) = N |A|$, and $L(J)$ is comparable to the expected number of eigenvalues in $J$ when $J$ lies in the bulk. The first main step of the proof is the following estimate, which holds for even and odd $N$. For all $\delta \in (0, 1/2]$ and $r_0 \ge 1$, there exists $C > 0$ such that if $r \le r_0$, $\overline J \subset \bulk_N(\delta)$, and $L(J) \le N^{1/8}$, then
\begin{equation}\label{e:overviewexact}
\Pgoe_N(J)^2 = \Dgue_N(J) \Big( \det \Amat_N(J) + O\big( e^{C L(J)} N^{-1/2} \big) \Big) , \qquad |\det \Amat_N(J)| \le C e^{C L(J)} .
\end{equation}

The estimate \eqref{e:overviewexact} is the combination of an exact identity, \zcref{c:exact}, with \zcref{p:endpoint}. The exact identity expresses $\Pgoe_N(J)^2$, for every $N \ge 2$, as $\Dgue_N(J)$ multiplied by the determinant of an explicit $(r + 2) \times (r + 2)$ matrix. This matrix contains $\Amat_N(J)$ as a principal submatrix, and we will not need its other entries in the overview (they are given in \zcref{d:Msharp}). The square appears because $\Pgoe_N(J)$ is a Pfaffian, by de Bruijn's integration formula \cite{deBruijn1955}, and the square of a Pfaffian is a determinant. The Pfaffian involves the sign function on the complement of $J$. After an integration by parts (\zcref{l:derivativetrick}), the sign function appears only in $r$ rank-one operators, one for each interval of $J$, and these produce the endpoint matrix. The two remaining rows of the $(r + 2) \times (r + 2)$ matrix involve the Hermite functions $h_{N - 1}$ and $h_N$ of the highest degrees. Using the bound \eqref{e:overviewinverse} below, we show that all entries of the matrix are $O(e^{C L(J)})$, and that these two rows are within $O(e^{C L(J)} N^{-1/2})$ of the corresponding rows of the identity matrix. Since the determinant is linear in each row, this gives \eqref{e:overviewexact}.

For even $N$, an identity of this type was found by Tracy and Widom \cite[(33)]{TracyWidom1996}, who manipulated Fredholm determinants of operators with matrix kernels and justified these manipulations in detail for the sine kernel. We derive the exact identity differently, starting from de Bruijn's formula and using only finite-dimensional linear algebra. Our derivation applies to even and odd $N$ (this matters, since \zcref{t:main} concerns all $N$). For odd $N$, one of the two remaining rows involves a function whose smallness is not evident from its definition. To describe it, we define
\begin{equation}\label{e:overvieweps}
(\eps f)(x) = \frac12 \int_\R \sgn(x - z) f(z) \dd z , \qquad p_N(y) = \int_\R K_N(x, y) \dd x , \qquad x, y \in \R ,
\end{equation}
for all $N \ge 1$ and integrable functions $f$. Note that $\int_\R h_{N - 1}(x) \dd x \neq 0$ when $N$ is odd. For odd $N$, we define
\begin{equation}\label{e:overviewomega}
c_N = \Big( \frac N2 \Big)^{-1/4} \Big( \int_\R h_{N - 1}(x) \dd x \Big)^{-1} , \qquad \omega_N = \Big( \frac N2 \Big)^{1/4} \eps h_N + c_N .
\end{equation}
In the bulk, the two terms of $\omega_N$ are of order one. However, we show in \zcref{l:parity}(b) that their sum is exactly a multiple of $1 - p_N$:
\begin{equation}\label{e:overviewodd}
\omega_N = c_N (1 - p_N) , \qquad N \text{ odd} .
\end{equation}
We also show that $p_N(y) = 1 + O(N^{-1/2})$ in the bulk, so \eqref{e:overviewodd} shows that $\omega_N = O(N^{-1/2})$ there. Landon and Sosoe observed a cancellation of this kind through an asymptotic analysis of integrals of Hermite functions \cite[Lemma~8.14 and Appendix~D]{LandonSosoe2022}.

\paragraph{Step 2: The sign function becomes local.}
Observe that the entries of $\Amat_N(J)$ involve the functions $g_{a_i}$ only on $J$. The second step, which has no counterpart in the proof for the GUE, is the observation that $g_t$ is concentrated near $t$, and that it can be computed from the kernel $K_N$ in the bulk. We set $\varphi_N = (N/2)^{1/4} h_N$ and $\psi_N = (N/2)^{1/4} h_{N - 1}$. We prove in \zcref{l:signkernel} the exact identity
\begin{equation}\label{e:overviewstep}
(K_N \eps_t)(y) = (\eps \varphi_N)(y) \, (\eps \psi_N)(y) + \int_t^y K_N(x, y) \dd x , \qquad t, y \in \R .
\end{equation}
In the bulk, the product $(\eps \varphi_N)(\eps \psi_N)$ is $O(N^{-1/2})$. Indeed, one of its factors is a constant multiple of $p_N$, namely $\eps \varphi_N = - c_N p_N$ for odd $N$ by \eqref{e:overviewodd}, and $\eps \psi_N = - p_N / \int_\R \varphi_N$ for even $N$ by \zcref{l:parity}(a). This factor is bounded in the bulk, since $p_N = 1 + O(N^{-1/2})$ there, and $c_N$ and $(\int_\R \varphi_N)^{-1}$ are bounded. By parity, the other factor equals $\int_0^y \varphi_N$ or $\int_0^y \psi_N$, which is the integral of an oscillating function and is $O(N^{-1/2})$ in the bulk. The integral of the kernel between $t$ and $y$ is close to $\sgn(y - t)/2$ when $\sqrt N |y - t|$ is large. In \zcref{c:gt}, we deduce from \eqref{e:overviewstep} that
\begin{equation}\label{e:overviewgbound}
|g_t(y)| \le \frac{C}{1 + \sqrt N |y - t|} + C N^{-1/4} , \qquad t, y \in \bulk_N(\delta) .
\end{equation}
The term $C N^{-1/4}$ comes from the oscillation of the kernel far from the diagonal and is not optimal. Further, we define
\begin{equation}\label{e:overviewgsine}
g_\sine(u) = - \sgn(u) \int_{|u|}^\infty S(v) \dd v , \qquad u \in \R ,
\end{equation}
and we show that
\begin{equation}\label{e:overviewglocal}
\Big| g_t\Big( t + \frac{u}{\varrho_N(t)} \Big) - g_\sine(u) \Big| \le C (1 + |u|) N^{-1/2} , \qquad t \in \bulk_N(\delta) , \ 0 < |u| \le N^{1/4} ,
\end{equation}
where we write $\varrho_N(t) = \sqrt{2N - t^2}$; see \zcref{f:gt}. The functions $p_N$, $\eps \varphi_N$, and $\eps \psi_N$ are defined by integrals over $\R$, but by parity each of them equals an explicit constant plus an integral over $[0, y]$ of Hermite functions, and the constants are given by exact formulas for $\int_\R h_n$. In particular, our arguments use the Hermite functions only in the bulk, and we will not need any estimates for them near the spectral edges $\pm \sqrt{2N}$.

\begin{figure}[t]
\centering
\includegraphics[width=\textwidth]{figures/gt.pdf}
\caption{The function $g_t = K_N \eps_t - \eps_t$ for $N = 60$ and $t = 2$. Panel (a) shows $K_N \eps_t$ and the step function $\eps_t$; the dotted lines are at $\pm \sqrt{2N}$. In the bulk, their difference is concentrated near $t$, up to small oscillations, while near and beyond the spectral edges it is of order one. Panel (b) shows $g_t(t + u/\varrho_N(t))$, where $\varrho_N(t) = \sqrt{2N - t^2}$, together with its limit $g_\sine(u)$.}
\label{f:gt}
\end{figure}

The bound \eqref{e:overviewgbound} is what makes the GOE hole probability factorize. We now assume, in addition to the assumptions of \eqref{e:overviewexact}, that the distance between any two intervals of $J$ is at least $1/\sqrt 2$, which corresponds to the distance $N^{-1/2}$ in the original coordinates $\lambda$, and that $L(J) \le \gamma \sqrt{\log N}$ for a fixed $\gamma > 0$. Then $g_{a_i}(y) = O(N^{-1/4})$ for $y$ in the intervals $(a_j, b_j)$ with $j \neq i$, by \eqref{e:overviewgbound}. Further, $|K_N(x, y)| \le C$ for $x$ and $y$ in different intervals, while $K_N(x, x)$ is of order $\sqrt N$ in the bulk. Using these bounds and the bound \eqref{e:overviewinverse} below, we show in \zcref{l:offdiagonal} that the off-diagonal entries of $\Amat_N(J)$ are $O(e^{C L(J)} N^{-1/4})$, and that each diagonal entry $\Amat_N(J)_{jj}$ differs by $O(e^{C L(J)} N^{-1/2})$ from the endpoint matrix of the single interval $J_j = (a_j, b_j)$. The latter is a $1 \times 1$ matrix, which we identify with its entry. We write each interval $J_j$ in the form of Step 4 below, with $x_0 = a_j$ and $s = 2 \varrho_N(a_j) (b_j - a_j) \le C L(J)$. The lower bound $\det \Amat_N(J_j) \ge c e^{-C L(J)}$ from Step 4 then shows that
\begin{equation}\label{e:overviewAfactor}
\det \Amat_N(J) = \big( 1 + O(e^{C L(J)} N^{-1/4}) \big) \prod_{j = 1}^r \Amat_N\big( (a_j, b_j) \big) .
\end{equation}
Next, the operator $K_J$ differs from the direct sum of the operators $K_{J_j}$ by the operator whose kernel equals $K_N(x, y)$ for $x$ and $y$ in different intervals of $J$ and vanishes otherwise. Its Hilbert--Schmidt norm is bounded by $C |J| \le C L(J) N^{-1/2}$. Together with \eqref{e:overviewinverse} below, a standard estimate for Fredholm determinants (\zcref{l:fischer}) gives
\begin{equation}\label{e:overviewDfactor}
\Dgue_N(J) = \big( 1 + O(e^{C L(J)} N^{-1}) \big) \prod_{j = 1}^r \Dgue_N(J_j) .
\end{equation}
To prove \eqref{e:overviewfactor}, we set $J = \sqrt{N/2} \, \bigcup_j A_j$, which satisfies $\overline J \subset \bulk_N(\delta/2)$ and $L(J) = N \sum_j |A_j| \le \gamma \sqrt{\log N} \le N^{1/8}$. We apply \eqref{e:overviewexact} to $J$ and to each $J_j$, and we insert \eqref{e:overviewAfactor} and \eqref{e:overviewDfactor}. Using the lower bound on $\det \Amat_N(J_j)$ once more and taking square roots, we obtain \eqref{e:overviewfactor}.

\paragraph{Step 3: An inverse bound from the uncertainty principle.}
The estimates in Steps 1 and 2 involve the norm of $(\Id - K_J)^{-1}$ on $L^2(J)$. Recall that the eigenvalues of $K_J$ lie in $[0, 1)$. Therefore, $\| (\Id - K_J)^{-1} \| \le \Dgue_N(J)^{-1}$. However, for an interval $J$, the right side grows like $e^{c L(J)^2}$, which is a power of $N$ when $L(J)$ is of order $\sqrt{\log N}$, and this would overwhelm the polynomially small errors. The third main step is the much better bound
\begin{equation}\label{e:overviewinverse}
\| (\Id - K_J)^{-1} \| \le C e^{C L(J)} , \qquad |J| \le 1 ,
\end{equation}
which we prove in \zcref{p:inverse} using a quantitative uncertainty principle of Nazarov \cite{Nazarov1993,Jaming2007}. To state it, we define the Fourier transform of $f \in L^1(\R) \cap L^2(\R)$ by
\begin{equation}
(\mathcal F f)(\xi) = \frac{1}{\sqrt{2\pi}} \int_\R f(x) e^{-i x \xi} \dd x , \qquad \xi \in \R ,
\end{equation}
and we extend $\mathcal F$ to a unitary operator on $L^2(\R)$ as usual. The principle states that for all measurable sets $A, \Sigma \subset \R$ of finite measure and all $f \in L^2(\R)$,
\begin{equation}\label{e:overviewnazarov}
\| f \|_2^2 \le C e^{C |A| |\Sigma|} \Big( \int_{\R \setminus A} |f(x)|^2 \dd x + \int_{\R \setminus \Sigma} |(\mathcal F f)(\xi)|^2 \dd \xi \Big) .
\end{equation}
We let $V_N$ denote the real span of $h_0, \dots, h_{N - 1}$. Each $h_n$ is an eigenfunction of $\mathcal F$, so $\mathcal F$ maps $V_N$ to the complex span of $h_0, \dots, h_{N - 1}$, and we show that for every $f$ in this complex span, the integral of $|f|^2$ over the complement of $[-8 \sqrt N, 8 \sqrt N]$ is bounded by $e^{-N} \| f \|_2^2$. Applying \eqref{e:overviewnazarov} with $A = J$ and $\Sigma = [-8 \sqrt N, 8 \sqrt N]$ to $f \in V_N$, and using $|J| |\Sigma| = 8 \sqrt 2 \, L(J)$, we find that
\begin{equation}
\| f \|_2^2 \le C e^{C L(J)} \big( \| f \|_2^2 - \| \1_J f \|_2^2 + e^{-N} \| f \|_2^2 \big) , \qquad f \in V_N .
\end{equation}
Since $|J| \le 1$, we have $L(J) \le \sqrt{2N}$, so the term $e^{-N}$ is negligible compared to $e^{-C L(J)}$, and we find $\| \1_J f \|_2^2 \le (1 - c e^{-C L(J)}) \| f \|_2^2$ for all $f \in V_N$. Since $\| K_J \|$ is the supremum of $\| \1_J f \|_2^2 / \| f \|_2^2$ over nonzero $f \in V_N$, this implies \eqref{e:overviewinverse}. For a bounded interval $J' \subset \R$, we let $S_{J'}$ denote the integral operator on $L^2(J')$ with kernel $S(u - v)$. Recall that $S(u - v)$ is the kernel of the orthogonal projection onto the space of functions in $L^2(\R)$ whose Fourier transform vanishes outside $[-1, 1]$. The same argument, with $V_N$ replaced by this space, $A = J'$, and $\Sigma = [-1, 1]$, so that $|A| |\Sigma| = 2 |J'|$ and there is no tail term, shows that $\| (\Id - S_{J'})^{-1} \| \le C e^{C |J'|}$.

\paragraph{Step 4: Comparison with the sine kernel.}
We now explain the proof of the one-interval estimate \eqref{e:overviewonehole}. We fix a point $y_0 \in [-2 + \delta, 2 - \delta]$ in eigenvalue coordinates and a number $0 < s \le \gamma \sqrt{\log N}$, and we define
\begin{equation}
x_0 = \sqrt{N/2} \, y_0 , \qquad \varrho = \varrho_N(x_0) , \qquad J = \Big( x_0 , x_0 + \frac{s}{2 \varrho} \Big) .
\end{equation}
Then $\varrho = \sqrt{N/2} \, \sigma(y_0)$, so $J = \sqrt{N/2} \, (y_0, y_0 + s/(N \sigma(y_0)))$, and $L(J) = s / \sigma(y_0) \le C s$. We apply \eqref{e:overviewexact} to $J$ with $r = 1$ and $\delta/4$ in place of $\delta$, so that the error term there is $O(e^{C s} N^{-1/2})$. We show in \zcref{l:localsine}, using the Plancherel--Rotach asymptotics of the Hermite functions \cite{PlancherelRotach1929,Szego1975}, that
\begin{equation}\label{e:overviewlocalsine}
\Big| \frac{1}{\varrho} K_N\Big( x_0 + \frac{u}{\varrho} , x_0 + \frac{v}{\varrho} \Big) - S(u - v) \Big| \le C N^{-1/2} , \qquad |u|, |v| \le 2 N^{1/4} .
\end{equation}
This shows that after the change of variables $x = x_0 + u/\varrho$, the operator $K_J$ is close to $S_{(0, s/2)}$. For $\ell > 0$, we let $\Amat_\sine(\ell)$ denote the number given by the definitions of Step 1 with $r = 1$, after $K_N$, $J$, and $g_{a_1}(x_0 + \cdot / \varrho)$ are replaced by the integral operator with kernel $S(u - v)$ on $L^2(\R)$, the interval $(0, \ell)$, and $g_\sine$. So $\Amat_\sine(s/2)$ is the sine-kernel analogue of $\det \Amat_N(J)$. Using \eqref{e:overviewlocalsine}, the convergence of $g_{x_0}(x_0 + u/\varrho)$ to $g_\sine(u)$ from Step 2, and the bounds on $(\Id - K_J)^{-1}$ and $(\Id - S_{(0, s/2)})^{-1}$ from Step 3, we show that
\begin{equation}
\begin{gathered}
\Dgue_N(J) = \det(\Id - S_{(0, s/2)}) \big( 1 + O(e^{C s} N^{-1/2}) \big) , \\
\det \Amat_N(J) = \Amat_\sine(s/2) + O(e^{C s} N^{-1/2}) .
\end{gathered}
\end{equation}
We translate $(0, s/2)$ to $(-s/4, s/4)$. On this symmetric interval, the operator with kernel $S(u - v)$ preserves the even and the odd functions. Its restrictions to these two subspaces are unitarily equivalent to $S_{+, s/4}$ and $S_{-, s/4}$, respectively. Therefore $\det(\Id - S_{(0, s/2)}) = D_+(s/4) D_-(s/4)$. Further, by an identity of Tracy and Widom \cite[(19), (25), (30)]{TracyWidom1996}, we have $\Amat_\sine(s/2) = D_+(s/4)/D_-(s/4)$. By \eqref{e:q1asymp} and the positivity and continuity of $D_\pm$, this ratio is bounded below by $c e^{-s/4}$ for all $s > 0$. Hence $\det \Amat_N(J) \ge c e^{-C L(J)}$, and the error term in \eqref{e:overviewexact} is small compared to $\det \Amat_N(J)$, where we used $s \le \gamma \sqrt{\log N}$. Combining these estimates, we have
\begin{equation}
\begin{split}
\Pgoe_N(J)^2 & = D_+(s/4) D_-(s/4) \cdot \frac{D_+(s/4)}{D_-(s/4)} \big( 1 + O(e^{C s} N^{-1/2}) \big) \\
& = q_1(s)^2 \big( 1 + O(e^{C s} N^{-1/2}) \big) .
\end{split}
\end{equation}
By the description of $J$ above, $\Pgoe_N(J)$ is the left side of \eqref{e:overviewonehole} with $y = y_0$, so this is \eqref{e:overviewonehole}.

Finally, we prove the upper bound \eqref{e:overviewquasi}, which requires no separation assumption. Let $J_j = \sqrt{N/2} \, A_j$ for $1 \le j \le r$, and let $J = \bigcup_j J_j$. Then $\overline J \subset \bulk_N(\delta/2)$ and $L(J) \le \gamma \sqrt{\log N} \le N^{1/8}$. By \eqref{e:overviewexact} with $\delta/2$ in place of $\delta$, we have $\Pgoe_N(J)^2 \le C e^{C L(J)} \Dgue_N(J)$. Next, the operator $\Id - K_J$ is positive definite, since $0 \le K_J < \Id$, so an inequality of Fischer type for Fredholm determinants of such operators gives $\Dgue_N(J) \le \prod_j \Dgue_N(J_j)$. Finally, we write $J_j$ in the form of the first part of this step, with $y_0$ the left endpoint of $A_j$ and $s = N \sigma(y_0) |A_j| \le 2 \gamma \sqrt{\log N}$. Using \eqref{e:overviewexact} for the single interval $J_j$ and the lower bound $\det \Amat_N(J_j) \ge c e^{-C L(J_j)}$, we find $\Dgue_N(J_j) \le C e^{C L(J_j)} \Pgoe_N(J_j)^2$. Combining these bounds, taking square roots, and using $L(J) \le \gamma \sqrt{\log N}$, we obtain \eqref{e:overviewquasi}.

\section{Organization of the paper}\label{s:organization}

\zcref{s:hermite} collects classical facts about Hermite functions and proves the exact identities, including \eqref{e:overviewodd} and \eqref{e:overviewstep}. \zcref{s:exact} proves the exact identity for $\Pgoe_N(J)^2$ described in Step 1. \zcref{s:bulk} contains the asymptotic estimates for Hermite functions in the bulk, and \zcref{s:inverse} proves the inverse bound \eqref{e:overviewinverse}. \zcref{s:endpoint} reduces this identity to the endpoint matrix and shows that the endpoint matrix of separated intervals is nearly diagonal. \zcref{s:oneinterval} proves the one-interval estimate \eqref{e:overviewonehole}, and \zcref{s:severalintervals} proves the estimates \eqref{e:overviewquasi} and \eqref{e:overviewfactor}. \zcref{s:firstmoment} computes the first moment \eqref{e:overviewfirst}, and \zcref{s:factorial} completes the proofs of \zcref{t:main}, \zcref{c:ordered}, and \zcref{r:universality}. \zcref{s:poisson} contains the criterion for Poisson convergence. The main logical dependencies between the sections are shown in \zcref{f:structure}.

\begin{figure}[t]
\centering
\input{figures/structure}
\caption{The main logical dependencies between the sections. An arrow from one section to another means that the second section uses results of the first. The classical facts and exact identities of \zcref{s:hermite} are used throughout.}
\label{f:structure}
\end{figure}