GOE bulk universality for fixed-degree random regular graphs

Proof. Every simple graph has the same number (d!)^n of stub realizations, which proves the conditional uniformity. A prescribed loop is one pair of stubs at a vertex. A prescribed double edge is an unordered pair of pairs between two distinct vertices, using distinct stubs. Let Z count all occurring prescribed defects of these two kinds. The expected counts tend to

respectively. Indeed, a prescribed family of e disjoint pairs occurs with probability

for fixed e.

For each fixed j, ordered j-tuples of defects on disjoint vertex sets give the leading contribution to \EE(Z)_j. Compatible overlapping tuples contribute o(1). To see this, regard their prescribed pairs as a multigraph with v vertices and e edges. Every vertex has degree at least two. A connected component all of whose degrees equal two can only be a single prescribed loop or a single prescribed double edge. Thus a genuinely overlapping component has v<e. There are O(n^v) choices and each has probability O(n^{-e}). There are only finitely many overlap patterns for fixed d,j. Consequently

The finite inclusion-exclusion bounds for \PP(Z=0), followed by increasing their order, give the asserted limit. Absence of these defects is exactly simplicity.

For a triangle through a fixed vertex, summing over the other two vertices and their stub choices gives O(n^2)\,O(n^{-3})=O(n^{-1}) in the matching model. Conditioning on simplicity changes this bound by a bounded factor. Finally, in any simple d-regular graph a uniform endpoint is a uniform vertex. Each previously chosen endpoint has only O_d(1) vertices equal or adjacent to it. A union bound over the fixed number of sampled endpoints proves the last assertion.□

Proof. We use the fixed-degree results of Huang and Yau [18, Theorems 1.2, 1.4, and 4.2], in the precise version arXiv:2102.00963v3. Their matrix is A/\sqrt{d-1}; hence their bulk coordinate is E/\sqrt{d-1}\in(-2,2). Rescaling the spectral parameter and the resolvent by this fixed factor gives our normalization.

For clarity, the relevant error and geometry can be checked directly in their Definition 4.1 and Equations (4.6)–(4.7). Choose fixed parameters

and neighborhood radius r_{\mathrm{HY}}=\mathfrak r\log_{d-1}n, rounded to an integer. Their Definition 1.1 and Proposition 2.1 give a tree-like event with probability tending to one, on which all but a polynomially small fraction of vertices have a tree neighborhood of radius (c/4)\log_{d-1}n. This is more than enough to contain every neighborhood needed for the resolvent approximation between two vertices of the same star. The tree extension with the semicircle weight in their Equation (2.12) is then the ordinary regular-tree Green function.

In a fixed bulk interval, their error is bounded by a logarithmic factor times

The factor S dominates every fixed power of \log n. After decreasing c_1 and increasing B, this gives the stated bounds and exceptional-set size. The same source provides polynomial rigidity and polylogarithmic complete delocalization. Its model and the quoted statements have no even-size restriction.

For the final assertion, apply the trace estimate at \zeta=E+\mathrm iS^2/n. Each eigenvalue with |\lambda_i-E|\le S^2/n contributes at least n/(2S^2) to \Im(\lambda_i-\zeta)^{-1}. The normalized imaginary trace is bounded, so the count is O(S^2).□

Proof. A Schur complement at the root of a forward tree gives q=(\zeta-(d-1)q)^{-1} in the convention (\zeta-A)^{-1}. Attaching d forward branches yields the displayed formula for -m_d. The boundary value \pi^{-1}\Im m_d(t+\mathrm i0) is the Kesten–McKay density above; the branch has the required mass at infinity and no additional pole.

The resolvent equation at the center and symmetry among its neighbors give the adjacent entry. Applying the same equation to a neighbor column gives, for the common distinct-neighbor entry g_2,

Taking boundary imaginary parts proves the covariance formulas. The tree formula is analytic across each side of a sufficiently small bulk neighborhood, so these limits have a polynomial rate. Substituting its boundary value into k_0=(1+Em_*)/d gives the quadratic identity.

To compute the potential, integrate its analytic derivative -m_d(\zeta)=q/(1-q^2) after changing from \zeta to q. Matching at infinity gives

At the bulk boundary,

These identities give (2.3) after substituting a=d(1-p^2)/(d-1). The boundary logarithmic integral exists because \rho_d is bounded near E. Finally,

by boundedness of the density near E, the substitution t-E=\eta u, and the integrable tail of \log(1+u^{-2}).□

Proof. Work first in the configuration model, with an ordering of the d half-edges at each vertex. Prescribe the vertex labels and the ordered list of desired neighbor labels at every vertex. For a fixed ordered tuple type T_0, there are at most \exp(nH(T_0)) such arrays, by the multinomial bound.

Let m_{jk}=ndQ_{jk} be the number of half-edges whose vertex has color j and whose prescribed neighbor has color k. Compatibility requires m_{jk}=m_{kj}, and m_{jj} must be even. The probability that a uniform matching respects all these prescriptions is

with zero-size factors omitted. Stirling’s bounds give its logarithm as

Thus the expected number of labelings of ordered type T_0 has the claimed exponential bound with H(T_0) in place of H(T).

Symmetrization increases entropy, so H(T_0)\le H(T); it leaves Q unchanged. There are at most \exp(O(M^{d+1}\log n)) possible ordered tuple types. Summing over those with a given symmetrization therefore changes only the sublinear error. Markov’s inequality, followed by a union bound over all types, gives (3.1) with a slightly larger sublinear allowance. The simplicity probability is bounded away from zero by Lemma 2.1, so the statement transfers to uniform simple graphs. The symmetrized histograms no longer depend on the half-edge ordering.□

Proof. Quantize each coordinate into intervals of width b=n^{-\gamma}, clipping all values outside a range of order [-n^\gamma,n^\gamma] into its two end bins. The resulting vertex alphabet has M=n^{2\ell\gamma+o(1)} elements, so Lemma 3.2 applies. Let \mathcal L be the random quantized labeling of F, and let T,Q be its random histograms. Entropy of the type index, followed by entropy bounded by logarithmic type size, gives

Here concavity was used only for the positive star term.

We next replace \EE H(Q) by H(\EE Q). Two independent edge locations have Gaussian observations whose off-diagonal covariance block has average Frobenius norm at most

Indeed endpoints from different sampled edges are independent uniform vertices, and \Tr\Gamma^2\le Cn/\eta^2. Each individual edge covariance has smallest eigenvalue at least c_1\eta, since \|A\|\le d. After whitening the two marginal blocks, Gaussian total variation from the product law is bounded by the off-block norm when that norm is small, and by one otherwise. The small-norm bound follows, for example, by differentiating the Gaussian density along the straight line joining the two covariance matrices. Averaging yields

where the marginal laws on the right are at the fixed sampled locations. Deterministic shifts and subsequent quantization do not increase this distance.

Consequently each cell frequency of Q has variance O(n^{-1/2}\eta^{-2}). There are M^2 edge cells, so

The finite-alphabet entropy continuity bound H_{\rm bin}(h)+h\log(M^2) at total variation distance h now proves

We record the differential-entropy passage carefully. Every fixed location covariance, including the complete array covariance, has spectrum in [c_1\eta,C/\eta]; every coordinate of the deterministic shift is n^{o(1)}. Unbounded binning and the entropy bound on a cube give

The displayed error accounts for clipping. It is in fact superpolynomially small: one-coordinate Gaussian tails beyond n^\gamma, together with their quantized entropy contributions, are superpolynomially small uniformly in the shifts, because \alpha<\gamma. The loss for the full array is at most the sum of the coordinate losses. Translation by the deterministic array \epsilon V does not change its differential entropy.

For either fixed-dimensional local observation X=F_{\rm star} or F_{\rm edge}, we also have

To justify the two-sided estimate, fill each bin uniformly. Bin averaging a Gaussian density changes it in total variation by at most Cb\eta^{-1/2}, by integrating its absolute first derivatives; the same bound holds for a mixture. Clipping adds only a superpolynomial tail error. Both resulting densities have polynomial supremum and second-moment bounds. Such bounds convert polynomial total variation into polynomial entropy error: subtract the common minimum of the two densities. The remaining densities have equal mass h; after division by h, their supremum and second moment are bounded by a polynomial divided by h, so their entropies have absolute value O(1+\log(n/h)). The entropy errors in mixing them with the common part are at most the binary entropy. This proves (3.8), and will also apply to fixed-dimensional projections below.

The laws \EE T,\EE Q are exactly the binned local laws. Substitute (3.7) and (3.8) into (3.6). The mesh terms cancel, since (d+1)\ell-(d/2)(2\ell)=\ell. This proves (3.4).

Finally the Gaussian determinant formula gives

The test function in this sum has Lipschitz constant O(\eta^{-1}). Rigidity therefore gives polynomial error in replacing the eigenvalues by Kesten–McKay quantiles. The quantile sum differs from its integral by O((n\eta)^{-1}), by summing the lengths of the quantile intervals. An omitted or extra bounded eigenvalue costs O(\log n/n). Lemma 2.4 gives the remaining regularization error O(\eta). This proves (3.5).□

Proof. The absolute Jacobian determinant of (Z,X)\mapsto(\mathcal B(Z,X),R_{\rm res}(Z,X)) is (\sqrt d\,a^{-(d-1)/2})^\ell. For each color, the residual variance of Y, averaged over centers, is at most \eta/(\pi\rho), since

The eigenvector equations give a residual O(n^{-1+o(1)}) for V, uniformly over locations and choices: |\lambda_{i(j)}-E|\le S^2/n. More precisely, put \Delta=\diag(\lambda_{i(j)}-E:1\le j\le\ell). At a uniform center v, the deterministic residual is \epsilon\Delta V(v), whose second-moment matrix is \epsilon^2\Delta^2 by eigenvector orthogonality. The Gaussian colors are independent and centered conditionally on the location, so cross terms vanish and

The relative correction inside the logarithmic entropy bound satisfies

The maximum-entropy bound by covariance and the fact that conditioning lowers entropy thus give

By the typical-star estimate at E+i\eta and the tree covariance in Lemma 2.4, the covariance of the noise \mathcal B(Y_{\rm star}) differs from \Id_D by a polynomially small error outside a polynomially small fraction of centers. Gaussian total variation after translation, followed by the entropy-continuity argument in Lemma 3.3, allows this noise to be replaced by G at polynomial entropy cost.

Similarly, the edge coordinates have the same entropy, up to polynomial error, as the observations of \epsilon B_V+G in directions w,v_i. The deterministic discrepancy in the neighbor coordinate is O(n^{-1+o(1)}), because the neighbor average is p times the center up to that error. The map from orthonormal coordinates on U_i to these two edge coordinates has absolute determinant (1-p^2)^{\ell/2}. It follows that the left side of (3.4) is at most the left side of (3.9) plus

Combine (3.4), (3.5), and the potential identity (2.3). Their constants agree exactly, proving (3.9).

For the remaining assertions, put \Lambda=\diag(\lambda_{i(1)},\ldots,\lambda_{i(\ell)}). Orthogonality of the normalized eigenvectors, the eigenvector equations, and uniform permutation of neighbors give the exact second-moment identities

The mean is zero because these eigenvectors are orthogonal to constants. These identities imply \EE B_VB_V^\top=\Id_D+O(n^{-1+o(1)}), using \Lambda=E\Id_\ell+O(n^{-1+o(1)}). Finally an actual uniform oriented edge is invariant under reversal. Its two coordinates differ from the w,v_i observations of B_V by O(n^{-1+o(1)}). This error remains of that order in every fixed moment because all coordinates are n^{o(1)}. The coordinate map is invertible since |p|<1, so reversal is precisely the asserted orthogonal swap on U_i.□

Proof. We use the weighted-frame geometry of Backhausz and Szegedy [3, Lemma 11.3], with the numbering of the published version, and then prove the tensor strictness needed here. In the spatial plane spanned by w,v_i, define the orthonormal bisecting directions

The swap fixes e_i^+, negates e_i^-, and hence interchanges the orthonormal directions s_i,t_i. Their squared central components are

Writing q_0=\sqrt{1-p^2}, the strict bulk condition is

For d\ge3, the two squared central components therefore strictly bracket 1/d. The weights

are positive, sum to one, and satisfy

Indeed its central component is one by the choice of weights. Its mixed block vanishes because \sum_iPe_i=0. The remaining block is a scalar identity on {\bf1}^{\perp}; its trace determines that scalar to be one.

Tensor (3.13) with the color identity, and use it in the first slot of \|H\|^2. In each summand project the other j-1 slots onto U_i. This can only decrease the value. After projection, swap invariance makes the squared norms associated with s_i and t_i equal. Since their weights sum to one, the resulting value is \frac12\sum_i\|P_i^{\otimes j}H\|^2. This proves nonnegativity of the left side of (3.12).

Suppose equality holds. Positivity of both weights implies that every discarded squared norm vanishes. If R_{i,s} and R_{i,t} are the projections onto s_i\otimes\R^\ell and t_i\otimes\R^\ell, then

and likewise for t; superscripts specify tensor slots. Adding the two identities and using symmetry gives

In particular each projection may be moved from one slot to another when acting on H. For i\ne k, the availability of three slots now gives

Thus [P_i,P_k] annihilates the first slot of H.

Let a_i' be the unit vector proportional to Pe_i. The spatial part of this commutator is

For d\ge3, it is invertible on \spanop(Pe_i,Pe_k) and vanishes on its orthogonal complement. The simultaneous kernel of all these commutators is therefore the central direction, tensored with the color space. Symmetry leaves H supported in (w\otimes\R^\ell)^{\otimes j}. But its projected swap invariance would also support it in (v_i\otimes\R^\ell)^{\otimes j}; the two spatial lines are distinct, so H=0. The continuous quadratic form is consequently strictly positive on the unit sphere of the finite-dimensional constraint space. Compactness proves (3.12).□

Proof. For y\in D, define the Hermite tensors \operatorname{He}_j(y) by the generating identity

Set H_j=\EE\operatorname{He}_j(B_n). We prove by induction that

This suffices because Hermite polynomials form a basis of the ordinary polynomials. The first- and second-moment hypotheses give the cases j=1,2.

Fix j\ge3, assuming the preceding orders have been proved. Choose \epsilon=n^{-\tau}, where \tau>0 satisfies

This choice may depend on j. Put u=\epsilon/\sqrt{1+\epsilon^2}. Dividing every entropy argument in the hypothesis by \sqrt{1+\epsilon^2} leaves that inequality unchanged, because the dimensions cancel. The normalized random vector is

Its second-moment matrix is \Id_D+u^2O(n^{-c_*}). Consequently (3.11), with this second-moment error included, turns the hypothesis into

We claim that

and that the analogous formula for P_iZ_u has P_i^{\otimes j}H_j in place of H_j. Let f_u denote the density ratio of Z_u to \varphi_D. The coefficient of u^k in its Taylor expansion is

This follows by integrating against the generating exponential: the resulting transform is \EE\exp(u\ip{B_n}{t}-u^2\|t\|^2/2). Pairing two generating exponentials under Gaussian measure gives Hermite orthogonality and the identity

The Taylor remainder through degree j is O(u^{j+1}n^{o(1)}) in L^3(\varphi_D). Indeed, at a fixed source value b_0, the ratio at parameter s is

For 0\le s\le u its required derivatives are bounded by polynomials in \|y\|,\|b_0\| times C\exp(Cs\|b_0\|\|y\|). Since \|b_0\|\le n^{o(1)} and u=n^{-\tau+o(1)}, their L^3(\varphi_D) norms are n^{o(1)}, uniformly over the source values. Taylor’s theorem and averaging prove the bound.

For a density ratio 1+h, the pointwise remainder after subtracting h+h^2/2 from (1+h)\log(1+h) is bounded by C|h|^3 for the entire range h\ge-1. As \int h\,\dd\varphi_D=0,

All lower H_k, 1\le k<j, are polynomially small by induction, whereas every fixed H_k is at most n^{o(1)}. The inequalities defining \tau make the lower-degree part, together with the Taylor remainder, of size u^j O(n^{-\delta}) in L^3 for some \delta>0. In particular \|f_u-1\|_{L^3}=u^j n^{o(1)}, so the cubic entropy error divided by u^{2j} is O(u^j n^{o(1)}), again polynomially small. Orthogonality now proves (3.17). The same reasoning applies to each projection.

Substitute these expansions into (3.16) and divide by u^{2j}/(2j!). The choice of \tau gives

Hermite tensors commute with orthogonal projections and are equivariant under orthogonal maps, so

in the Gaussian space U_i. Since \operatorname{He}_j is a fixed polynomial, the assumed moment symmetries give the swap constraints of Lemma 3.5 with polynomially small errors. In fixed finite dimensions the distance to the common kernel of a collection of linear maps is bounded by a constant times the sizes of their images. Project H_j onto that kernel. This changes it by a polynomially small amount, and changes the quadratic expression in (3.18) by a polynomially small amount, since \|H_j\|=n^{o(1)}. Lemma 3.5 therefore proves (3.15) and closes the induction. Each step uses only the common constants in the hypotheses, proving uniformity.□

Proof of Proposition 3.1. Condition on a graph in the common events of Lemmas 3.2 and 3.4 and Proposition 2.2. Apply Lemma 3.6, with q replaced by \max(q,2), to B_n=B_V. Complete delocalization gives its deterministic subpolynomial bound. Lemma 3.4 supplies its first two moments and all fixed-order projected moment symmetries. Averaging (3.9) over i gives the criterion’s entropy hypothesis. All these inputs are uniform over the selected eigenindices. Hence every moment of B_V of degree at most q differs polynomially little from the corresponding standard Gaussian moment in D.

The observations of that Gaussian in directions w,v_i have covariance C_\ell. The uniformly O(n^{-1+o(1)}) eigenvector residual transfers the moment estimate to the actual edge values (V(v),V(v')). For k conditionally independent edge locations, a coordinate monomial factors across those locations. The one-edge estimates give the asserted joint moment estimate, with independent edge samples and, equivalently, independent N(0,C_k) vectors across the \ell eigenindices.

The common graph events do not depend on the selected eigenindices: Lemma 3.2 applies simultaneously to every quantized labeling, and the spectral and eigenvector residual bounds are uniform throughout the band. The Gaussian criterion preserves this uniformity at every fixed moment order. This proves the proposition.□

Proof. It suffices to prove the statement for k\ge2, since the smaller matrices are compressions. Throughout this proof k,s,y are fixed. The exponents in polynomial errors may depend on fixed moment orders.

Initial estimates. Put z_0=s+\mathrm i S. Proposition 2.2 and Lemma 2.4 give, for some c>0,

with probability tending to one. To justify the off-block entries in the second assertion, write R=(A-E-z_0/n)^{-1}. For independent uniform vertices v,w, the resolvent identity gives

Endpoints from distinct edge samples have this independent uniform distribution. Markov’s inequality makes their resolvent entries O(S^{-c}), after decreasing c. Within a single sampled edge the typical-star estimate applies, with adjacent entry (1+E m_*)/d. The same reasoning gives the last assertion of (4.3). These estimates also imply that the number of indices in |x_i|\le S^2 is O(S^2)=n^{o(1)}.

Concentration under a trace bound. Descend from S to y by halving the height, shortening the last step if necessary, and keep real part s. Let \mathcal Z_n denote these O(\log S) arguments. For every fixed C<\infty, we claim that, outside an event with limiting upper probability O_{k,D,C}(y^{-D}), simultaneously for z\in\mathcal Z_n,

Here c_2>0 can be chosen independently of D,k,s,C, with constants in the error allowed to depend on the fixed parameters.

We compare z with z_0. Once S exceeds a fixed multiple of |s|+y+1, for |x|>S^2 we have

Thus the scalar and matrix contributions outside the band are O(S^{-1}) on the event in (4.3). For the matrix bound, use that each residue b_i b_i^{\mathsf T} is positive semidefinite.

Inside the band, let

If the scalar trace at z is bounded, then

Replace the marks for this calculation by independent Gaussians of covariance C_k. For every fixed integer j, each entry of

then has 2j-th absolute moment at most C_{j,k,C}(\Im z)^{-j}. In the moment expansion a spectral index occurring only once has zero expectation. For every remaining partition of the indices, bounded Gaussian moments and (4.6) give the stated bound.

The same estimate, with an additional polynomially small error, holds conditionally on each graph satisfying the required fixed-order moment estimates. There are only n^{o(1)} terms in the expansion, and the weights are bounded because \Im z\ge y\ge1. The event |m_n(z)|\le C depends only on the graph, so it causes no difficulty in this conditional calculation. Choose

Markov’s inequality at threshold (\Im z)^{-c_2} gives a failure probability bounded by a constant times (\Im z)^{-j(1-2c_2)}, plus a polynomially small error in n. Summing over the geometric sequence of heights is bounded by C_{j,k,C}y^{-j/2}+o(1). Taking j sufficiently large proves (4.4). The graph exceptional events and the failure of (4.3) have probability tending to zero. Factors from S and the number of steps are n^{o(1)}, so the polynomial moment errors remain negligible.

A switch gives the scalar equation. The graph/list bijection below uses the measure-preserving switching mechanism of Bauerschmidt et al. [5, Section 7.2, Propositions 7.3–7.4]. Figure 1 illustrates the two-edge operation.

With probability 1-O_k(n^{-1}), the 2k sampled endpoints are distinct and the only graph edges among them are the sampled pairs. This follows directly from fixed degree and conditional independence of the edge samples. Restrict temporarily to this event. Exchange the ending vertices of the first two pairs and update their entries in the edge list. The operation is a bijection of the restricted graph/list space, and all its elements have the same probability. In particular the switched data, denoted by hats, satisfy (4.4) with the same failure bound. The restriction itself is preserved by the switch.

The adjacency update has bounded rank. The eigenvalue counting functions before and after the update consequently differ by a bounded amount. Integrating their difference against the derivative of (t-z)^{-1} gives

This calculation is in the microscopic spectral variable; equivalently, m_n(z)=\Tr(n(A-E)-z)^{-1}.

Write the first two oriented pairs as (v_1,w_1),(v_2,w_2), and put

The switch adds F B F^{\mathsf T}, where

Let U=F^{\mathsf T}(A-E-z/n)^{-1}F, and define \widehat U using the switched graph but the same old columns F. The resolvent identity yields

Suppose now that m=m_n(z) lies in a fixed bounded range. Equations (4.4) and (4.7), applied to both data sets with a sufficiently large trace threshold, give

The sign in the second matrix comes from expressing the updated edge list in the old columns: the right difference column is negated. The off-diagonal entries of (4.8) now imply

As a quadratic in m, its exact left side has simple conjugate roots in the strict bulk. By Lemma 2.4, the upper root is m_*. The lower root stays a fixed distance from every m with \Im m\ge0. Hence, throughout the bounded range,

Closing the descent. There is no a priori trace assumption at a lower step. Instead, if \eta/2\le\eta'\le\eta, the resolvent sum gives

Indeed |x-s-\mathrm i\eta'|\ge(\eta'/\eta)|x-s-\mathrm i\eta| term by term. Choose a fixed trace range containing the initial value and every possible next-step value in (4.12) when the preceding value is within 1 of m_*. Enlarge it once more for the switched trace using (4.7). Finally take the terminal height y sufficiently large that (4.11) returns a value within 1 of m_*. Induction down the heights proves the required trace bounds and the endpoint estimates at every step.

Combining the two concentration failure bounds, and decreasing the final exponent from c_2 to a smaller \xi>0, absorbs the fixed constants into the threshold y^{-\xi}. This proves (4.2). All constants are independent of s: the only additional restriction was that S eventually exceed a fixed multiple of |s|+y+1, which holds before the limit for every fixed s,y.□

Proof.

Tightness. For fixed L,y>0, positivity gives

Proposition 4.1, with y a sufficiently large fixed number, therefore proves tightness of the local counts. It also proves compact-open tightness of the analytic functions. The choice of y is made separately for each probability tolerance: choose it so that the bound C_{k,D}y^{-D} is below that tolerance, and then take n large. A single height is not required to handle all tolerances. Finitely many initial sizes do not affect tightness. Here is a direct way to see the needed compact bounds. For a scalar resolvent sum f with nonnegative residues, and a compact set \mathcal K in the upper half-plane,

Both inequalities follow by comparing the corresponding kernels on the real line. Apply them to m_n, and to the real quadratic forms of M_{n,k}; polarization controls the matrix entries. Thus bounded values at \mathrm i y give uniform bounds on each fixed compact, and Cauchy’s estimates give equicontinuity on smaller compacts.

We may consequently extract the asserted joint limit by tightness, first on growing compact sets and finitely many channel dimensions, and then diagonally. Local vague limits with locally bounded counts are point measures with integer multiplicities. Limits of the analytic functions have nonnegative imaginary parts in the same scalar or matrix sense.

Conditional Gaussian marks on compact intervals. Choose nested compact position intervals exhausting \R, with boundaries almost surely uncharged by \nu. Such deterministic boundaries can be chosen because the set of points with positive probability of carrying an atom is countable. On each interval, truncate its point count at a fixed integer. Order the prelimit indices using graph data alone, resolving ties by the chosen eigenbasis ordering. The positions and the corresponding marks are tight: their number is bounded, and the conditional second-moment estimates of Proposition 3.1 control the marks. Extract these marked lists together with the preceding limit.

For completeness, fixed-order moment matching suffices here even though the selected indices depend on the graph. For any fixed number of indices, the approximation is uniform over their choices. By a diagonal selection of the graph exceptional events, one may require that an increasing finite collection of moment orders match with errors tending to zero, on events of probability tending to one. More explicitly, for each fixed list length and k, choose integers J_n\to\infty sufficiently slowly that all orders at most J_n have errors at most J_n^{-1} on a common graph event of probability at least 1-J_n^{-1}. This uses only finitely many fixed-order assertions at each stage, not a uniform growing-order estimate. If the conditional mark laws failed to approach their product Gaussian law in a metric for weak convergence, one could select a deterministic sequence of graphs and index lists with this failure while every fixed moment converged. Tightness and moment determinacy of the Gaussian law contradict that selection. Thus conditional weak convergence is uniform over these bounded lists.

Bounded Lipschitz tests, first with the list length truncated, now show that the same joint limit is obtained by assigning independent Gaussian marks to the prelimit positions independently of the graph. Letting the truncation increase identifies the conditional marking law given the full unmarked process \nu. The covariance in 2k endpoint coordinates is C_k. At a repeated limiting position the list can be retained on an extension of the probability space, with a separate independent mark for each multiplicity; its summed outer product does not require an intrinsic ordering of the tied points. Diagonal extraction in the intervals and in k gives consistent marks

for all finite channel dimensions.

Removing imaginary tails. The convergence of compact marked sums alone does not yet identify the full resolvents. For fixed z in the upper half-plane and large L, the scalar kernels satisfy, when |x|>L,

The constants can be uniform for z in a fixed compact set. For matrix sums use positivity of the residues and bound the norm by the trace of the corresponding positive imaginary part. Proposition 4.1 shows that these traces at \mathrm i L are tight as L tends to infinity, after taking the upper limit in n. The prefactor L^{-1} in (4.17) therefore removes the tails in probability. Compact marked sums already converge, so

One may establish these identities first at countably many arguments and then use positivity and continuity. The imaginary sums and the sums of absolute values of the difference kernels are finite. In particular no imaginary spectral mass, including a linear Herglotz term, is lost at infinity. Real additive constants are not excluded by this argument and will be identified below.

Uniform bounds after the microscopic limit. At every fixed s+\mathrm i y, Proposition 4.1 passes to the limit by the portmanteau inequality for the open event where its error exceeds the stated threshold. Its constants remain independent of the fixed real part s. Decrease \xi if necessary so that \xi<1.

Fix 0<\delta<1. In a padded region

take a unit grid. There are O(R^3) grid points. For each fixed dyadic R, this is a finite collection of arguments in the exact limiting law. The failure probability from (4.2) is at most C R^{3-\delta D}. Choose D sufficiently large and apply Borel–Cantelli over dyadic R. The kernel identities in (4.18) give

Imaginary parts at arguments within bounded distance and comparable height are comparable by their positive kernels. The error in filling the unit grid is consequently O(1/u). Padding allows all sufficiently large R, not just dyadic values. This proves (4.13). The limit in n preceded the growing grid argument; no finite-n uniform estimate on expanding real windows has been used. A countable set of \delta’s implies all the stated fixed-\delta conclusions, by slightly decreasing \delta when necessary.

From the Poisson transform to counts. Equation (4.13) first gives \nu([-R,R])=O(R), by evaluating \Im m(\mathrm i R). For -R\le v\le w\le R, integrate \pi^{-1}\Im m(e+\mathrm i\sqrt R) over e\in[v,w]. The result is

To compare it with the sharp count, put

The number of points within distance R^{3/4} of either endpoint is O(R^{3/4}): use (4.13) at that endpoint and height R^{3/4}. For points of [-2R,2R] outside these two strips,

Their total error is O(R^{3/4}), since the number of points is O(R). Finally, for |x|>2R, P_{v,w}(x)\le C R^{3/2}/x^2; the linear count bound and dyadic summation make this contribution O(\sqrt R). These estimates are uniform in v,w, and prove (4.14) for any sufficiently small deterministic

The estimates also show that endpoint conventions have no effect on the stated error bound.

At this point the limiting positions have a deterministic density with a power-saving count error, and all finite spectral residues have been identified. It remains to determine the real constants that the compact convergence did not see.

Scalar principal values and their constant. Partial summation in (4.14) gives, almost surely,

Only large R is involved, so possible points at zero do not enter these expressions. To see the first assertion, apply the count discrepancy separately on the positive and negative half-lines; their leading density terms cancel in the symmetric sum, and integration of O(t^{1-\chi})/t^2 leaves O(R^{-\chi}). The second assertion follows directly from the linear count bound.

Thus

exists locally uniformly and is analytic. Comparison with the constant-density measure gives

Indeed, on each half-line the counting error is O(1+|v|^{1-\chi}). Partial summation bounds its contribution by a constant times

whereas the symmetric transform of constant density is \mathrm i\pi\rho. By (4.18), m-S_\nu is a constant function, and its imaginary part is zero. Equations (4.13) and (4.20) identify that real constant as m_{\rm r}.

Marked principal values and their constant. Fix k and condition on \nu. Set

These matrices are independent and centered, with bounded second moments depending only on k,d,E. The series \sum_j Z_j/(x_j-\mathrm i), ordered by increasing |x_j|, converges almost surely entrywise by square summability and independence. Also, conditionally almost surely,

because its conditional expectation is finite by (4.19). Consequently the corresponding difference series between z and \mathrm i converges absolutely and locally uniformly. Together with the scalar principal value this proves local uniform convergence of

The linear count bound also gives

Thus every entry of the centered marked series at \mathrm i R has conditional second moment O(R^{-1}), with a sample-dependent constant. It follows that

The limiting version of Proposition 4.1 gives M_k(\mathrm i R)\to K_k(m_*) in probability. On the other hand (4.18) shows that M_k-S_{k,\nu} is almost surely constant in z. That constant must therefore equal

almost surely. This proves both identities in (4.15).

Finally, (4.14) supplies infinitely many points, and C_k is positive definite. Any 2k independently Gaussian marks span \R^{2k} almost surely. Their positive contributions to the imaginary-part sum in (4.18) make \Im M_k(z) positive definite for every z in the upper half-plane on the same event. Countability of the channel dimensions completes the simultaneous assertion.□

Proof. First suppose that O is a permutation matrix. In a finite graph, let W=W_{n,k} contain the starting endpoint columns followed by the ending endpoint columns of the k sampled oriented edges. Except with probability O_k(n^{-1}), all endpoints are distinct and their induced subgraph is exactly the marked matching. On this event,

permutes the pairing between the two endpoint lists. The resulting graph is simple and d-regular. Its marked list has the same restriction, and the inverse permutation recovers the input. Every graph with its ordered list has the same joint weight, so this bijection preserves the restricted law.

Write R(z)=(A-E-z/n)^{-1}. The determinant and inverse identities for finite-rank perturbations give

Indeed m_n=n^{-1}\Tr R, and differentiation in z supplies exactly this factor n^{-1}. The formulas also follow without requiring invertibility of the prelimit compression by using (\Id+M_{n,k}B_O)^{-1}M_{n,k}. In the limit, the positive definite imaginary part already established makes the inverse formula regular. Compact-open convergence of analytic functions also gives convergence of their derivatives. We may therefore pass to the limit, proving the proposition for permutation matrices.

Next let U\in O(k), and set Q_U=\diag(U,U). The matrices C_k and K_k(m_{\rm r}) commute with Q_U. The Gaussian representation (4.15) consequently shows that

preserves the joint law. This assertion follows even conditionally on \nu, but we will need switching invariance only unconditionally. Since Q_UD_OQ_U^\top=D_{UOU^\top} and Q_UB_OQ_U^\top=B_{UOU^\top}, the inverse and determinant formulas give

Thus the law is invariant under every orthogonal conjugate of a permutation switch.

Finally, the switching maps form an action: applying \mathcal T_O and then \mathcal T_{O'} gives \mathcal T_{OO'}. To check this directly, observe that

Successive inverse updates therefore give the asserted second component. The determinant ratios for the two updates multiply, which gives the first component after logarithmic differentiation. A transposition is a hyperplane reflection. Its orthogonal conjugates are all hyperplane reflections, and their finite products generate O(k). The proposition follows.□

Proof. The centered matrix in (5.16) is a sum of independent self-adjoint matrices

where v_\alpha is Gaussian of covariance P. Its coefficients may have either sign. Gaussian fourth moments give

Truncate each original Gaussian column on the event that all its entries have absolute value at most \log r, and center the truncated summand. There are only polynomially many entries, so this alters the matrices with superpolynomially small probability. Gaussian tails also make the resulting centering correction superpolynomially small.

Each truncated summand has norm at most

If X is a self-adjoint random matrix, then \EE[(X-\EE X)^2]=\EE[X^2]-(\EE X)^2\le\EE[X^2]. Also truncating v_\alpha v_\alpha^\top by an event can only decrease the expectation of its square in positive semidefinite order. Thus (5.18), including for negative \beta_{\alpha\alpha}, bounds the matrix variance parameter by

The bounded self-adjoint matrix Bernstein inequality [24, Theorem 1.4] gives, for a sum of independent centered self-adjoint matrices of norm at most b and variance parameter \sigma_{\rm mat}^2, a two-sided tail bounded by

Its dimension factor is polynomial in r. Taking a sufficiently large power of \log r in t proves (5.16).

For the cross block, independence of distinct Gaussian rows and direct second-moment calculation give

Equation (5.17) follows from (5.15) and Markov’s inequality, with one additional logarithmic factor.□

Proof. Conjugating the matrix to be inverted by U_F gives

The second term has norm at most \delta_F, by Lemma 5.3, the bound |s_1|=O(T), and (5.15). Since \delta_F=o(1), a Neumann expansion proves (5.19). In particular, Q_F=U_F^{-1}G_FU_F^{-1} has bounded norm. By (5.10),

This proves the assertion about the far trace.

For the Schur correction, replacing G_F by D_F^{-1} costs at most

in the norm weighted on both sides by U_I. It remains to estimate A'_{IF}D_F^{-1}A'_{FI}. Condition further on the inner Gaussian row block \mathcal G_I, and set

The columns of A'_{IF}, indexed by j\in F, are then independent centered Gaussian vectors a_j with covariance Q. Consequently

Also \EE[Q_{aa}^2\mid\nu] =s_2^2+2\Tr\beta^4=O(T^2). Since Q is positive semidefinite, |Q_{ab}|^2\le Q_{aa}Q_{bb}; Cauchy–Schwarz bounds the remaining second moments by O(T^2). Multiplying entry a,b by (U_I)_{aa}(U_I)_{bb}, summing its second moment over a,b, and using (5.10) and (5.15), we obtain

Markov’s inequality therefore bounds this norm by T\operatorname{polylog}(r)(K^{-\chi}+K^{-1/2}) with probability tending to one.

Combining this estimate with (5.21) gives

with probability tending to one. After multiplication by T^5, all terms tend to zero by a fixed power of r. For example, the three potentially slow terms have exponents

respectively. They are negative under (5.4); the terms T^6K^{-1/2} and T^6(1+T)/K have still ample margins. This proves (5.20).□

Proof. Every diagonal entry of U_I(D_I-ih)U_I has modulus at least 1/\sqrt2. The added term has norm at most

with probability tending to one. Hence

The identity G_I(ih)^*G_I(ih)=h^{-1}\Im G_I(ih) now gives

Alternatively the same mixed-weight estimates follow from the preceding weighted inverse and \|U_I\|\le h^{-1/2}. Use the exact twice-iterated resolvent identity

The unweighted middle resolvent in the last term has norm at most (\Im z)^{-1}. Conjugating the identity by U_I^{-1}, the three terms are bounded by

respectively. This proves (5.22).

The Schur complement formula for the full inverse has upper-left block \widetilde G_I=(G_I^{-1}-\Sigma)^{-1}. Set R_I=U_I^{-1}G_IU_I^{-1} and S_I=U_I\Sigma U_I. Then

By (5.22) and (5.20), a Neumann expansion gives

Multiplying the last bound by \Tr U_I^2=O(\log r) shows that their traces differ by o_{\PP}(1).

The full block-inverse trace is

The far trace is o_{\PP}(1) by Lemma 5.4. For the last term recall Q_F=U_F^{-1}G_FU_F^{-1}. It has bounded norm, and

No positivity is asserted for this complex square. The ordinary operator-norm and Hilbert–Schmidt trace inequality instead gives

This proves (5.23). All estimates used the same Gaussian matrices, so taking a finite union gives the joint assertion.□

Proof. Write

where the g_\alpha are independent standard q_*-dimensional Gaussian vectors. The covariance identity

shows that gg^\top-\Id and a raw GOE matrix have the same entry covariances. Their means are both zero. Thus replace the summands one at a time by \beta_{\alpha\alpha}\mathcal W_\alpha, with independent raw GOE matrices \mathcal W_\alpha. Negative coefficients are allowed: both covariances are multiplied by \beta_{\alpha\alpha}^2.

For either kind of summand X,

For the rank-one term this follows from Gaussian moments of \|g\|^2; for the GOE term the Hilbert–Schmidt norm gives a sufficient bound. Fix spectral parameters z_1,\ldots,z_b\in\C_+, and a bounded smooth test \Phi on the corresponding real and imaginary trace coordinates, with bounded derivatives. For any real symmetric matrix A, its resolvents satisfy \|(A-z_\ell)^{-1}\|\le(\Im z_\ell)^{-1}. Differentiating the resolvent in a real symmetric direction \Delta gives, for j=1,2,3, a trace derivative bounded by C_{z,j}q_*\|\Delta\|^j. The chain rule therefore bounds the third directional derivative of the composed test by

This estimate holds at every intermediate real symmetric matrix, independently of eigenvalue gaps.

Taylor expansion to second order cancels the means and covariances in each replacement. Summing the third-order remainders gives total error

Finally \sum_\alpha\beta_{\alpha\alpha}\mathcal W_\alpha has exactly the law \sqrt{s_2}\,\mathcal W, proving the lemma.□

Completion of Proposition 5.2. The coefficient bounds were proved after (5.13). That equation, Lemma 5.5, and Lemma 5.6 prove the conditional comparison with \mathcal S_r in the proposition. Each comparison concerns bounded tests, so its conditional error can be averaged by bounded convergence. By Proposition 5.1, the unconditional law of \mathcal S_r is the law of m for every r. This proves (5.9); bounded smooth tests suffice to determine finite-dimensional weak convergence.

We finish by explaining why these trace limits imply convergence of point measures, without assuming a uniform bound on their total number of points. Let \nu^{(r)} be the eigenvalue counting measure of H_r, and define

Its total mass is \Im\Tr(H_r-i)^{-1}, whose laws are tight by (5.9). Regarding \sigma_r as a finite positive measure on the one-point compactification \overline{\R}=\R\cup\{\infty\}, its laws are therefore tight. For z=e+i\eta\in\C_+, the kernel

extends continuously to \overline{\R}, with k_z(\infty)=\eta. Its integral against \sigma_r is \Im\Tr(H_r-z)^{-1}. It follows that in any subsequential limit \sigma, the joint law of these integrals at a countable dense set of z’s is that of \Im m(z). Continuity extends this conclusion to all z.

These integrals determine the de-weighted restriction of \sigma to \R. Indeed, for f\in C_c(\R), integrate the kernel against \pi^{-1}f(e)\dd e and let \eta\downarrow0. On \R the resulting kernel tends to (1+x^2)f(x), by the Poisson approximate identity. The integrated kernels are bounded uniformly over x\in\overline{\R} and 0<\eta\le1: near the compact support this is the Poisson mass bound, and away from it the factor 1+x^2 cancels the quadratic denominator. At infinity their value is \eta\pi^{-1}\int f(e)\dd e, which tends to zero. Dominated convergence is thus applicable to each finite measure \sigma, including any mass at infinity.

By (4.15), applying the same inversion to \Im m(z) gives \int f\,\dd\nu. Taking a countable determining family of compactly supported tests, the joint distribution of these integrals is therefore exactly the law of \nu. The map from \sigma to its de-weighted restriction is continuous for local vague convergence: for a compactly supported f, the function (1+x^2)f(x), extended by zero at infinity, is continuous on \overline{\R}. Every subsequential local vague limit of \nu^{(r)} consequently has the law of \nu. This proves the point-process assertion.□

Proof. The marginal statement is Theorem 2.2 of Landon et al. [21], with Definition 2.1 and Equations (2.5)–(2.9) there, specialized to energy zero and equal densities \rho_{\mathrm{fc},t}(0)=\rho_{\mathrm{sc}}(0)=1/\pi. The energy lies in the required interior window for any fixed interior fraction. That theorem gives a power-decaying bound for (). To pass to the last assertion, multiply each integral by

This is the Jacobian and ordered-tuple multiplicity for the unit-mass marginal convention, and it tends to \pi^a for fixed a.□

Proof. First, \|V_*\| is bounded by (6.7), since T=o(K). We verify a slightly enlarged regularity domain,

so that subsequent dilation by c_r=1+o(1) will be harmless. The normalized resolvent trace is

Writing e=q_*\Re\zeta/\alpha_*-s_1 and u=q_*\Im\zeta/\alpha_*, the enlarged domain has

The lower height exceeds a fixed positive power of K. The uniform complex resolvent estimate controls the Poisson sum down to heights of order T^{1/2}, while the count discrepancy controls truncation at K. Equation (4.13), with a radius comparable to K and any fixed exponent smaller than 50\tau, therefore gives

uniformly in (6.9); here we used the imaginary part of (4.15). For u/K sufficiently small, the omitted terms with |x_j|>K sum to O(u/K), uniformly for |e|=o(K), by the linear count bound following from (4.14). Hence the retained sum has a positive lower bound. For u/K bounded below, every retained point has |x_j-e|\le2K eventually, so

Its upper bound follows by positivity from the full sum. This proves regularity on the enlarged domain.

With m=q_*, the relevant powers are

In particular,

The inequalities m^{\tau-1}\le g\le m^{-\tau} and G\le m^{-\tau} follow as well. Our bound on \tau makes all these exponent inequalities strict.

It remains to determine the free-convolution density. At \zeta=0, Equation (6.3) becomes

For |w-i|\le1/4, the argument z=\alpha_*s_2w-s_1 has \Im z\asymp T and |z|=O(T). Equations (4.13) and (4.15) imply

uniformly in this disk. The real constant m_{\rm r} has been subtracted in this formula. Furthermore,

The first bound is partial summation of (4.14); the second uses its linear count bound and |z|=o(K). Thus F_r(w)=i+o(1) uniformly on |w-i|\le1/4.

Rouché’s theorem applied to w-F_r(w) and w-i shows that (6.10) has exactly one zero, counted with multiplicity, in this disk for all large r. The zero w_r is therefore simple, and w_r\to i. The implicit function theorem continues it to a solution near \zeta=0 of the full fixed-point equation, still with positive imaginary part. This solution is the physical branch. Indeed, every upper-half-plane solution w, at \Im\zeta>0, satisfies

If two distinct such solutions existed, subtracting their equations and dividing by their difference would give

Cauchy–Schwarz and (6.11) make the absolute value of the right-hand side strictly less than one. This contradiction proves uniqueness and identifies the analytic continuation. Its boundary value gives \rho_{\mathrm{fc},t_*}(0)=\Im w_r/\pi\to1/\pi.

Finally, simultaneous replacement of V_*,t_* by c_rV_*,c_r^2t_* dilates the free-convolution law by c_r. For example its transform is c_r^{-1}m_{\mathrm{fc},t_*}(\zeta/c_r), as follows directly from (6.3). Its density at zero is exactly

Since c_r\to1, the enlarged regularity domain proves (6.1) for c_rV_* with the original g,G. The strict power margins preserve (6.4) for c_r^2t_*, completing the proof.□

Proof. We bound the GOE count by twice a Gaussian unitary count plus one, using the Gaussian decimation identity of Forrester and Rains [15, arXiv version 2, Theorem 5.2 and Equation (5.9)]. We derive that identity in our normalization, then use the unitary ensemble’s projection kernel to bound its count moments.

Gaussian decimation. Write \mathrm O_m for the decreasingly ordered points with density proportional to

They have the law of the eigenvalues of \sqrt m\,W_m/\sqrt2. Thus the count \mathcal G_m(J) is the \mathrm O_m count in

Superpose independent \mathrm O_m and \mathrm O_{m+1}, and order the resulting points as z_1>\cdots>z_{2m+1}. Apart from a constant and the common Gaussian factor, their density is the polynomial

where each sublist is in the inherited order. The parity–Vandermonde factorization of Forrester and Rains [15, arXiv version 2, Lemma 4.1 and Equation (4.2)] is

If z_i=z_j with i,j of the same parity, terms having both indices in one sublist vanish. The other terms cancel in pairs under exchanging membership of i,j: the sign ratio is (-1)^{j-i-1}=-1, since each intervening index belongs to exactly one sublist. Hence every factor on the right divides P. Both sides have degree m^2, so the quotient is constant. It is positive on strictly decreasing z, proving (6.13).

Fix the even points y_i=z_{2i}. The odd points interlace them, one in each of

To integrate their Vandermonde times their Gaussian factors, use the monic polynomials

Writing the Vandermonde as the determinant of these polynomials reduces the integral, by multilinearity, to a determinant of integrals over the displayed intervals. Replace each row by the sum of it and all preceding rows. For column j\ge1, the first m rows become

and the last row is zero. In column zero the last row is \int_\R b_0. Expansion along that last row leaves a Vandermonde in the y_i times \prod_i b_0(y_i). More explicitly, the ascending-degree polynomial columns give signs (-1)^{m(m+1)/2} for the odd-point Vandermonde and (-1)^{m(m-1)/2} for the even-point Vandermonde, while expansion along the last row contributes (-1)^m. Their product is one. Thus, if x_0>y_1>x_1>\cdots>y_m>x_m, the integral just evaluated is exactly

Together with the even-point factors already present in (6.13), this proves that the even points have density proportional to

Projection kernel and count moments. The even-point density now gives a direct bound on count moments. Let \psi_j(x) be the degree-j orthonormal polynomial for the weight e^{-x^2}, multiplied by e^{-x^2/2}. The kernel

is the rank-m orthogonal projection in L^2(\R) onto their span. The symmetric density in (6.14) is (m!)^{-1}\det(K_m(y_i,y_j))_{i,j=1}^m. Indeed the determinant is a squared weighted Vandermonde, and orthonormality gives its integral m!. More generally its order-a factorial moment density is \det(K_m(y_i,y_j))_{i,j=1}^a. For completeness, expand a size-a determinant by Cauchy–Binet as the sum of squared minors indexed by a basis functions. Integrating its last variable uses their orthogonality, and each minor of size a-1 occurs m-a+1 times. Iterating this identity proves the factorial density formula with its stated normalization.

We need a bound uniform in both m and location. The operator -\partial_x^2+x^2 has eigenvalue 2j+1 on \psi_j. This follows by conjugating it by e^{-x^2/2}: on the polynomial factor it acts as -p''+2xp'+p, preserves successive degree spaces, is symmetric for the Gaussian inner product, and has leading coefficient multiplier 2j+1. Consequently every f in the range of K_m satisfies

The elementary bound |f(x)|^2\le2\|f'\|_2\|f\|_2, obtained by integrating the derivative of |f|^2, gives

If N_m^{\mathrm{even}} is the number of even points in J_m, the Gram determinant inequality yields, for every a\ge0,

For a>m the left side is zero. Hence, for every v>0,

In any interval the total superposition count is at most twice its even count plus one. The original \mathrm O_m count is no greater than the superposition count. Thus

which proves the assertion.□

Proof of Theorem 1.1. Start with any sequence of admissible graph sizes tending to infinity. Proposition 4.2 supplies a further subsequence on which the unscaled microscopic measures \nu_n converge in law to a process \nu of the form used above. The graph point measure in the theorem is \mathcal D_\rho\nu_n. Bounded continuity of its Laplace functional and (6.19) show that its Laplace expectation, along this further subsequence, tends to the Laplace expectation of \mathcal S_{\mathrm{GOE}}. The reference GOE expectations at dimension n tend to the same value, along all sizes. Every original subsequence therefore has a further subsequence along which their difference tends to zero, which proves the required full-sequence limit. All spectral estimates and negligible exclusions used to construct the subsequential law were events in the original simple-graph model; no conditioning on connectedness, bipartiteness, or any additional graph property enters the conclusion.□