Fix an integer $dge 3$. Let $G_{n,d}$ be the uniform random graph
among all simple labelled $d$-regular graphs on the vertex set
${1,...,n}$: no loops, no repeated edges, every vertex of degree
$d$, all such graphs equally likely. An integer $n$ is admissible
when this set is nonempty, that is, when $nd$ is even and $dle n-1$.
Write $A_{n,d}$ for its adjacency matrix: $A_{ij}=1$ if $i$ and $j$ are
joined, $A_{ij}=0$ otherwise. Each row has exactly $d$ ones, so $A_{n,d}$
is a very sparse matrix, and it is real symmetric. Let
[
  lambda_1ge lambda_2ge...ge lambda_n
]
be its eigenvalues, counted with multiplicity. They are real and lie in
$[-d,d]$.

A dense $ntimes n$ Gaussian matrix, by contrast, has all $tfrac12 n(n+1)$
entries random. Its eigenvalues are completely different objects at the
macroscopic scale. And yet, as we now explain, near any fixed interior
energy the two matrices produce the same local eigenvalue
statistics. That is the surprise.


Kesten (1959) and McKay (1981) computed the limiting eigenvalue density
of $G_{n,d}$: the empirical measure
$frac1nsum_{i=1}^ndelta_{lambda_i}$ converges weakly, in probability,
to the Kesten--McKay density

  rho_d(t);=;frac{d,√(4(d-1)-t^2)}{2pi,(d^2-t^2)}
  qquadfor |t|<2√(d-1),

and $rho_d(t)=0$ otherwise. It integrates to $1$. Its support is the
compact interval $[-2√(d-1),,2√(d-1)]$, the spectrum of the
$d$-regular tree.

For example, when $d=3$ the support is $[-2sqrt2,2sqrt2]$ and
$rho_3(0)=3sqrt8/(2picdot9)approx 0.150$, so the mean gap between
neighbouring eigenvalues near $0$ is about $1/0.150approx 6.7$ in the
variable $lambda$.

So near a fixed energy $E$ with $|E|<2√(d-1)$ there are about
$nrho_d(E)$ eigenvalues in a window of fixed length, hence gaps of
order $1/n$. Macroscopic statements such as (eq:km-density) cannot
distinguish two ensembles whose gaps differ by a factor of $2$; they
only see averages over windows containing millions of eigenvalues. The
interesting question is what happens at scale $1/n$, where individual
eigenvalues are visible.


Fix once and for all an energy
[
  Ein(-2√(d-1),,2√(d-1)),
]
and write $rho=rho_d(E)>0$. Two positive numbers measure everything:
the point $E$ and the density $rho$. The two-step zoom is
[
  lambda_i;⟼; x_i:=n(lambda_i-E)
  qquad(spectral coordinates),
]
and then, finally, we multiply by $rho$ to make the density at the
origin equal to $1$:
[
  Xi_{n,d,E};:=;sum_{i=1}^n
    delta_{,rho n(lambda_i-E)}.
]
A point measure is a locally finite sum of unit point masses, one
per eigenvalue. The prefactor $rho n$ is chosen so that near $E$ the
mean number of points per unit length of $Xi_{n,d,E}$ tends to $1$:
neighbouring eigenvalues sit about distance one apart in the new
coordinate.

A point measure can be tested against a continuous function $h$ with
compact support: the pairing is $int h,ddXi=sum_ih(points of Xi)$.
Two random point measures have the same law exactly when
[
  h;⟼;EEexp(-int h,ddXi)
]
is the same function for every bounded continuous compactly supported
$hge0$. These are the Laplace functionals; they are the natural
notion of convergence of point processes, and they are exactly what
Theorem~thm:main compares.

The comparison object is the Gaussian Orthogonal Ensemble. Let $H_n$ be
a real symmetric random matrix whose entries on and above the diagonal
are independent centred Gaussians, with variance $1/n$ off the diagonal
and $2/n$ on it. Its eigenvalues $mu_i$ follow the semicircle density
$frac1{2pi}√(4-t^2)$ on $[-2,2]$, whose value at the origin is
$1/pi$. Following the same zoom-and-unfold recipe, define
[
  Xi_n^{GOE}:=sum_{i=1}^ndelta_{(n/pi)mu_i}.
]


[Main result]
For every fixed integer $dge3$, every fixed energy $E$ in the open
interval $(-2√(d-1),2√(d-1))$, and every nonnegative continuous
compactly supported $h:RtoR$,
[
  |
  EEexp(-int h,ddXi_{n,d,E})
  -
  EEexp(-int h,ddXi_n^{GOE})
  |→ 0
]
as $nto∞$ through all admissible sizes.


Because the GOE point processes converge in law (Valk'o--Vir'ag), the
theorem identifies the full limiting local point process of the graph:
all finite joint distributions of eigenvalue counts in windows of
size comparable to the typical spacing, not only the one-point density,
not only energy-averaged correlations.

[In plain words]
Take a random $d$-regular graph on a huge number $n$ of vertices. Stand
at an energy $E$ inside the spectrum and zoom in on the eigenvalues
there until the typical distance between neighbouring eigenvalues is
$1$. The random picture of points you see has the same statistics as the
analogous picture for a dense random Gaussian matrix. This is true for
every interval, jointly, and it is happening although the two matrices
have nothing superficially in common: one has only $dn/2$ random bits,
the other has about $n^2/2$.


Two features of the statement deserve emphasis.


item No conditioning. The theorem makes no restriction to
connected or non-bipartite graphs. The relevant local information is
present in the unconditioned ensemble, and events of vanishing
probability are discarded freely.
item Fixed degree, fixed energy. Both $d$ and $E$ are held fixed
while $nto∞$. Degree $d=2$ is genuinely different: a $2$-regular
graph is a disjoint union of cycles, all eigenvalues are explicit
cosines, and the local statistics are not the Gaussian ones. Thus
$dge3$ is the honest hypothesis, and ``$d$ grows with $n$'' is a
different (earlier) theorem.



A random $d$-regular graph looks locally like the infinite $d$-regular
tree: for a fixed radius, the ball around a typical vertex is a tree
of degree $d$. Trees are not Gaussian; their spectral measure is
$rho_d$, not the semicircle. And yet Theorem~thm:main says the
microscopic eigenvalue statistics agree with those of a dense
Gaussian matrix. Locality of the graph and universality of the
eigenvalues coexist.

The conjecture in this form was made by Jakobson, Miller, Rivin and
Rudnick (1999), who also gave numerical evidence; they phrased it in
terms of nearest-neighbour spacings, the older and less complete
formulation of the same physics. What was known before this result:


item Macroscopic law. Kesten, McKay: (eq:km-density).
item Local laws. Bauerschmidt--Huang--Yau and Huang--Yau proved
the local Kesten--McKay law, delocalization and rigidity for fixed
$dge 3$; these are the spectral inputs used below. They are estimates
for smoothed eigenvalue counts (spectral windows of width at
least $n^{-1+o(1)}$), much coarser than individual spacings.
item Growing degree. Bauerschmidt--Huang--Knowles--Yau obtained
bulk gap universality for polynomial degrees; Bourgade--Huang obtained
fixed-energy point-process universality for $(log n)^{24}ll dle n^{1/2}$
and stated the fixed-degree problem as their Conjecture 2.11.
item Edges. Huang--McKenzie--Yau proved optimal rigidity and
GOE edge universality for every fixed $dge3$.
item Gaussian waves. Backhausz--Szegedy proved qualitative
Gaussian limits for ``almost eigenvectors''; the entropy framework of
Backhausz--Szegedy--Trees and Backhausz--Bordenave--Szegedy is the
engine we upgrade to a quantitative, variance-one statement.


Two new mechanisms are needed for fixed $d$, and they are the heart of
the proof.

First, weak convergence of eigenvector statistics is not enough: we need
quantitative Gaussian moments, polynomially accurate, uniformly
for all eigenvectors in a spectral band. The entropy argument delivers
exactly this, with an error $n^{-c}$ on an event of probability tending
to one, uniformly over graph-dependent choices of eigenvectors.

Second, the local laws know only about sums of resolvents over many
eigenvalues. The bridge to a single energy is a two-edge switch:
an exact symmetry of the uniform graph ensemble that turns a purely
spectral identity into a self-consistency equation for the resolvent at
a fixed microscopic height, which we then drive down by a halving
argument.

Once the microscopic marked law exists, entropy is not needed again:
we use the switch symmetry together with the Gaussian marks to
insert a genuine Gaussian matrix perturbation, and the classical
Dyson Brownian motion universality finishes the job.


Everything rests on four propositions, proved in order, plus a
collection of imported theorems (local laws, edge universality,
Dyson Brownian motion universality, semicircle edge limits).


item the entropy argument: Gaussian moments at sampled
edges. For fixed $d,E$ and any fixed number of eigenvectors in a
microscopic band, the joint values of those eigenvectors at random edge
endpoints are Gaussian to polynomial precision.
Mechanism: a colored-star counting inequality becomes a
differential-entropy inequality after smoothing by a Gaussian field;
then a strict tensor inequality turns the small entropy deficit into
moment matching, by induction over Hermite degree.

item the microscopic-limit step: endpoint resolvents and the
microscopic limit. The moments upgrade to control of resolvent
compressions at sampled edges, uniformly down to fixed microscopic
height. A two-edge switch gives the scalar self-consistency equation
that pins the trace. Then a subsequential limit of the microscopic
process exists and comes with Gaussian endpoint marks.

item the switching step: continuous switching and
Gaussian insertion. The discrete switch symmetry, combined with
orthogonal invariance of the Gaussian marks, becomes an exact
continuous symmetry of the limiting law. Composing many small
orthogonal switches produces, in the limit, a diagonal matrix plus a
GOE perturbation with the same point process as the microscopic law.

item the Dyson Brownian motion step: identification. The diagonal
matrices satisfy the hypotheses of Dyson Brownian motion universality
at a single energy, after an exact adjustment of the free-convolution
density. A uniform exponential bound on GOE counts upgrades correlation
comparisons to Laplace functionals. That is Theorem~thm:main.


The logical shape is therefore:
[
local laws⇒Gaussian moments⇒
microscopic limit⇒GOE noise⇒
GOE statistics.
]
No step uses a quantitatively strong statement about a spectral window
of microscopic width on the original graph; the microscopic information
is produced by the entropy step and then transported by symmetry and
Brownian motion.


Throughout, $dge3$ and $E$ with $|E|<2√(d-1)$ are fixed, and every
constant may depend on $d$ and $E$ without comment. We write $c,C$ for
positive constants whose values may change from one occurrence to the
next, and $o(1)$ for a quantity tending to zero as $nto∞$. On
events of probability tending to one we write $n^{o(1)}$ for a
deterministic subpolynomial bound. The notation $O_{PP}(b)$ means a
random quantity whose ratio to $b$ is bounded in probability. All
eigenvectors are real eigenbasis vectors of $A_{n,d}$, normalized to
have squared norm $n$; the choice inside an eigenspace is arbitrary
and measurable.


This section fixes notation, records the three counting facts about the
uniform regular-graph model that we use, states precisely the spectral
estimates imported from Huang--Yau, and derives the tree identities
that control the eigenvector covariance and the entropy cost. Only the
spectral estimates are quoted; everything else here is proved.

[In plain words]
Three pieces of background are fixed here once and for all. First, the
uniform $d$-regular graph is studied through the configuration model:
$d$ labelled stubs at every vertex, matched uniformly at random, then
conditioned on producing a simple graph; the conditioning is uniform,
with probability bounded below. Second, the sharp local spectral
estimates of Huang--Yau give rigidity, delocalization, and a resolvent
approximation by the infinite tree value at every height above $S/n$.
Third, on the infinite $d$-regular tree, boundary values of the resolvent
give the exact covariance of eigenvector values at nearby vertices; this
is what makes the Gaussian endpoint marks appear later.



For $Imzeta>0$ put
[
  m_d(zeta)=int_{R}(rho_d(t))/(t-zeta),dd t,
]
the Stieltjes transform of the Kesten--McKay density. It is analytic in
the upper half-plane with $Im m_d>0$. The boundary value from above
exists and has imaginary part $pirho_d$ on the open support. Unless a
subscript says otherwise we set, once and for all,

  rho=rho_d(E),qquad p=frac Ed,qquad
  a=(d)/(d-1)(1-p^2),qquad
  m_*=m_d(E+mathrm i0)=m_rmi0)=m_{rm r}+mathrm irmipirho .

Then $rho>0$, $a>0$, $pin(-1,1)$, and $m_{rm r}$ is a real number.
Note that

  1-p^2=(d^2-E^2)/(d^2),qquad
  a=(d^2-E^2)/(d(d-1))=(d)/(d-1)(1-p^2).

Beware that $a$ and $1-p^2$ differ by the factor $d/(d-1)$; both forms
are used below and they are not interchangeable.

Let $u_1,...,u_n$ be a real orthonormal eigenbasis of $A=A_{n,d}$,
normalized by $|u_i|^2=n$, so $Au_i=lambda_iu_i$ and
$< u_i,u_j>=ndelta_{ij}$. Define the microscopic
coordinates and the empirical spectral measure of the zoomed picture:


  x_i=n(lambda_i-E),qquad
  nu_n=sum_{i=1}^ndelta_{x_i},qquad
  m_n(z)&=sum_{i=1}^nfrac1{x_i-z}
        =frac1nTr R_n(z),
  R_n(z)&=(A-E-frac zn)^{-1},quad Im z>0 .


Thus $m_n$ is the Stieltjes transform of $nu_n$ in the coordinates
$x_i$, and $nu_n$ is recovered from $m_n$ by the Stieltjes inversion
formula. The final unfolded coordinates in
Theorem~thm:main are obtained from the $x_i$ by one more
multiplication by $rho$; this final scaling is a deterministic
continuous operation and will be applied only at the very end.


For $kge1$ define

  J_k=0&I_k I_k&0,qquad
  C_k=I_{2k}+pJ_k,qquad
  K_k(m)=mC_k+d^{-1}J_k .

The ordering of the $2k$ coordinates is always: the $k$ starting
vertices of $k$ sampled oriented edges, then their $k$ ending vertices.
This ordering matters and is kept fixed. The eigenvalues of $J_k$ are
$+1$ and $-1$, each with multiplicity $k$, so the eigenvalues of $C_k$
are $1+p$ and $1-p$: both positive because $|p|<1$. In particular
$C_k$ is positive definite and $C_k^{1/2}$ is well defined and bounded
with bounded inverse, uniformly in $k$.

The matrix $C_k$ is the covariance we will find for the Gaussian marks:
the value of an eigenvector at a sampled vertex has variance one after
normalization, and the two ends of a sampled edge have covariance $p$.


We take graphs from the configuration (pairing) model and condition on
simplicity. The facts we need are elementary and are proved here rather
than quoted.


In the configuration model with $d$ labelled half-edges, or stubs, at
each of $n$ labelled vertices, and with a uniformly random matching of
the $nd$ stubs, conditioning the matching to produce a simple graph
gives exactly the uniform simple labelled $d$-regular graph $G_{n,d}$.
Moreover, along all admissible $n$,
[
  PP(the matching produces a simple graph)
   → exp(-(d^2-1)/(4))>0 .
]
In a uniform simple $d$-regular graph: a fixed vertex lies in a triangle
with probability $O(n^{-1})$; and, for any fixed $k$, $k$
independently sampled uniform oriented edges have distinct endpoints
and induce exactly the sampled matching, with probability
$1-O_k(n^{-1})$.


Proof. 
Each simple $d$-regular graph corresponds to exactly $(d!)^n$ matchings
(one for each ordering of the stubs at the $n$ vertices), so conditioning
on simplicity is uniform.

For the asymptotic probability of being simple, count the expected
number of defects. A loop is a matched pair of stubs at the same
vertex; there are $binom d2$ candidate pairs at each of $n$ vertices,
each realized with probability $1/(nd-1)$. So the expected number of
loops is
[
  n,(d(d-1))/(2)cdotfrac1{nd-1}→ (d-1)/(2).
]
A double edge is a pair of matched stub pairs between the same
two distinct vertices. Choose the two vertices, two stubs at each, and a
pairing of the four stubs into two pairs: the expected number of double
edges is
[
  binom n2cdot ((binom d2)^2cdot 2)/((nd-1)(nd-3))
  → ((d-1)^2)/(4).
]
A prescribed family of $e$ disjoint stub pairs is realized with
probability $prod_{j=0}^{e-1}(nd-1-2j)^{-1}sim(nd)^{-e}$, for fixed
$e$. Hence, for fixed $j$, the $j$-th factorial moment of the total
defect count $Z$ (loops plus double edges) is asymptotic to
$((d-1)/(2)+((d-1)^2)/(4))^j$: for disjoint defect
tuples this is immediate from the product formula, and for tuples that
share vertices one checks that they contribute $o(1)$. Indeed, view the
prescribed pairs of an overlapping tuple as a multigraph with $v$
vertices and $e$ edges. Every vertex has degree at least $2$. A
connected component with all degrees exactly $2$ is a single prescribed
loop or a single prescribed double edge; a genuinely overlapping
component has $v<e$, there are $O(n^v)$ choices, and each occurs with
probability $O(n^{-e})$. Since there are only finitely many overlap
patterns for fixed $d,j$, the contribution is $o(1)$. Since $Z$
converges in distribution to a Poisson variable with parameter
$(d^2-1)/(4)$, and $Z=0$ is exactly simplicity, we get the stated
limit; in particular it is positive.

For a triangle through a fixed vertex, sum over the other two vertices
and their stub choices: $O(n^2)cdot O(n^{-3})=O(n^{-1})$ in the
matching model, and conditioning on simplicity multiplies this by a
bounded factor.

Finally, in a simple $d$-regular graph a uniform vertex is uniform. Each
previously chosen endpoint has only $O_d(1)$ vertices at distance at
most one from it, and a uniform oriented edge has a uniform starting
vertex. A union bound over the fixed number of sampled edges gives the
assertion.



Put

  S=exp(√(log n)),

a quantity that is $n^{o(1)}$ but larger than any fixed power of
$log n$. We import the following deliberately weak package from
Huang--Yau.

[Spectral input]
There exist constants $c_0,c_1>0$, $B<∞$, and an open interval
$I_E$ containing $E$ whose closure lies in
$(-2√(d-1),2√(d-1))$, such that the following hold with
probability tending to one, along all admissible $n$.

item After deleting the eigenvalue attached to the constant vector, the
eigenvalues are within $n^{-c_0}$ of their Kesten--McKay classical
locations; and every eigenvector entry is bounded by a deterministic
$n^{o(1)}$.
item Uniformly for $Rezetain I_E$ and $S/nleImzetale 1$,
[
  |frac1nTr(A-zeta)^{-1}-m_d(zeta)|
  le (log n)^B(n^{-c_1}+(nImzeta)^{-c_1}).
]
Except possibly at $n^{1-c_1}$ centres, the resolvent entries at a
centre and at its $d$ neighbours agree with the corresponding infinite
$d$-regular-tree entries up to the same bound.

In particular $#{i:|x_i|le S^2}=O(S^2)$ with probability tending to
one.


Proof. 
This is a restatement in our normalization of Theorems 1.2, 1.4 and 4.2
of Huang--Yau, in the arXiv version 2102.00963v3. Their matrix is
$A/√(d-1)$, so their bulk energy is $E/√(d-1)in(-2,2)$;
rescaling the spectral parameter and the resolvent by the fixed factor
$√(d-1)$ gives our normalization.

We keep track of why their error has our shape. Choose fixed
$0<c<1$, $a_{rm HY}ge12$, $b_{rm HY}ge25a_{rm HY}$ and
$0<mathfrak r<c/32$,frakr<c/32$, and let $r_{rm HY}$ be
$mathfrak rfrakrlog_{d-1}n$ rounded to an integer. Their Definition 1.1
and Proposition 2.1 produce an event of probability tending to one on
which all but a polynomially small fraction of vertices have a tree
neighbourhood of radius $(c/4)log_{d-1}n$; that radius is much larger
than $r_{rm HY}$, and in particular larger than any neighbourhood
appearing below. On that event their resolvent comparison, with their
error (4.6)--(4.7), bounds the discrepancy by a logarithmic factor
times
[
  (d-1)^{-r_{rm HY}}+(nImzeta)^{-1/2}+(nImzeta)^{-2/3}.
]
Since $S$ dominates every fixed power of $log n$, and after decreasing
$c_1$ and increasing $B$, this is the stated bound, uniform in the
indicated range. The same two statements supply polynomial rigidity and
polylogarithmic complete delocalization, which is the bound on
eigenvector entries.

For the count, apply the trace estimate at $zeta=E+mathrm irmi S^2/n$.
Every eigenvalue with $|lambda_i-E|le S^2/n$ contributes at least
$n/(2S^2)$ to the imaginary part of $(x_i-zeta/n)^{-1}$; the
normalized imaginary trace is bounded by the estimator above, and the
count is $O(S^2)$.



Proposition~prop:spectral-input describes events in the
unconditioned ensemble of simple graph. We never condition the final
statement on connectedness, bipartiteness, or a tree-neighbourhood
property. All our estimates hold on events of probability tending to
one, and the observables in Theorem~thm:main are bounded, so
discarding vanishing events cannot change a limit.



The Gaussian-wave covariance on a regular tree is the object that
records how eigenvector values at nearby vertices are correlated. We
derive it, together with the logarithmic potential that appears in the
entropy computation, so that no normalization is hidden.

Fix a locally finite graph with adjacency matrix $A$ and a point $E$.
For vertices $v,w$ define the covariance kernel from the boundary
values of $G^zeta=(A-zeta)^{-1}$:

  K(v,w)=lim_{eta↓0}frac{Im< e_v,(A-E-ieta)^{-1}e_w>}{pirho}.

The kernel in the numerator, $eta/((A-E)^2+eta^2)$, is
positive semidefinite, so $K$ is a covariance kernel.


Let $q(zeta)$ be the solution of
[
  zeta=q^{-1}+(d-1)q
]
that tends to $0$ as $zetato∞$ in the upper half-plane. Then

  -m_d(zeta)=frac1{zeta-dq(zeta)} .

The covariance kernel (eq:cov-definition) of the infinite
$d$-regular tree is: variance $1$ at every vertex, covariance $p$ between
a centre and a neighbour, and covariance
[
  (E^2-d)/(d(d-1))
]
between two distinct neighbours of a centre. Moreover

  k_0=(1+Em_*)/(d)
  quadsatisfiesquad k_0^2-k_0-m_*^2=0 .

Finally, with $U(E)=intlog|t-E|rho_d(t),dd t$,

  -U(E)=-frac12log d+(d-1)/(2)log a-frac d4log(1-p^2),

and replacing $E$ by $E+mathrm irmieta$ changes the real part of the
logarithmic potential by $O(eta)$.


Proof. 
The Green function. Let $o$ be the root of the full tree and let
$F(zeta)$ be the diagonal entry of $(zeta-A_{rm branch})^{-1}$ at the
root of a rooted tree in which the root has $d-1$ children, recursively.
Each child subtree is again such a branch, so a Schur complement at the
root gives
[
  F=frac1{zeta-(d-1)F},
  qquadi.e.qquad zeta=F^{-1}+(d-1)F,
]
and $Fto0$ at infinity; thus $F=q$. At the root of the full tree, a
Schur complement over the $d$ branches gives
[
  bigl< e_o,(zeta-A)^{-1}e_obigr>=frac1{zeta-dq}.
]
The left side is $intrho_d(t)/(zeta-t)dd t=-m_d(zeta)$, which proves
(eq:tree-green). Its boundary imaginary part is
$pi^{-1}Im m_d(t+mathrm i0)=rmi0)=rho_d(t)$: this follows by checking that
the boundary function is nonnegative, compactly supported on
$[-2√(d-1),2√(d-1)]$, and has total mass $1$, the latter from
$m_d(zeta)sim-1/zeta$ at infinity and a contour argument.

Covariance. Write $G=(A-zeta)^{-1}$ on the tree and abbreviate
$g=G_{oo}=m_d$, $g_1=G_{ox}$ for a neighbour $x$, $g_2=G_{xy}$ for two
distinct neighbours of $o$; these are independent of the chosen
vertices by homogeneity. The resolvent equation $(zeta-A)G=I$ at the
root gives
[
  zeta g-1=d g_1 ,
  qquadsoqquad g_1=(1+zeta m_d(zeta))/(d),
]
and applied to a neighbour column gives
[
  zeta g_1=g+(d-1)g_2 .
]
Taking $zeta↓ E$, the imaginary parts of $g=pirho$ and of
$g_1=ppirho$ give, after dividing by $pirho$, the variance $1$ and
the centre-neighbour covariance $p$; then
[
  (d-1)Im g_2=Im(zeta g_1)-Im g
  =(E^2)/(d)pirho-pirho,
]
i.e. $Im g_2/(pirho)=(E^2-d)/(d(d-1))$, which is the stated
distinct-neighbour covariance. By (eq:tree-green) the function
$q$ and hence $m_d$ are analytic across each side of a small complex
neighbourhood of $E$, so all these boundary values are attained with
polynomial rate in the imaginary part.

The quadratic equation. Attach the value $q$ at $zeta=E+i0$ to
the identity $zeta-dq=-1/m$; eliminating $q$ using
$zeta q-(d-1)q^2=1$ gives the quadratic
[
  (d^2-E^2)m^2+E(d-2)m+(d-1)=0quadat m=m_* .
]
The same computation with $m$ arbitrary is exactly
$k_0^2-k_0-m^2=0$ for $k_0=(1+Em)/d$, since
[
  ((1+Em)/(d))^2-(1+Em)/(d)-m^2
  =((E^2-d^2)m^2+E(2-d)m+(1-d))/(d^2),
]
which is a negative multiple of the left side above. This is
(eq:k0). In the strict bulk the quadratic has two real roots, one
of which is $m_*$ with $Im m_*=+pirho$.

The potential. Differentiate the analytic function
$-m_d(zeta)=q/(1-q^2)$ and integrate along a path with
$zeta=q^{-1}+(d-1)q$:
[
  Reintlog(zeta-t),rho_d(t),dd t
  =-log|q|-(d-2)/(2)log|1-q^2| ,
]
matched at infinity using $m_d(zeta)sim-1/zeta$. At the boundary,
[
  |q|=(d-1)^{-1/2},qquad
  |1-q^2|=frac{√(d^2-E^2)}{d-1} ,
]
by evaluating the quadratic at $zeta=E$. Substituting into
(eq:tree-potential) and using $a=d(1-p^2)/(d-1)$ in the form
$1-p^2=(d^2-E^2)/d^2$ gives the formula. (As a check at $d=3$, $E=0$:
$|q|=1/sqrt2$, $|1-q^2|=3/2$, and the right side is
$log(sqrt3/2)$.)

Finally, $tmapstolog|t-E|$ is bounded near $E$ because $rho_d$ is,
and
[
  intlog|1-(mathrm irmieta)/(t-E)|rho_d(t),dd t
  =frac12intlog(1+(eta^2)/((t-E)^2))rho_d(t),dd t
  =O(eta),
]
by the substitution $t-E=eta u$ and the integrability of
$log(1+u^{-2})$ at infinity.



The first of the four steps gives a fully quantitative answer to a
qualitative question: do eigenvector values look Gaussian? The answer is
yes, with polynomial accuracy, uniformly over all choices of eigenvectors
in a microscopic band around $E$. This is the only place where the graph
structure (as opposed to imported spectral input) is used in the first
half of the proof.

Why do we need such precision? A fixed moment of a fixed number of
eigenvectors is later summed over eigenindices, and the number of
indices in a band of width $S^2/n$ is $n^{o(1)}$; a fixed negative power
$n^{-c}$ is more than enough to beat any subpolynomial factor. What
would not be enough is abstract weak convergence, which carries no rate.

The strategy is the entropy method of Backhausz--Szegedy, made
quantitative. A counting inequality for colourings of the graph
(the coloured-star count) bounds the number of labelings with
given local statistics. After a Gaussian smoothing, which turns the
counting bound into a differential-entropy inequality, the deficit
between the entropy of the star observation and half that of the edge
observation is small. The projection step removes the
one nearly deterministic relation (the eigenvector equation), and
the tensor step shows that this small entropy deficit
forces all Hermite moments of low degree to be small. That is Gaussian
behaviour.


[Uniform endpoint moments]
Fix positive integers $k,ell$ and a nonnegative integer $q$. There are
constants $c,C>0$, and events of graphs with probability tending to one,
with the following property. Choose distinct indices
$i(1),...,i(ell)$ satisfying $|x_{i(j)}|le S^2$, by any rule that
depends only on the graph. Independently of that choice, sample $k$
uniform oriented edges $(v_1,v'_1),...,(v_k,v'_k)$. Set
[
  b_j=(u_{i(j)}(v_1),...,u_{i(j)}(v_k),
           u_{i(j)}(v'_1),...,u_{i(j)}(v'_k))^{mathsf T}sfT
  inR^{2k},
  qquad j=1,...,ell .
]
Then for every monomial $P$ of total degree at most $q$ in the
coordinates of $b_1,...,b_ell$,

  |EE[P(b_1,...,b_ell)mid A]
        -EE P(Z_1,...,Z_ell)|le Cn^{-c},

where $Z_1,...,Z_ell$ are independent centred Gaussian vectors with
covariance $C_k$ from (eq:channel-matrices).


The constants and the graph events are uniform over the choices of
eigenindices. They depend on $d,E,k,ell,q$ but not on $n$; no assertion
is made for dimensions or moment orders that grow with $n$. The
covariance $C_k$ is the one predicted by the tree: variance $1$, and
covariance $p$ across each sampled edge, with zero covariance between
different edges.

[In plain words]
Choose finitely many eigenvectors whose eigenvalues sit within
$S^2/n$ of $E$. Look at the values of these eigenvectors at the two
ends of a few random edges. Those values behave, to accuracy $n^{-c}$
and with probability $1-o(1)$, like independent Gaussian vectors whose
covariance matrix is the tree covariance. ``Being Gaussian'' here means
every fixed polynomial moment matches, which is enough for everything
we do later.



Write
[
  V(v)=(u_{i(1)}(v),...,u_{i(ell)}(v))inR^{ell}
]
for the vector of the chosen eigenvectors at a vertex. A star
observation is the value at a uniform vertex together with the values
at its $d$ neighbours in uniformly random order; an edge
observation is the pair of values at the two ends of a uniform oriented
edge. All randomness in an observation includes the choice of location;
the graph is conditioned on.

We need two notions of entropy. For a random variable on a finite set we
write $H$ for the discrete entropy; for a random vector with a density
we write $mathcal H$calH$ for the differential entropy; both use natural
logarithms. The only fact about $H$ we use besides its definition is
that the number of labelings with a given empirical distribution of
local observations is at most exponential in $n$ times the entropy of
that distribution; the only facts about $mathcal H$calH$ are the Gaussian
maximum-entropy bound and the continuity estimates for binning, both
spelled out below.

Positive polynomial error exponents may decrease from one occurrence to
the next; this is harmless because only finitely many steps occur.


The first ingredient does not involve eigenvectors at all. Give the
vertices of the configuration model labels from an alphabet of size $M$.
Let $T$ be the distribution of the ordered star observation, symmetrized
over the $d!$ orders of the neighbours, and let $Q$ be the distribution
of the oriented edge observation. Then $T$ and $Q$ are determined by the
labeling, and they are the only local information the counting
inequality sees. The following is the coloured-star count of
Backhausz--Szegedy, with the dependence on a growing alphabet kept
explicit.

[Uniform type count]
Suppose $M=M_n$ and, for some fixed $kappa>0$,

  M^{d+1}log nle n^{1-kappa}.

With probability tending to one, simultaneously for every pair of
histograms $T,Q$, the number of vertex labelings giving these histograms
is at most

  expBigl{nH(T)-(nd)/(2)H(Q)+O(n^{1-c})Bigr},

where $c>0$ and the constant in the error depend only on $d$ and
$kappa$.


Proof. 
Work in the configuration model with the $d$ stubs at each vertex
ordered. Fix an ordered star type $T_0$, i.e. a distribution on
the ordered observation $(z,x_1,...,x_d)$. Prescribe for every vertex
its label and the labels of its $d$ ordered stubs' desired partners.
For this fixed $T_0$ there are at most $exp(nH(T_0))$ such arrays,
because there are $n$ vertices and the empirical distribution of their
type is $T_0$.

Now fix also $Q$, and let $m_{jk}=nd,Q_{jk}$ be the number of stubs at
vertices of label $j$ whose desired partner has label $k$. A compatible
matching requires $m_{jk}=m_{kj}$ for all $j,k$, and requires $m_{jj}$
to be even for each $j$. Given these conditions, the number of uniform
matchings of the $nd$ stubs that realize the prescription is
[
  prod_{j<k}m_{jk}!;prod_j (m_{jj}-1)!! ,
]
the first factor choosing the bijection between the $j$-to-$k$ stubs and
their $k$-to-$j$ partners and the second pairing the $j$-to-$j$ stubs
among themselves. Dividing by the total number $(nd-1)!!$ of matchings
and expanding the factorials with Stirling's formula, the logarithm of
the probability is
[
  -(nd)/(2)H(Q)+O(M^2log n).
]
(The computation: each $mlog m$ term with $m=ndQ_{jk}$ contributes
$Q_{jk}log Q_{jk}$ times $nd$, the $-tfrac12 ndlog(nd)$ parts of the
numerators cancel the denominator, and the residual terms are
$O(M^2log n)$.)

Thus the expected number of labelings of ordered type $T_0$ and edge
type $Q$ is at most $exp{nH(T_0)-tfrac{nd}2H(Q)+O(M^2log n)}$.
Symmetrizing over neighbour orders can only increase the entropy of the
star type, so $H(T_0)le H(T)$; it does not touch $Q$. The number of
ordered types compatible with a given symmetrized $T$ is at most
$exp(O(M^{d+1}log n))$, because a type is a probability
vector on $M^{d+1}$ ordered observations. Summing over these, and using
(eq:type-count-hypothesis), changes only the sublinear error term.
Markov's inequality together with a union bound over all pairs
$(T,Q)$ upgrades the statement about expectations to a statement
holding for all pairs simultaneously with probability tending to one:
there are only $exp(O(M^{d+1}log n))$ pairs, and each has negligible
tail.

Finally, by Lemma~lem:model the probability that a uniform
matching is simple is bounded below by a positive constant, so
restricting to simple graphs multiplies the count and the probability by
bounded factors. The statement does not refer to the stub order, so it
passes to the uniform simple graph.



Counting alone sees only the discrete labels of a colouring. To connect
it to eigenvectors we smooth by a Gaussian field whose covariance is the
tree covariance at height $eta$, so that the labelings of the smoothed
field approximate a genuine continuous law.

Choose small fixed parameters $gamma,alpha>0$ with

  2ell(d+1)gamma<1,qquad
  4ellgamma+2alpha<frac14,qquad
  0<alpha<gamma,

and take $alpha$ also smaller than the rigidity exponent $c_0$ of
Proposition~prop:spectral-input. Put

  eta=n^{-alpha},

and, conditionally on the graph, let
$Y=(Y^{(1)},...,Y^{(ell)})$ be $ell$ independent Gaussian fields on
the vertices, each with the covariance

  Gamma=(1)/(pirho),(eta)/((A-E)^2+eta^2).

This is exactly the covariance kernel (eq:cov-definition) at
height $eta$, so its entries approximate the tree covariance as
$eta↓0$. For $0leεle1$ define the smoothed field
[
  F=ε V+Y ,
]
where $V$ is the eigenvector field $V(v)=(u_{i(1)}(v),...,u_{i(ell)}(v))$.
The parameter $eta$ regularizes the almost-deterministic eigenvector
relation; the independent parameter $ε$ measures how much of the
eigenvector signal is kept. We keep the entropy inequality uniformly in
$ε$ and later choose a small power of $n^{-1}$ when extracting
each fixed moment.

[Star--edge entropy]
On events of graphs of probability tending to one, uniformly over all
choices of $V$ and all $0leεle1$,

  mathcal H(F_calH(F_{rm star})-frac d2mathcal H(F_calH(F_{rm edge})
   ge frac1nmathcal H(Y_calH(Y_{rm global})-O(n^{-c}),

where $Y_{rm global}$ is the full $nell$-dimensional Gaussian array.
Moreover

  frac1nmathcal H(Y_calH(Y_{rm global})
  =ell(frac12logfrac{2pi{rm e}eta}{pirho}-U(E))
   +O(n^{-c}),

with $U(E)=intlog|t-E|rho_d(t),dd t$ from
Lemma~lem:tree-identities.


Proof. 
Step 1: quantize and count. Quantize every coordinate of $F$ into
intervals of width $b=n^{-gamma}$, clipping all values outside
$[-n^{gamma},n^{gamma}]$ into the two end bins. The resulting vertex
alphabet has

  M=(2n^{2gamma})^{ell}=n^{2ellgamma+o(1)}quadelements,

so by $2ell(d+1)gamma<1$ the hypothesis
(eq:type-count-hypothesis) holds with room to spare. Let
$mathcal L$calL$ be the resulting random quantized labeling of the vertices,
and let $T,Q$ be its star and edge histograms. Entropy of the type
index. The total number of labelings is $exp(nH(mathcal L))$,calL))$, the
number of possible type indices is at most
$exp(O(M^{d+1}log n))$, and for a given type index the count is at
most (eq:ent-type-count). Hence, by concavity of $H$ in its first
argument,

  frac1nH(mathcal L)calL)
  le H(EE T)-frac d2,EE H(Q)+O(n^{-c}) .


Step 2: concentrate the edge histogram. Recall that two distinct
uniform oriented edges have independent uniform endpoints. The
covariance of the Gaussian observations at the endpoints of two edges
therefore has, after averaging over locations, an off-block Frobenius
norm of at most
[
  C((TrGamma^2)/(n^2))^{1/2}
  le frac{C}{√(neta)},
]
because $TrGamma^2le Cn/eta$, and the smallest eigenvalue of a
single-edge covariance is at least $ceta$ (use $|A|le d$). After
whitening the two marginal blocks, the total-variation distance from the
product law is at most $|{rm off,block}|/lambda_{min}le
Cn^{-1/2}eta^{-3/2}$; for reference this follows by differentiating the
Gaussian density along the segment joining the two covariance matrices.
Consequently, for two independent uniform oriented edges and any cell of
the quantization, the covariance of the two cell indicators is at most
that distance. The edge histogram $Q$ is an average over the $nd$
oriented edges, and the joint law of two distinct oriented edges is
exactly that of two independent uniform oriented edges; hence every cell
frequency of $Q$ has variance at most
[
  frac1{nd}+C,n^{-1/2}eta^{-2}le C,n^{-1/2+2alpha},
]
the diagonal term being smaller. There are $M^2=n^{4ellgamma+o(1)}$
cells, so Cauchy--Schwarz bounds the mean $ell^1$ deviation
[
  EE|Q-EE Q|_1
  le (M^2cdot Cn^{-1/2+2alpha})^{1/2}
  =n^{2ellgamma-1/4+alpha+o(1)}=O(n^{-c}),
]
because $4ellgamma+2alpha<tfrac14$. For two probability vectors on a
set of size $M^2$ at $ell^1$ distance $h$, the entropies differ by at
most $hlog(M^2)+H_{rm bin}(h)$, where $H_{rm bin}$ is the binary
entropy; therefore

  EE H(Q)=H(EE Q)+O(n^{-c}).


Step 3: from discrete to differential entropy. Every marginal
covariance occurring here, including that of the whole array, has
spectrum in $[ceta,C/eta]$, because $Gamma$ does and because a
deterministic shift does not change covariance. Unbounded binning of a
Gaussian vector and the entropy bound on a cube give
[
  frac1nH(mathcal L)calL) ge frac1nmathcal H(Y_calH(Y_{rm global})
  -elllog b-O(n^{-c});
]
the clipping error is in fact superpolynomially small, since
$alpha<gamma$ makes the Gaussian tails beyond $n^{gamma}$
superpolynomially unlikely and their quantized entropy contributions
superpolynomially small, uniformly in the shift. For either local
observation, of dimension $ell(d+1)$ for the star and $2ell$ for the
edge, we also have the two-sided expansion

  H(binned X)=mathcal H(X)-(calH(X)-(dim X)log b+O(n^{-c}).

Here is the justification, since everything downstream rests on it.
Bin averaging changes a Gaussian density, or a mixture of Gaussians, in
total variation by at most $Cbeta^{-1/2}$, by integrating its absolute
first derivatives along each coordinate. Clipping adds only a
superpolynomial tail error. Both resulting densities have polynomial
supremum and second-moment bounds, and such bounds convert polynomial
total variation into polynomial entropy error: subtract the common
minimum of the two densities; the remainders have equal mass $h$ and,
after division by $h$, entropy bounded by $O(1+log(n/h))$; mixing the
parts with a common piece costs at most the binary entropy.

The laws $EE T$ and $EE Q$ are exactly the binned local laws.
Substituting (eq:ent-edge-concentration) and
(eq:ent-binning) into (eq:ent-quantized), the mesh terms
cancel because
$(d+1)ell-frac d2(2ell)=ell$, and the shift $ε V$ leaves the
differential entropy unchanged; since $F_{rm global}$ is a
deterministic shift of $Y_{rm global}$ we obtain
(eq:ent-star-edge).

Step 4: the global entropy. The determinant formula for a
Gaussian gives
[
  frac1nmathcal H(Y_calH(Y_{rm global})
  =(ell)/(2)logfrac{2pi{rm e}eta}{pirho}
   -(ell)/(n)sum_{i=1}^nlog|lambda_i-E-mathrm irmieta|,
]
because the eigenvalues of $Gamma$ are
$etabigl/(pirho,((lambda_i-E)^2+eta^2))$, each with
multiplicity $ell$. The function $tmapstolog|t-E-mathrm irmieta|$ is
Lipschitz with constant at most $eta^{-1}$. By rigidity we may replace
the eigenvalues by Kesten--McKay quantiles at a total error of at most
$ncdoteta^{-1}cdot n^{-c_0}$; the quantile sum differs from its
integral against $rho_d$ by $O((neta)^{-1})$, by summing the
lengths of the quantile intervals; an omitted or extra bounded
eigenvalue costs $O(log n/n)$. The choice $alpha<c_0$ makes the first
error $o(n)$. Lemma~lem:tree-identities supplies the remaining
regularization error $O(eta)$ in turning $log|t-E-mathrm irmieta|$ into
$log|t-E|$. This proves (eq:ent-global-entropy).



The entropy inequality (eq:ent-star-edge) still contains one
nearly deterministic direction. An eigenvector satisfies
$sum_{ysim v}u(y)=lambda u(v)$, so at a star with centre $v$ the sum
of the $d$ neighbour values is nearly $E$ times the centre value; this
is not a random relation, and a Gaussian model fitted to the raw star
would spend a whole direction on it. We remove that direction, so that
the surviving Gaussian has identity covariance, which is what the
Hermite expansion in the tensor step needs.

Write a star as $(Z,X)$ with $ZinR^{ell}$ (the $ell$ colour values
at the centre) and $XinR^{d} ⊗ R^{ell}$ (the colour values at
the $d$ neighbours). Let $P$ be the orthogonal projection of $R^d$ onto
${bf 1}^{perp}$, the hyperplane orthogonal to the all-ones vector, and
define the linear map to the ``clean'' coordinates, and the residual,
by
[
  mathcal B(Z,X)=(Z,(PX)/(calB(Z,X)=(Z,(PX)/(sqrt a)),qquad
  R_{rm res}(Z,X)=sum_{i=1}^dX_i-EZ .
]
Everything is tensored with the identity on the $ell$ colours, so
$D:=(R ⊕ {bf 1}^{perp}) ⊗ R^{ell}$ is the space in which
$mathcal B$calB$ takes its values. Let $w$ denote the unit vector along the
central direction $R$ inside $R ⊕ {bf 1}^{perp}$, and put
[
  v_i=p,w+sqrt a,Pe_i,qquad
  U_i=spanop(w,v_i) ⊗ R^{ell}quad(i=1,...,d),
]
where $e_i$ is the $i$-th standard basis vector of $R^d$. Since
$|Pe_i|^2=(d-1)/d$ and $a(d-1)/d=1-p^2$ by (eq:a-identity), we
get $|v_i|=1$ and $ip{w}{v_i}=p$. These two directions are the two
observable coordinates of the $i$-th edge: the centre and the $i$-th
neighbour. Their Gram matrix is
$1&pp&1$, with determinant $1-p^2>0$, so
they are linearly independent and the coordinate map on $U_i$ is
invertible.

Let $P_i$ be the orthogonal projection of $D$ onto $U_i$, put
$B_V=mathcal B(V_calB(V_{rm star})$ where $V_{rm star}$ is the true
eigenvector star, and let $G$ be an independent standard Gaussian vector
in $D$. All entropies on $D$ and on $U_i$ are computed with their
Euclidean volume.

[Projected-star entropy]
On the graph events of Lemma~ent:star-edge, uniformly for
$0leεle1$ and for each $i$,

  mathcal H(εcalH(ε B_V+G)-frac d2mathcal H(P_i(εcalH(P_i(ε
  B_V+G)) ge -O(n^{-c}).

Moreover $B_V$ has mean zero and second-moment matrix
$Id_D+O(n^{-1+o(1)})$, and every fixed moment of $P_iB_V$ is invariant,
up to $O(n^{-1+o(1)})$, under the orthogonal map of $U_i$ that
interchanges $w$ and $v_i$ and acts trivially on the colours.


Proof. 
The Jacobian. Describe the map
$(Z,X)mapsto(mathcal B(Z,X),R_calB(Z,X),R_{rm res}(Z,X))$ in the orthonormal
coordinates of $R^{d+1}=R ⊕ R{bf 1} ⊕ {bf 1}^{perp}$, one
colour at a time. On the central line it is the identity. On the
all-ones line, $R_{rm res}$ reads $sqrt d,t-EZ$ where $X=t{bf
1}/sqrt d+...$, so the diagonal entry is $sqrt d$. On each of the
$d-1$ remaining directions, $mathcal B$calB$ divides by $sqrt a$, so the
diagonal entry is $a^{-1/2}$. The map is triangular in these
coordinates, hence its absolute determinant, per colour, is
$sqrt d,a^{-(d-1)/2}$, and for $ell$ colours it is
$(sqrt d,a^{-(d-1)/2})^{ell}$.

The residual is small. For the noise, note the exact identity
[
  (A-E)Gamma(A-E)
  =(eta)/(pirho),((A-E)^2)/((A-E)^2+eta^2)
  le(eta)/(pirho),Id .
]
For the eigenvector part, the eigenvector equation gives
[
  sum_{i=1}^d u(y_i)=lambda,u(v)
  quadsoquad
  R_{rm res}(V_{rm star})(v)=(lambda-E)V(v),
]
and on the graph event $|lambda-E|le S^2/n$ for every selected index.
Thus, putting $Delta=diag(lambda_{i(1)}-E,...,lambda_{i(ell)}-E)$,
the deterministic residual of $F_{rm star}=ε V_{rm star}+Y_{rm star}$
at a uniform centre has, conditionally on the location, the covariance
matrix at most

  Cov(R_{rm res}(F_{rm star}))
  le(eta)/(pirho)Id_{ell}+ε^2Delta^2
  le((eta)/(pirho)+(S^4)/(n^2))Id_{ell},

because the Gaussian colours are independent and centred, so the cross
terms vanish, and because $|Delta|le S^2/n$ entrywise. The maximum
entropy bound by the covariance, together with the fact that
conditioning decreases entropy, therefore gives
[
  mathcal H(R_calH(R_{rm res}(F_{rm star})mid mathcal B(F_calB(F_{rm star}))
  le(ell)/(2)logfrac{2pi{rm e}eta}{pirho}+O(n^{-c}).
]
The relative correction inside the logarithm is
[
  log(1+(pirho S^4)/(n^2eta))
  =O(n^{-2+alpha+o(1)}),
]
which is polynomial because $S^4=n^{o(1)}$ and $eta=n^{-alpha}$ with
$alpha<1$.

The noise is standard. By the typical-star estimate of
Proposition~prop:spectral-input at the height $E+mathrm irmieta$,
together with the tree covariance of Lemma~lem:tree-identities,
the covariance matrix of the smoothed star noise
$mathcal B(Y_calB(Y_{rm star})$ differs from $Id_D$ by a polynomial error
outside a polynomially small fraction of centres. Two centred Gaussian
vectors whose covariance matrices differ by $delta$ in operator norm
have laws at total variation distance at most $Cdelta$ after
whitening; hence, using the entropy continuity argument of
Lemma~ent:star-edge, we may replace the actual noise by the
standard Gaussian $G$ at a polynomial entropy cost. The same argument
applies to the edge observation: its two coordinates in $U_i$ differ
from the observations of $ε B_V+G$ along $w$ and $v_i$ by
$O(n^{-1+o(1)})$, since the mean neighbour value is $p$ times the
centre value up to the eigenvector residual. The map from the
orthonormal coordinates on $U_i$ to the two edge coordinates has
absolute determinant $(1-p^2)^{ell/2}$.

Bookkeeping: the constants cancel. For a linear bijection $Phi$
with Jacobian $J_Phi$
[
  mathcal H(calH(Phi(X))=mathcal H(X)+calH(X)+log|det J_Phi| ,
]
since the density of $Phi(X)$ is the density of $X$ divided by
$|det J_Phi|$. Applied to the triangular map
$Phi=(Z,X)mapsto(mathcal B(Z,X),R_calB(Z,X),R_{rm res}(Z,X))$, whose determinant
per colour is $sqrt d,a^{-(d-1)/2}$, and combining with the
conditional bound on the residual, we get
[
  mathcal H(F_calH(F_{rm star})
  lemathcal H(calH(mathcal B(F_calB(F_{rm star}))
   +(ell)/(2)logfrac{2pi{rm e}eta}{pirho}
   -ell(frac12log d-(d-1)/(2)log a)+O(n^{-c}).
]
The edge observation is the image of $P_iF_{rm star}$ under the
coordinate map of $U_i$, which has determinant
$(1-p^2)^{ell/2}$, so
[
  mathcal H(F_calH(F_{rm edge})
  =mathcal H(P_iF_calH(P_iF_{rm star})+(ell)/(2)log(1-p^2)+O(n^{-c}).
]
Substituting both formulas into (eq:ent-star-edge), and using
(eq:ent-global-entropy) for the right-hand side, the terms with
$log(2pi{rm e}eta/pirho)$ cancel, and the remaining constant is
[
  -ell U(E)+ell(frac12log d-(d-1)/(2)log a)
  +(dell)/(4)log(1-p^2),
]
which is exactly zero by the potential identity
(eq:tree-potential). Hence
[
  mathcal H(calH(mathcal B(F_calB(F_{rm star}))
  -frac d2,mathcal H(P_iF_calH(P_iF_{rm star}) ge -O(n^{-c}).
]
Replacing the smoothed eigenvector star $F_{rm star}$ by
$ε B_V+G$ costs polynomial entropy in exactly the same way as in
the noise-replacement step above: the covariance of
$mathcal B(F_calB(F_{rm star})$ is $Id_D+O(n^{-1+o(1)})$ uniformly in
$ε$, and the entropy is a continuous function of the covariance
on the fixed finite-dimensional space $D$. This proves
(eq:ent-projected-entropy).

Second moments and the reversal symmetry. Put
$Λ=diag(lambda_{i(1)},...,lambda_{i(ell)})$. The exact
eigenvector relations, orthogonality of distinct eigenvectors, and
uniform permutation of the $d$ neighbours give the exact second-moment
identities
[
  EE ZZ^{mathsf T}sfT=EE X_iX_i^{mathsf T}sfT=Id_{ell},qquad
  EE ZX_i^{mathsf T}sfT=(Λ)/(d),qquad
  EE X_iX_j^{mathsf T}sfT=frac{Λ^2-d,Id_{ell}}{d(d-1)}
  quad(ine j).
]
The mean is zero because every selected eigenvector is orthogonal to the
constant vector. Since $Λ=E,Id_ell+O(n^{-1+o(1)})$ on the graph
event, these identities give, using $1-(E^2-d)/(d(d-1))=a$,
[
  EE B_VB_V^{mathsf T}sfT=Id_D+O(n^{-1+o(1)}).
]
Finally, a uniform oriented edge is invariant in law under reversal. Its
two coordinate vectors differ from the observations of $B_V$ along $w$
and $v_i$ by $O(n^{-1+o(1)})$, and since every coordinate of $V$ is
$n^{o(1)}$ that error stays of order $n^{-1+o(1)}$ in every fixed
moment. The coordinate map is invertible because $|p|<1$, so reversal is
exactly the orthogonal map of $U_i$ that interchanges $w$ and $v_i$ and
acts trivially on the colours. This proves the invariance assertion.



We have the inequality (eq:ent-projected-entropy), the approximate
reversal symmetry, and the second moments. What is still missing is the
inference from a small entropy deficit to small moments.
Two facts do the work: a formal identity expressing the entropy deficit
as a difference of relative entropies, and a strict inequality for
symmetric tensors that forbids a nonzero first Hermite coefficient.

Write $varphi_L$ for the standard Gaussian density on a Euclidean space
$L$, and for a random vector $Z$ with density $f$ on $L$ define the
relative entropy
[
  Ent(Law(Z)midvarphi_L)=int_L flog(f)/(varphi_L) .
]
Relative entropy is nonnegative, and it vanishes only for $Z$ standard
Gaussian.


Let $Z$ be a centred random vector in $D$ with identity covariance
matrix and finite differential entropies for $Z$ and for each $P_iZ$.
Then

  mathcal H(Z)-calH(Z)-frac12sum_{i=1}^dmathcal H(P_iZ)calH(P_iZ)
  =-Ent(Law(Z)midvarphi_D)
   +frac12sum_{i=1}^dEnt(Law(P_iZ)midvarphi_{U_i}).



Proof. 
For any $X$ in a Euclidean space $L$ with a density,
$mathcal H(X)=calH(X)=mathcal H(calH(varphi_L)-Ent(Law(X)midvarphi_L)$; this is
the definition of relative entropy together with $intvarphilogvarphi
=-mathcal H(calH(varphi)$. Also $mathcal H(calH(varphi_L)=(dim L)/(2)log(2pi{rm e})$.
Substituting into the left side of (eq:ent-deficit-sign), the
Gaussian entropies contribute
$(dim D)/(2)-frac12sum_i(dim U_i)/(2)$ times
$log(2pi{rm e})$, and this vanishes because
$dim D=dell$ while $dim U_i=2ell$ and there are $d$ terms. What
remains is exactly the right side.


Thus (eq:ent-projected-entropy) says that the full relative
entropy is at most half the sum of the projected relative entropies,
up to $O(n^{-c})$. The next lemma says that a distribution which is
symmetric under the swap of a centre and a neighbour and has variance
one cannot have most of its relative entropy invisible to all the
projections: the small deficit forces Gaussianity.

For an integer $jge1$ use the Euclidean norm on symmetric tensors in
$D^{ ⊗  j}$, and let $sigma_i$ be the orthogonal map of $U_i$ that
interchanges $w$ and $v_i$ and fixes the colour space.

[Tensor coercivity]
For every fixed $jge3$ there is $c_j>0$ such that every symmetric
$j$-tensor $H$ on $D$ satisfying
[
  sigma_i^{ ⊗  j}P_i^{ ⊗  j}H=P_i^{ ⊗  j}H
  qquadfor all i=1,...,d
]
obeys

  |H|^2-frac12sum_{i=1}^d|P_i^{ ⊗  j}H|^2 ge c_j|H|^2 .



Proof. 
A weighted frame. Work in the two-dimensional spatial plane
$spanop(w,v_i)$, which is the same abstract plane for every $i$ since
$|v_i|=1$ and $ip{w}{v_i}=p$. In an orthonormal basis with
$w=(1,0)$ and $v_i=(p,√(1-p^2))$, define the orthonormal
bisecting directions
[
  e_i^+=frac{w+v_i}{√(2(1+p))},qquad
  e_i^-=frac{w-v_i}{√(2(1-p))},qquad
  s_i=(e_i^++e_i^-)/(sqrt2),qquad
  t_i=(e_i^+-e_i^-)/(sqrt2).
]
The swap fixes $e_i^+$, negates $e_i^-$, and hence interchanges $s_i$
and $t_i$. Writing $q_0=√(1-p^2)$, a direct computation gives
[
  ip{w}{s_i}^2=(1+q_0)/(2),qquad
  ip{w}{t_i}^2=(1-q_0)/(2),qquad
  ip{Pe_i,s_i}^2=(1-q_0)/(2),qquad
  ip{Pe_i,t_i}^2=(1+q_0)/(2).
]
The strict bulk condition is
[
  p^2<(4(d-1))/(d^2)
  quad⇔quad
  q_0>(d-2)/(d),
]
with equality exactly at the spectral edges $|E|=2√(d-1)$; this is
the only place where the open bulk hypothesis and $dge3$ are used. In
particular the two numbers
[
  b_s=(q_0-(d-2)/d)/(2q_0),qquad
  b_t=(q_0+(d-2)/d)/(2q_0)
]
are positive and sum to one. With these weights the frame identity

  sum_{i=1}^d(b_s,s_is_i^{mathsf T}sfT+b_t,t_it_i^{mathsf T}sfT)
  =Id_{R ⊕ {bf 1}^{perp}}

holds in $R ⊕ {bf 1}^{perp}$, tensored with the identity on the
colours. To verify it, expand both sides in the orthogonal decomposition
$R ⊕ {bf 1} ⊕ {bf 1}^{perp}$: the central component of the
left side is
$d(b_sfrac{1+q_0}2+b_tfrac{1-q_0}2)=1$; the mixed
$R$--${bf 1}^{perp}$ block vanishes because $s_i$ and $t_i$ have
parallel central and $Pe_i$ components with $i$-independent
coefficients and $sum_iPe_i=0$; and the ${bf 1}^{perp}$ block is
$cP$ with
$c=b_s(1-q_0)/(2)+b_t(1+q_0)/(2)=1$, because
$sum_i(Pe_i)(Pe_i)^{mathsf T}sfT=P$.

Nonnegativity. Let $R^{(1)}_{i,s}$ denote the projection of the
first tensor slot onto $s_i$, tensored with the identity on the other
slots, and similarly for $t_i$. Inserting (eq:ent-frame) in the
first slot of $ip{H}{H}$ gives the exact identity
[
  |H|^2
  =sum_{i=1}^d(b_sbigl|R^{(1)}_{i,s}Hbigr|^2
                  +b_tbigl|R^{(1)}_{i,t}Hbigr|^2).
]
Since $s_i,t_iin U_i$, each norm decreases when the remaining slots are
projected onto $U_i$:
[
  bigl|R^{(1)}_{i,s}Hbigr|^2
  gebigl|R^{(1)}_{i,s}P_i^{(2)}... P_i^{(j)}Hbigr|^2,
]
and likewise for $t_i$. On $P_i^{ ⊗  j}H$ the swap invariance makes
the two norms equal, because $sigma_i$ interchanges $s_i$ and $t_i$ and
commutes with $P_i$. Since $s_i,t_i$ are orthonormal and span the
spatial part of $U_i$,
[
  bigl|R^{(1)}_{i,s}P_i^{ ⊗  j}Hbigr|^2
  +bigl|R^{(1)}_{i,t}P_i^{ ⊗  j}Hbigr|^2
  =bigl|P_i^{ ⊗  j}Hbigr|^2 ,
]
so each of the two is $frac12|P_i^{ ⊗  j}H|^2$. Using
$b_s+b_t=1$ now gives
[
  |H|^2 ge frac12sum_{i=1}^dbigl|P_i^{ ⊗  j}Hbigr|^2 ,
]
which is the nonnegativity half of (eq:ent-tensor).

Where equality is possible. Suppose equality holds. Then every
inequality above is an equality, so for every $i$ and every slot
$age2$ the discarded part vanishes:
$R^{(1)}_{i,s}H=R^{(1)}_{i,s}P_i^{(2)}... P_i^{(j)}H$ and likewise
for $t$. Adding the two identities over the orthonormal pair
$s_i,t_i$, and using symmetry of $H$ to move a projection from slot one
to any other slot, gives

  P_i^{(a)}H=P_i^{ ⊗  j}Hqquadfor every slot a.

Now take $ine k$ and three distinct slots, which exist because $jge3$:
[
  P_i^{(1)}P_k^{(1)}H=P_i^{(1)}P_k^{(2)}H
  =P_k^{(2)}P_i^{(3)}H=P_i^{(3)}P_k^{(1)}H=P_k^{(1)}P_i^{(1)}H .
]
Thus the commutator $[P_i,P_k]$ annihilates the first slot of $H$, for
all pairs $ine k$.

Let $a_i'in{bf 1}^{perp}$ be the unit vector along $Pe_i$. On
$R ⊕ {bf 1}^{perp}$ the spatial part of this commutator is
[
  ip{a_i'}{a_k'}(a_i'(a_k')^{mathsf T}sfT-a_k'(a_i')^{mathsf T}sfT),
  qquad ip{a_i'}{a_k'}=-frac1{d-1}ne0  (dge3),
]
which is an invertible rotation generator on
$spanop(Pe_i,Pe_k)$ and vanishes on its orthogonal complement. As
$i,k$ range, these generate all rotations of ${bf 1}^{perp}$, so
their common kernel in $R ⊕ {bf 1}^{perp}$ is the central line
$R w$. Hence $H$ is supported in $(w ⊗ R^{ell})^{ ⊗  j}$.

But $win U_i$, so $P_i^{ ⊗  j}H=H$, and the invariance
$sigma_i^{ ⊗  j}P_i^{ ⊗  j}H=P_i^{ ⊗  j}H$ reads
$sigma_i^{ ⊗  j}H=H$, i.e.
$v_i^{ ⊗  j} ⊗  h=w^{ ⊗  j} ⊗  h$ for the
corresponding symmetric $j$-tensor $h$ on $R^ell$. Since $w$ and
$v_i$ are linearly independent ($|p|<1$), this forces $h=0$, hence
$H=0$.

The quadratic form in (eq:ent-tensor) is continuous on the
finite-dimensional linear subspace defined by the swap constraints, and
it is strictly positive on the unit sphere of that subspace. By
compactness it is bounded below by $c_j|H|^2$ there.



We now extract Gaussian moments from the three ingredients above. The
spaces $D$, $U_i$, projections $P_i$ and swaps $sigma_i$ keep the
geometry of the projection step.

[Entropy criterion for Gaussian moments]
Fix an integer $qge2$. Let $B_n$ be random vectors in $D$ with
$|B_n|le K_n$ almost surely, where $K_n=n^{o(1)}$ is deterministic.
Suppose that for some $c_*>0$
[
  |EE B_n|+bigl|EE B_nB_n^{mathsf T}sfT-Id_Dbigr|=O(n^{-c_*}),
]
and that for every $1le jle q$ and every $i$
[
  bigl|(sigma_i^{ ⊗  j}-Id),EE(P_iB_n)^{ ⊗  j}bigr|
  =O(n^{-c_*}).
]
Let $G$ be an independent standard Gaussian in $D$ and suppose in
addition that, uniformly for $0leεle1$,
[
  mathcal H(εcalH(ε B_n+G)-frac12sum_{i=1}^d
    mathcal H(P_i(εcalH(P_i(ε B_n+G)) ge -O(n^{-c_*}).
]
Then every fixed polynomial moment of $B_n$ of degree at most $q$
differs from its standard Gaussian value by $O(n^{-c})$ for some $c>0$.
The conclusion is uniform over all families satisfying the hypotheses
with the same $c_*$, the same error constants, and the same bound
$K_n$.


Proof. 
Hermite tensors. For $yin D$ define symmetric tensors
$He_j(y)$ by the generating identity
[
  exp(ip{y}{t}-(|t|^2)/(2))
  =sum_{j=0}^∞frac{ip{He_j(y)}{t^{ ⊗  j}}}{j!},
  qquad tin D .
]
Set $H_j=EEHe_j(B_n)$. We show by induction on $j$ that

  |H_j|=O(n^{-c_j})qquad(1le jle q),

which is enough because Hermite polynomials form a basis of the ordinary
polynomials. The first- and second-moment hypotheses give $j=1,2$.

The induction step. Fix $jge3$ and assume the preceding orders
proved. Choose
[
  ε=n^{-tau},qquad
  2jtau<frac{c_*}2,qquad (j-k)tau<frac{c_k}2quad(1le k<j),
]
and put $u=ε/√(1+ε^2)$. Dividing every entropy
argument in the hypothesis by $√(1+ε^2)$ leaves that
inequality unchanged, because the dimensions cancel; the normalized
vector is
[
  Z_u=uB_n+√(1-u^2),G,
]
whose second-moment matrix is $Id_D+u^2O(n^{-c_*})$. Feeding this into
(eq:ent-deficit-sign), and absorbing the second-moment error into
the Gaussian entropies, gives

  Ent(Law(Z_u)midvarphi_D)
  -frac12sum_{i=1}^dEnt(Law(P_iZ_u)midvarphi_{U_i})
  le O(n^{-c_*}).


The relative entropy is essentially the $j$-th Hermite moment.
Let $f_u$ be the density ratio of $Z_u$ to $varphi_D$. The coefficient
of $u^k$ in the Taylor expansion of $f_u$ is
[
  frac1{k!}ip{H_k}{He_k(y)} .
]
Indeed, integrating the generating exponential against the law of
$B_n$ gives
$EEexp(uip{B_n}{t}-u^2|t|^2/2)$, and pairing two generating
exponentials under Gaussian measure gives Hermite orthogonality
together with
[
  Bigl|frac{ip{H_k}{He_k}}{k!}Bigr|_{L^2(varphi_D)}^2
  =(|H_k|^2)/(k!) .
]
The Taylor remainder through degree $j$ is
$O(u^{j+1}n^{o(1)})$ in $L^3(varphi_D)$. At a fixed source value $b_0$
the ratio at parameter $sin[0,u]$ is
[
  (1-s^2)^{-dim D/2}
  exp((|y|^2)/(2)-(|y-sb_0|^2)/(2(1-s^2))),
]
whose required derivatives are bounded by polynomials in
$|y|,|b_0|$ times $Cexp(Cs|b_0||y|)$; since
$|b_0|le K_n=n^{o(1)}$ and $u=n^{-tau+o(1)}$, their
$L^3(varphi_D)$ norms are $n^{o(1)}$, uniformly over the source
values. Consequently
[
  int_D|f_u-1|^3,ddvarphi_Dle Cu^3n^{o(1)} ,
]
and, since $|(1+h)log(1+h)-h-h^2/2|le C|h|^3$ for $|h|letfrac12$,
while $int(f_u-1)ddvarphi_D=0$, we get

  Ent(Law(Z_u)midvarphi_D)
  =frac12int_D(f_u-1)^2ddvarphi_D+O(u^3n^{o(1)}).


Extracting the $j$-th Hermite coefficient. Orthogonality of the
Hermite tensors gives, for the squared norm above,
[
  int_D(f_u-1)^2ddvarphi_D
  =u^{2j}(|H_j|^2)/(j!)
   +sum_{kne j}u^{2k}(|H_k|^2)/(k!) .
]
The terms with $k<j$ are $n^{-2ktau}O(n^{-2c_k})$, hence polynomially
small. For the whole sum we have a cruder bound that suffices: the
left side is a $chi^2$ divergence, and conditionally on $B_n=b_0$ the
density ratio is
$rho_{b_0,u}(y)=(1-u^2)^{-dim D/2}
 e^{< y,mu>-|mu|^2/2}
 e^{-u^2|y|^2/(2(1-u^2))}$
with $mu=ub_0/(1-u^2)$; a Gaussian integral then gives
$chi^2=O(u^2(dim D+|b_0|^2))$ for small $u$, whose mean is
$O(u^2)$, because $EE|B_n|^2=dim D+o(1)$. Therefore

  Ent(Law(Z_u)midvarphi_D)
  =frac{u^{2j}}{2}(|H_j|^2)/(j!)+O(n^{-2tau+o(1)}),

and the same identity holds on $U_i$ with $H_j$ replaced by
$P_i^{ ⊗  j}H_j$.

The quadratic-form inequality. Substitute
(eq:ent-relative-expansion) and its projected analogue into
(eq:ent-relative-deficit), multiply by $2j!/u^{2j}$, and use
$u^{-2j}=n^{2jtau+o(1)}$. The result is

  |H_j|^2-frac12sum_{i=1}^dbigl|P_i^{ ⊗  j}H_jbigr|^2
  le n^{-c_*+2jtau+o(1)}le n^{-delta}

for some $delta>0$, because $2jtau<c_*/2$.

Returning to the kernel. Let $K$ be the common kernel, i.e. the
closed subspace of symmetric $j$-tensors satisfying
$sigma_i^{ ⊗  j}P_i^{ ⊗  j}H=P_i^{ ⊗  j}H$ for every $i$.
Write $L_i=(sigma_i^{ ⊗  j}-Id)P_i^{ ⊗  j}$, so that
$L_iH_j$ has norm $O(n^{-c_*})$ by hypothesis, and note that
$K=∩_iker L_i$. The equality analysis of
Lemma~ent:tensor shows precisely that no nonzero tensor in
$K^{perp}$ is annihilated by all the $L_i$; since $D^{ ⊗  j}$ is
finite-dimensional and the $L_i$ are fixed, there is $c=c(d,E,j)>0$
with

  sum_{i=1}^d|L_iH|^2 ge c;dist(H,K)^2

for every symmetric $H$. Hence $dist(H_j,K)=O(n^{-c_*})$.
Let $H$ be the orthogonal projection of $H_j$ onto $K$ and $H'$ the
complement. The quadratic form
$Q(X)=|X|^2-frac12sum_i|P_i^{ ⊗  j}X|^2$ is bounded on
$D^{ ⊗  j}$, so
[
  |Q(H_j)-Q(H)|le C|H_j+H|,|H'|le n^{-c_*+o(1)},
]
where $|H_j|=n^{o(1)}$ was used. By
(eq:ent-quadratic-small), $Q(H)le n^{-delta'}$ for some fixed
$delta'>0$, and by Lemma~ent:tensor,
$Q(H)ge c_j|H|^2$. Therefore
[
  |H_j|le|H|+|H'|=O(n^{-delta'/2})+O(n^{-c_*})
  =O(n^{-c_j})
]
with a new fixed $c_j>0$. This closes the induction.

Polynomial moments. Every monomial of degree at most $q$ equals
$ip{T}{y^{ ⊗  j}}$ for some fixed symmetric tensor $T$ and some
$jle q$, and $y^{ ⊗  j}$ differs from
$He_j(y)$ by lower-order polynomial terms that are
already controlled. Hence
$EEip{T}{B_n^{ ⊗  j}}=ip{T}{H_j}+o(n^{-c})$ equals its Gaussian
value in every fixed degree. The constants in the bounds above depend
only on $d,E,j,c_*$ and the fixed error constants; the fact that the
choice of $ε$ depends on $j$ is harmless because only finitely
many $jle q$ occur. This proves the lemma.


Proof. [Proof of Proposition~prop:edge-moments]
Condition on a graph in the common event of Lemma~ent:type-count
and Proposition~prop:spectral-input, and let $B_n=B_V$, the
projected star field of the projection step. Complete
delocalization gives the deterministic bound $|B_n|le n^{o(1)}$.
Lemma~ent:projected-entropy supplies the entropy hypothesis, the
first two moments $EE B_n=0$ and
$EE B_nB_n^{mathsf T}sfT=Id_D+O(n^{-1+o(1)})$, and the moment
symmetries $(sigma_i^{ ⊗  j}-Id)EE(P_iB_n)^{ ⊗  j}=O(n^{-1+o(1)})$
for every fixed $j$, all uniformly over the selected eigenindices and
in $0leεle1$. Averaging (eq:ent-projected-entropy) over
$i$ gives the criterion's entropy hypothesis, with the harmless factor
$d$.
Lemma~ent:gaussian-criterion, applied with $q$ replaced by
$max(q,2)$, now gives that every moment of $B_n$ of degree at most $q$
differs from the corresponding standard Gaussian moment in $D$ by
$O(n^{-c})$, uniformly over the choices.


From the projected star to the edge values. The observations of
the standard Gaussian $G$ along $w$ and $v_i$ form a $2ell$-vector
with covariance $C_ell$, because $G$ has covariance $Id_D$ in the
$dell$ dimensions of $D$ and $ip{w}{v_i}=p$. It remains to compare
these coordinates with the actual edge values
$(V(v),V(v'))inR^{ell} ⊕ R^{ell}$ of the coloured-star count. Write $B_V=(Z,frac1{sqrt a}PX)$ as in
the projection step. Then
[
  ip{v_i}{B_V}
  =pip{w}{B_V}+ip{Pe_i}{PX}
  =pZ+X_i-frac1dsum_{j=1}^dX_j ,
]
and $sum_jX_j=R_{rm res}+EZ$ by definition of the residual, so
[
  ip{v_i}{B_V}
  =X_i+(p-frac Ed)Z-frac1dR_{rm res} .
]
Here $p=E/d$ exactly, and $|R_{rm res}|=O(n^{-1+o(1)})$
because $|lambda-E|le S^2/n$ and $|V(v)|=n^{o(1)}$. Hence the edge
value vector equals the $(w,v_i)$-observations of $B_V$ up to
$O(n^{-1+o(1)})$ in every coordinate, and since all coordinates are at
most $n^{o(1)}$ in absolute value, every fixed coordinate monomial
changes by $O(n^{-1+o(1)})$. This is the only place where the exact
identity $p=E/d$ enters, and it is
essential.

Several edges. For $k$ conditionally independent uniform
oriented edges, a coordinate monomial factors over the edges, because
the edges are chosen independently and the eigenvectors are fixed once
the graph is. The one-edge estimate, applied with $ell$ channels,
therefore gives the product estimate with an error of order
$k,n^{-c+o(1)}$, and $k$ is fixed. This concludes the proof of
Proposition~prop:edge-moments.


The section has produced Gaussian moments at one spectral
energy $E$, for a fixed number of eigenvectors, fixed number of edges
and fixed moment order, with a uniform error $n^{-c}$ over graph
dependent choices of eigenvectors. Nothing here is uniform in the
moment order or the number of edges, and nothing here applies to a
window of energies; both limitations are consistent with our later use
and are the reason the proof strategy lets $nto∞$ before any of
those parameters grow.



We now have Gaussian moments at sampled edges. Two obstacles remain.
First, all of the entropy argument concerns a graph on which we
already know the local law and rigidity; those are statements about
smoothed eigenvalue counts, and at the end we must pass to an
unsmoothed statement at the single energy $E$. Second, moments of
eigenvector entries are not yet a point process. This section solves
both: it constructs the microscopic limit at $E$ and identifies its
unmarked marginal with the Kesten--McKay density in the correct
microscopic coordinates.

The only tool is an exact symmetry of the uniform graph ensemble, the
two-edge switch, which we turn into a self-consistency identity for the
resolvent at a fixed height. It is worth saying what that identity says
in one line:

  average over edges of  G_{vw}
   = (1+zeta m_n(zeta))/(d) + small,

where $m_n(zeta)=frac1nTr(A-zeta)^{-1}$. The right side is the
Kesten--McKay (tree) formula with the actual $m_n$ inserted
rather than the deterministic $m_d$. Such an identity, valid at every
fixed height, is exactly what a stochastic Airy-type argument normally
needs and normally lacks; it is available here because the switch is an
exact graph operation, not an approximation.


Fix four distinct vertices $v,w,x,y$ such that $(v,x)$ and $(w,y)$ are
edges and $(v,w)$, $(x,y)$ are not. The switch deletes the two
edges and inserts $(v,w)$ and $(x,y)$. The result is again a simple
$d$-regular graph on the same vertex set: degrees are preserved, the new
edges join previously non-adjacent pairs, and the deleted edges are
distinct. Moreover the switch is invertible, the inverse switch being
the same operation applied to the new edge pair; hence the switch is an
involution of the finite set of simple $d$-regular graphs, and a uniform
graph pushed forward by the switch is again uniform. In particular, for
every function $F$ of the graph alone,

  EE F(A)=EE F(A+Z),

where $Z$ is the matrix of the switch. If the switch is not defined at a
sampled quadruple, both sides are interpreted by applying the identity
to the identity operation; the undefined event has probability bounded
away from one and its contribution is controlled by the same error
estimates as below. Exactness holds on the event that the sampled pair
really is switchable.

The defining convenience of the switch is that $Z$ has rank two in the
following explicit sense. Put
[
  u=e_v-e_y,qquad v=e_w-e_x,qquad
  Z=uv^{mathsf T}sfT+vu^{mathsf T}sfT.
]
Expanding,
[
  Z=e_ve_w^{mathsf T}sfT+e_we_v^{mathsf T}sfT+e_xe_y^{mathsf T}sfT
    +e_ye_x^{mathsf T}sfT
    -e_ve_x^{mathsf T}sfT-e_xe_v^{mathsf T}sfT
    -e_we_y^{mathsf T}sfT-e_ye_w^{mathsf T}sfT,
]
so $Z$ deletes $(v,x)$ and $(w,y)$ and inserts $(v,w)$ and $(x,y)$,
as required. Because $Z$ has rank at most two, the resolvent equation
has an explicit elementary solution.

[Rank-two resolvent identity]
Let $W=A-zeta$ with $Imzeta>0$, $G=W^{-1}$, $G'=(W+Z)^{-1}$, and
[
  A_{uv}=v^{mathsf T}sfTGu,qquad A_{uu}=u^{mathsf T}sfTGu,qquad
  A_{vv}=v^{mathsf T}sfTGv,qquad
  D=(1+A_{uv})^2-A_{uu}A_{vv}.
]
Then $Dneq0$ and

  G'_{xy}
  &=G_{xy}-D^{-1}Bigl[(G_{xv}-G_{xy})
      [(1+A_{uv})(G_{yw}-G_{yx})-A_{vv}(G_{yv}-G_{yy})]
      notag
  &hspace{22mm}+(G_{xw}-G_{xx})
      [(1+A_{uv})(G_{yv}-G_{yy})-A_{uu}(G_{yw}-G_{yx})]
      Bigr],

and the same formula with $(x,y)$ replaced by $(v,w)$ throughout gives
$G'_{vw}$.


Proof. 
From $(W+uv^{mathsf T}sfT+vu^{mathsf T}sfT)G'=I$ we get
[
  G'=G-Guv^{mathsf T}sfTG'-Gvu^{mathsf T}sfTG' .
]
Set $X=v^{mathsf T}sfTG'$ and $Y=u^{mathsf T}sfTG'$, row vectors. Multiplying
the displayed identity on the left by $v^{mathsf T}sfT$ and by
$u^{mathsf T}sfT$ gives the linear system
[
  1+A_{uv}&A_{vv} A_{uu}&1+A_{uv}
  X Y
  =v^{mathsf T}sfTG u^{mathsf T}sfTG,
]
whose determinant is $D$. Since $Im G$ is negative definite and $Z$ is
symmetric with zero exponential type, $W+Z$ has its spectrum in the open
upper half-plane, so $G'$ exists; consistency of the system gives
$Dne0$ (otherwise $W+Z$ would have a null vector with real
eigenvalue). Cramer's rule gives
[
  X=D^{-1}[(1+A_{uv})v^{mathsf T}sfTG-A_{vv}u^{mathsf T}sfTG],quad
  Y=D^{-1}[(1+A_{uv})u^{mathsf T}sfTG-A_{uu}v^{mathsf T}sfTG].
]
Now take the $(x,y)$ entry of $G'=G-Guv^{mathsf T}sfTG'-Gvu^{mathsf T}sfTG'$:
[
  G'_{xy}=G_{xy}-(Gu)_x(e_y^{mathsf T}sfTv^{mathsf T}sfTG')
                  -(Gv)_x(e_y^{mathsf T}sfTu^{mathsf T}sfTG'),
]
and $e_y^{mathsf T}sfTv^{mathsf T}sfTG'=(G'v)_y=(X)_y$,
$e_y^{mathsf T}sfTu^{mathsf T}sfTG'=(Y)_y$, while
$(Gu)_x=G_{xv}-G_{xy}$ and $(Gv)_x=G_{xw}-G_{xx}$ because
$u=e_v-e_y$ and $v=e_w-e_x$. Substituting the expressions for $X,Y$
proves (eq:switch-entries).



The local law already controls the trace $m_n$ and the adjacent resolvent
entry along an edge. What it does not produce is a closed equation that
would pin down $m_n$ at a single microscopic energy: those estimates are
smoothed, so they see averages over windows and cannot resolve individual
spacings. A two-edge switch fills that gap. It is an exact symmetry of
the uniform graph ensemble, and Woodbury's formula converts the graph
identity into a matrix relation among resolvent entries at the sampled
endpoints. The next subsection extracts the scalar equation from that
relation.


Conditionally on the graph, sample an infinite sequence of independent
uniform oriented edges. For fixed $kge1$, let $W_{n,k}$ be the $ntimes2k$
matrix whose columns are the coordinate vectors of the first $k$
starting vertices followed by their $k$ ending vertices. Define

 M_{n,k}(z)=W_{n,k}^{mathsf T}sfT(A-E-frac zn)^{-1}W_{n,k}
 =sum_ifrac{b_ib_i^{mathsf T}sfT}{x_i-z},qquad
 b_i=W_{n,k}^{mathsf T}sfTu_i,quad Im z>0 .

The normalization $|u_i|^2=n$ accounts for the missing factor $n^{-1}$.
We show that, at any fixed microscopic height, the scalar $m_n$ and the
matrix $M_{n,k}$ concentrate on their tree values $m_*$ and $K_k(m_*)$.

[Endpoint control]
There is a deterministic $xi>0$, depending only on $d,E$, such that for
every fixed $kge1$ and $D>0$ there are constants $C_{k,D}<∞$ and
$y_{k,D}<∞$ with the following property. For every fixed $sinR$
and every fixed $yge y_{k,D}$,

 limsup_{nto∞}PP(
 |m_n(s+mathrm iy)-m_*|rmiy)-m_*|
 +|M_{n,k}(s+mathrm iy)-K_k(m_*)rmiy)-K_k(m_*)|>y^{-xi})le C_{k,D},y^{-D}.

The constants and the threshold do not depend on $s$, and $s+mathrm iy$rmiy$
is fixed before the limit.


Proof. 
It suffices to treat $kge2$, since smaller matrices are compressions;
throughout $k,s,y$ are fixed and error exponents may decrease from one
occurrence to the next.

Initial estimates. At $z_0=s+mathrm iS$,rmiS$,
Proposition~prop:spectral-input and Lemma~lem:tree-identities
give, for some $c>0$,

 m_n(z_0)=m_*+O(S^{-c}),quad M_{n,k}(z_0)=K_k(m_*)+O(S^{-c}),quad
 Im m_n(mathrm iS^2)+rmiS^2)+|Im M_{n,k}(mathrm iS^2)rmiS^2)|=O(1).

For the off-block entries, independent uniform vertices $v,w$ satisfy
$EE_{v,w}|R_{vw}|^2=n^{-2}Tr R^*R=Im m_n(z_0)/S$, so Markov bounds
them by $O(S^{-c})$. Endpoints of distinct sampled edges are independent
uniform vertices, and within one edge the typical-star estimate gives the
adjacent entry $(1+Em_*)/d$. Also $#{|x_i|le S^2}=O(S^2)=n^{o(1)}$.

Concentration under a trace bound. Descend from height $S$ to $y$
by halving the imaginary part at fixed real part $s$; let $mathcal Z_n$calZ_n$
be the $O(log S)$ dyadic arguments. For every fixed $C<∞$, outside
an event of limiting probability $O_{k,D,C}(y^{-D})$, simultaneously for
$zinmathcal Z_n$,calZ_n$,

 |m_n(z)|le C ⇒ M_{n,k}(z)=K_k(m_n(z))+O((Im z)^{-c_2}),

with $c_2>0$ independent of $D,k,s,C$. To see this, compare with $z_0$:
for $|x|>S^2$ and large $S$,
$|(x-z)^{-1}-(x-z_0)^{-1}|le CSx^{-2}le CS^{-1}Im(x-mathrm iS^2)^rmiS^2)^{-1}$,
so the outside-band parts are $O(S^{-1})$ (positive residues make the
matrix bound immediate). Inside the band, with
$w_i(z)=(x_i-z)^{-1}-(x_i-z_0)^{-1}$, a bounded trace gives
$sum_{|x_i|le S^2}|w_i(z)|^2le 2Im m_n(z)/Im z+2Im m_n(z_0)/Sle C'/Im z$.
Replacing the marks by independent $N(0,C_k)$ vectors, each entry of
$sum(b_ib_i^{mathsf T}sfT-C_k)w_i(z)$ has $2j$-th absolute moment
$O_{j,k,C}((Im z)^{-j})$ (any spectral index occurring once has zero
expectation; the remaining partitions are bounded by Gaussian moments and
the weight squares). Conditionally on a graph satisfying the fixed-order
moment estimates, the same holds up to polynomial error, with only
$n^{o(1)}$ terms. Pick $0<c_2<min(c,1/4)$; Markov at threshold
$(Im z)^{-c_2}$ bounds the failure by $(Im z)^{-j(1-2c_2)}$, and summing
over the dyadic heights gives $C_{j,k,C}y^{-j/2}+o(1)$, so large $j$
proves the claim.

A switch gives the scalar equation. Apply the cross-switch that
exchanges the ending vertices of the first two sampled edges: the pairs
$(v_1,w_1)$ and $(v_2,w_2)$ become $(v_1,w_2)$ and $(v_2,w_1)$. This is
a different involution from Lemma~lem:switch-identity, but the
Woodbury computation below applies to it identically. With probability
$1-O_k(n^{-1})$ the $2k$ endpoints are distinct and induce exactly the
sampled matching; on that event the switch is a measure-preserving
involution of the restricted graph/list space. The rank-two update changes the eigenvalue counting
function by a bounded amount, so

 |widehatwidehat m_n(z)-m_n(z)|m_n(z)-m_n(z)|le C/Im z .

Write the first two pairs as $(v_1,w_1),(v_2,w_2)$ and set
$f_1=(e_{v_1}-e_{v_2})/sqrt2$, $f_2=(e_{w_1}-e_{w_2})/sqrt2$,
$F=(f_1,f_2)$, $B=0&-2-2&0$; the switch adds
$FBF^{mathsf T}sfT$, and the resolvent identity reads
$(I+UB)widehatwidehat U=U$U=U$ with $U=F^{mathsf T}sfT(A-E-z/n)^{-1}F$. Assuming
$|m_n(z)|le C$ and applying (eq:lim-concentration-implication) and
(eq:lim-rank-trace) to both data sets,

 U=m&qq&m+O((Im z)^{-c_2}),quad
 widehatwidehat U=m&-qU=m&-q-q&m+O((Im z)^{-c_2}),
 quad q=(1+Em)/(d),quad m=m_n(z),

the sign change in $widehatwidehat U$U$ coming from expressing the updated pairs
in the old columns. The off-diagonal of $(I+UB)widehatwidehat U=U$U=U$ gives

 q^2-q-m^2=O((Im z)^{-c_2}).

As a quadratic in $m$ this has simple conjugate roots in the strict bulk;
by Lemma~lem:tree-identities the upper root is $m_*$ and the lower
stays a fixed distance from ${Im mge0}$. Hence, in the bounded range,

 m_n(z)=m_*+O((Im z)^{-c_2}).


Closing the descent. No trace bound is assumed at lower steps.
Instead, for $eta/2leeta'leeta$, term by term
$|x-s-mathrm irmieta'|ge(eta'/eta)|x-s-mathrm irmieta|$, so

 |m_n(s+mathrm irmieta')-m_n(s+mathrm irmieta)|
 le(eta-eta')/(eta')Im m_n(s+mathrm irmieta)leIm m_n(s+mathrm irmieta).

Choose a fixed trace range containing the initial value and every
next-step value of (eq:lim-halving-bound) whenever the preceding
value is within $1$ of $m_*$, enlarged once more for the switched trace
via (eq:lim-rank-trace); take $y$ large enough that
(eq:lim-scalar-stability) returns within $1$ of $m_*$. Induction
down the heights proves the trace bounds and the endpoint estimates at
every step. Absorbing fixed constants into a smaller $xi$ proves
(eq:lim-endpoint-control); the only $s$-dependence is the harmless
requirement $Sge C(|s|+y+1)$, which holds before the limit.



For a locally finite point measure $nu=sum_jdelta_{x_j}$, write $PV$
for the symmetric truncation $lim_{Rto∞}sum_{|x_j|le R}$.

[Microscopic limits]
From every sequence of admissible sizes tending to infinity one can
extract a subsequence along which
[
 (nu_n,m_n,(M_{n,k})_{kge1})⇒(nu,m,(M_k)_{kge1})
]
jointly, in the local vague topology for point measures and the
compact-open topology for analytic functions. The limit $nu=sum_jdelta_{x_j}$
is locally finite, with multiplicities permitted, and there are
deterministic $xi,chi>0$ such that the following hold.

For every fixed $0<delta<1$, almost surely for all sufficiently large $R$,

 m(e+mathrm iu)=m_*+O(u^rmiu)=m_*+O(u^{-xi}),qquad |e|le R, R^deltale ule R^2,

and

 nu([v,w])=rho(w-v)+O(R^{1-chi}),qquad -Rle vle wle R, Rto∞ .

On a possible extension of the probability space, for every fixed $k$
there are standard Gaussian vectors $g_jinR^{2k}$, independent over
$j$ conditionally on $nu$, such that

 m(z)=m_{rm r}+PVsum_jfrac1{x_j-z},qquad
 M_k(z)=K_k(m_{rm r})+PVsum_jfrac{C_k^{1/2}g_jg_j^{mathsf T}sfTC_k^{1/2}}{x_j-z}.

These representations are consistent as $k$ grows by keeping the first
$k$ channels, a separate mark is assigned to each occurrence of a
repeated point, the principal values converge locally uniformly on $C_+$,
and $Im M_k(z)$ is positive definite there, almost surely simultaneously
for all finite $k$ and all $zinC_+$.


Proof. 
Tightness. For fixed $L,y>0$, positivity gives
$nu_n([-L,L])le(L^2+y^2)y^{-1}Im m_n(mathrm iy)$;rmiy)$; using
Proposition~prop:endpoint-control with large fixed $y$ (chosen per
tolerance) gives tightness of the local counts. For compact-open
tightness, fix a compact $mathcal KcalKsubsetC_+$ and $y>0$. For a resolvent
sum $f$ with nonnegative residues, comparing kernels on the real line
gives $Im f(z)le C_{mathcal K,y}calK,yIm f(mathrm iy)$rmiy)$ and
$|f(z)-f(mathrm iy)|rmiy)|le C_{mathcal K,y}calK,yIm f(mathrm iy)$rmiy)$ for
$zinmathcal K$;calK$; applying these to $m_n$ and, by polarization, to the
real quadratic forms of $M_{n,k}$ bounds each fixed compact, and Cauchy's
estimates give equicontinuity on smaller compacts. Extracting diagonally
over growing compacts and channel dimensions gives the joint limit;
local vague limits with locally bounded counts are point measures with
integer multiplicities, and the analytic limits have nonnegative
imaginary parts.

Conditional Gaussian marks on compact intervals. Choose nested
compact position intervals exhausting $R$ with deterministic boundaries
almost surely uncharged by $nu$ (possible since only countably many
points have positive probability of being atoms). On each interval
truncate the point count at a fixed integer, order the prelimit indices
using graph data alone, and extract positions and marks. By
Proposition~prop:edge-moments every fixed joint moment of the marks
is polynomially close to that of independent $N(0,C_k)$ vectors, uniformly
over the selected indices. Diagonal selection of the graph exceptional
events (choose $J_nto∞$ slowly enough that all orders at most $J_n$
match to $J_n^{-1}$ on a common event of probability $ge1-J_n^{-1}$; only
finitely many fixed-order assertions enter at each stage) upgrades
fixed-moment matching to conditional weak convergence to the product
Gaussian law, uniformly over bounded lists: a list failing to converge
would admit a deterministic sequence of graphs and index lists failing
while all fixed moments converged, contradicting tightness and Gaussian
moment determinacy. Passing to the truncation limit identifies the marking
law given $nu$, with covariance $C_k$ over the $2k$ endpoints; on an
extension of the probability space one assigns an independent mark
$b_j=C_k^{1/2}g_j$ to each occurrence of a repeated point.

Removing imaginary tails. For fixed $zinC_+$ and large $L$, the
scalar kernels satisfy, for $|x|>L$,
$Im(x-z)^{-1}+|(x-z)^{-1}-(x-mathrm i)^rmi)^{-1}|le C_zL^{-1}Im(x-mathrm iL)^rmiL)^{-1}$;
for matrix sums use positivity of the residues and bound the norm by the
trace of the positive imaginary part. By Proposition~prop:endpoint-control
these traces at $mathrm iL$rmiL$ are tight as $Lto∞$ after the upper
limit in $n$; the factor $L^{-1}$ removes the tails in probability, and
compact marked sums converge, so

 
 Im m(z)&=sum_jImfrac1{x_j-z},&
 m(z)-m(mathrm i)&=rmi)&=sum_j(frac1{x_j-z}-frac1{x_j-mathrm i}rmi),
 Im M_k(z)&=sum_jb_jb_j^{mathsf T}sfTImfrac1{x_j-z},&
 M_k(z)-M_k(mathrm i)&=rmi)&=sum_jb_jb_j^{mathsf T}sfT(frac1{x_j-z}-frac1{x_j-mathrm i}rmi).
 

No imaginary spectral mass is lost at infinity; real additive constants
are not yet excluded and are fixed below.

Uniform bounds after the limit. Pass
Proposition~prop:endpoint-control to the limit by the portmanteau
inequality on the open event of large error (its constants are independent
of $s$), and slightly decrease $xi<1$. Fix $0<delta<1$ and take a unit
grid on the padded region $|e|le2R$, $tfrac12R^deltale ule2R^2$,
giving $O(R^3)$ arguments. For fixed dyadic $R$ this is a finite
collection in the exact limit; the failure probability is at most
$CR^{3-delta D}$, and Borel--Cantelli over dyadic $R$ with large $D$
controls it. The kernel identities give $|m'(z)|leIm m(z)/Im z$, and
imaginary parts at comparable height and bounded distance are comparable,
so filling the grid costs $O(1/u)$; padding reaches all large $R$. A
countable set of $delta$'s gives every fixed-$delta$ conclusion. The
limit in $n$ preceded this growing-grid argument.

From smoothing to counts. Equation~(eq:mesoscopic-limit) gives
$nu([-R,R])=O(R)$. Integrating $pi^{-1}Im m(e+mathrm irmisqrt R)$ over
$ein[v,w]$ yields $rho(w-v)+O(R^{1-xi/2})$. Compare with the sharp
count using $P_{v,w}(x)=pi^{-1}int_v^w(sqrt R)/((x-e)^2+R)dd e$:
points within $R^{3/4}$ of an endpoint number $O(R^{3/4})$; on
$[-2R,2R]$ off the two strips, $|P_{v,w}(x)-mathbf 1_bf1_{[v,w]}(x)|=O(R^{-1/4})$,
for total error $O(R^{3/4})$; for $|x|>2R$ the bound
$P_{v,w}le CR^{3/2}/x^2$ with the linear count contributes
$O(sqrt R)$. This gives (eq:discrepancy) for any
$0<chi<min(1/4,xi/2)$.

Scalar principal values. Partial summation in (eq:discrepancy)
gives $PVsum_{|x_j|>R}x_j^{-1}=O(R^{-chi})$ and
$sum_{|x_j|>R}x_j^{-2}=O(R^{-1})$ (the leading density terms cancel in
the symmetric sum on the two half-lines). Hence
$S_nu(z)=PVsum_j(x_j-z)^{-1}$ is analytic and locally uniformly
convergent, and comparing with the constant density gives
$S_nu(mathrm iR)=rmiR)=mathrm irmipirho+o(1)$. By
(eq:lim-difference-identities), $m-S_nu$ is constant with zero
imaginary part; the two displays identify it as $m_{rm r}$.

Marked principal values. Fix $k$, condition on $nu$, and set
$Z_j=b_jb_j^{mathsf T}sfT-C_k$, independent centered with bounded second
moments. Ordered by increasing $|x_j|$, $sum_jZ_j/(x_j-mathrm i)$rmi)$
converges almost surely entrywise, and
$sum_{|x_j|>1}|Z_j|x_j^{-2}<∞$ conditionally almost surely, so
the difference series between $z$ and $mathrm i$rmi$ converges absolutely
and locally uniformly; with the scalar PV this gives local uniform
convergence of $S_{k,nu}(z)=PVsum_jb_jb_j^{mathsf T}sfT(x_j-z)^{-1}$.
The linear count gives $sum_j|x_j-mathrm iR|^rmiR|^{-2}=O(R^{-1})$, so each
entry of the centered marked series at $mathrm iR$rmiR$ has conditional
second moment $O(R^{-1})$, whence
$S_{k,nu}(mathrm iR)rmiR)tomathrm irmipirho,C_k$ in probability. Since
$M_k-S_{k,nu}$ is almost surely constant in $z$ and
$M_k(mathrm iR)rmiR)to K_k(m_*)$ in probability, that constant equals
$K_k(m_*)-mathrm irmipirho,C_k=K_k(m_{rm r})$.

Finally, (eq:discrepancy) supplies infinitely many points and $C_k$
is positive definite; $2k$ independent Gaussian marks span $R^{2k}$
almost surely, so the positive contributions to $Im M_k(z)$ make it
positive definite for every $zinC_+$, simultaneously over a countable
set of $k$.



We work with a microscopic subsequential law $nu$ supplied by
Proposition~prop:microscopic-limit. Its endpoint marks are exactly
Gaussian, while genuine graph switches preserve its law. Combining the two
gives exact orthogonal switching symmetries; composing many small
orthogonal switches then inserts a Gaussian matrix perturbation.

[In plain words]
The microscopic-limit step leaves us with the microscopic limiting law and
its Gaussian endpoint marks. Those marks give a symmetry that a general
law would not have: conjugating $M_k$ by an orthogonal $Q_U$ preserves
the law. Composing many small conjugations then produces an honest GOE
perturbation: the natural graph trace $mathcal S_r(z)$calS_r(z)$ is approximated,
in law and up to $o(1)$, by $m_{rm r} + Tr(H_r-z)^{-1}$, where $H_r$
is a fixed diagonal plus independent GOE noise.



Fix $kge2$, write $M=M_k$, and for $Oin O(k)$ set
$D_O=I_k&00&O$ and
$B_O=D_OJ_kD_O^{mathsf T}sfT-J_k$, which is symmetric. For an analytic
symmetric matrix function $M$ with positive definite imaginary part, put

 mathcal T_O(m,M)=(m-∂_zcalT_O(m,M)=(m-∂_zlogdet(I_{2k}+B_OM),,
 D_O^{mathsf T}sfT(M^{-1}+B_O)^{-1}D_O).

This is well defined: $M$ is invertible and
$Im(M^{-1})=-M^{-*}(Im M)M^{-1}<0$, so $M^{-1}+B_O$ is invertible and
its inverse has positive definite imaginary part. Only the logarithmic
derivative is used, so no branch of the logarithm is chosen.

[Continuous switching symmetry]
For every fixed $kge2$ and every $Oin O(k)$, the joint law of $(m,M_k)$
under $nu$ is invariant under $mathcal T_O$.calT_O$.


Proof. 
Let $O$ be a permutation matrix. On the event, of probability
$1-O_k(n^{-1})$, that the $2k$ sampled endpoints are distinct and induce
exactly the sampled matching, the map $Amapsto A+WB_OW^{mathsf T}sfT$,
$Wmapsto WD_O$ permutes the pairing between the two endpoint lists, is a
bijection of the restricted graph/list space preserving all weights, and
is inverted by the inverse permutation. The determinant and inverse
identities for finite-rank perturbations give
$widehatwidehat m_n=m_n-∂_zm_n=m_n-∂_zlogdet(I_{2k}+B_OM_{n,k})$ and
$widehatwidehat M_M_{n,k}=D_O^{mathsf T}sfT(M_{n,k}^{-1}+B_O)^{-1}D_O$ (the factor
$n^{-1}$ in $m_n=n^{-1}Tr R$ supplies the derivative). Compact-open
convergence of analytic functions gives convergence of derivatives, and
positive-definite imaginary parts make the inverse regular, so passing to
the limit proves the claim for permutation matrices.

For $Uin O(k)$ put $Q_U=diag(U,U)$. As $C_k$ and $K_k(m_{rm r})$
commute with $Q_U$, the Gaussian representation
(eq:marked-representation) shows that
$mathcal R_U(m,M)=(m,Q_U^calR_U(m,M)=(m,Q_U^{mathsf T}sfTMQ_U)$ preserves the joint law. Since
$Q_UD_OQ_U^{mathsf T}sfT=D_{UOU^{mathsf T}sfT}$ and
$Q_UB_OQ_U^{mathsf T}sfT=B_{UOU^{mathsf T}sfT}$, we get
$mathcal R_U^calR_U^{-1}mathcal T_OcalT_Omathcal R_U=calR_U=mathcal T_calT_{UOU^{mathsf T}sfT}$;
so the law is invariant under every orthogonal conjugate of a permutation
switch. The maps satisfy $B_{OO'}=B_O+D_OB_{O'}D_O^{mathsf T}sfT$ and
$D_{OO'}=D_OD_{O'}$, so $mathcal T_OcalT_Omathcal T_calT_{O'}=mathcal T_calT_{OO'}$
(the inverse updates compose and the determinant ratios multiply before
differentiation). A transposition is a hyperplane reflection; its
orthogonal conjugates are all hyperplane reflections, and finite products
generate $O(k)$.



Fix a discrepancy exponent $0<chi<1/4$ from (eq:discrepancy). Use
$r$ planar rotations through a small angle $theta$, with aggregate
variance of order $T=rtheta^2$. A cutoff $K$ marks the finite block that
will receive GOE noise, and a much larger temporary cutoff $L$ first turns
the principal-value resolvent into a finite matrix. Choose once and for all

 0<tau<(chi)/(100000),quad k=2r,quad K=r^{1/100},quad T=r^{tau},
quad theta=√(T/r),quad L=r^{100/chi},quad rinN .

For each $r$ let $O=O_r$ be the direct sum of $r$ planar rotations of angle
$theta$, write $B=B_{O_r}$ and $M^0=K_k(m_{rm r})$. The real symmetric
matrix

 Gamma_r=C_k^{1/2}(B^{-1}+M^0)^{-1}C_k^{1/2}

is well defined for large $r$: the singular values of $B$ are all
$2|sin(theta/2)|$ and $|M^0|$ is bounded independently of $r$. Let
$beta$ be a diagonal matrix of eigenvalues of $Gamma_r$ and set

 s_1=Trbeta,qquad s_2=Trbeta^2 .

Index the limiting positions with multiplicity, put
$I={j:|x_j|le K}$, $q_*=|I|$ and $D_I=diag(x_j:jin I)$. Conditionally
on $nu$, let $mathcal W$calW$ be a raw GOE matrix of size $q_*$ (variances $1$
off the diagonal and $2$ on it), independent of $nu$ apart from its
dimension, and set

 H_r=D_I+s_1Id_{q_*}+√(s_2),mathcal WcalW .


[Gaussian insertion]
The coefficients satisfy

 |beta|le Ctheta,qquad cTle s_2le CT,qquad |s_1|le CT

with deterministic positive constants. As $rto∞$,

 m_{rm r}+Tr(H_r-z)^{-1} ⇒ m(z)

jointly at every fixed finite list of $zinC_+$, and the unscaled
eigenvalue point measure of $H_r$ converges locally vaguely in law to
$nu$. More precisely, generating the Gaussian marks of $M_k$
conditionally on $nu$ separately for each $r$, and writing
$mathcal S_r(z)=m(z)-∂_zcalS_r(z)=m(z)-∂_zlogdet(I_{2k}+BM_k(z))$, expectations of
every bounded smooth function with bounded derivatives of finitely many
real and imaginary trace values differ by $o(1)$ when
$m_{rm r}+Tr(H_r-z)^{-1}$ is replaced by $mathcal S_r(z)$.calS_r(z)$.


First cut at $L$. Elementary consequences of (eq:discrepancy)
are

 nu([-R,R])=2rho R+O(R^{1-chi}),quad
 |PVsum_{|x_j|>R}frac1{x_j}|=O(R^{-chi}),quad
 sum_{|x_j|>R}frac1{x_j^2}=O(R^{-1}),

the last two by partial summation on the two half-lines. Cut
[
  M_L(z)=M^0+sum_{|x_j|le L}C_k^{1/2}g_jg_j^top C_k^{1/2}(x_j-z)^{-1}.
]
For fixed $zinC_+$,

 |M-M_L|+|M'-M_L'|=O_{PP}(r^2L^{-chi}),qquad |M'|=O_{PP}(r^2),

since the mean tail is $O(L^{-chi})C_k$ (and $O(L^{-1})C_k$ for the
derivative), each centered entry of the Gaussian tail has variance
$O(L^{-1})$, and there are $(4r)^2$ entries. Keeping only $|x_j|le r^4$
(of which there are $sim2rho r^4$) and using Gaussian fourth moments
gives

 Im M(z)ge r^{-5}Id_{4r},qquad Im M_L(z)ge r^{-5}Id_{4r}

with probability tending to one (using $Lgg r^4$). As $B^{-1}$ is real
symmetric, this bounds $|(B^{-1}+M)^{-1}|$ and $|(B^{-1}+M_L)^{-1}|$ by
$r^5$, so cutting both scalar and marked sums in
$mathcal S_r(z)=m(z)-calS_r(z)=m(z)-Tr(B^{-1}+M(z))^{-1}M'(z)$ costs $o_{PP}(1)$: the
marked error is $O_{PP}(rcdot r^5cdot r^{-98})+O_{PP}(rcdot r^{10}
cdot r^{-98}cdot r^2)=O_{PP}(r^{-85})$, and the scalar tail is
$O(L^{-chi})$. Diagonalizing $Gamma_r$ and rotating the independent
Gaussian row vectors in channel space gives a matrix $mathcal G$calG$ of
independent standard Gaussian entries with $4r$ columns, and the
finite-rank inverse identity (Woodbury) yields

 mathcal S_r(z)=m_calS_r(z)=m_{rm r}+Tr(D+A'-z)^{-1}+o_{PP}(1),qquad
 A'=mathcal GcalGbetamathcal G^calG^top,quad D=diag(x_j:|x_j|le L).

The coefficient bounds follow from
$(B^{-1}+M^0)^{-1}=(Id+BM^0)^{-1}B=B+O(theta^2)$: every singular value
is comparable to $theta$, so $|beta|=O(theta)$ and $s_2asymp
rtheta^2=T$; and $s_1=Tr(C_kB)+O(rtheta^2)$ with
$Tr(C_kB)=2p,Tr(O-Id_k)=O(rtheta^2)$.

Removing the outer positions. Let $F={j:K<|x_j|le L}$ and put
$h=T^2$, $U_I=(|D_I|+h)^{-1/2}$, $U_F=|D_F|^{-1/2}$. Partial summation
using (eq:sw-tail-sums) gives

 Tr U_I^2+Tr U_F^2=O(log r),qquad |U_I^2|le h^{-1},quad
 |U_F^2|le K^{-1}.

In the estimates below, $polylog(r)$ denotes a fixed power
of $log r$ whose exponent may rise from one occurrence to the next.

[Weighted block estimates]
For $mathcal AcalAin{I,F}$ put $a_I=h$, $a_F=K$. With probability tending
to one,

 bigl|U_{mathcal A}calA(A'_{mathcal AcalAmathcal A}calA-s_1Id)U_{mathcal A}calAbigr|
 lepolylog(r)(sqrt{T/a_{mathcal A}calA}+theta),

and, for $P_0=U_IA'_{IF}U_F$,

 |P_0|_{HS}^2le Tpolylog(r).



Proof. 
The centered matrix in (eq:sw-bernstein) is a sum of independent
self-adjoint $sum_{alpha=1}^{4r}beta_{alphaalpha}(v_alpha
v_alpha^top-P)$ with $P=U_{mathcal A}calA^2$ and $v_alpha$ Gaussian of
covariance $P$. Gaussian fourth moments give
$EE[(v_alpha v_alpha^top)^2]=(Tr P)P+2P^2$. Truncating each Gaussian
column at $|entry|lelog r$ (a superpolynomially small event, with
superpolynomially small centering correction) leaves summands of norm at
most $C|beta|(log r)^2Tr Plethetapolylog(r)$, and
bounds the matrix variance parameter by $Cs_2(Tr P)|P|le
(T/a_{mathcal A}calA)polylog(r)$. The bounded matrix Bernstein
inequality gives a two-sided tail $2|mathcal A|calA|exp(-ct^2/
(sigma_{rm mat}^2+bt))$, and a large power of $log r$ in $t$
proves (eq:sw-bernstein). For the cross block, independence of
distinct Gaussian rows gives $EE|P_0|_{HS}^2=s_2(Tr U_I^2)(Tr U_F^2)$,
so Markov with (eq:sw-weight-bounds) gives (eq:sw-cross-block).


Fix $zinC_+$ and write $G_F=(D_F+A'_{FF}-z)^{-1}$,
$G_I=(D_I+A'_{II}-z)^{-1}$.

[Far block and Schur correction]
Let $delta_F=polylog(r)((1+T)/K+√(T/K)+theta)$.
Then $delta_F=o(1)$ and, with probability tending to one,

 |U_F^{-1}(G_F-D_F^{-1})U_F^{-1}|ledelta_F=o(1),qquad
 Tr G_F=o_{PP}(1),

and, for $Σ=A'_{IF}G_FA'_{FI}$,

 |U_IΣ U_I|=o_{PP}(T^{-5}).



Proof. 
Conjugating by $U_F$,
$U_F(D_F+A'_{FF}-z)U_F=sign(D_F)+U_F(A'_{FF}-z)U_F$, and the
second term has norm at most $delta_F=o(1)$ by
Lemma~sw:weighted-blocks, $|s_1|=O(T)$ and (eq:sw-weight-bounds);
a Neumann expansion proves (eq:sw-far-resolvent), so
$Q_F=U_F^{-1}G_FU_F^{-1}$ has bounded norm. Also
$|Tr(G_F-D_F^{-1})|ledelta_FTr U_F^2=o(1)$, so $Tr G_F=o_{PP}(1)$.

For $Σ$, replacing $G_F$ by $D_F^{-1}$ costs at most
$|P_0|_{HS}^2delta_Fle Tpolylog(r)delta_F$ in the
$U_I$-weighted norm. For $A'_{IF}D_F^{-1}A'_{FI}$, condition on the inner
block and set $Q=mathcal G_IcalG_Ibeta^2mathcal G_I^calG_I^top$; the columns $a_j$
of $A'_{IF}$, $jin F$, are independent centered Gaussian with covariance
$Q$, so $EE[sum_{jin F}a_ja_j^top/x_jmidcdot]=Qsum_{jin F}x_j^{-1}$
and the entrywise variance is $(Q_{aa}Q_{bb}+Q_{ab}^2)sum_{jin F}x_j^{-2}$.
With $EE[Q_{aa}^2]=s_2^2+2Trbeta^4=O(T^2)$ and $|Q_{ab}|^2le Q_{aa}Q_{bb}$,
(eq:sw-tail-sums) and (eq:sw-weight-bounds) give
$EE|U_IA'_{IF}D_F^{-1}A'_{FI}U_I|_{HS}^2le C T^2(log r)^2
(K^{-2chi}+K^{-1})$, so by Markov this norm is
$Tpolylog(r)(K^{-chi}+K^{-1/2})$ with probability tending to
one. Combining and multiplying by $T^5$, the exponents
$6tau-chi/100$, $tfrac{13}2tau-tfrac1{200}$ and
$tfrac{13}2tau-tfrac12$ are all negative under (eq:sw-parameters),
giving (eq:sw-schur-small).


[Inner resolvent and trace cut]
For every fixed $zinC_+$,

 |U_I^{-1}G_I(z)U_I^{-1}|=O_{PP}(h),

and

 Tr(D+A'-z)^{-1}=Tr(D_I+A'_{II}-z)^{-1}+o_{PP}(1),

jointly for any fixed finite list of $z$'s.


Proof. 
Every diagonal entry of $U_I(D_I-mathrm ih)U_I$rmih)U_I$ has modulus at least
$1/sqrt2$, and the added term has norm $o(1)$ by
Lemma~sw:weighted-blocks, $|s_1|=O(T)$ and (eq:sw-weight-bounds).
Hence $|U_I^{-1}G_I(mathrm ih)U_I^rmih)U_I^{-1}|le C$, and
$G_I(mathrm ih)^*G_I(rmih)^*G_I(mathrm ih)=h^rmih)=h^{-1}Im G_I(mathrm ih)$rmih)$ gives
$|G_I(mathrm ih)U_I^rmih)U_I^{-1}|le Ch^{-1/2}$. The twice-iterated resolvent
identity $G_I(z)=G_I(mathrm ih)+(z-rmih)+(z-mathrm ih)G_I(rmih)G_I(mathrm ih)^2+rmih)^2+
(z-mathrm ih)^2G_I(rmih)^2G_I(mathrm ih)G_I(z)G_I(rmih)G_I(z)G_I(mathrm ih)$,rmih)$, conjugated by
$U_I^{-1}$, bounds the three terms by $C$, $C|z-mathrm ih|/h$rmih|/h$ and
$C|z-mathrm ih|^2/(hrmih|^2/(hIm z)$, proving (eq:sw-inner-resolvent).

For (eq:sw-second-cut), the upper-left block of the full inverse is
$widetilde G_I=(G_I^{-1}-Σ)^{-1}$. With $R_I=U_I^{-1}G_IU_I^{-1}$ and
$S_I=U_IΣ U_I$, a Neumann expansion using
(eq:sw-inner-resolvent) and (eq:sw-schur-small) gives
$|U_I^{-1}(widetilde G_I-G_I)U_I^{-1}|le O(h^2)o(T^{-5})=o(T^{-1})$, so
$Trwidetilde G_I-Tr G_I=o_{PP}(1)$. The full block-inverse trace is
$Trwidetilde G_I+Tr G_F+Tr(widetilde G_IA'_{IF}G_F^2A'_{FI})$; the
middle term is $o_{PP}(1)$, and with
$Q_FU_F^2Q_F=U_F^{-1}G_F^2U_F^{-1}$ of norm $O(K^{-1})$ the last term is at
most $|U_I^{-1}widetilde G_IU_I^{-1}||P_0|_{HS}^2|Q_FU_F^2Q_F|
=polylog(r)O(T^3/K)=o_{PP}(1)$.


[Gaussian replacement]
Conditional on almost every $nu$, bounded smooth tests with bounded
derivatives of finitely many trace values have expectations differing by
$o(1)$ between $D_I+A'_{II}$ and
$H_r=D_I+s_1Id_{q_*}+√(s_2),mathcal W$.calW$.


Proof. 
Write $A'_{II}-s_1Id_{q_*}=sum_{alpha=1}^{4r}beta_{alphaalpha}
(g_alpha g_alpha^top-Id_{q_*})$ with independent standard Gaussian
$g_alpha$. Since $Cov(g_ag_b,g_cg_d)=delta_{ac}delta_{bd}+delta_{ad}
delta_{bc}$, each $g_alpha g_alpha^top-Id$ and a raw GOE matrix have
equal mean (zero) and entry covariances; replace the summands one at a time
by $beta_{alphaalpha}mathcal W_calW_alpha$ with independent raw GOE
$mathcal W_calW_alpha$ (negative coefficients allowed, both covariances scaling
by $beta_{alphaalpha}^2$). For either kind of summand $X$,
$EE|X|^3le Ctheta^3q_*^3$ (Gaussian moments of $|g|^2$; Hilbert--
Schmidt norm for GOE). Differentiating the resolvent in a symmetric
direction $Delta$ gives trace derivatives bounded by $C_{z,j}q_*|Delta|^j$,
so for bounded smooth $Phi$ the third directional derivative is at most
$C_{z,Phi}q_*^3|Delta|^3$ at every real symmetric intermediate matrix.
Taylor expansion to second order cancels means and covariances in each
replacement, and the summed third-order remainder is
$C_{z,Phi}q_*^6rtheta^3=O(r^{-0.44+(3/2)tau})=o(1)$. Finally
$sum_alphabeta_{alphaalpha}mathcal W_calW_alpha$ has the law
$√(s_2),mathcal W$.calW$.


Proof. [Completion of Proposition~prop:gaussian-insertion]
The coefficient bounds were proved after (eq:sw-finite-matrix). That
equation, Lemma~sw:trace-cut and Lemma~sw:lindeberg prove the
conditional comparison with $mathcal S_r$.calS_r$. Bounded convergence averages
the conditional error, and by Proposition~prop:continuous-switching
the unconditional law of $mathcal S_r$calS_r$ is the law of $m$ for every $r$;
this proves (eq:sw-trace-limit).

For the point measures, let $nu^{(r)}$ be the eigenvalue counting measure
of $H_r$ and $sigma_r(dd x)=nu^{(r)}(dd x)/(1+x^2)$; its total mass
$ImTr(H_r-mathrm i)^rmi)^{-1}$ is tight by (eq:sw-trace-limit), so
$sigma_r$ is tight on the one-point compactification $overlineR$. For
$z=e+mathrm irmieta$ the kernel $k_z(x)=(1+x^2)eta((x-e)^2+eta^2)^{-1}$
extends continuously to $overlineR$ with $k_z(∞)=eta$, and its
integral against $sigma_r$ is $ImTr(H_r-z)^{-1}$; so any subsequential
limit $sigma$ has the law of $Im m$ at a countable dense set of $z$,
whence at all $z$ by continuity. For $fin C_c(R)$, integrating $k_z$
against $pi^{-1}f(e)dd e$ and letting $eta↓0$ gives $(1+x^2)f(x)$
on $R$ (Poisson identity), with uniformly bounded integrals and value
$etapi^{-1}int fto0$ at infinity, so dominated convergence applies;
applying the same inversion to $Im m(z)$ gives $int f,ddnu$. Every
subsequential local vague limit of $nu^{(r)}$ therefore has the law of
$nu$.



Proposition~prop:gaussian-insertion approximates every microscopic
subsequential law by the eigenvalues of a finite diagonal matrix plus GOE
noise. We now apply fixed-energy universality to these finite matrices.

[In plain words]
The switching step reduces the microscopic limit to a finite
matrix $H_r=$ diagonal $+$ GOE noise. To identify its eigenvalue law, we
apply fixed-energy Dyson Brownian motion universality
(Landon--Sosoe--Yau): a diagonal matrix perturbed by a small GOE noise
has, in the bulk, the same local eigenvalue statistics as pure GOE,
provided the free-convolution comparison density at the origin is
$1/pi$. Lemma~dbm:regularity-density verifies the regularity
hypotheses, and the $c_r$ dilation adjusts that density to exactly
$1/pi$.



For a real diagonal $V=diag(v_1,...,v_m)$ write
$m_V(zeta)=m^{-1}sum_j(v_j-zeta)^{-1}$. For $g,G>0$ call $V$
$(g,G)$-regular if, for fixed $c,C,C_V>0$,

 cleIm m_V(e+mathrm irmieta)le Cquad (|e|le G, gleetale10),
 qquad |V|le m^{C_V},

with $g,G$ possibly $m$-dependent and, for a fixed $delta>0$,

 m^{delta-1}le gle m^{-delta},qquad Gle m^{-delta}.

Let $W_m$ be standard normalized GOE (off-diagonal variance $1/m$,
diagonal $2/m$). The free-convolution comparison density associated with
$V+sqrt t,W_m$ is described by

 m_{{rm fc},t}(zeta)=m_V(zeta+t,m_{{rm fc},t}(zeta)),qquad
 rho_{{rm fc},t}(e)=pi^{-1}lim_{eta↓0}Im m_{{rm fc},t}(e+mathrm irmieta),

using the upper-half-plane solution. Write $p_H^{(m)}$ for the unit-mass
symmetrized joint eigenvalue density of $H$, and $p_H^{(a)}$ for its
unit-mass $a$-point marginal.

[Fixed-energy DBM input]
Suppose a sequence of deterministic diagonal $V$ of dimension $m$
satisfies (eq:dbm-regularity) and (eq:dbm-parameter-domain)
with fixed constants. Fix $sigma>0$ and

 g,m^{sigma}le tle m^{-sigma}G^2,

and assume $rho_{{rm fc},t}(0)=1/pi$. Then for every fixed $age1$ and
$Oin C_c^∞(R^a)$ the difference of
$int O(boldsymbolalpha),p^{(a)}_{V+sqrt tW_m}(piboldsymbolalpha/m)ddboldsymbolalpha$
and $int O(boldsymbolalpha),p^{(a)}_{W_m}(piboldsymbolalpha/m)ddboldsymbolalpha$
tends to zero. Consequently the expected sums of $O$ over ordered tuples
of distinct points of the two eigenvalue processes, both rescaled by
$m/pi$, differ by $o(1)$.


Proof. 
This is Theorem~2.2 of Landon--Sosoe--Yau (with Definition~2.1 and
Equations~(2.5)--(2.9) there), specialized to energy zero and equal
densities. For the last assertion multiply each integral by
$(m)_a/(m/pi)^atopi^a$, the Jacobian and ordered-tuple factor for the
unit-mass convention.



Fix a microscopic subsequential law $nu$ from
Proposition~prop:microscopic-limit and the parameters and
conditional matrices of Proposition~prop:gaussian-insertion:
$K=r^{1/100}$, $T=r^{tau}$ with $0<tau<chi/100000$,
$I={j:|x_j|le K}$, $q_*=|I|$, and

 H_r=D_I+s_1Id+√(s_2),mathcal W_r,calW_r,qquad |s_1|le CT, cTle s_2le CT,

with $mathcal W_r$calW_r$ raw GOE of dimension $q_*$. By (eq:discrepancy),

 q_*=2rho K+O(K^{1-chi}),qquad T=q_*^{100tau+o(1)}.

Let $alpha_*=pirho$. Scaling by $alpha_*/q_*$,

 (alpha_*)/(q_*)H_rstackrel{law}{=}V_*+√(t_*),W_{q_*},
 qquad V_*=(alpha_*)/(q_*)(D_I+s_1Id),quad t_*=(alpha_*^2s_2)/(q_*).



For almost every $nu$, $V_*$ is eventually regular with
$g=T^{1/2}/q_*$, $G=q_*^{-1/8}$, satisfying
(eq:dbm-parameter-domain) and (eq:dbm-time-window) with
$m=q_*$, $delta=sigma=tau$. Moreover
$rho_{{rm fc},t_*}(0)to1/pi$; with $c_r=pirho_{{rm fc},t_*}(0)$ the
data $c_rV_*$, $c_r^2t_*$ satisfy Theorem~thm:dbm-input, including
the exact density condition.


Proof. 
$|V_*|$ is bounded since $T=o(K)$. Apply a slightly enlarged domain
$|Rezeta|le2G$, $g/2leImzetale20$; the normalized resolvent trace is
$m_{V_*}(zeta)=alpha_*^{-1}sum_{jin I}(x_j-(q_*zeta/alpha_*-s_1))^{-1}$.
Writing $e=q_*Rezeta/alpha_*-s_1$, $u=q_*Imzeta/alpha_*$ gives

 |e|=O(K^{7/8}+T)=o(K),qquad cT^{1/2}le ule CK.

Equation~(eq:mesoscopic-limit) (radius comparable to $K$, fixed
exponent below $50tau$) and the imaginary part of
(eq:marked-representation) give
$sum_j u((x_j-e)^2+u^2)^{-1}=pirho+o(1)$ uniformly in
(eq:dbm-original-domain); terms with $|x_j|>K$ sum to $O(u/K)$, and
the retained sum has positive lower bound $ge q_*u/(4K^2+u^2)ge c>0$ in
the complementary regime. This gives regularity. The powers are
$g=m^{-1+50tau+o(1)}$, $t_*=m^{-1+100tau+o(1)}$, $G=m^{-1/8}$, so
$t_*/(gm^tau)=m^{49tau+o(1)}to∞$ and
$t_*/(m^{-tau}G^2)=m^{-3/4+101tau+o(1)}to0$, verifying the time window
and domain for small $tau$.

For the density, $zeta=0$ in (eq:dbm-free-convolution) reads

 w=F_r(w),qquad F_r(w)=alpha_*^{-1}sum_{jin I}(x_j-(alpha_*s_2w-s_1))^{-1}.

On $|w-mathrm i|rmi|le1/4$ the argument $z=alpha_*s_2w-s_1$ has
$Im zasymp T$, $|z|=O(T)$; (eq:mesoscopic-limit) and
(eq:marked-representation) give
$PVsum_j(x_j-z)^{-1}=mathrm irmipirho+o(1)$, and
$PVsum_{|x_j|>K}(x_j-z)^{-1}=O(K^{-chi}+T/K)$. Hence $F_r(w)=mathrm i+o(1)$rmi+o(1)$
uniformly; Rouch'e for $w-F_r(w)$ against $w-mathrm i$rmi$ gives a unique
simple zero $w_rtomathrm i$rmi$ in the disk, continued by the implicit
function theorem to a positive-imaginary solution near $zeta=0$. That this
is the physical branch follows from the identity

 frac1{q_*}sum_{jin I}|(V_*)_{jj}-zeta-t_*w|^{-2}
 =(Im w)/(Imzeta+t_*Im w)<frac1{t_*}

for every upper-half-plane solution; two distinct solutions would give,
after subtracting and dividing, $1$ equal to a quantity of absolute value
strictly less than one by Cauchy--Schwarz and
(eq:dbm-uniqueness-identity). Its boundary value gives
$rho_{{rm fc},t_*}(0)=Im w_r/pito1/pi$.

Finally, replacing $V_*,t_*$ by $c_rV_*,c_r^2t_*$ dilates the
free-convolution law by $c_r$, with transform
$c_r^{-1}m_{{rm fc},t_*}(zeta/c_r)$ and density
$c_r^{-1}rho_{{rm fc},t_*}(0)=1/pi$ at zero; since $c_rto1$ the enlarged
domain gives (eq:dbm-regularity) and the strict power margins
preserve (eq:dbm-time-window) for $c_r^2t_*$.


By (eq:dbm-normalized-matrix) the rescaled eigenvalues supplied by
Theorem~thm:dbm-input are

 (q_*)/(pi)cdot(c_ralpha_*)/(q_*)lambda_j(H_r)=c_rrho,lambda_j(H_r),

with $c_r=1$ at exceptional initial values.


Let $mathcal G_m=calG_m=sum_{j=1}^mdelta_{(m/pi)lambda_j(W_m)}$. These
converge vaguely in law along all dimensions to a common locally finite
process $mathcal S_calS_{GOE}$: the eigenvalues of $sqrt m,W_m$ have
density $propto e^{-frac14sum y_j^2}prod_{i<j}|y_i-y_j|$, and
Theorem~1 of Valk'o--Vir'ag (for $beta=1$, center zero) gives vague
convergence after multiplying by $2sqrt m$; dividing by $2pi$ gives
$mathcal G_m$.calG_m$.


For every compact $JsubsetR$ and every $b>0$,
$sup_{mge1}EEexp(b,mathcal G_m(J))<∞$.calG_m(J))<∞$.


Proof. 
Bound the GOE count by twice a Gaussian-unitary count plus one via the
Gaussian decimation identity of Forrester--Rains. Let $mathrm O_m$rmO_m$ be the
decreasingly ordered points with density
$proptoDelta(y)prod_i b_0(y_i)$, $b_0=e^{-y^2/2}$, the law of
$sqrt m,W_m/sqrt2$; then $mathcal G_m(J)$calG_m(J)$ is the $mathrm O_m$rmO_m$ count in
$J_m=pi J/√(2m)$, so $|J_m|=O(m^{-1/2})$. Superposing independent
$mathrm O_m$rmO_m$ and $mathrm O_rmO_{m+1}$ gives density
$P(z)=sum_{|S|=m}Delta(z_S)Delta(z_{S^c})$, which by the
parity--Vandermonde factorization equals

 P(z)=c_m,Delta(z_1,z_3,...,z_{2m+1}),Delta(z_2,z_4,...,z_{2m}),

with $c_m>0$: every factor on the right divides $P$ (same-parity terms
vanish at coincidences; the others cancel in pairs with sign ratio
$(-1)^{j-i-1}=-1$), both sides have degree $m^2$, and the quotient is
positive. Integrating the odd points over their interlacing intervals with
the monic polynomials $p_j=(-1)^jb_0^{(j)}/b_0$ (each row replaced by its
partial row sums, last row $int b_0$, expansion along it) gives that the
even points have density

 proptoDelta(y)^2prod_i e^{-y_i^2}.


This density equals $(m!)^{-1}det(K_m(y_i,y_j))$ for the
projection kernel $K_m=sum_{j=0}^{m-1}psi_jpsi_j$ of the orthonormal
polynomials for $e^{-x^2}$, with order-$a$ factorial moment density
$det(K_m(y_i,y_j))_{i,j=1}^a$. The operator $-∂_x^2+x^2$
has eigenvalue $2j+1$ on $psi_j$, so $|f'|_2^2le(2m-1)|f|_2^2$ on the
range of $K_m$, and $|f(x)|^2le2|f'|_2|f|_2$ gives
$K_m(x,x)le2√(2m)$. Then for every $age0$,
$EE(N_m^{even})_ale(2√(2m),|J_m|)^ale C_J^a$, so
$EE e^{vN_m^{even}}leexp(C_J(e^v-1))$. As the
superposition count is at most twice its even count plus one and dominates
the $mathrm O_m$rmO_m$ count,
$EE e^{bmathcal G_m(J)}calG_m(J)le e^bEE e^{2bN_m^{even}}$.



Let $nu^{(r)}=sum_jdelta_{lambda_j(H_r)}$; by
Proposition~prop:gaussian-insertion,

 nu^{(r)}⇒nu

unconditionally in local vague topology. Fix $0le hin C_c^∞(R)$,
$f=1-e^{-h}$, and for a finite point measure $zeta=sum_idelta_{y_i}$ set
$A_0=1$, $A_a=(a!)^{-1}sum_{i_1,...,i_a distinct}prod_{j=1}^a f(y_{i_j})$,
$B_ell=sum_{a=0}^{ell}(-1)^aA_a$. The Bonferroni inequalities give

 B_{2ell+1}(zeta)leexp(-int h,ddzeta)le B_{2ell}(zeta),

seen by Bernoulli variables of success probabilities $f(y_i)$. Conditional
on almost every $nu$, let $zeta_r$ have coordinates $c_rrho,lambda_j(H_r)$;
Lemma~dbm:regularity-density and Theorem~thm:dbm-input with
test $prod_j f(alpha_j)$ give, for each fixed $ell$,
$EE[B_ell(zeta_r)midnu]-EE B_ell(mathcal G_calG_{q_*})to0$. To close
the brackets uniformly in $q_*$, take $Jsupsetsupp h$; since
$A_a(mathcal G_m)calG_m)lebinom{mathcal G_m(J)}calG_m(J)ale2^{-a}3^{mathcal G_m(J)}calG_m(J)$
(from $sum_a2^abinom Na=3^N$), Lemma~lem:goe-counts gives

 sup_mEE[B_{2ell}(mathcal G_m)-B_calG_m)-B_{2ell+1}(mathcal G_m)]calG_m)]le C_J2^{-(2ell+1)}.

Using (eq:dbm-bonferroni) for both models, sending $rto∞$ at
fixed $ell$ and then $ellto∞$, gives

 EEBigl[exp(-int h,ddzeta_r)midnuBigr]to
 EEexp(-int h,ddmathcal S_calS_{GOE}).

Both sides lie in $[0,1]$, so bounded convergence removes the conditioning.
For $c>0$ let $mathcal D_ccalD_czeta$ scale positions by $c$; this map is
continuous in $(c,zeta)$ for $c$ near a positive number and $zeta$
locally vaguely convergent. Since $c_rto1$ almost surely,
$zeta_r=mathcal D_calD_{c_rrho}nu^{(r)}⇒mathcal D_calD_rhonu$, and
bounded continuity of the Laplace functional of a compactly supported
continuous nonnegative test gives

 EEexp(-int h,ddmathcal D_calD_rhonu)
 =EEexp(-int h,ddmathcal S_calS_{GOE})

for $0le hin C_c^∞$, extended to $0le hin C_c$ by convolution with
a nonnegative approximate identity and bounded convergence.

Proof. [Proof of Theorem~thm:main]
Start with any sequence of admissible sizes; Proposition~prop:microscopic-limit
supplies a further subsequence on which $nu_n⇒nu$ as above. The
graph point measure is $mathcal D_calD_rhonu_n$, and by bounded continuity of
its Laplace functional and (eq:dbm-limit-identification) its Laplace
expectation along this subsequence tends to that of $mathcal S_calS_{GOE}$;
the reference GOE expectations at dimension $n$ tend to the same value along
all sizes. Every original subsequence thus has a further subsequence on
which the difference tends to zero, proving the full-sequence limit. All
excluded events had probability tending to zero in the original
simple-graph model; no conditioning on connectedness or bipartiteness enters.