\chapter{Pair statistics of the Gaussian orthogonal ensemble}\label{s:goe}

In this appendix, we prove \zcref{p:goepairs}. Throughout, we fix $\alpha \in (0, 1/4)$, $p > 1$, and a function $f$ as in \zcref{p:goepairs}, together with the bulk interval $I$ and the constant $C_f$ given there.

We follow the approach of Feng, Li, and Yao \cite{feng2024small}, who treat the Gaussian symplectic ensemble in detail and sketch the orthogonal case. The correlation functions of the GOE are Pfaffians of a $2 \times 2$ matrix kernel built from Hermite functions. The mean of the statistic in \eqref{e:prelim-goepairs} is an integral of the two-point function, and the variance is controlled by the two-, three-, and four-point functions. However, the GOE analogues of the kernel bounds in \cite[Lemma~3.2]{feng2024small} are asserted in \cite[Section~7.1]{feng2024small} without proof, and the expansion of the two-point function that we need is not stated there, so we give complete proofs. They rest on elementary estimates for Hermite functions, proved in \zcref{s:goe-hermite}, and on the density of states, which we deduce from eigenvalue rigidity in \zcref{l:goe-density}. The main step is the exact identity for the slope of the two-point function at the diagonal in \zcref{l:goe-two}; together with the density of states, it determines the two-point function at separations below $N^{-1}$. The most technical part of the argument is the decorrelation estimate for two separated pairs of nearby points in \zcref{l:goe-corr}, which controls the variance.

In this appendix only, we abbreviate $\lambda_i = \lambda_i(G)$. For $\delta \in (0, 1)$, we set $I_\delta = [-2 + \delta, 2 - \delta]$. Constants $C, c > 0$ and thresholds $N_0$ may depend on $\delta$, $\alpha$, $p$, $I$, and $f$, but not on $N$.

\section{Hermite functions}\label{s:goe-hermite}

Let $\mathrm{He}_n$ denote the monic Hermite polynomials orthogonal with respect to the weight $e^{-w^2/2}$. They are determined by the generating function
\begin{equation}\label{e:goe-genfun}
\sum_{n = 0}^\infty \frac{\mathrm{He}_n(w) \zeta^n}{n!} = e^{w \zeta - \zeta^2/2}, \qquad w, \zeta \in \R.
\end{equation}
For integers $n \ge 0$, we define the Hermite functions
\begin{equation}\label{e:goe-hermite-def}
\mathsf h_n(w) = (2 \pi)^{-1/4} (n!)^{-1/2} \mathrm{He}_n(w) e^{-w^2/4}, \qquad w \in \R,
\end{equation}
and we set $\mathsf h_{-1} = 0$. The functions $\mathsf h_n$ form an orthonormal family in $L^2(\R)$ (see \cite[Section~3.2]{anderson2010introduction}), and $\mathsf h_n(-w) = (-1)^n \mathsf h_n(w)$. Differentiating \eqref{e:goe-genfun} in $w$ and in $\zeta$, we find $\mathrm{He}_n' = n \mathrm{He}_{n - 1}$ and $w \mathrm{He}_n = \mathrm{He}_{n + 1} + n \mathrm{He}_{n - 1}$. These two relations imply $\mathrm{He}_n'' - w \mathrm{He}_n' + n \mathrm{He}_n = 0$. We set
\begin{equation}\label{e:goe-Q-def}
\mathsf e_n(w) = n + \frac12 - \frac{w^2}{4}, \qquad n \ge 0, \ w \in \R.
\end{equation}
In terms of the Hermite functions, the three relations above become
\begin{equation}\label{e:goe-hermite-id}
w \mathsf h_n = \sqrt{n + 1} \, \mathsf h_{n + 1} + \sqrt n \, \mathsf h_{n - 1}, \qquad \mathsf h_n' = \frac{\sqrt n \, \mathsf h_{n - 1} - \sqrt{n + 1} \, \mathsf h_{n + 1}}{2}, \qquad \mathsf h_n'' = - \mathsf e_n \mathsf h_n,
\end{equation}
for all $n \ge 0$. For $w, w' \in \R$ and $j \ge 0$, the first identity in \eqref{e:goe-hermite-id} gives
\begin{align*}
(w - w') \mathsf h_j(w) \mathsf h_j(w') &= \sqrt{j + 1} \, \big( \mathsf h_{j + 1}(w) \mathsf h_j(w') - \mathsf h_j(w) \mathsf h_{j + 1}(w') \big) \\
&\quad - \sqrt j \, \big( \mathsf h_j(w) \mathsf h_{j - 1}(w') - \mathsf h_{j - 1}(w) \mathsf h_j(w') \big).
\end{align*}
Summing over $0 \le j \le n - 1$, we see that the right side telescopes, and we obtain the Christoffel--Darboux formula
\begin{equation}\label{e:goe-cd}
\sum_{j = 0}^{n - 1} \mathsf h_j(w) \mathsf h_j(w') = \sqrt n \, \frac{\mathsf h_n(w) \mathsf h_{n - 1}(w') - \mathsf h_{n - 1}(w) \mathsf h_n(w')}{w - w'}, \qquad n \ge 1, \ w \neq w'.
\end{equation}
Letting $w' \to w$ in \eqref{e:goe-cd}, we also have
\begin{equation}\label{e:goe-cd-diag}
\sum_{j = 0}^{n - 1} \mathsf h_j(w)^2 = \sqrt n \big( \mathsf h_n'(w) \mathsf h_{n - 1}(w) - \mathsf h_{n - 1}'(w) \mathsf h_n(w) \big).
\end{equation}

For integers $n \ge 0$ and $w \in \R$, we define
\begin{equation}\label{e:goe-F-def}
\mathsf F_n(w) = \frac12 \int_\R \sgn(w - w') \mathsf h_n(w') \dd w',
\end{equation}
and we set $\mathsf F_{-1} = 0$. Then $\mathsf F_n' = \mathsf h_n$. Since $\mathsf h_n'$ is integrable and $\mathsf h_n(w) \to 0$ as $|w| \to \infty$, we have
\begin{equation*}
\frac12 \int_\R \sgn(w - w') \mathsf h_n'(w') \dd w' = \mathsf h_n(w), \qquad w \in \R.
\end{equation*}
Inserting the second identity in \eqref{e:goe-hermite-id} into the left side, we obtain
\begin{equation}\label{e:goe-F-id}
\mathsf h_n = \frac{\sqrt n \, \mathsf F_{n - 1} - \sqrt{n + 1} \, \mathsf F_{n + 1}}{2}, \qquad n \ge 0.
\end{equation}

We will need the following estimates.

\begin{lemma}\label{l:goe-hermite}
Let $\delta \in (0, 1]$. There exists $C > 0$, depending only on $\delta$, such that the following statements hold for every integer $n \ge 1$.
\begin{enumerate}[label=(\alph*)]
\item For $|w| \le (2 - \delta) \sqrt n$ and $k \in \{0, 1, 2, 3\}$, we have $|\mathsf h_n^{(k)}(w)| \le C n^{k/2 - 1/4}$.
\item For $|w|, |w'| \le (2 - \delta) \sqrt n$, we have $\big| \int_w^{w'} \mathsf h_n \big| \le C n^{-3/4}$.
\item If $n$ is even, then $C^{-1} n^{-1/4} \le \int_\R \mathsf h_n \le C n^{-1/4}$. If $n$ is odd, then $\big| \int_0^\infty \mathsf h_n \big| \le C n^{-1/4}$.
\item For $|w| \le (2 - \delta) \sqrt n$, we have $|\mathsf F_n(w)| \le C n^{-1/4}$.
\end{enumerate}
\end{lemma}

\begin{proof}
We begin with (a) for $k \in \{0, 1\}$, using the argument of Sonine and P\'olya (see \cite[Section~7.31]{szego1975orthogonal}). On the interval $[0, (4n + 2)^{1/2})$, on which $\mathsf e_n > 0$, we define
\begin{equation*}
\Theta(w) = \mathsf h_n(w)^2 + \frac{\mathsf h_n'(w)^2}{\mathsf e_n(w)}.
\end{equation*}
Using $\mathsf h_n'' = - \mathsf e_n \mathsf h_n$ and $\mathsf e_n' = - w/2$, we find
\begin{equation*}
\Theta' = 2 \mathsf h_n \mathsf h_n' + \frac{2 \mathsf h_n' \mathsf h_n''}{\mathsf e_n} - \frac{\mathsf e_n' (\mathsf h_n')^2}{\mathsf e_n^2} = \frac{w (\mathsf h_n')^2}{2 \mathsf e_n^2} \ge 0.
\end{equation*}
Hence $\Theta$ is nondecreasing on $[0, (4n + 2)^{1/2})$.

Set $w_1 = (2 - \delta) \sqrt n$ and $w_2 = (2 - \delta/2) \sqrt n$. Then $w_2 < (4n + 2)^{1/2}$, and
\begin{equation*}
\mathsf e_n(w) \ge n \big( 1 - (1 - \delta/4)^2 \big) \ge \frac{\delta n}{4}, \qquad w \in [w_1, w_2].
\end{equation*}

Let $0 \le w \le w_1$. Since $\Theta$ is nondecreasing, $\Theta(w) \le \Theta(w')$ for all $w' \in [w_1, w_2]$. Averaging over $w'$, we obtain
\begin{align*}
\Theta(w) \le \frac{1}{w_2 - w_1} \int_{w_1}^{w_2} \Theta(w') \dd w' &\le \frac{2}{\delta \sqrt n} \Big( \int_\R \mathsf h_n^2 + \frac{4}{\delta n} \int_\R (\mathsf h_n')^2 \Big) \\
&= \frac{2}{\delta \sqrt n} \Big( 1 + \frac{2n + 1}{\delta n} \Big) \le \frac{8}{\delta^2 \sqrt n},
\end{align*}
where in the equality we used the identity
\begin{equation*}
\int_\R (\mathsf h_n')^2 = \frac{2n + 1}{4},
\end{equation*}
which follows from the second identity in \eqref{e:goe-hermite-id} and orthonormality. Consequently,
\begin{equation*}
\mathsf h_n(w)^2 \le \Theta(w) \le \frac{8}{\delta^2 \sqrt n}, \qquad \mathsf h_n'(w)^2 \le \mathsf e_n(w) \Theta(w) \le \frac{8 (n + 1)}{\delta^2 \sqrt n}.
\end{equation*}
The same bounds hold for $-w_1 \le w \le 0$ by the parity of $\mathsf h_n$, which proves (a) for $k \in \{0, 1\}$.

The cases $k \in \{2, 3\}$ follow from the identity $\mathsf h_n'' = - \mathsf e_n \mathsf h_n$ and its consequence
\begin{equation}\label{e:goe-hermite-third}
\mathsf h_n''' = (w/2) \mathsf h_n - \mathsf e_n \mathsf h_n',
\end{equation}
together with $|\mathsf e_n(w)| \le n + 1$ and $|w| \le 2 \sqrt n$.

Next, we prove (b). We have
\begin{equation*}
\mathsf e_n(w) \ge n \big( 1 - (1 - \delta/2)^2 \big) \ge \frac{\delta n}{2}, \qquad |w| \le (2 - \delta) \sqrt n.
\end{equation*}
On this interval, $\mathsf h_n = - \mathsf h_n''/\mathsf e_n$, and two integrations by parts give
\begin{equation*}
\int_w^{w'} \mathsf h_n = - \Big[ \frac{\mathsf h_n'}{\mathsf e_n} - \Big( \frac{1}{\mathsf e_n} \Big)' \mathsf h_n \Big]_w^{w'} - \int_w^{w'} \Big( \frac{1}{\mathsf e_n} \Big)'' \mathsf h_n.
\end{equation*}
We have
\begin{equation*}
\Big( \frac{1}{\mathsf e_n} \Big)' = \frac{w}{2 \mathsf e_n^2}, \qquad \Big( \frac{1}{\mathsf e_n} \Big)'' = \frac{1}{2 \mathsf e_n^2} + \frac{w^2}{2 \mathsf e_n^3},
\end{equation*}
so $1/\mathsf e_n$, $(1/\mathsf e_n)'$, and $(1/\mathsf e_n)''$ are bounded by $C n^{-1}$, $C n^{-3/2}$, and $C n^{-2}$ on the interval. Using (a) and $|w' - w| \le 4 \sqrt n$, we find
\begin{equation*}
\Big| \int_w^{w'} \mathsf h_n \Big| \le C \big( n^{-1} n^{1/4} + n^{-3/2} n^{-1/4} + \sqrt n \, n^{-2} n^{-1/4} \big) \le C n^{-3/4}.
\end{equation*}

We now prove (c). In the computations below, the exchange of sums and integrals is justified by absolute convergence.

Suppose first that $n = 2l$ with $l \ge 1$. Multiplying \eqref{e:goe-genfun} by $e^{-w^2/4}$ and integrating over $w \in \R$, we have
\begin{equation*}
\sum_{n' = 0}^\infty \frac{\zeta^{n'}}{n'!} \int_\R e^{-w^2/4} \mathrm{He}_{n'}(w) \dd w = \int_\R e^{-(w - 2\zeta)^2/4 + \zeta^2/2} \dd w = 2 \sqrt \pi \, e^{\zeta^2/2}.
\end{equation*}
Comparing the coefficients of $\zeta^{2l}$, we obtain
\begin{equation*}
\int_\R e^{-w^2/4} \mathrm{He}_{2l}(w) \dd w = \frac{2 \sqrt \pi \, (2l)!}{2^l l!},
\end{equation*}
and therefore
\begin{equation*}
\int_\R \mathsf h_{2l} = \frac{2 \sqrt \pi}{(2 \pi)^{1/4}} \cdot \frac{\sqrt{(2l)!}}{2^l l!} = \frac{2 \sqrt \pi}{(2 \pi)^{1/4}} \Big( 4^{-l} \binom{2l}{l} \Big)^{1/2}.
\end{equation*}
The claim for even $n$ now follows from the elementary bounds
\begin{equation*}
(4l)^{-1/2} \le 4^{-l} \binom{2l}{l} \le (\pi l)^{-1/2}.
\end{equation*}

Suppose next that $n = 2l + 1$ with $l \ge 1$. For $\zeta \in \R$, we have
\begin{equation*}
\int_0^\infty e^{-w^2/4} e^{w \zeta - \zeta^2/2} \dd w = e^{\zeta^2/2} \int_{-2\zeta}^\infty e^{-v^2/4} \dd v = \sqrt \pi \, e^{\zeta^2/2} + 2\zeta \int_0^1 e^{\zeta^2 (1/2 - v^2)} \dd v.
\end{equation*}
The first term on the right side is even in $\zeta$. Let
\begin{equation*}
\mathsf j_l = \int_0^1 (1 - 2 v^2)^l \dd v.
\end{equation*}
Comparing the coefficients of $\zeta^{2l + 1}$, we obtain
\begin{equation*}
\int_0^\infty \mathsf h_{2l + 1} = \frac{2^{1 - l} \sqrt{(2l + 1)!}}{(2 \pi)^{1/4} l!} \mathsf j_l = \frac{2 \sqrt{2l + 1}}{(2 \pi)^{1/4}} \Big( 4^{-l} \binom{2l}{l} \Big)^{1/2} \mathsf j_l.
\end{equation*}

To bound $\mathsf j_l$, we use the inequalities
\begin{align*}
0 \le 1 - 2v^2 &\le e^{-2 v^2}, & v &\in [0, 2^{-1/2}], \\
0 \le 2v^2 - 1 &\le 1 - 2(1 - v) \le e^{-2 (1 - v)}, & v &\in [2^{-1/2}, 1].
\end{align*}
They give
\begin{equation*}
|\mathsf j_l| \le \int_0^\infty e^{-2 l v^2} \dd v + \int_0^\infty e^{-2 l v} \dd v \le 2 l^{-1/2},
\end{equation*}
and therefore
\begin{equation*}
\Big| \int_0^\infty \mathsf h_{2l + 1} \Big| \le C l^{1/2} l^{-1/4} l^{-1/2} \le C n^{-1/4}.
\end{equation*}
The case $n = 1$ is trivial.

Finally, we prove (d). When $n$ is even, $\mathsf h_n$ is even, and when $n$ is odd, $\mathsf h_n$ is odd. Therefore
\begin{equation*}
\mathsf F_n(w) = \int_0^w \mathsf h_n \quad (n \text{ even}), \qquad \mathsf F_n(w) = \int_0^w \mathsf h_n - \int_0^\infty \mathsf h_n \quad (n \text{ odd}),
\end{equation*}
and (b) with $w' = 0$ and (c) give the claim.
\end{proof}

\section{The Pfaffian formula}\label{s:goe-pfaffian}

The entries $G_{ab}$ with $a \le b$ are independent centered Gaussian random variables with variance $2/N$ for $a = b$ and $1/N$ for $a < b$, so $G$ has density proportional to $\exp(-N \Tr A^2/4)$ with respect to Lebesgue measure on the space of $N \times N$ real symmetric matrices $A$. Therefore, the eigenvalues of $G$, in uniformly random order, have the joint density
\begin{equation}\label{e:goe-density}
\wp_N(x_1, \dots, x_N) = \frac{1}{\mathcal Z_N} \prod_{1 \le i < j \le N} |x_i - x_j| \prod_{i = 1}^N e^{-N x_i^2/4}
\end{equation}
on $\R^N$, where $\mathcal Z_N$ is a normalizing constant (see \cite[Section~2.5]{anderson2010introduction}; this is \cite[(1.1)]{feng2024small} for the orthogonal ensemble). For $n \in [N]$, we define the $n$-point correlation function
\begin{equation}\label{e:goe-corr-def}
R_{n, N}(x_1, \dots, x_n) = \frac{N!}{(N - n)!} \int_{\R^{N - n}} \wp_N(x_1, \dots, x_N) \dd x_{n + 1} \cdots \dd x_N.
\end{equation}
It is nonnegative and symmetric in its arguments. Since $\wp_N$ is symmetric, we have
\begin{equation}\label{e:goe-corr-identity}
\E \sum_{(i_1, \dots, i_n)} \mathcal F(\lambda_{i_1}, \dots, \lambda_{i_n}) = \int_{\R^n} \mathcal F(x_1, \dots, x_n) R_{n, N}(x_1, \dots, x_n) \dd x_1 \cdots \dd x_n
\end{equation}
for every measurable $\mathcal F \colon \R^n \to [0, \infty)$, where the sum runs over $n$-tuples of pairwise distinct indices in $[N]$. By linearity, \eqref{e:goe-corr-identity} also holds for every bounded measurable $\mathcal F$ with compact support.

We recall the definition of the Pfaffian. For an antisymmetric $2n \times 2n$ matrix $\mathsf W$, we set
\begin{equation}\label{e:goe-pf-def}
\operatorname{Pf} \mathsf W = \sum_{\mathfrak m} \sgn(\mathfrak m) \prod_{\{\sigma, \sigma'\} \in \mathfrak m, \, \sigma < \sigma'} \mathsf W_{\sigma \sigma'},
\end{equation}
where the sum runs over the perfect matchings $\mathfrak m$ of $[2n]$, and $\sgn(\mathfrak m)$ is the sign of the permutation $(\sigma_1, \sigma_1', \dots, \sigma_n, \sigma_n')$ of $[2n]$ obtained by listing the pairs $\{\sigma_r, \sigma_r'\}$ of $\mathfrak m$ with $\sigma_r < \sigma_r'$ (this sign does not depend on the order of the pairs). We use three properties. First, $\operatorname{Pf}(B^T \mathsf W B) = \det(B) \operatorname{Pf}(\mathsf W)$ for every $2n \times 2n$ matrix $B$ \cite[Section~2.2]{feng2024small}. Second, let $1 \le l < n$ and set $P_1 = [2l]$ and $P_2 = [2n] \setminus [2l]$. Let $\mathsf W_{P_1 P_1}$ and $\mathsf W_{P_2 P_2}$ denote the corresponding principal submatrices. Then the terms in \eqref{e:goe-pf-def} indexed by the matchings $\mathfrak m$ such that every pair of $\mathfrak m$ is a subset of $P_1$ or of $P_2$ add up to $\operatorname{Pf}(\mathsf W_{P_1 P_1}) \operatorname{Pf}(\mathsf W_{P_2 P_2})$. Third, the number of perfect matchings of $[2n]$ is $(2n - 1)!!$, so that $|\operatorname{Pf} \mathsf W| \le (2n - 1)!! \max_{\sigma, \sigma'} |\mathsf W_{\sigma \sigma'}|^n$.

We now define the matrix kernel of the GOE, following \cite[(7.2)]{feng2024small} up to the names of the functions.

\begin{definition}\label{d:goe-kernel}
Let $N \ge 2$. For $x, y \in \R$, we define
\begin{equation*}
\mathsf K_N(x, y) = \sqrt N \sum_{j = 0}^{N - 1} \mathsf h_j(\sqrt N x) \mathsf h_j(\sqrt N y), \qquad \mathsf P(x, y) = \frac{N}{2} \mathsf h_{N - 1}(\sqrt N x) \mathsf F_N(\sqrt N y),
\end{equation*}
and $\mathsf S = \mathsf K_N + \mathsf P$. We define $\mathsf q(x) = 0$ if $N$ is even and
\begin{equation*}
\mathsf q(x) = \sqrt N \, \mathsf h_{N - 1}(\sqrt N x) \Big( \int_\R \mathsf h_{N - 1} \Big)^{-1}
\end{equation*}
if $N$ is odd, where $\int_\R \mathsf h_{N - 1} > 0$ by the proof of \zcref{l:goe-hermite}(c). We further set $\mathsf Q(x) = (1/2) \int_\R \sgn(x - x') \mathsf q(x') \dd x'$ and
\begin{equation*}
\mathsf T(x, y) = \mathsf S(x, y) + \mathsf q(x), \qquad \mathsf V(x, y) = - \partial_y \mathsf S(x, y),
\end{equation*}
\begin{equation*}
\mathsf I(x, y) = \frac12 \int_\R \sgn(x - x') \mathsf S(x', y) \dd x' - \frac12 \sgn(x - y) + \mathsf Q(x) - \mathsf Q(y).
\end{equation*}
Finally, we define the $2 \times 2$ matrix
\begin{equation*}
\mathsf A(x, y) = \begin{pmatrix} \mathsf I(x, y) & \mathsf T(y, x) \\ - \mathsf T(x, y) & - \mathsf V(x, y) \end{pmatrix}.
\end{equation*}
\end{definition}

We suppress the dependence on $N$ in $\mathsf P$, $\mathsf S$, $\mathsf q$, $\mathsf Q$, $\mathsf T$, $\mathsf V$, $\mathsf I$, and $\mathsf A$. (The function $\mathsf K_N$ is the correlation kernel of the eigenvalues of the Gaussian unitary ensemble normalized so that $\E |G_{ab}|^2 = 1/N$ for all $a, b$, which is the normalization of \zcref{s:hermitian} \cite[Section~3.2]{anderson2010introduction}, and $\mathsf P$ and $\mathsf q$ are the corrections for the orthogonal ensemble.) The functions $\mathsf K_N$, $\mathsf P$, $\mathsf S$, $\mathsf q$, $\mathsf T$, and $\mathsf V$ are smooth, since $\mathsf F_N' = \mathsf h_N$. Observe that $\mathsf S(\cdot, y)$ and $\mathsf q$ are integrable, so $\mathsf I$ and $\mathsf Q$ are well defined, and that $\mathsf Q' = \mathsf q$.

\begin{lemma}\label{l:goe-pfaffian}
Let $N \ge 2$ and $n \in [N]$. For all pairwise distinct $x_1, \dots, x_n \in \R$,
\begin{equation}\label{e:goe-pf}
R_{n, N}(x_1, \dots, x_n) = \operatorname{Pf} \big[ \mathsf A(x_i, x_j) \big]_{i, j = 1}^n,
\end{equation}
where the $2n \times 2n$ matrix on the right has the $2 \times 2$ block $\mathsf A(x_i, x_j)$ in block position $(i, j)$.
\end{lemma}

This is the classical Pfaffian, or quaternion determinant, formula of Mehta for the correlation functions of the GOE \cite{mehta2004random}. In the present normalization, it is \cite[(7.1)--(7.2)]{feng2024small}. In the notation of \cite{feng2024small}, the matrix size is $n$ instead of $N$, the functions $\mathsf T$, $\mathsf V$, $\mathsf I$, $\mathsf q$, and $\mathsf Q$ are denoted $S^{(1)}_n + \alpha_n$, $V^{(1)}_n$, $J^{(1)}_n$, $\alpha_n$, and $\gamma_n$, and $\mathsf A(x, y)$ is the matrix $J K^{(1)}_n(x, y)$; these symbols are used only in this paragraph. For even $N$, the matrix kernel is derived by Tracy and Widom \cite[(9.1) and Section~10]{tracy1998correlation}, and for odd $N$ it is given by Adler, Forrester, Nagao, and van Moerbeke \cite[(2.11) and (4.3)]{adler2000classical}. These references use the weight $e^{-x^2/2}$, which corresponds to the eigenvalues of $G$ multiplied by $(N/2)^{1/2}$. By \zcref{l:goe-identities}(a) below, the matrix on the right side of \eqref{e:goe-pf} is antisymmetric.

We will use the following identities for the kernel.

\begin{lemma}\label{l:goe-identities}
Let $N \ge 2$ and $x, y \in \R$.
\begin{enumerate}[label=(\alph*)]
\item We have $\mathsf V(x, y) = - \mathsf V(y, x)$ and $\mathsf I(x, y) = - \mathsf I(y, x)$. In particular, $\mathsf A(x, y)^T = - \mathsf A(y, x)$.
\item We have
\begin{equation}\label{e:goe-I-formula}
\mathsf I(x, y) = \int_y^x \mathsf T(x', y) \dd x' - \frac12 \sgn(x - y).
\end{equation}
Consequently, $\partial_x \mathsf I(x, y) = \mathsf T(x, y)$ and $\partial_y \mathsf I(x, y) = - \mathsf T(y, x)$ for $x \neq y$.
\item We have $R_{1, N}(x) = \mathsf T(x, x)$, and if $x \neq y$, then
\begin{equation}\label{e:goe-R2}
R_{2, N}(x, y) = \mathsf T(x, x) \mathsf T(y, y) - \mathsf T(x, y) \mathsf T(y, x) + \mathsf I(x, y) \mathsf V(x, y).
\end{equation}
\end{enumerate}
\end{lemma}

\begin{proof}
Throughout the proof, we write $w = \sqrt N x$ and $w' = \sqrt N y$. We begin with $\mathsf V$. The second identity in \eqref{e:goe-hermite-id} gives
\begin{equation*}
\sum_{j = 0}^{N - 1} \mathsf h_j(w) \mathsf h_j'(w') = \frac12 \sum_{j = 1}^{N - 1} \sqrt j \big( \mathsf h_j(w) \mathsf h_{j - 1}(w') - \mathsf h_{j - 1}(w) \mathsf h_j(w') \big) - \frac{\sqrt N}{2} \mathsf h_{N - 1}(w) \mathsf h_N(w').
\end{equation*}
Using $\mathsf F_N' = \mathsf h_N$, we have
\begin{equation*}
\partial_y \mathsf P(x, y) = \frac{N^{3/2}}{2} \mathsf h_{N - 1}(w) \mathsf h_N(w'),
\end{equation*}
and therefore
\begin{equation}\label{e:goe-V-antisym}
\begin{split}
\mathsf V(x, y) &= - N \Big( \sum_{j = 0}^{N - 1} \mathsf h_j(w) \mathsf h_j'(w') + \frac{\sqrt N}{2} \mathsf h_{N - 1}(w) \mathsf h_N(w') \Big) \\
&= - \frac N2 \sum_{j = 1}^{N - 1} \sqrt j \big( \mathsf h_j(w) \mathsf h_{j - 1}(w') - \mathsf h_{j - 1}(w) \mathsf h_j(w') \big),
\end{split}
\end{equation}
which is antisymmetric in $(x, y)$.

Next, we consider $\mathsf I$. Substituting $\sqrt N x'$ as the variable of integration, we find
\begin{equation*}
\frac12 \int_\R \sgn(x - x') \sqrt N \mathsf h_j(\sqrt N x') \dd x' = \mathsf F_j(w),
\end{equation*}
and therefore
\begin{equation*}
\frac12 \int_\R \sgn(x - x') \mathsf S(x', y) \dd x' = \sum_{j = 0}^{N - 1} \mathsf F_j(w) \mathsf h_j(w') + \frac{\sqrt N}{2} \mathsf F_{N - 1}(w) \mathsf F_N(w').
\end{equation*}
Inserting \eqref{e:goe-F-id} for $\mathsf h_j(w')$ and arguing as for \eqref{e:goe-V-antisym}, we find
\begin{equation}\label{e:goe-I-antisym}
\frac12 \int_\R \sgn(x - x') \mathsf S(x', y) \dd x' = \frac12 \sum_{j = 1}^{N - 1} \sqrt j \big( \mathsf F_j(w) \mathsf F_{j - 1}(w') - \mathsf F_{j - 1}(w) \mathsf F_j(w') \big).
\end{equation}
Let $\mathsf I_0(x, y)$ denote the right side of \eqref{e:goe-I-antisym} plus $\mathsf Q(x) - \mathsf Q(y)$, so that
\begin{equation*}
\mathsf I(x, y) = \mathsf I_0(x, y) - \frac12 \sgn(x - y).
\end{equation*}
The function $\mathsf I_0$ is continuous and antisymmetric, and the function $(x, y) \mapsto \sgn(x - y)$ is antisymmetric. Therefore $\mathsf I$ is antisymmetric, which completes the proof of (a), since the identity $\mathsf A(x, y)^T = - \mathsf A(y, x)$ follows from the antisymmetry of $\mathsf I$ and $\mathsf V$.

We now prove (b). Since $\mathsf I_0$ is continuous and antisymmetric, we have $\mathsf I_0(y, y) = 0$. Further, the fundamental theorem of calculus and $\mathsf Q' = \mathsf q$ give
\begin{equation*}
\partial_x \mathsf I_0(x, y) = \mathsf S(x, y) + \mathsf q(x) = \mathsf T(x, y), \qquad x \in \R.
\end{equation*}
Hence
\begin{equation*}
\mathsf I_0(x, y) = \int_y^x \mathsf T(x', y) \dd x',
\end{equation*}
which is \eqref{e:goe-I-formula}. The formula for $\partial_x \mathsf I$ follows at once, and the formula for $\partial_y \mathsf I$ follows from $\mathsf I(x, y) = - \mathsf I(y, x)$.

It remains to prove (c). By (a), we have $\mathsf I(x, x) = 0$ and $\mathsf V(x, x) = 0$. Hence $\mathsf A(x, x)$ has the entries $\mathsf T(x, x)$ and $- \mathsf T(x, x)$ off the diagonal and zeros on the diagonal, and \eqref{e:goe-pf} with $n = 1$ gives $R_{1, N}(x) = \mathsf T(x, x)$. For $n = 2$, the definition of $\mathsf A$ and (a) show that the matrix $\mathsf W$ in \eqref{e:goe-pf} is
\begin{equation*}
\mathsf W = \begin{pmatrix} 0 & \mathsf T(x, x) & \mathsf I(x, y) & \mathsf T(y, x) \\ - \mathsf T(x, x) & 0 & - \mathsf T(x, y) & - \mathsf V(x, y) \\ - \mathsf I(x, y) & \mathsf T(x, y) & 0 & \mathsf T(y, y) \\ - \mathsf T(y, x) & \mathsf V(x, y) & - \mathsf T(y, y) & 0 \end{pmatrix}.
\end{equation*}
Therefore
\begin{equation*}
R_{2, N}(x, y) = \mathsf W_{12} \mathsf W_{34} - \mathsf W_{13} \mathsf W_{24} + \mathsf W_{14} \mathsf W_{23},
\end{equation*}
which is \eqref{e:goe-R2}.
\end{proof}

\section{Kernel estimates}\label{s:goe-kernel}

For functions of two real variables, we write $\partial_1$ and $\partial_2$ for the partial derivatives in the first and second variables. We first bound the kernel $\mathsf K_N$ and its derivatives in the bulk.

\begin{lemma}\label{l:goe-gue}
Let $\delta \in (0, 1)$. There exist $C > 0$ and $N_0$ such that the following statements hold for all $N \ge N_0$ and $x, y \in I_\delta$.
\begin{enumerate}[label=(\alph*)]
\item For all integers $k_1, k_2 \ge 0$ with $k_1 + k_2 \le 3$, we have $|\partial_x^{k_1} \partial_y^{k_2} \mathsf K_N(x, y)| \le C N^{1 + k_1 + k_2}$.
\item If $x \neq y$, then $|\mathsf K_N(x, y)| \le C |x - y|^{-1}$ and $|\partial_y \mathsf K_N(x, y)| \le C (N |x - y|^{-1} + |x - y|^{-2})$.
\item The derivative of the function $x \mapsto \mathsf K_N(x, x)$ is bounded by $C N$ on $I_\delta$.
\item We have
\begin{equation}\label{e:goe-K-second}
\partial_y^2 \mathsf K_N(x, y) \big|_{y = x} = - \frac N3 \Big( N \Big( 1 - \frac{x^2}{4} \Big) + \frac12 \Big) \mathsf K_N(x, x) - \frac{N^2}{3} \mathsf h_N(\sqrt N x) \mathsf h_{N - 1}'(\sqrt N x).
\end{equation}
\end{enumerate}
\end{lemma}

\begin{proof}
Let $w = \sqrt N x$ and $w' = \sqrt N y$, and define
\begin{equation*}
\mathsf K^\circ(w, w') = \sum_{j = 0}^{N - 1} \mathsf h_j(w) \mathsf h_j(w'), \qquad \text{so that} \quad \mathsf K_N(x, y) = \sqrt N \mathsf K^\circ(w, w').
\end{equation*}
We have $|w|, |w'| \le (2 - \delta) \sqrt N$, and there exists $N_0$ such that $(2 - \delta) \sqrt N \le (2 - \delta/2) \sqrt n$ for all $N \ge N_0$ and $n \in \{N - 1, N\}$. Hence \zcref{l:goe-hermite}(a) with $\delta/2$ in place of $\delta$ gives
\begin{equation}\label{e:goe-hermite-N}
|\mathsf h_n^{(k)}(w)| + |\mathsf h_n^{(k)}(w')| \le C N^{k/2 - 1/4}, \qquad n \in \{N - 1, N\}, \ k \in \{0, 1, 2, 3\}.
\end{equation}

We begin with (a). Using \eqref{e:goe-cd-diag} with $n = N$ and \eqref{e:goe-hermite-N}, we obtain
\begin{equation*}
\sum_{j = 0}^{N - 1} \mathsf h_j(w)^2 = \sqrt N \big( \mathsf h_N'(w) \mathsf h_{N - 1}(w) - \mathsf h_{N - 1}'(w) \mathsf h_N(w) \big) \le C \sqrt N.
\end{equation*}
Squaring the second identity in \eqref{e:goe-hermite-id} and using the Cauchy--Schwarz inequality, we find
\begin{equation*}
\sum_{j = 0}^{N - 1} \mathsf h_j'(w)^2 \le \frac12 \sum_{j = 0}^{N - 1} \big( j \mathsf h_{j - 1}(w)^2 + (j + 1) \mathsf h_{j + 1}(w)^2 \big) \le N \sum_{j = 0}^{N} \mathsf h_j(w)^2 \le C N^{3/2},
\end{equation*}
where in the last step we used \eqref{e:goe-hermite-N} for $\mathsf h_N(w)$. Next, for $j \le N - 1$, \eqref{e:goe-hermite-id} and \eqref{e:goe-hermite-third} give
\begin{align*}
|\mathsf h_j''(w)| &= |\mathsf e_j(w)| |\mathsf h_j(w)| \le 2N |\mathsf h_j(w)|, \\
|\mathsf h_j'''(w)| &\le \frac{|w|}{2} |\mathsf h_j(w)| + 2N |\mathsf h_j'(w)| \le \sqrt N |\mathsf h_j(w)| + 2N |\mathsf h_j'(w)|.
\end{align*}
Combining the last three bounds, we find
\begin{equation}\label{e:goe-hermite-sums}
\sum_{j = 0}^{N - 1} \mathsf h_j^{(k)}(w)^2 \le C N^{k + 1/2}, \qquad k \in \{0, 1, 2, 3\},
\end{equation}
and the same bound holds at $w'$. The Cauchy--Schwarz inequality and \eqref{e:goe-hermite-sums} now give
\begin{align*}
|\partial_x^{k_1} \partial_y^{k_2} \mathsf K_N(x, y)| &= N^{(1 + k_1 + k_2)/2} \Big| \sum_{j = 0}^{N - 1} \mathsf h_j^{(k_1)}(w) \mathsf h_j^{(k_2)}(w') \Big| \\
&\le N^{(1 + k_1 + k_2)/2} \big( C N^{k_1 + 1/2} \big)^{1/2} \big( C N^{k_2 + 1/2} \big)^{1/2} = C N^{1 + k_1 + k_2}.
\end{align*}

To prove (b), we write \eqref{e:goe-cd} in the form
\begin{equation*}
\mathsf K_N(x, y) = \sqrt N \, \frac{\mathsf Z(w, w')}{x - y}, \qquad \mathsf Z(w, w') = \mathsf h_N(w) \mathsf h_{N - 1}(w') - \mathsf h_{N - 1}(w) \mathsf h_N(w').
\end{equation*}
By \eqref{e:goe-hermite-N}, we have
\begin{equation*}
|\mathsf Z(w, w')| \le C N^{-1/2}, \qquad |\partial_{w'} \mathsf Z(w, w')| \le C.
\end{equation*}
The bounds in (b) then follow from the last two displays and the identity
\begin{equation*}
\partial_y \mathsf K_N(x, y) = \frac{N \partial_{w'} \mathsf Z(w, w')}{x - y} + \frac{\sqrt N \mathsf Z(w, w')}{(x - y)^2}.
\end{equation*}

Next, we prove (c). The second identity in \eqref{e:goe-hermite-id} and a telescoping sum give
\begin{equation*}
\frac{\dd}{\dd x} \mathsf K_N(x, x) = 2N \sum_{j = 0}^{N - 1} \mathsf h_j(w) \mathsf h_j'(w) = - N^{3/2} \mathsf h_{N - 1}(w) \mathsf h_N(w),
\end{equation*}
and the right side is bounded by $C N$ by \eqref{e:goe-hermite-N}.

It remains to prove (d), which holds for all $x \in \R$. We fix $w$ and set $\mathsf z(w') = \mathsf Z(w, w')$, so that $\mathsf z$ is entire and $\mathsf z(w) = 0$. By \eqref{e:goe-cd} and the Taylor expansion of $\mathsf z$ at $w$, we have
\begin{equation*}
\mathsf K^\circ(w, w') = - \sqrt N \, \frac{\mathsf z(w')}{w' - w} = - \sqrt N \sum_{l = 1}^\infty \frac{\mathsf z^{(l)}(w)}{l!} (w' - w)^{l - 1}, \qquad w' \neq w.
\end{equation*}
Hence
\begin{equation*}
\mathsf K^\circ(w, w) = - \sqrt N \mathsf z'(w), \qquad \partial_{w'}^2 \mathsf K^\circ(w, w') \big|_{w' = w} = - \frac{\sqrt N}{3} \mathsf z'''(w).
\end{equation*}

To compute $\mathsf z'''$, we use \eqref{e:goe-hermite-third} with $n \in \{N - 1, N\}$, together with $\mathsf e_{N - 1} = \mathsf e_N - 1$. We obtain
\begin{align*}
\mathsf z'''(w) &= - \mathsf e_{N - 1}(w) \mathsf h_N(w) \mathsf h_{N - 1}'(w) + \mathsf e_N(w) \mathsf h_{N - 1}(w) \mathsf h_N'(w) \\
&= - \mathsf e_N(w) \mathsf z'(w) + \mathsf h_N(w) \mathsf h_{N - 1}'(w).
\end{align*}
Combining the last two displays, we obtain
\begin{equation*}
\partial_{w'}^2 \mathsf K^\circ(w, w') \big|_{w' = w} = - \frac{\mathsf e_N(w)}{3} \mathsf K^\circ(w, w) - \frac{\sqrt N}{3} \mathsf h_N(w) \mathsf h_{N - 1}'(w).
\end{equation*}

Finally, inserting the relations
\begin{equation*}
\partial_y^2 \mathsf K_N(x, y) = N^{3/2} \partial_{w'}^2 \mathsf K^\circ(w, w'), \quad \mathsf K^\circ(w, w) = \frac{\mathsf K_N(x, x)}{\sqrt N}, \quad \mathsf e_N(\sqrt N x) = N \Big( 1 - \frac{x^2}{4} \Big) + \frac12
\end{equation*}
into the preceding formula for $\partial_{w'}^2 \mathsf K^\circ(w, w') |_{w' = w}$, we obtain \eqref{e:goe-K-second}.
\end{proof}

We now turn to the functions $\mathsf T$, $\mathsf V$, and $\mathsf I$ of \zcref{d:goe-kernel}.

\begin{lemma}\label{l:goe-kernel}
Let $\delta \in (0, 1)$. There exist $C > 0$ and $N_0$ such that the following statements hold for all $N \ge N_0$ and $x, y \in I_\delta$.
\begin{enumerate}[label=(\alph*)]
\item For all integers $k_1, k_2 \ge 0$ with $k_1 + k_2 \le 2$, we have $|\partial_x^{k_1} \partial_y^{k_2} \mathsf T(x, y)| \le C N^{1 + k_1 + k_2}$ and $|\partial_x^{k_1} \partial_y^{k_2} \mathsf V(x, y)| \le C N^{2 + k_1 + k_2}$.
\item If $x \neq y$, then $|\mathsf T(x, y)| \le C (|x - y|^{-1} + N^{1/2})$ and $|\mathsf V(x, y)| \le C (N |x - y|^{-1} + |x - y|^{-2} + N)$.
\item We have $|\mathsf I(x, y)| \le C (1 + \log(1 + N |x - y|))$.
\item We have $|\mathsf T(x, x) - \mathsf K_N(x, x)| \le C N^{1/2}$, and the derivative of the function $x \mapsto \mathsf T(x, x)$ is bounded by $C N^{3/2}$ on $I_\delta$.
\item We have
\begin{equation}\label{e:goe-slope-raw}
\Big| \partial_2 \mathsf V(x, x) - \frac N3 \Big( N \Big( 1 - \frac{x^2}{4} \Big) + \frac12 \Big) \mathsf K_N(x, x) \Big| \le C N^2.
\end{equation}
\end{enumerate}
\end{lemma}

\begin{proof}
We write $w = \sqrt N x$ and $w' = \sqrt N y$. Let $N_0$ be as in the proof of \zcref{l:goe-gue}, so that \eqref{e:goe-hermite-N} holds. Since $|w'| \le (2 - \delta/2) \sqrt N$, \zcref{l:goe-hermite}(d) and (c), with $\delta/2$ in place of $\delta$, give
\begin{equation*}
|\mathsf F_N(w')| \le C N^{-1/4}, \qquad \int_\R \mathsf h_{N - 1} \ge c N^{-1/4} \quad \text{if } N \text{ is odd}.
\end{equation*}

Recall that $\mathsf T = \mathsf K_N + \mathsf P + \mathsf q$ and $\mathsf V = - \partial_y \mathsf K_N - \partial_y \mathsf P$. We first bound the corrections $\mathsf P$ and $\mathsf q$. For integers $k_1 \ge 0$ and $k_2 \ge 1$ with $k_1 + k_2 \le 3$, we have
\begin{align*}
\partial_x^{k_1} \mathsf P(x, y) &= \frac{N^{1 + k_1/2}}{2} \mathsf h_{N - 1}^{(k_1)}(w) \mathsf F_N(w'), \\
\partial_x^{k_1} \partial_y^{k_2} \mathsf P(x, y) &= \frac{N^{1 + (k_1 + k_2)/2}}{2} \mathsf h_{N - 1}^{(k_1)}(w) \mathsf h_N^{(k_2 - 1)}(w'),
\end{align*}
and, for odd $N$,
\begin{equation*}
\mathsf q^{(k_1)}(x) = N^{(1 + k_1)/2} \mathsf h_{N - 1}^{(k_1)}(w) \Big( \int_\R \mathsf h_{N - 1} \Big)^{-1}.
\end{equation*}
Therefore
\begin{equation}\label{e:goe-corrections}
|\partial_x^{k_1} \mathsf P(x, y)| + |\mathsf q^{(k_1)}(x)| \le C N^{k_1 + 1/2}, \qquad |\partial_x^{k_1} \partial_y^{k_2} \mathsf P(x, y)| \le C N^{k_1 + k_2},
\end{equation}
where we used \eqref{e:goe-hermite-N} and the bounds on $\mathsf F_N(w')$ and $\int_\R \mathsf h_{N - 1}$. The bounds in (a) follow from \eqref{e:goe-corrections} and \zcref{l:goe-gue}(a), and the bounds in (b) follow from \eqref{e:goe-corrections} with $k_1 + k_2 \le 1$ and \zcref{l:goe-gue}(b).

To prove (c), we use \eqref{e:goe-I-formula} and bound the contributions of $\mathsf K_N$, $\mathsf P$, and $\mathsf q$ to the integral there separately. By \zcref{l:goe-gue}(a) and (b), we have $|\mathsf K_N(x', y)| \le C \min(N, |x' - y|^{-1})$ for $x' \in I_\delta$, so
\begin{equation*}
\Big| \int_y^x \mathsf K_N(x', y) \dd x' \Big| \le C \int_0^{|x - y|} \min(N, v^{-1}) \dd v \le C (1 + \log(1 + N |x - y|)).
\end{equation*}
Next, substituting $\sqrt N x'$ as the variable of integration, we find
\begin{equation*}
\int_y^x \mathsf P(x', y) \dd x' = \frac{\sqrt N}{2} \mathsf F_N(w') \int_{w'}^{w} \mathsf h_{N - 1}, \qquad \int_y^x \mathsf q(x') \dd x' = \Big( \int_\R \mathsf h_{N - 1} \Big)^{-1} \int_{w'}^w \mathsf h_{N - 1},
\end{equation*}
where the second identity is for odd $N$. Since $|w|, |w'| \le (2 - \delta/2) \sqrt{N - 1}$, \zcref{l:goe-hermite}(b) with $\delta/2$ in place of $\delta$ gives
\begin{equation*}
\Big| \int_{w'}^w \mathsf h_{N - 1} \Big| \le C N^{-3/4},
\end{equation*}
and therefore the integrals of $\mathsf P(\cdot, y)$ and $\mathsf q$ from $y$ to $x$ are bounded in absolute value by $C N^{-1/2}$. Together with \eqref{e:goe-I-formula}, this completes the proof of (c).

We now prove (d). We have
\begin{align*}
\mathsf T(x, x) - \mathsf K_N(x, x) &= \mathsf P(x, x) + \mathsf q(x), \\
\frac{\dd}{\dd x} \mathsf T(x, x) &= \frac{\dd}{\dd x} \mathsf K_N(x, x) + \partial_1 \mathsf P(x, x) + \partial_2 \mathsf P(x, x) + \mathsf q'(x),
\end{align*}
and \eqref{e:goe-corrections} and \zcref{l:goe-gue}(c) bound the right sides by $C N^{1/2}$ and $C N^{3/2}$, respectively.

Finally, we prove (e). We have
\begin{equation*}
\partial_2 \mathsf V(x, x) = - \partial_y^2 \mathsf K_N(x, y) \big|_{y = x} - \partial_y^2 \mathsf P(x, y) \big|_{y = x},
\end{equation*}
and the second term is bounded by $C N^2$ by \eqref{e:goe-corrections}. By \eqref{e:goe-K-second}, the first term differs from the main term in \eqref{e:goe-slope-raw} by $N^2 |\mathsf h_N(w) \mathsf h_{N - 1}'(w)|/3$, which is bounded by $C N^2$ by \eqref{e:goe-hermite-N}.
\end{proof}

We now deduce the density of states from eigenvalue rigidity, using the bound on the derivative of $x \mapsto \mathsf T(x, x)$ in \zcref{l:goe-kernel}(d) to pass from intervals to points. Recall that $c(E) = (\pi^2/6) \rhosc(E)^3$, and note that $1 - x^2/4 = \pi^2 \rhosc(x)^2$ for $|x| \le 2$.

\begin{lemma}\label{l:goe-density}
Let $\delta \in (0, 1)$. There exist $C > 0$ and $N_0$ such that for all $N \ge N_0$ and $x \in I_\delta$,
\begin{equation}\label{e:goe-density-bounds}
\begin{gathered}
|R_{1, N}(x) - N \rhosc(x)| + |\mathsf K_N(x, x) - N \rhosc(x)| \le C N^{7/8}, \\
\Big| \frac12 \partial_2 \mathsf V(x, x) - N^3 c(x) \Big| \le C N^{23/8}.
\end{gathered}
\end{equation}
\end{lemma}

\begin{proof}
Set $\eps = 1/8$, and choose $c_* \in (0, 1/4)$ such that
\begin{equation*}
\int_{-2}^{-2 + \delta/4} \rhosc > 2 c_*.
\end{equation*}
For $N \ge c_*^{-1}$, the definition of the classical locations and the symmetry $\gamma_{N - i} = - \gamma_i$ then give
\begin{equation*}
\gamma_i < -2 + \frac\delta4 \quad (i \le c_* N + 1), \qquad \gamma_i > 2 - \frac\delta4 \quad (i \ge (1 - c_*) N - 1).
\end{equation*}
Let $\mathcal E_N$ be the event that $|\lambda_i - \gamma_i| \le N^{-1 + \eps}$ for all $i$ such that $c_* N \le i \le (1 - c_*) N$. By \zcref{p:regularity}(b) for $G$, we have $\Prob(\mathcal E_N^c) \le N^{-2}$ for $N \ge N_0$.

We first show that for all $y_1 < y_2$ such that $[y_1, y_2] \subset I_{\delta/2}$,
\begin{equation}\label{e:goe-counting}
\Big| \E \# \{ i : \lambda_i \in [y_1, y_2] \} - N \int_{y_1}^{y_2} \rhosc \Big| \le C N^{\eps}.
\end{equation}
Suppose that $\mathcal E_N$ holds. Then, for $i < c_* N$, we have
\begin{equation*}
\lambda_i \le \lambda_{\lceil c_* N \rceil} \le \gamma_{\lceil c_* N \rceil} + N^{-1 + \eps} < y_1,
\end{equation*}
provided $N \ge N_0$, and similarly $\lambda_i > y_2$ for $i > (1 - c_*) N$.

For the remaining indices, we have the implications
\begin{equation*}
\begin{gathered}
\lambda_i \in [y_1, y_2] \implies \gamma_i \in [y_1 - N^{-1 + \eps}, y_2 + N^{-1 + \eps}], \\
\gamma_i \in [y_1 + N^{-1 + \eps}, y_2 - N^{-1 + \eps}] \implies \lambda_i \in [y_1, y_2].
\end{gathered}
\end{equation*}
Since $\rhosc \le 1/\pi$, consecutive classical locations are at distance at least $\pi/N$. Hence each of the intervals $[y_1 - N^{-1 + \eps}, y_1 + N^{-1 + \eps}]$ and $[y_2 - N^{-1 + \eps}, y_2 + N^{-1 + \eps}]$ contains at most $C N^{\eps}$ classical locations, and $[y_1, y_2]$ contains $N \int_{y_1}^{y_2} \rhosc + O(1)$ of them. Therefore, on $\mathcal E_N$,
\begin{equation*}
\# \{ i : \lambda_i \in [y_1, y_2] \} = N \int_{y_1}^{y_2} \rhosc + O(N^{\eps}).
\end{equation*}
On $\mathcal E_N^c$, the left side is bounded by $N$, and $N \Prob(\mathcal E_N^c) \le N^{-1}$. Taking expectations, we obtain \eqref{e:goe-counting}.

Now let $x \in I_\delta$, so that $[x, x + N^{-3/4}] \subset I_{\delta/2}$. Using \eqref{e:goe-corr-identity} with $n = 1$, we have
\begin{equation*}
\E \# \{ i : \lambda_i \in [x, x + N^{-3/4}] \} = \int_x^{x + N^{-3/4}} R_{1, N}(x') \dd x'.
\end{equation*}
Note that $R_{1, N}(x') = \mathsf T(x', x')$ by \zcref{l:goe-identities}(c), so \zcref{l:goe-kernel}(d) on $I_{\delta/2}$ gives
\begin{equation*}
|R_{1, N}(x') - R_{1, N}(x)| \le C N^{3/2} N^{-3/4} = C N^{3/4}, \qquad x' \in [x, x + N^{-3/4}].
\end{equation*}
Using also \eqref{e:goe-counting} and the fact that $\rhosc$ is Lipschitz on $I_{\delta/2}$, we find
\begin{align*}
|R_{1, N}(x) - N \rhosc(x)| &\le N^{3/4} \Big| \int_x^{x + N^{-3/4}} R_{1, N} - N \int_x^{x + N^{-3/4}} \rhosc \Big| + C N^{3/4} + C N^{1/4} \\
&\le C \big( N^{\eps + 3/4} + N^{3/4} \big) \le C N^{7/8}.
\end{align*}
By \zcref{l:goe-kernel}(d), we have $|\mathsf K_N(x, x) - R_{1, N}(x)| \le C N^{1/2}$, which gives the first bound in \eqref{e:goe-density-bounds}. Inserting it into \eqref{e:goe-slope-raw} and using $1 - x^2/4 = \pi^2 \rhosc(x)^2$, we obtain
\begin{equation*}
\frac12 \partial_2 \mathsf V(x, x) = \Big( \frac{\pi^2}{6} N^2 \rhosc(x)^2 + \frac N{12} \Big) \big( N \rhosc(x) + O(N^{7/8}) \big) + O(N^2) = N^3 c(x) + O(N^{23/8}),
\end{equation*}
which is the second bound in \eqref{e:goe-density-bounds}.
\end{proof}

\section{Correlation functions}\label{s:goe-correlations}

We now estimate the correlation functions. In the next lemma, part (a) determines the mean of the statistic in \eqref{e:prelim-goepairs}, and the cruder bound in part (b) will be used for the tails.

\begin{lemma}\label{l:goe-two}
Let $\delta \in (0, 1)$. There exist $C > 0$ and $N_0$ such that the following statements hold for all $N \ge N_0$.
\begin{enumerate}[label=(\alph*)]
\item For all $x, y \in I_\delta$ such that $0 < y - x \le N^{-1}$,
\begin{equation}\label{e:goe-two-expansion}
\begin{gathered}
\big| R_{2, N}(x, y) - N^3 c(x) (y - x) \big| \le C \big( N^{23/8} (y - x) + N^4 (y - x)^2 \big), \\
|R_{2, N}(x, y)| \le C N^3 (y - x).
\end{gathered}
\end{equation}
\item For all $x, y \in I_\delta$ such that $|x - y| \ge N^{-1}$, we have $|R_{2, N}(x, y)| \le C N^2 \log N$.
\end{enumerate}
\end{lemma}

\begin{proof}
We begin with (a). We set $y_v = x + v$ for $v \in (0, y - x]$, and we define
\begin{equation*}
\mathsf r(v) = R_{2, N}(x, y_v) = \mathsf T(x, x) \mathsf T(y_v, y_v) - \mathsf T(x, y_v) \mathsf T(y_v, x) + \mathsf I(x, y_v) \mathsf V(x, y_v),
\end{equation*}
where the second equality is \eqref{e:goe-R2}.
By \zcref{l:goe-identities}(b) and the antisymmetry of $\mathsf I$, we have
\begin{equation*}
\mathsf I(x, y') = \frac12 - \int_x^{y'} \mathsf T(x', x) \dd x', \qquad \partial_2 \mathsf I(x, y') = - \mathsf T(y', x), \qquad y' > x,
\end{equation*}
so $\mathsf r$ is smooth on $(0, y - x]$. Differentiating, we find
\begin{align*}
\mathsf r'(v) &= \mathsf T(x, x) \big( \partial_1 \mathsf T + \partial_2 \mathsf T \big)(y_v, y_v) - \partial_2 \mathsf T(x, y_v) \mathsf T(y_v, x) - \mathsf T(x, y_v) \partial_1 \mathsf T(y_v, x) \\
&\quad - \mathsf T(y_v, x) \mathsf V(x, y_v) + \mathsf I(x, y_v) \partial_2 \mathsf V(x, y_v).
\end{align*}

We now let $v \downarrow 0$. Then $\mathsf I(x, y_v) \to 1/2$, and $\mathsf V(x, x) = 0$ by \zcref{l:goe-identities}(a). Since $\mathsf T$ and $\mathsf V$ are smooth, the first three terms in $\mathsf r'(v)$ cancel in the limit, and the fourth term tends to $- \mathsf T(x, x) \mathsf V(x, x) = 0$. Hence $\mathsf r(v) \to 0$ and
\begin{equation}\label{e:goe-slope}
\lim_{v \downarrow 0} \mathsf r'(v) = \frac12 \partial_2 \mathsf V(x, x).
\end{equation}

Next, we bound $\mathsf r'$ and $\mathsf r''$. Since $v \le N^{-1}$, \zcref{l:goe-kernel}(c) gives $|\mathsf I(x, y_v)| \le C$. The derivatives of $\mathsf T$ and $\mathsf V$ of order at most two are bounded on $I_\delta^2$ by \zcref{l:goe-kernel}(a), and $\partial_2^2 \mathsf I(x, y') = - \partial_1 \mathsf T(y', x)$ for $y' > x$. Differentiating $\mathsf r'$ once more and inserting these bounds, we obtain
\begin{equation*}
|\mathsf r'(v)| \le C N^3, \qquad |\mathsf r''(v)| \le C N^4, \qquad v \in (0, y - x].
\end{equation*}
For example, the last term of $\mathsf r'$ contributes $|\mathsf I \, \partial_2^2 \mathsf V| + |\partial_2 \mathsf I \, \partial_2 \mathsf V| \le C N^4$ to $\mathsf r''$.

Since $\mathsf r(v) \to 0$ as $v \downarrow 0$, we have
\begin{equation*}
R_{2, N}(x, y) = \int_0^{y - x} \mathsf r'(v) \dd v = \frac{y - x}{2} \partial_2 \mathsf V(x, x) + \int_0^{y - x} \Big( \mathsf r'(v) - \frac12 \partial_2 \mathsf V(x, x) \Big) \dd v,
\end{equation*}
where the last integral is bounded by $C N^4 (y - x)^2$ by \eqref{e:goe-slope} and the bound on $\mathsf r''$. Combining this with \zcref{l:goe-density}, we obtain the first bound in \eqref{e:goe-two-expansion}. The second bound follows from the first identity in the last display and $|\mathsf r'| \le C N^3$.

We now prove (b), using \eqref{e:goe-R2} and \zcref{l:goe-kernel}. Since $|x - y| \ge N^{-1}$, parts (a), (b), and (c) of that lemma give $|\mathsf T| \le C N$ on $I_\delta^2$ and
\begin{equation*}
|\mathsf V(x, y)| \le C (N^2 + N^2 + N) \le C N^2, \qquad |\mathsf I(x, y)| \le C (1 + \log(1 + 4N)) \le C \log N,
\end{equation*}
where we also used $|x - y| \le 4$. Inserting these bounds into \eqref{e:goe-R2}, we obtain $|R_{2, N}(x, y)| \le C N^2 \log N$.
\end{proof}

To bound the higher correlation functions, we normalize the entries of the matrix in \eqref{e:goe-pf}. Let $\mathsf D_N = \diag(1, N^{-1})$. We let $B_n$ denote the block diagonal matrix with $n$ diagonal blocks equal to $\mathsf D_N$. Then $\det B_n = N^{-n}$, and the first property of the Pfaffian gives
\begin{equation}\label{e:goe-pf-normalized}
R_{n, N}(x_1, \dots, x_n) = N^n \operatorname{Pf} \big[ \mathsf D_N \mathsf A(x_i, x_j) \mathsf D_N \big]_{i, j = 1}^n,
\end{equation}
where
\begin{equation*}
\mathsf D_N \mathsf A(x, y) \mathsf D_N = \begin{pmatrix} \mathsf I(x, y) & N^{-1} \mathsf T(y, x) \\ - N^{-1} \mathsf T(x, y) & - N^{-2} \mathsf V(x, y) \end{pmatrix}.
\end{equation*}
For the matrix in \eqref{e:goe-pf-normalized}, we call row $2i - 1$ the first row of the point $x_i$ and row $2i$ its second row, and similarly for columns. By \zcref{l:goe-kernel}(a), all entries of $\mathsf D_N \mathsf A(x, y) \mathsf D_N$ other than $\mathsf I(x, y)$ are bounded by a constant for $x, y \in I_\delta$. For $\varrho > 0$, we set
\begin{equation}\label{e:goe-Delta}
\Delta_0(\varrho) = (N \varrho)^{-1} + N^{-1/2}.
\end{equation}

\begin{lemma}\label{l:goe-corr}
Let $\delta \in (0, 1)$. There exist $C > 0$ and $N_0$ such that the following statements hold for all $N \ge N_0$.
\begin{enumerate}[label=(\alph*)]
\item For all pairwise distinct $x_1, x_2, x_3 \in I_\delta$ such that $|x_i - x_j| \le N^{-1}$ for all $i, j$, we have $|R_{3, N}(x_1, x_2, x_3)| \le C N^3$.
\item For all pairwise distinct $x_1, x_2, x_3, x_4 \in I_\delta$, we have $|R_{4, N}(x_1, x_2, x_3, x_4)| \le C N^4 (\log N)^2$.
\item Let $x_1, x_2, x_3, x_4 \in I_\delta$ with $x_1 < x_2 \le x_1 + N^{-1}$ and $x_3 < x_4 \le x_3 + N^{-1}$, and set $\varrho = |x_3 - x_1|$. Suppose that $\varrho \ge \max(4(x_2 - x_1 + x_4 - x_3), N^{-1})$. Then
\begin{equation}\label{e:goe-decorrelation}
\big| R_{4, N}(x_1, x_2, x_3, x_4) - R_{2, N}(x_1, x_2) R_{2, N}(x_3, x_4) \big| \le C N^4 (\log N) \Delta_0(\varrho).
\end{equation}
\end{enumerate}
\end{lemma}

\begin{proof}
We begin with (a). The points are pairwise within distance $N^{-1}$, so \zcref{l:goe-kernel}(c) shows that the entries $\mathsf I(x_i, x_j)$ are bounded. Hence all entries of the $6 \times 6$ matrix in \eqref{e:goe-pf-normalized} are bounded, and the third property of the Pfaffian gives $|R_{3, N}| \le 15 C^3 N^3$.

Next, we prove (b). The entries $\mathsf I(x_i, x_j)$ of the $8 \times 8$ matrix in \eqref{e:goe-pf-normalized} are bounded by $C \log N$ by \zcref{l:goe-kernel}(c), and all other entries are bounded. The entries $\mathsf I(x_i, x_j)$ occupy the positions where a first row meets a first column. There are four first rows, so each perfect matching of $[8]$ contains at most two pairs $\{\sigma, \sigma'\}$ for which $\sigma$ and $\sigma'$ are first rows. Hence each term in \eqref{e:goe-pf-def} is bounded by $C (\log N)^2$, and $|R_{4, N}| \le 105 \, C N^4 (\log N)^2$.

We now prove (c), following the use of congruence transformations in \cite[Section~5.2]{feng2024small}. Let $\mathsf W$ denote the $8 \times 8$ matrix in \eqref{e:goe-pf-normalized} for the points $x_1, x_2, x_3, x_4$. We define $\mathsf W'$ as the matrix obtained from $\mathsf W$ by subtracting row $1$ from row $3$ and row $5$ from row $7$, followed by the same operations on the columns. Thus row $3$ of $\mathsf W'$ is the difference of the first rows of $x_2$ and $x_1$, and row $7$ is the difference of the first rows of $x_4$ and $x_3$. We have $\mathsf W' = B^T \mathsf W B$ for a matrix $B$ with $\det B = 1$, so
\begin{equation*}
\operatorname{Pf} \mathsf W' = \operatorname{Pf} \mathsf W = N^{-4} R_{4, N}(x_1, x_2, x_3, x_4).
\end{equation*}
Let $P_1 = \{1, 2, 3, 4\}$ and $P_2 = \{5, 6, 7, 8\}$. The transformation does not mix $P_1$ and $P_2$, so $\mathsf W'_{P_1 P_1}$ is the congruence transform of $\mathsf W_{P_1 P_1}$ by a matrix of determinant one, and \eqref{e:goe-pf-normalized} with $n = 2$ gives
\begin{equation*}
\operatorname{Pf} \mathsf W'_{P_1 P_1} = \operatorname{Pf} \mathsf W_{P_1 P_1} = N^{-2} R_{2, N}(x_1, x_2), \qquad \operatorname{Pf} \mathsf W'_{P_2 P_2} = N^{-2} R_{2, N}(x_3, x_4).
\end{equation*}

We next bound the entries of $\mathsf W'$. Observe that for $x' \in [x_1, x_2]$ and $j \in \{3, 4\}$,
\begin{equation*}
|x' - x_j| \ge \varrho - (x_2 - x_1) - (x_4 - x_3) \ge \frac{\varrho}{2},
\end{equation*}
and the same holds for $x' \in [x_3, x_4]$ and $j \in \{1, 2\}$. Each entry of $\mathsf W'_{P_1 P_1}$ is a signed sum of at most four entries of $\mathsf W_{P_1 P_1}$, and $|x_1 - x_2| \le N^{-1}$, so the entries of $\mathsf W'_{P_1 P_1}$ are bounded by a constant by \zcref{l:goe-kernel}(a) and (c). The same holds for $\mathsf W'_{P_2 P_2}$.

Now consider an entry $\mathsf W'_{\sigma \sigma'}$ with $\sigma \in P_1$ and $\sigma' \in P_2$. When $\sigma$ and $\sigma'$ are second rows, we have $\mathsf W'_{\sigma \sigma'} = - N^{-2} \mathsf V(x_i, x_j)$ for some $i \in \{1, 2\}$ and $j \in \{3, 4\}$. When exactly one of $\sigma$ and $\sigma'$ is a first row, the entry $\mathsf W'_{\sigma \sigma'}$ is of the form $\pm N^{-1} \mathsf T(x_i, x_j)$ or $\pm N^{-1} \mathsf T(x_j, x_i)$ with $i \in \{1, 2\}$ and $j \in \{3, 4\}$, or a difference of two such entries. By \zcref{l:goe-kernel}(b) and $N \varrho \ge 1$,
\begin{align*}
N^{-1} |\mathsf T(x_i, x_j)| + N^{-1} |\mathsf T(x_j, x_i)| &\le C \big( (N \varrho)^{-1} + N^{-1/2} \big) \le C \Delta_0(\varrho), \\
N^{-2} |\mathsf V(x_i, x_j)| &\le C \big( (N \varrho)^{-1} + (N \varrho)^{-2} + N^{-1} \big) \le C \Delta_0(\varrho)
\end{align*}
for $i \in \{1, 2\}$ and $j \in \{3, 4\}$, so $|\mathsf W'_{\sigma \sigma'}| \le C \Delta_0(\varrho)$ in these cases.

It remains to consider the first rows $\sigma \in \{1, 3\}$ and $\sigma' \in \{5, 7\}$, and there are four cases. First, the entry $\mathsf W'_{15} = \mathsf I(x_1, x_3)$ is not small, but it is bounded by $C \log N$ by \zcref{l:goe-kernel}(c). Second, since $x_3 \notin [x_1, x_2]$, \eqref{e:goe-I-formula} gives
\begin{equation*}
\mathsf W'_{35} = \mathsf I(x_2, x_3) - \mathsf I(x_1, x_3) = \int_{x_1}^{x_2} \mathsf T(x', x_3) \dd x',
\end{equation*}
so \zcref{l:goe-kernel}(b) and the assumption that $x_2 - x_1 \le N^{-1}$ give
\begin{equation*}
|\mathsf W'_{35}| \le C (x_2 - x_1) (\varrho^{-1} + N^{1/2}) \le C \Delta_0(\varrho).
\end{equation*}
Third, by the antisymmetry of $\mathsf I$ and \eqref{e:goe-I-formula},
\begin{equation*}
\mathsf W'_{17} = \mathsf I(x_1, x_4) - \mathsf I(x_1, x_3) = - \int_{x_3}^{x_4} \mathsf T(x', x_1) \dd x',
\end{equation*}
which is bounded in the same way. Fourth, the same argument gives
\begin{equation*}
\mathsf W'_{37} = \int_{x_1}^{x_2} \big( \mathsf T(x', x_4) - \mathsf T(x', x_3) \big) \dd x',
\end{equation*}
which is also bounded by $C \Delta_0(\varrho)$.

Finally, we expand $\operatorname{Pf} \mathsf W'$ according to \eqref{e:goe-pf-def}. By the second property of the Pfaffian, the matchings $\mathfrak m$ such that every pair of $\mathfrak m$ is a subset of $P_1$ or of $P_2$ contribute $\operatorname{Pf} \mathsf W'_{P_1 P_1} \operatorname{Pf} \mathsf W'_{P_2 P_2}$. Since $|P_1| = 4$ is even, every other matching contains an even number of pairs joining $P_1$ and $P_2$, hence at least two. Among these pairs, only the pair $\{1, 5\}$, which occurs at most once, carries an entry that is not bounded by $C \Delta_0(\varrho)$. Hence the corresponding term contains at least one factor bounded by $C \Delta_0(\varrho)$, at most one factor bounded by $C \log N$, and other factors bounded by $C$. Using $\Delta_0(\varrho) \le 2$, we conclude that
\begin{equation*}
\big| N^{-4} R_{4, N}(x_1, x_2, x_3, x_4) - N^{-4} R_{2, N}(x_1, x_2) R_{2, N}(x_3, x_4) \big| \le 105 \, C (\log N) \Delta_0(\varrho),
\end{equation*}
which is \eqref{e:goe-decorrelation}.
\end{proof}

\section{Mean, variance, and tails}\label{s:goe-proof}

We now prove \zcref{p:goepairs}. Recall from \eqref{e:prelim-scales} that $\ell_\alpha = N^{-1 - \alpha}$ and $M_\alpha = N^{1 - 2\alpha} = N^3 \ell_\alpha^2$. We choose $\delta \in (0, 1/4)$ such that $I \subset I_{4 \delta}$. For a bounded measurable function $g \colon [0, \infty) \times \R \to \R$, we define
\begin{equation}\label{e:goe-stat}
\mathcal X_\alpha(G; g) = \sum_{1 \le i < j \le N} g \Big( \frac{\lambda_j - \lambda_i}{\ell_\alpha}, \frac{\lambda_i + \lambda_j}{2} \Big), \qquad \mathcal L(g) = \int_{-2}^2 \int_0^\infty c(E) \, r g(r, E) \dd r \dd E,
\end{equation}
where $\mathcal L(g)$ is defined whenever the integral converges absolutely. The statistic $\mathcal X_\alpha(G; g)$ extends the pair statistic \eqref{e:comparison-pair} of $G$ to functions $g$ that are not products of a function of $E$ and a function of $r$. We also set
\begin{equation*}
g^\sharp(x, y) = g \Big( \frac{|y - x|}{\ell_\alpha}, \frac{x + y}{2} \Big), \qquad x, y \in \R,
\end{equation*}
so that
\begin{equation*}
\mathcal X_\alpha(G; g) = \frac12 \sum_{i \neq j} g^\sharp(\lambda_i, \lambda_j).
\end{equation*}

Throughout, we consider functions $g$ that vanish for $E \notin I$. If such an $g$ is nonnegative or vanishes for all sufficiently large $r$, then \eqref{e:goe-corr-identity} with $n = 2$ and the substitution $y = x + v$ give
\begin{equation}\label{e:goe-first-moment}
\begin{split}
\E \mathcal X_\alpha(G; g) &= \int_{x < y} g \Big( \frac{y - x}{\ell_\alpha}, \frac{x + y}{2} \Big) R_{2, N}(x, y) \dd x \dd y \\
&= \int_0^\infty \int_\R g \Big( \frac v{\ell_\alpha}, x + \frac v2 \Big) R_{2, N}(x, x + v) \dd x \dd v.
\end{split}
\end{equation}
The integrand in \eqref{e:goe-first-moment} vanishes unless $x + v/2 \in I \subset I_{4 \delta}$, and in that case $x, x + v \in I_{3 \delta}$ whenever $v \le \delta$.

\paragraph{The mean.} Let $r_0 \ge 1$. Let $g$ be bounded and measurable with $g(r, E) = 0$ for $r > r_0$ or $E \notin I$. We claim that
\begin{equation}\label{e:goe-mean}
\lim_{N \to \infty} \frac{\E \mathcal X_\alpha(G; g)}{M_\alpha} = \mathcal L(g).
\end{equation}
Let $N$ be large enough that $r_0 \ell_\alpha \le N^{-1} \le \delta$. The integrand in \eqref{e:goe-first-moment} vanishes unless $v \le r_0 \ell_\alpha$, and in that case $x, x + v \in I_{3 \delta}$. Using \zcref{l:goe-two}(a) with $3 \delta$ in place of $\delta$, we obtain
\begin{align*}
\E \mathcal X_\alpha(G; g) &= N^3 \int_0^{r_0 \ell_\alpha} \int_\R g \Big( \frac v{\ell_\alpha}, x + \frac v2 \Big) c(x) v \dd x \dd v + O \Big( \int_0^{r_0 \ell_\alpha} \big( N^{23/8} v + N^4 v^2 \big) \dd v \Big) \\
&= N^3 \ell_\alpha^2 \int_0^{r_0} \int_I g(r, E) \, c \Big( E - \frac{\ell_\alpha r}{2} \Big) r \dd E \dd r + O \big( N^{23/8} \ell_\alpha^2 + N^4 \ell_\alpha^3 \big),
\end{align*}
where the implicit constants depend on $r_0$ and $\|g\|_\infty$, and in the second line we substituted $v = \ell_\alpha r$ and $x = E - \ell_\alpha r/2$. Since the function $E \mapsto c(E)$ is uniformly continuous on $[-2, 2]$, we have $c(E - \ell_\alpha r/2) - c(E) \to 0$ uniformly in $E \in I$ and $r \in [0, r_0]$. Dividing by $M_\alpha = N^3 \ell_\alpha^2$ and using $N^{23/8} \ell_\alpha^2 = N^{-1/8} M_\alpha$ and $N^4 \ell_\alpha^3 = N^{-\alpha} M_\alpha$, we obtain \eqref{e:goe-mean}.

\paragraph{The variance.} Let $r_0$ and $g$ be as in the previous paragraph. We claim that
\begin{equation}\label{e:goe-variance}
\lim_{N \to \infty} \frac{\Var \mathcal X_\alpha(G; g)}{M_\alpha^2} = 0.
\end{equation}
Expanding the square of $\sum_{i \neq j} g^\sharp(\lambda_i, \lambda_j)$ according to the number of indices that the two pairs share, and using \eqref{e:goe-corr-identity} for $n \in \{2, 3, 4\}$, we find
\begin{align}
\Var \mathcal X_\alpha(G; g) &= \frac12 \int g^\sharp(x_1, x_2)^2 R_{2, N}(x_1, x_2) + \int g^\sharp(x_1, x_2) g^\sharp(x_1, x_3) R_{3, N}(x_1, x_2, x_3) \notag \\
&\quad + \frac14 \int g^\sharp(x_1, x_2) g^\sharp(x_3, x_4) \big( R_{4, N}(x_1, x_2, x_3, x_4) - R_{2, N}(x_1, x_2) R_{2, N}(x_3, x_4) \big), \label{e:goe-var-expansion}
\end{align}
where the integrals are over $\R^2$, $\R^3$, and $\R^4$, and where we used
\begin{equation*}
\big( \E \mathcal X_\alpha(G; g) \big)^2 = \frac14 \int g^\sharp(x_1, x_2) g^\sharp(x_3, x_4) R_{2, N}(x_1, x_2) R_{2, N}(x_3, x_4) \dd x_1 \cdots \dd x_4.
\end{equation*}
Let $N$ be large enough that $8 r_0 \ell_\alpha \le N^{-1} \le \delta$. We bound the three terms in \eqref{e:goe-var-expansion} in turn, with constants depending on $r_0$ and $\|g\|_\infty$.

Since $R_{2, N} \ge 0$, the first term in \eqref{e:goe-var-expansion} is bounded by $\|g\|_\infty \E \mathcal X_\alpha(G; |g|)$, which is $O(M_\alpha)$ by \eqref{e:goe-mean} applied to $|g|$. In the second term, the integrand vanishes unless $|x_2 - x_1|, |x_3 - x_1| \le r_0 \ell_\alpha$ and $x_1, x_2, x_3 \in I_{3 \delta}$. The three points are then pairwise within distance $2 r_0 \ell_\alpha \le N^{-1}$, and \zcref{l:goe-corr}(a) bounds the second term by
\begin{equation*}
C N^3 \cdot |I_{3\delta}| (2 r_0 \ell_\alpha)^2 \le C M_\alpha.
\end{equation*}

We now turn to the third term. The integrand is invariant under exchanging $x_1$ and $x_2$ and under exchanging $x_3$ and $x_4$, so we may replace the factor $1/4$ by $1$ and restrict the integral to the region where $x_1 < x_2$ and $x_3 < x_4$. On this region, the integrand vanishes unless $x_2 - x_1 \le r_0 \ell_\alpha$, $x_4 - x_3 \le r_0 \ell_\alpha$, and all four points lie in $I_{3 \delta}$. Set $\varrho_0 = N^{-1/2}$ and $\varrho = |x_3 - x_1|$. We first consider the part of the region where $\varrho \le \varrho_0$. There we use \zcref{l:goe-corr}(b) and the bound $|R_{2, N}(x_1, x_2)| \le C N^3 (x_2 - x_1) \le C N^2$ from \zcref{l:goe-two}(a). This part of the support has measure at most $|I_{3\delta}| (r_0 \ell_\alpha)^2 \cdot 2 \varrho_0$, so its contribution is bounded by
\begin{equation*}
C N^4 (\log N)^2 \ell_\alpha^2 \varrho_0.
\end{equation*}
In the remaining part, where $\varrho > \varrho_0$, we have
\begin{equation*}
\varrho \ge N^{-1} \ge 8 r_0 \ell_\alpha \ge 4 (x_2 - x_1 + x_4 - x_3),
\end{equation*}
so \zcref{l:goe-corr}(c) applies. The contribution of this part is bounded by
\begin{align*}
& C N^4 (\log N) (r_0 \ell_\alpha)^2 \int_{I_{3\delta}} \int_{\{ x_3 : \varrho_0 < |x_3 - x_1| \le 4 \}} \Delta_0(|x_3 - x_1|) \dd x_3 \dd x_1 \\
& \qquad \le C N^4 \ell_\alpha^2 (\log N) \big( N^{-1} \log N + N^{-1/2} \big).
\end{align*}
Collecting the bounds and using $N^4 \ell_\alpha^2/M_\alpha^2 = N^{-2} \ell_\alpha^{-2} = N^{2 \alpha}$, we obtain
\begin{equation*}
\frac{\Var \mathcal X_\alpha(G; g)}{M_\alpha^2} \le C \big( M_\alpha^{-1} + N^{2\alpha} (\log N)^2 N^{-1/2} \big) = C \big( N^{2\alpha - 1} + (\log N)^2 N^{2\alpha - 1/2} \big).
\end{equation*}
The right side tends to zero because $\alpha < 1/4$, which proves \eqref{e:goe-variance}. This is the only place in the proof where the assumption $\alpha < 1/4$ is used.

\paragraph{The tails.} Let $r_0 \ge 1$. Let $g$ be measurable with $|g(r, E)| \le C_f (1 + r^2)^{-p}$ for all $(r, E)$ and $g(r, E) = 0$ for $r < r_0$ or $E \notin I$. We claim that
\begin{equation}\label{e:goe-tail}
\limsup_{N \to \infty} \frac{\E \mathcal X_\alpha(G; |g|)}{M_\alpha} \le C C_f \int_{r_0}^\infty r (1 + r^2)^{-p} \dd r.
\end{equation}
To prove this, we let $N \ge \delta^{-1}$ and split the integral over $v$ in \eqref{e:goe-first-moment}, with $|g|$ in place of $g$, into three ranges. In the range $r_0 \ell_\alpha \le v \le N^{-1}$, \zcref{l:goe-two}(a) gives $|R_{2, N}(x, x + v)| \le C N^3 v$, and the substitution $v = \ell_\alpha r$ bounds the contribution of this range by
\begin{equation*}
C C_f |I| \int_{r_0 \ell_\alpha}^{N^{-1}} \Big( 1 + \frac{v^2}{\ell_\alpha^2} \Big)^{-p} N^3 v \dd v \le C C_f M_\alpha \int_{r_0}^\infty r (1 + r^2)^{-p} \dd r.
\end{equation*}
In the range $N^{-1} \le v \le \delta$, we use $(1 + v^2/\ell_\alpha^2)^{-p} \le (\ell_\alpha/v)^{2p}$ and \zcref{l:goe-two}(b) on $I_{3 \delta}$. The contribution of this range is bounded by
\begin{align*}
C C_f N^2 (\log N) \int_{N^{-1}}^\delta \Big( \frac{\ell_\alpha}{v} \Big)^{2p} \dd v &\le C N^{2p + 1} \ell_\alpha^{2p} \log N \\
&= C N^{1 - 2\alpha p} \log N = C M_\alpha N^{-2\alpha (p - 1)} \log N.
\end{align*}
The range $v > \delta$ contributes the expectation of the sum of $|g((\lambda_j - \lambda_i)/\ell_\alpha, (\lambda_i + \lambda_j)/2)|$ over the pairs $i < j$ such that $\lambda_j - \lambda_i > \delta$. There are at most $N^2$ such pairs, and each term is bounded by $C_f (\ell_\alpha/\delta)^{2p}$, so the contribution of this range is bounded by $C N^2 \ell_\alpha^{2p} = C N^{2 - 2p - 2\alpha p}$. Since $p > 1$, the contributions of the last two ranges are $o(M_\alpha)$, which proves \eqref{e:goe-tail}.

\paragraph{Conclusion.} Let $r_0 \ge 1$. Set $g_1(r, E) = \1_{r \le r_0} f(r, E)$ and $g_2 = f - g_1$. Then $g_1$ satisfies the hypotheses of \eqref{e:goe-mean} and \eqref{e:goe-variance}, and $g_2$ satisfies the hypotheses of \eqref{e:goe-tail}. By the Cauchy--Schwarz inequality, \eqref{e:goe-mean}, and \eqref{e:goe-variance},
\begin{equation*}
\E \Big| \frac{\mathcal X_\alpha(G; g_1)}{M_\alpha} - \mathcal L(g_1) \Big| \le \Big( \frac{\Var \mathcal X_\alpha(G; g_1)}{M_\alpha^2} + \Big( \frac{\E \mathcal X_\alpha(G; g_1)}{M_\alpha} - \mathcal L(g_1) \Big)^2 \Big)^{1/2} \longrightarrow 0.
\end{equation*}
Since $\mathcal X_\alpha(G; f) = \mathcal X_\alpha(G; g_1) + \mathcal X_\alpha(G; g_2)$ and $\mathcal L(f) = \mathcal L(g_1) + \mathcal L(g_2)$, it follows from \eqref{e:goe-tail} that
\begin{equation*}
\limsup_{N \to \infty} \E \Big| \frac{\mathcal X_\alpha(G; f)}{M_\alpha} - \mathcal L(f) \Big| \le \limsup_{N \to \infty} \frac{\E \mathcal X_\alpha(G; |g_2|)}{M_\alpha} + |\mathcal L(g_2)| \le C C_f \int_{r_0}^\infty r (1 + r^2)^{-p} \dd r,
\end{equation*}
where in the last inequality we also used
\begin{equation*}
|\mathcal L(g_2)| \le C_f \Big( \sup_{E \in [-2, 2]} c(E) \Big) |I| \int_{r_0}^\infty r (1 + r^2)^{-p} \dd r.
\end{equation*}
Since $p > 1$, the right side tends to zero as $r_0 \to \infty$. Hence $\mathcal X_\alpha(G; f)/M_\alpha \to \mathcal L(f)$ in $L^1$, and therefore in probability by Markov's inequality. Since $\mathcal X_\alpha(G; f)/M_\alpha$ and $\mathcal L(f)$ are the two sides of \eqref{e:prelim-goepairs}, this completes the proof of \zcref{p:goepairs}.