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 1/2 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);=;d,√(4(d-1)-t^2)/2pi,(d^2-t^2)
   for |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/(2pi·9)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)
  (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.
 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.
 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:
 Macroscopic law. Kesten, McKay: (eq:km-density).
 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.
 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.
 Edges. Huang--McKenzie--Yau proved optimal rigidity and
GOE edge universality for every fixed dge3.
 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).
 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 ineuality becomes a
differential-entropy ineuality after smoothing by a Gaussian field;
then a strict tensor ineuality turns the small entropy deficit into
moment matching, by induction over Hermite degree.
 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.
 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.
 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), p=frac Ed,
  a=(d)/(d-1)(1-p^2),
  m_*=m_d(E+i0)=m_rmi0)=m_{r}+irmipirho .

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

  1-p^2=(d^2-E^2)/(d^2),
  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),
  nu_n=sum_{i=1}^ndelta_{x_i},
  m_n(z)&=sum_{i=1}^nfrac1{x_i-z}
        =frac1nTr R_n(z),
  R_n(z)&=(A-E-frac zn)^{-1}, 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,
  C_k=I_{2k}+pJ_k,
  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)·frac1{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 n2· ((binom d2)^2· 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)· 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.
 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)}.
 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_{HY}ge12, b_{HY}ge25a_{HY} and
0<r<c/32,r<c/32, and let r_{HY} be
rlog_{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_{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_{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+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+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_{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},
   i.e. 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+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 ,
   so  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},
  |1-q^2|=√(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, t|->log|t-E| is bounded near E because rho_d is,
and
[
  intlog|1-(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 ineuality 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 ineuality, 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))^{T}
  inR^{2k},
   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 H 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 H 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
ineuality 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 -1/2 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 ineuality 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,
  4ellgamma+2alpha<frac14,
  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 ineuality 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,

  H(F_{star})-frac d2H(F_{edge})
   ge frac1nH(Y_{global})-O(n^{-c}),

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

  frac1nH(Y_{global})
  =ell(frac12logfrac{2pi{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
L 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(L)), 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(L)
  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 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 |{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^2· Cn^{-1/2+2alpha})^{1/2}
  =n^{2ellgamma-1/4+alpha+o(1)}=O(n^{-c}),
]
because 4ellgamma+2alpha<1/4. 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_{bin}(h), where H_{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(L) ge frac1nH(Y_{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)=H(X)-(H(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_{global} is a
deterministic shift of Y_{global} we obtain
(eq:ent-star-edge).

Step 4: the global entropy. The determinant formula for a
Gaussian gives
[
  frac1nH(Y_{global})
  =(ell)/(2)logfrac{2pi{e}eta}{pirho}
   -(ell)/(n)sum_{i=1}^nlog|lambda_i-E-irmieta|,
]
because the eigenvalues of Gamma are
etabigl/(pirho,((lambda_i-E)^2+eta^2)), each with
multiplicity ell. The function t|->log|t-E-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
n·eta^{-1}· 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-irmieta| into
log|t-E|. This proves (eq:ent-global-entropy).



The entropy ineuality (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
{ 1}^{perp}, the hyperplane orthogonal to the all-ones vector, and
define the linear map to the ``clean'' coordinates, and the residual,
by
[
  B(Z,X)=(Z,(PX)/(B(Z,X)=(Z,(PX)/(sqrt a)),
  R_{res}(Z,X)=sum_{i=1}^dX_i-EZ .
]
Everything is tensored with the identity on the ell colours, so
D:=(R ⊕ { 1}^{perp}) ⊗ R^{ell} is the space in which
B takes its values. Let w denote the unit vector along the
central direction R inside R ⊕ { 1}^{perp}, and put
[
  v_i=p,w+sqrt a,Pe_i,
  U_i=spanop(w,v_i) ⊗ R^{ell}(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=B(V_{star}) where V_{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,

  H(εcalH(ε B_V+G)-frac d2H(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)|->(B(Z,X),R_calB(Z,X),R_{res}(Z,X)) in the orthonormal
coordinates of R^{d+1}=R ⊕ R{ 1} ⊕ { 1}^{perp}, one
colour at a time. On the central line it is the identity. On the
all-ones line, R_{res} reads sqrt d,t-EZ where X=t{
1}/sqrt d+..., so the diagonal entry is sqrt d. On each of the
d-1 remaining directions, B 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_{res}(V_{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_{star}=ε V_{star}+Y_{star}
at a uniform centre has, conditionally on the location, the covariance
matrix at most

  Cov(R_{res}(F_{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
[
  H(R_{res}(F_{star})mid B(F_{star}))
  le(ell)/(2)logfrac{2pi{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+irmieta,
together with the tree covariance of Lemma~lem:tree-identities,
the covariance matrix of the smoothed star noise
B(Y_{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
[
  H( H(X)+H(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)|->(B(Z,X),R_calB(Z,X),R_{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
[
  H(F_{star})
  leH( B(F_calB(F_{star}))
   +(ell)/(2)logfrac{2pi{e}eta}{pirho}
   -ell(frac12log d-(d-1)/(2)log a)+O(n^{-c}).
]
The edge observation is the image of P_iF_{star} under the
coordinate map of U_i, which has determinant
(1-p^2)^{ell/2}, so
[
  H(F_{edge})
  =H(P_iF_{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{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
[
  H( B(F_calB(F_{star}))
  -frac d2,H(P_iF_{star}) ge -O(n^{-c}).
]
Replacing the smoothed eigenvector star F_{star} by
ε B_V+G costs polynomial entropy in exactly the same way as in
the noise-replacement step above: the covariance of
B(F_{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^{T} X_iX_i^{T}{ell},
  EE ZX_i^{T}Λ)/(d),
  EE X_iX_j^{T}{Λ^2-d,Id_{ell}}{d(d-1)}
  (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^{T}+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 ineuality (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 ineuality for
symmetric tensors that forbids a nonzero first Hermite coefficient.

Write phi_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)midphi_L)=int_L flog(f)/(phi_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

  H(Z)-H(Z)-frac12sum_{i=1}^dH(P_iZ)
  =-Ent(Law(Z)midphi_D)
   +frac12sum_{i=1}^dEnt(Law(P_iZ)midphi_{U_i}).



Proof. 
For any X in a Euclidean space L with a density,
H(X)= H(H(phi_L)-Ent(Law(X)midphi_L); this is
the definition of relative entropy together with intphilogphi
=-H( Also H( L)/(2)log(2pi{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{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
   for 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^+=w+v_i/√(2(1+p)),
  e_i^-=w-v_i/√(2(1-p)),
  s_i=(e_i^++e_i^-)/(sqrt2),
  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),
  ip{w}{t_i}^2=(1-q_0)/(2),
  ip{Pe_i,s_i}^2=(1-q_0)/(2),
  ip{Pe_i,t_i}^2=(1+q_0)/(2).
]
The strict bulk condition is
[
  p^2<(4(d-1))/(d^2)
  ⇔
  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),
  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^{T}+b_t,t_it_i^{T}
  =Id_{R ⊕ { 1}^{perp}}

holds in R ⊕ { 1}^{perp}, tensored with the identity on the
colours. To verify it, expand both sides in the orthogonal decomposition
R ⊕ { 1} ⊕ { 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--{ 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 { 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)^{T}

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
ineuality 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}H for 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{ 1}^{perp} be the unit vector along Pe_i. On
R ⊕ { 1}^{perp} the spatial part of this commutator is
[
  ip{a_i'}{a_k'}(a_i'(a_k')^{T}-a_k'(a_i')^{T},
   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 { 1}^{perp}, so
their common kernel in R ⊕ { 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^{T}-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,
[
  H(εcalH(ε B_n+G)-frac12sum_{i=1}^d
    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!},
   tin D .
]
Set H_j=EEHe_j(B_n). We show by induction on j that

  |H_j|=O(n^{-c_j})(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},
  2jtau<frac{c_*}2, (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
ineuality 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)midphi_D)
  -frac12sum_{i=1}^dEnt(Law(P_iZ_u)midphi_{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 phi_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(phi_D)}^2
  =(|H_k|^2)/(k!) .
]
The Taylor remainder through degree j is
O(u^{j+1}n^{o(1)}) in L^3(phi_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(phi_D) norms are n^{o(1)}, uniformly over the source
values. Consequently
[
  int_D|f_u-1|^3,ddphi_Dle Cu^3n^{o(1)} ,
]
and, since |(1+h)log(1+h)-h-h^2/2|le C|h|^3 for |h|le1/2,
while int(f_u-1)ddphi_D=0, we get

  Ent(Law(Z_u)midphi_D)
  =frac12int_D(f_u-1)^2ddphi_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)^2ddphi_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)midphi_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 ineuality. 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^{T}+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_{res}+EZ by definition of the residual, so
[
  ip{v_i}{B_V}
  =X_i+(p-frac Ed)Z-frac1dR_{res} .
]
Here p=E/d exactly, and |R_{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, v=e_w-e_x,
  Z=uv^{T}+vu^{T}
]
Expanding,
[
  Z=e_ve_w^{T}+e_we_v^{T}+e_xe_y^{T}
    +e_ye_x^{T}
    -e_ve_x^{T}-e_xe_v^{T}
    -e_we_y^{T}-e_ye_w^{T},
]
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^{T}, A_{uu}=u^{T},
  A_{vv}=v^{T},
  D=(1+A_{uv})^2-A_{uu}A_{vv}.
]
Then Dne0 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^{T}+vu^{T}'=I we get
[
  G'=G-Guv^{T}'-Gvu^{T}' .
]
Set X=v^{T}' and Y=u^{T}', row vectors. Multiplying
the displayed identity on the left by v^{T} and by
u^{T} gives the linear system
[
  1+A_{uv}&A_{vv} A_{uu}&1+A_{uv}
  X Y
  =v^{T} u^{T},
]
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^{T}-A_{vv}u^{T}],
  Y=D^{-1}[(1+A_{uv})u^{T}-A_{uu}v^{T}].
]
Now take the (x,y) entry of G'=G-Guv^{T}'-Gvu^{T}':
[
  G'_{xy}=G_{xy}-(Gu)_x(e_y^{T}{T}')
                  -(Gv)_x(e_y^{T}{T}'),
]
and e_y^{T}{T}'=(G'v)_y=(X)_y,
e_y^{T}{T}'=(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}^{T}-E-frac zn)^{-1}W_{n,k}
 =sum_ifrac{b_ib_i^{T}}{x_i-z},
 b_i=W_{n,k}^{T}, 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+iy)-m_*|rmiy)-m_*|
 +|M_{n,k}(s+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+iyrmiy
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+iS,rmiS,
Proposition~prop:spectral-input and Lemma~lem:tree-identities
give, for some c>0,

 m_n(z_0)=m_*+O(S^{-c}), M_{n,k}(z_0)=K_k(m_*)+O(S^{-c}),
 Im m_n(iS^2)+rmiS^2)+|Im M_{n,k}(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 Z_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
zinZ_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-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^{T}-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

 |mhat_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^{T}, and the resolvent identity reads
(I+UB)Uhat=UU=U with U=F^{T}-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}),
 Uhat=m&-qU=m&-q-q&m+O((Im z)^{-c_2}),
  q=(1+Em)/(d), m=m_n(z),

the sign change in UUhat coming from expressing the updated pairs
in the old columns. The off-diagonal of (I+UB)Uhat=UU=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-irmieta'|ge(eta'/eta)|x-s-irmieta|, so

 |m_n(s+irmieta')-m_n(s+irmieta)|
 le(eta-eta')/(eta')Im m_n(s+irmieta)leIm m_n(s+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+iu)=m_*+O(u^rmiu)=m_*+O(u^{-xi}), |e|le R, R^deltale ule R^2,

and

 nu([v,w])=rho(w-v)+O(R^{1-chi}), -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_{r}+PVsum_jfrac1{x_j-z},
 M_k(z)=K_k(m_{r})+PVsum_jfrac{C_k^{1/2}g_jg_j^{T}{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(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 K+ and y>0. For a resolvent
sum f with nonnegative residues, comparing kernels on the real line
gives Im f(z)le C_{K,y}K,yIm f(iy)rmiy) and
|f(z)-f(iy)|rmiy)|le C_{K,y}K,yIm f(iy)rmiy) for
zinK; 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-i)^rmi)^{-1}|le C_zL^{-1}Im(x-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 iLrmiL 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(i)&=rmi)&=sum_j(frac1{x_j-z}-frac1{x_j-i}rmi),
 Im M_k(z)&=sum_jb_jb_j^{T}{x_j-z},&
 M_k(z)-M_k(i)&=rmi)&=sum_jb_jb_j^{T}{x_j-z}-frac1{x_j-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
ineuality 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, 1/2R^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+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)-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(iR)=rmiR)=irmipirho+o(1). By
(eq:lim-difference-identities), m-S_nu is constant with zero
imaginary part; the two displays identify it as m_{r}.

Marked principal values. Fix k, condition on nu, and set
Z_j=b_jb_j^{T}-C_k, independent centered with bounded second
moments. Ordered by increasing |x_j|, sum_jZ_j/(x_j-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 irmi converges absolutely
and locally uniformly; with the scalar PV this gives local uniform
convergence of S_{k,nu}(z)=PVsum_jb_jb_j^{T}-z)^{-1}.
The linear count gives sum_j|x_j-iR|^rmiR|^{-2}=O(R^{-1}), so each
entry of the centered marked series at iRrmiR has conditional
second moment O(R^{-1}), whence
S_{k,nu}(iR)rmiR)toirmipirho,C_k in probability. Since
M_k-S_{k,nu} is almost surely constant in z and
M_k(iR)rmiR)to K_k(m_*) in probability, that constant equals
K_k(m_*)-irmipirho,C_k=K_k(m_{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 S_r(z) is approximated,
in law and up to o(1), by m_{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^{T}-J_k, which is symmetric. For an analytic
symmetric matrix function M with positive definite imaginary part, put

 T_O(m,M)=(m-∂_zcalT_O(m,M)=(m-∂_zlogdet(I_{2k}+B_OM),,
 D_O^{T}{-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 T_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 A|-> A+WB_OW^{T},
W|-> 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
mhat_n=m_n-∂_zm_n=m_n-∂_zlogdet(I_{2k}+B_OM_{n,k}) and
Mhat_M_{n,k}=D_O^{T}{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_{r})
commute with Q_U, the Gaussian representation
(eq:marked-representation) shows that
R_U(m,M)=(m,Q_U^R_U(m,M)=(m,Q_U^{T} preserves the joint law. Since
Q_UD_OQ_U^{T}{UOU^{T}} and
Q_UB_OQ_U^{T}{UOU^{T}}, we get
R_U^{-1}T_O R_U=R_U=T_{UOU^{T}};
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^{T} and
D_{OO'}=D_OD_{O'}, so T_O T_calT_{O'}=T_{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), k=2r, K=r^{1/100}, T=r^{tau},
 theta=√(T/r), L=r^{100/chi}, 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_{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, 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 W 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),W .


[Gaussian insertion]
The coefficients satisfy

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

with deterministic positive constants. As rto∞,

 m_{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
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_{r}+Tr(H_r-z)^{-1} is replaced by S_r(z).


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

 nu([-R,R])=2rho R+O(R^{1-chi}),
 |PVsum_{|x_j|>R}frac1{x_j}|=O(R^{-chi}),
 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}), |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}, 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
S_r(z)=m(z)-S_r(z)=m(z)-Tr(B^{-1}+M(z))^{-1}M'(z) costs o_{PP}(1): the
marked error is O_{PP}(r· r^5· r^{-98})+O_{PP}(r· r^{10}
· r^{-98}· 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 G of
independent standard Gaussian entries with 4r columns, and the
finite-rank inverse identity (Woodbury) yields

 S_r(z)=m_{r}+Tr(D+A'-z)^{-1}+o_{PP}(1),
 A'=G G^G^top, 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), |U_I^2|le h^{-1},
 |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 A{I,F} put a_I=h, a_F=K. With probability tending
to one,

 bigl|U_{A}'_{A A}A-s_1Id)U_{A}|
 lepolylog(r)(sqrt{T/a_{A}}+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_{A} 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_{A} The bounded matrix Bernstein
ineuality gives a two-sided tail 2|A|A|exp(-ct^2/
(sigma_{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),
 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=G_I G_I^G_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_jmid·]=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-1/{200} and
tfrac{13}2tau-1/2 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-ih)U_Irmih)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(ih)U_I^rmih)U_I^{-1}|le C, and
G_I(ih)^*G_I(rmih)^*G_I(ih)=h^rmih)=h^{-1}Im G_I(ih)rmih) gives
|G_I(ih)U_I^rmih)U_I^{-1}|le Ch^{-1/2}. The twice-iterated resolvent
identity G_I(z)=G_I(ih)+(z-rmih)+(z-ih)G_I(rmih)G_I(ih)^2+rmih)^2+
(z-ih)^2G_I(rmih)^2G_I(ih)G_I(z)G_I(rmih)G_I(z)G_I(ih),rmih), conjugated by
U_I^{-1}, bounds the three terms by C, C|z-ih|/hrmih|/h and
C|z-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),W.


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}W_ with independent raw GOE
W_ (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}W_ has the law
√(s_2),W.


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 S_r. Bounded convergence averages
the conditional error, and by Proposition~prop:continuous-switching
the unconditional law of S_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-i)^rmi)^{-1} is tight by (eq:sw-trace-limit), so
sigma_r is tight on the one-point compactification overlineR. For
z=e+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+irmieta)le Cquad (|e|le G, gleetale10),
  |V|le m^{C_V},

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

 m^{delta-1}le gle m^{-delta}, 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_{{fc},t}(zeta)=m_V(zeta+t,m_{{fc},t}(zeta)),
 rho_{{fc},t}(e)=pi^{-1}lim_{eta↓0}Im m_{{fc},t}(e+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_{{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),W_r, |s_1|le CT, cTle s_2le CT,

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

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

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

 (alpha_*)/(q_*)H_rstackrel{law}{=}V_*+√(t_*),W_{q_*},
  V_*=(alpha_*)/(q_*)(D_I+s_1Id), 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_{{fc},t_*}(0)to1/pi; with c_r=pirho_{{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), 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), F_r(w)=alpha_*^{-1}sum_{jin I}(x_j-(alpha_*s_2w-s_1))^{-1}.

On |w-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}=irmipirho+o(1), and
PVsum_{|x_j|>K}(x_j-z)^{-1}=O(K^{-chi}+T/K). Hence F_r(w)=i+o(1)rmi+o(1)
uniformly; Rouch'e for w-F_r(w) against w-irmi gives a unique
simple zero w_rtoirmi 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_{{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_{{fc},t_*}(zeta/c_r) and density
c_r^{-1}rho_{{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)·(c_ralpha_*)/(q_*)lambda_j(H_r)=c_rrho,lambda_j(H_r),

with c_r=1 at exceptional initial values.


Let G_m={j=1}^mdelta_{(m/pi)lambda_j(W_m)}. These
converge vaguely in law along all dimensions to a common locally finite
process S_{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
G_m.


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


Proof. 
Bound the GOE count by twice a Gaussian-unitary count plus one via the
Gaussian decimation identity of Forrester--Rains. Let O_mrmO_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 G_m(J) is the O_mrmO_m count in
J_m=pi J/√(2m), so |J_m|=O(m^{-1/2}). Superposing independent
O_mrmO_m and 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 O_mrmO_m count,
EE e^{bG_m(J)} 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 ineualities 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(G_{q_*})to0. To close
the brackets uniformly in q_*, take Jsupsetsupp h; since
A_a(G_m){G_m(J)}{-a}3^{G_m(J)}
(from sum_a2^abinom Na=3^N), Lemma~lem:goe-counts gives

 sup_mEE[B_{2ell}(G_m)-B_calG_m)-B_{2ell+1}(G_m)]G_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,ddS_{GOE}).

Both sides lie in [0,1], so bounded convergence removes the conditioning.
For c>0 let D_c 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=D_{c_rrho}nu^{(r)}⇒D_, and
bounded continuity of the Laplace functional of a compactly supported
continuous nonnegative test gives

 EEexp(-int h,ddD_
 =EEexp(-int h,ddS_{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 D_, and by bounded continuity of
its Laplace functional and (eq:dbm-limit-identification) its Laplace
expectation along this subsequence tends to that of S_{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.