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 \[\frac{d-1}{2}\quad\hbox{and}\quad\frac{(d-1)^2}{4},\] respectively. Indeed, a prescribed family of \(e\) disjoint pairs occurs with probability \[\frac{1}{(nd-1)(nd-3)\cdots(nd-2e+1)} \sim(nd)^{-e}\] for fixed \(e\).

For each fixed \(j\), ordered \(j\)-tuples of defects on disjoint vertex sets give the leading contribution to \(\mathbb E(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 \[\mathbb E(Z)_j\longrightarrow \left(\frac{d-1}{2}+\frac{(d-1)^2}{4}\right)^j.\] The finite inclusion-exclusion bounds for \(\mathbb P(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 (Huang and Yau 2024, 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 \[0<c<1,\qquad a_{\mathrm{HY}}\ge12,\qquad b_{\mathrm{HY}}\ge25a_{\mathrm{HY}},\qquad 0<\mathfrak r<c/32,\] 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 \[(d-1)^{-r_{\mathrm{HY}}} +(n\Im\zeta)^{-1/2} +(n\Im\zeta)^{-2/3}.\] 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\), \[m_d(\zeta)+(d-1)g_2 =\zeta\,\frac{1+\zeta m_d(\zeta)}d.\] 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 \[\Re\int\log(\zeta-t)\rho_d(t)\,\mathrm dt =-\log|q|-\frac{d-2}{2}\log|1-q^2|.\] At the bulk boundary, \[|q|=(d-1)^{-1/2},\qquad |1-q^2|=\frac{\sqrt{d^2-E^2}}{d-1}.\] These identities give (3) after substituting \(a=d(1-p^2)/(d-1)\). The boundary logarithmic integral exists because \(\rho_d\) is bounded near \(E\). Finally, \[\int \left(\log|t-E-\mathrm i\eta|-\log|t-E|\right)\rho_d(t)\,\mathrm dt =\frac12\int\log\!\left(1+\frac{\eta^2}{(t-E)^2}\right)\rho_d(t)\,\mathrm dt =O(\eta),\] 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 \[\frac{\displaystyle\prod_{j<k}m_{jk}! \prod_j(m_{jj}-1)!!}{(nd-1)!!},\] with zero-size factors omitted. Stirling’s bounds give its logarithm as \[-\frac{nd}{2}H(Q)+O(M^2\log n).\] 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 (4) with a slightly larger sublinear allowance. The simplicity probability is bounded away from zero by Lemma 2, 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 7 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 \[ \frac1nH(\mathcal L) \le H(\mathbb ET)-\frac d2\mathbb EH(Q)+O(n^{-c}). \tag{9}\] Here concavity was used only for the positive star term.

We next replace \(\mathbb EH(Q)\) by \(H(\mathbb EQ)\). Two independent edge locations have Gaussian observations whose off-diagonal covariance block has average Frobenius norm at most \[C\left(\frac{\mathop{\mathrm{Tr}}\Gamma^2}{n^2}\right)^{1/2} \le \frac{C}{\sqrt n\,\eta}.\] Indeed endpoints from different sampled edges are independent uniform vertices, and \(\mathop{\mathrm{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 \[\mathbb E_{\rm locations}\, \mathop{\mathrm{dist}}_{\mathrm{TV}}\bigl(\mathop{\mathrm{Law}}(Y_{\rm edge},Y_{\rm edge}'), \mathop{\mathrm{Law}}(Y_{\rm edge})\otimes \mathop{\mathrm{Law}}(Y_{\rm edge}')\bigr) \le \frac{C}{\sqrt n\,\eta^2},\] 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 \[\mathbb E\|Q-\mathbb EQ\|_1 \le n^{4\ell\gamma-1/4+\alpha+o(1)}=O(n^{-c}).\] The finite-alphabet entropy continuity bound \(H_{\rm bin}(h)+h\log(M^2)\) at total variation distance \(h\) now proves \[ \mathbb EH(Q)=H(\mathbb EQ)+O(n^{-c}). \tag{10}\]

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 \[\frac1nH(\mathcal L) \ge \frac1n\mathcal H(Y_{\rm global}) -\ell\log b-O(n^{-c}).\] 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 \[ H(\text{binned }X) =\mathcal H(X)-(\dim X)\log b+O(n^{-c}). \tag{11}\] 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 (11), and will also apply to fixed-dimensional projections below.

The laws \(\mathbb ET,\mathbb EQ\) are exactly the binned local laws. Substitute (10) and (11) into (9). The mesh terms cancel, since \((d+1)\ell-(d/2)(2\ell)=\ell\). This proves (7).

Finally the Gaussian determinant formula gives \[\frac1n\mathcal H(Y_{\rm global}) =\frac{\ell}{2}\log\frac{2\pi{\rm e}\eta}{\pi\rho} -\frac{\ell}{n}\sum_i\log|\lambda_i-E-i\eta|.\] 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 5 gives the remaining regularization error \(O(\eta)\). This proves (8). ◻

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 \[(A-E)\Gamma(A-E) =\frac{\eta}{\pi\rho} \frac{(A-E)^2}{(A-E)^2+\eta^2} \le\frac{\eta}{\pi\rho}I.\] 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=\mathop{\mathrm{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 \[ \mathop{\mathrm{Cov}}\bigl(R_{\rm res}(F_{\rm star})\bigr) \le \frac{\eta}{\pi\rho}I_\ell+\epsilon^2\Delta^2 \le \left(\frac{\eta}{\pi\rho}+\frac{S^4}{n^2}\right)I_\ell. \tag{13}\] The relative correction inside the logarithmic entropy bound satisfies \[\log\left(1+\frac{\pi\rho S^4}{n^2\eta}\right) =O\bigl(n^{-2+\alpha+o(1)}\bigr).\] The maximum-entropy bound by covariance and the fact that conditioning lowers entropy thus give \[\mathcal H(R_{\rm res}(F_{\rm star})\mid \mathcal B(F_{\rm star})) \le \frac{\ell}{2}\log\frac{2\pi{\rm e}\eta}{\pi\rho} +O(n^{-c}).\]

By the typical-star estimate at \(E+i\eta\) and the tree covariance in Lemma 5, the covariance of the noise \(\mathcal B(Y_{\rm star})\) differs from \(I_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 8, 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 (7) is at most the left side of (12) plus \[\ell\left( \frac12\log\frac{2\pi{\rm e}\eta}{\pi\rho} -\frac12\log d+\frac{d-1}{2}\log a -\frac d4\log(1-p^2)\right)+O(n^{-c}).\] Combine (7), (8), and the potential identity (3). Their constants agree exactly, proving (12).

For the remaining assertions, put \(\Lambda=\mathop{\mathrm{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 \[\mathbb EZZ^\top=\mathbb EX_iX_i^\top=I_\ell,\qquad \mathbb EZX_i^\top=\Lambda/d,\qquad \mathbb EX_iX_j^\top=\frac{\Lambda^2-dI_\ell}{d(d-1)} \quad(i\ne j).\] The mean is zero because these eigenvectors are orthogonal to constants. These identities imply \(\mathbb EB_VB_V^\top=I_D+O(n^{-1+o(1)})\), using \(\Lambda=EI_\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 (2019, 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 \[e_i^+=\frac{w+v_i}{\sqrt{2(1+p)}},\qquad e_i^-=\frac{w-v_i}{\sqrt{2(1-p)}},\qquad s_i=\frac{e_i^++e_i^-}{\sqrt2},\qquad t_i=\frac{e_i^+-e_i^-}{\sqrt2}.\] The swap fixes \(e_i^+\), negates \(e_i^-\), and hence interchanges the orthonormal directions \(s_i,t_i\). Their squared central components are \[\frac{1+\sqrt{1-p^2}}2,\qquad \frac{1-\sqrt{1-p^2}}2.\] Writing \(q_0=\sqrt{1-p^2}\), the strict bulk condition is \[p^2<\frac{4(d-1)}{d^2} \quad\Longleftrightarrow\quad q_0>\frac{d-2}{d}.\] For \(d\ge3\), the two squared central components therefore strictly bracket \(1/d\). The weights \[b_s=\frac{q_0-(d-2)/d}{2q_0},\qquad b_t=\frac{q_0+(d-2)/d}{2q_0}\] are positive, sum to one, and satisfy \[ \sum_{i=1}^d\bigl(b_s s_is_i^\top+b_t t_it_i^\top\bigr) =I_{\mathbb R\oplus{\bf1}^{\perp}}. \tag{16}\] 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 (16) 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 (15).

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\mathbb R^\ell\) and \(t_i\otimes\mathbb R^\ell\), then \[R_{i,s}^{(1)}H =R_{i,s}^{(1)}P_i^{(2)}\cdots P_i^{(j)}H,\] and likewise for \(t\); superscripts specify tensor slots. Adding the two identities and using symmetry gives \[ P_i^{(a)}H=P_i^{\otimes j}H \quad\text{for every slot }a. \tag{17}\] 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 \[\begin{align*} 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 . \end{align*}\] 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 \[\langle a_i',a_k'\rangle \bigl(a_i'(a_k')^\top-a_k'(a_i')^\top\bigr), \qquad \langle a_i',a_k'\rangle=-\frac1{d-1}.\] For \(d\ge3\), it is invertible on \(\mathop{\mathrm{span}}(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\mathbb R^\ell)^{\otimes j}\). But its projected swap invariance would also support it in \((v_i\otimes\mathbb 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 (15). ◻

Proof. For \(y\in D\), define the Hermite tensors \(\operatorname{He}_j(y)\) by the generating identity \[\exp\bigl(\langle y,t\rangle-\|t\|^2/2\bigr) =\sum_{j=0}^\infty \frac{\langle \operatorname{He}_j(y),t^{\otimes j}\rangle}{j!}, \qquad t\in D.\] Set \(H_j=\mathbb E\operatorname{He}_j(B_n)\). We prove by induction that \[ \|H_j\|=O(n^{-c_j}) \quad(1\le j\le q). \tag{18}\] 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 \[2j\tau<c_*/2,\qquad (j-k)\tau<c_k/2\quad(1\le k<j).\] 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 \[Z_u=uB_n+\sqrt{1-u^2}\,G.\] Its second-moment matrix is \(I_D+u^2O(n^{-c_*})\). Consequently (14), with this second-moment error included, turns the hypothesis into \[ \mathop{\mathrm{Ent}}(\mathop{\mathrm{Law}}(Z_u)\mid\varphi_D) -\frac12\sum_{i=1}^d \mathop{\mathrm{Ent}}(\mathop{\mathrm{Law}}(P_iZ_u)\mid\varphi_{U_i}) \le O(n^{-c_*}). \tag{19}\]

We claim that \[ \mathop{\mathrm{Ent}}(\mathop{\mathrm{Law}}(Z_u)\mid\varphi_D) =u^{2j}\left(\frac{\|H_j\|^2}{2j!}+O(n^{-c'})\right), \tag{20}\] 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 \[\frac1{k!}\langle H_k,\operatorname{He}_k(y)\rangle.\] This follows by integrating against the generating exponential: the resulting transform is \(\mathbb E\exp(u\langle B_n,t\rangle-u^2\|t\|^2/2)\). Pairing two generating exponentials under Gaussian measure gives Hermite orthogonality and the identity \[\left\|\frac{\langle H_k,\operatorname{He}_k\rangle}{k!}\right\|_{L^2(\varphi_D)}^2 =\frac{\|H_k\|^2}{k!}.\]

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 \[(1-s^2)^{-\dim D/2} \exp\left(\frac{\|y\|^2}{2} -\frac{\|y-sb_0\|^2}{2(1-s^2)}\right).\] 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\,\,\mathrm d\varphi_D=0\), \[\mathop{\mathrm{Ent}}((1+h)\varphi_D\mid\varphi_D) =\frac12\int h^2\,\,\mathrm d\varphi_D +O\left(\int|h|^3\,\,\mathrm d\varphi_D\right).\] 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 (20). The same reasoning applies to each projection.

Substitute these expansions into (19) and divide by \(u^{2j}/(2j!)\). The choice of \(\tau\) gives \[ \|H_j\|^2-\frac12\sum_{i=1}^d \|P_i^{\otimes j}H_j\|^2 \le O(n^{-c''}). \tag{21}\] Hermite tensors commute with orthogonal projections and are equivariant under orthogonal maps, so \[P_i^{\otimes j}H_j =\mathbb E\operatorname{He}_j(P_iB_n)\] in the Gaussian space \(U_i\). Since \(\operatorname{He}_j\) is a fixed polynomial, the assumed moment symmetries give the swap constraints of Lemma 10 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 (21) by a polynomially small amount, since \(\|H_j\|=n^{o(1)}\). Lemma 10 therefore proves (18) and closes the induction. Each step uses only the common constants in the hypotheses, proving uniformity. ◻

Proof of Proposition 6. Condition on a graph in the common events of Lemmas 7 and 9 and Proposition 3. Apply Lemma 11, with \(q\) replaced by \(\max(q,2)\), to \(B_n=B_V\). Complete delocalization gives its deterministic subpolynomial bound. Lemma 9 supplies its first two moments and all fixed-order projected moment symmetries. Averaging (12) 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 7 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.

Proof.

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, \[A\longmapsto A+WB_OW^\top,\qquad W\longmapsto WD_O\] 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 \[\begin{align*} \widehat m_n(z) &=m_n(z)-\partial_z\log\det(I_{2k}+B_OM_{n,k}(z)),\\ \widehat M_{n,k}(z) &=D_O^\top(M_{n,k}(z)^{-1}+B_O)^{-1}D_O. \end{align*}\] Indeed \(m_n=n^{-1}\mathop{\mathrm{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 \((I+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=\mathop{\mathrm{diag}}(U,U)\). The matrices \(C_k\) and \(K_k(m_{\rm r})\) commute with \(Q_U\). The Gaussian representation (36) consequently shows that \[\mathcal R_U(m,M)=(m,Q_U^\top M Q_U)\] 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 \[\mathcal R_U^{-1}\circ\mathcal T_O\circ\mathcal R_U =\mathcal T_{UOU^\top}.\] 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 \[ B_{OO'}=B_O+D_OB_{O'}D_O^\top,\qquad D_{OO'}=D_OD_{O'}. \tag{44}\] 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 (57) is a sum of independent self-adjoint matrices \[\sum_{\alpha=1}^{4r}\beta_{\alpha\alpha} (v_\alpha v_\alpha^\top-P),\qquad P=U_{\mathcal A}^2,\] where \(v_\alpha\) is Gaussian of covariance \(P\). Its coefficients may have either sign. Gaussian fourth moments give \[ \mathbb E[(v_\alpha v_\alpha^\top)^2] =(\mathop{\mathrm{Tr}}P)P+2P^2. \tag{59}\] 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 \[C\|\beta\|(\log r)^2\mathop{\mathrm{Tr}}P \le \theta\operatorname{polylog}(r).\] If \(X\) is a self-adjoint random matrix, then \(\mathbb E[(X-\mathbb EX)^2]=\mathbb E[X^2]-(\mathbb EX)^2\le\mathbb E[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 (59), including for negative \(\beta_{\alpha\alpha}\), bounds the matrix variance parameter by \[C s_2(\mathop{\mathrm{Tr}}P)\|P\| \le \frac{T}{a_{\mathcal A}}\operatorname{polylog}(r).\] The bounded self-adjoint matrix Bernstein inequality (Tropp 2012, 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 \[2|\mathcal A|\exp\left( -\frac{c t^2}{\sigma_{\rm mat}^2+bt}\right).\] Its dimension factor is polynomial in \(r\). Taking a sufficiently large power of \(\log r\) in \(t\) proves (57).

For the cross block, independence of distinct Gaussian rows and direct second-moment calculation give \[\mathbb E\|P_0\|_{\mathrm{HS}}^2 =s_2(\mathop{\mathrm{Tr}}U_I^2)(\mathop{\mathrm{Tr}}U_F^2).\] Equation (58) follows from (56) and Markov’s inequality, with one additional logarithmic factor. ◻

Proof. Conjugating the matrix to be inverted by \(U_F\) gives \[U_F(D_F+A'_{FF}-z)U_F =\operatorname{sign}(D_F)+U_F(A'_{FF}-z)U_F.\] The second term has norm at most \(\delta_F\), by Lemma 16, the bound \(|s_1|=O(T)\), and (56). Since \(\delta_F=o(1)\), a Neumann expansion proves (60). In particular, \(Q_F=U_F^{-1}G_FU_F^{-1}\) has bounded norm. By (51), \[\left|\sum_{j\in F}\frac1{x_j}\right|=O(K^{-\chi}), \qquad |\mathop{\mathrm{Tr}}(G_F-D_F^{-1})| \le\delta_F\mathop{\mathrm{Tr}}U_F^2=o(1).\] This proves the assertion about the far trace.

For the Schur correction, replacing \(G_F\) by \(D_F^{-1}\) costs at most \[ \|P_0\|_{\mathrm{HS}}^2\delta_F \le T\operatorname{polylog}(r)\delta_F \tag{62}\] 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 \[Q=\mathcal G_I\beta^2\mathcal G_I^\top.\] The columns of \(A'_{IF}\), indexed by \(j\in F\), are then independent centered Gaussian vectors \(a_j\) with covariance \(Q\). Consequently \[\begin{align*} \mathbb E\left[\sum_{j\in F}\frac{a_ja_j^\top}{x_j} \,\middle|\,\mathcal G_I,\nu\right] &=Q\sum_{j\in F}\frac1{x_j},\\ \mathop{\mathrm{Var}}\left(\sum_{j\in F}\frac{(a_j)_a(a_j)_b}{x_j} \,\middle|\,\mathcal G_I,\nu\right) &=(Q_{aa}Q_{bb}+Q_{ab}^2) \sum_{j\in F}\frac1{x_j^2}. \end{align*}\] Also \(\mathbb E[Q_{aa}^2\mid\nu] =s_2^2+2\mathop{\mathrm{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 (51) and (56), we obtain \[\mathbb E\left\| U_IA'_{IF}D_F^{-1}A'_{FI}U_I \right\|_{\mathrm{HS}}^2 \le C T^2(\log r)^2(K^{-2\chi}+K^{-1}).\] 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 (62) gives \[\|U_I\Sigma U_I\| \le \operatorname{polylog}(r)\,O\left( TK^{-\chi}+TK^{-1/2}+\frac{T(1+T)}K +T\sqrt{T/K}+T\theta\right)\] 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 \[6\tau-\frac{\chi}{100},\qquad \frac{13}{2}\tau-\frac1{200},\qquad \frac{13}{2}\tau-\frac12,\] respectively. They are negative under (45); the terms \(T^6K^{-1/2}\) and \(T^6(1+T)/K\) have still ample margins. This proves (61). ◻

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 \[\|U_IA'_{II}U_I\| \le C T/h+\operatorname{polylog}(r) \left(\sqrt{T/h}+\theta\right)=o(1)\] with probability tending to one. Hence \[\|U_I^{-1}G_I(ih)U_I^{-1}\|\le C.\] The identity \(G_I(ih)^*G_I(ih)=h^{-1}\Im G_I(ih)\) now gives \[\|G_I(ih)U_I^{-1}\|+\|U_I^{-1}G_I(ih)\| \le C h^{-1/2}.\] 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 \[\begin{align*} G_I(z) ={}&G_I(ih)+(z-ih)G_I(ih)^2\\ &+(z-ih)^2G_I(ih)G_I(z)G_I(ih). \end{align*}\] 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 \[C,\qquad C|z-ih|/h,\qquad C|z-ih|^2/(h\Im z),\] respectively. This proves (63).

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 \[U_I^{-1}\widetilde G_IU_I^{-1} =(R_I^{-1}-S_I)^{-1}.\] By (63) and (61), a Neumann expansion gives \[\|U_I^{-1}\widetilde G_IU_I^{-1}\|=O_{\mathbb P}(h),\qquad \|U_I^{-1}(\widetilde G_I-G_I)U_I^{-1}\| \le O_{\mathbb P}(h^2)o_{\mathbb P}(T^{-5})=o_{\mathbb P}(T^{-1}).\] Multiplying the last bound by \(\mathop{\mathrm{Tr}}U_I^2=O(\log r)\) shows that their traces differ by \(o_{\mathbb P}(1)\).

The full block-inverse trace is \[\mathop{\mathrm{Tr}}\widetilde G_I+\mathop{\mathrm{Tr}}G_F+ \mathop{\mathrm{Tr}}\bigl(\widetilde G_IA'_{IF}G_F^2A'_{FI}\bigr).\] The far trace is \(o_{\mathbb P}(1)\) by Lemma 17. For the last term recall \(Q_F=U_F^{-1}G_FU_F^{-1}\). It has bounded norm, and \[U_F^{-1}G_F^2U_F^{-1}=Q_FU_F^2Q_F,\qquad \|Q_FU_F^2Q_F\|=O_{\mathbb P}(K^{-1}).\] No positivity is asserted for this complex square. The ordinary operator-norm and Hilbert–Schmidt trace inequality instead gives \[\begin{align*} \left|\mathop{\mathrm{Tr}}\bigl(\widetilde G_IA'_{IF}G_F^2A'_{FI}\bigr)\right| &\le \|U_I^{-1}\widetilde G_IU_I^{-1}\|\, \|P_0\|_{\mathrm{HS}}^2\, \|Q_FU_F^2Q_F\|\\ &\le \operatorname{polylog}(r)\,O_{\mathbb P}(T^3/K)=o_{\mathbb P}(1). \end{align*}\] This proves (64). All estimates used the same Gaussian matrices, so taking a finite union gives the joint assertion. ◻

Proof. Write \[A'_{II}-s_1I_{q_*} =\sum_{\alpha=1}^{4r} \beta_{\alpha\alpha}(g_\alpha g_\alpha^\top-I_{q_*}),\] where the \(g_\alpha\) are independent standard \(q_*\)-dimensional Gaussian vectors. The covariance identity \[\mathop{\mathrm{Cov}}(g_ag_b,g_cg_d) =\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc}\] shows that \(gg^\top-I\) 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\), \[\mathbb E\|X\|^3\le C\theta^3q_*^3.\] 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\mathbb 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 \[C_{z,\Phi}q_*^3\|\Delta\|^3.\] 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 \[C_{z,\Phi}q_*^6 r\theta^3 =O\bigl(r^{-0.44+(3/2)\tau}\bigr)=o(1).\] Finally \(\sum_\alpha\beta_{\alpha\alpha}\mathcal W_\alpha\) has exactly the law \(\sqrt{s_2}\,\mathcal W\), proving the lemma. ◻

Completion of Proposition 15. The coefficient bounds were proved after (54). That equation, Lemma 18, and Lemma 19 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 14, the unconditional law of \(\mathcal S_r\) is the law of \(m\) for every \(r\). This proves (50); 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 \[\sigma_r(\,\mathrm dx)=\frac{\nu^{(r)}(\,\mathrm dx)}{1+x^2}.\] Its total mass is \(\Im\mathop{\mathrm{Tr}}(H_r-i)^{-1}\), whose laws are tight by (50). Regarding \(\sigma_r\) as a finite positive measure on the one-point compactification \(\overline{\mathbb R}=\mathbb R\cup\{\infty\}\), its laws are therefore tight. For \(z=e+i\eta\in\mathbb C_+\), the kernel \[k_z(x)=(1+x^2)\frac{\eta}{(x-e)^2+\eta^2}\] extends continuously to \(\overline{\mathbb R}\), with \(k_z(\infty)=\eta\). Its integral against \(\sigma_r\) is \(\Im\mathop{\mathrm{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 \(\mathbb R\). Indeed, for \(f\in C_c(\mathbb R)\), integrate the kernel against \(\pi^{-1}f(e)\,\mathrm de\) and let \(\eta\downarrow0\). On \(\mathbb 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{\mathbb 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)\,\mathrm de\), which tends to zero. Dominated convergence is thus applicable to each finite measure \(\sigma\), including any mass at infinity.

By (36), applying the same inversion to \(\Im m(z)\) gives \(\int f\,\,\mathrm d\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{\mathbb 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. (2019), 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 [eq:dbm-marginal-test]. To pass to the last assertion, multiply each integral by \[\frac{(m)_a}{(m/\pi)^a}, \qquad (m)_a=m(m-1)\cdots(m-a+1).\] 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 (71), since \(T=o(K)\). We verify a slightly enlarged regularity domain, \[|\Re\zeta|\le2G,\qquad g/2\le\Im\zeta\le20,\] so that subsequent dilation by \(c_r=1+o(1)\) will be harmless. The normalized resolvent trace is \[m_{V_*}(\zeta)= \alpha_*^{-1}\sum_{j\in I} \frac1{x_j-(q_*\zeta/\alpha_*-s_1)}.\] Writing \(e=q_*\Re\zeta/\alpha_*-s_1\) and \(u=q_*\Im\zeta/\alpha_*\), the enlarged domain has \[ |e|=O(K^{7/8}+T)=o(K),\qquad cT^{1/2}\le u\le CK. \tag{73}\] 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 (34), with a radius comparable to \(K\) and any fixed exponent smaller than \(50\tau\), therefore gives \[\sum_j\frac{u}{(x_j-e)^2+u^2} =\pi\rho+o(1)\] uniformly in (73); here we used the imaginary part of (36). 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 (35). Hence the retained sum has a positive lower bound. For \(u/K\) bounded below, every retained point has \(|x_j-e|\le2K\) eventually, so \[\sum_{j\in I}\frac{u}{(x_j-e)^2+u^2} \ge \frac{q_*u}{4K^2+u^2}\ge c>0.\] Its upper bound follows by positivity from the full sum. This proves regularity on the enlarged domain.

With \(m=q_*\), the relevant powers are \[g=m^{-1+50\tau+o(1)},\qquad t_*=m^{-1+100\tau+o(1)},\qquad G=m^{-1/8}.\] In particular, \[\frac{t_*}{g m^\tau}=m^{49\tau+o(1)}\longrightarrow\infty, \qquad \frac{t_*}{m^{-\tau}G^2} =m^{-3/4+101\tau+o(1)}\longrightarrow0.\] 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 (67) becomes \[ w=F_r(w),\qquad F_r(w)=\alpha_*^{-1}\sum_{j\in I} \frac1{x_j-(\alpha_*s_2w-s_1)}. \tag{74}\] For \(|w-i|\le1/4\), the argument \(z=\alpha_*s_2w-s_1\) has \(\Im z\asymp T\) and \(|z|=O(T)\). Equations (34) and (36) imply \[\mathop{\mathrm{PV}}\sum_j\frac1{x_j-z}=i\pi\rho+o(1)\] uniformly in this disk. The real constant \(m_{\rm r}\) has been subtracted in this formula. Furthermore, \[\mathop{\mathrm{PV}}\sum_{|x_j|>K}\frac1{x_j-z} =\mathop{\mathrm{PV}}\sum_{|x_j|>K}\frac1{x_j} +\sum_{|x_j|>K}\frac{z}{x_j(x_j-z)} =O(K^{-\chi}+T/K).\] The first bound is partial summation of (35); 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 (74) 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 \[ \frac1{q_*}\sum_{j\in I} |(V_*)_{jj}-\zeta-t_*w|^{-2} =\frac{\Im w}{\Im\zeta+t_*\Im w}<\frac1{t_*}. \tag{75}\] If two distinct such solutions existed, subtracting their equations and dividing by their difference would give \[1=\frac{t_*}{q_*}\sum_{j\in I} \frac1{((V_*)_{jj}-\zeta-t_*w_1) ((V_*)_{jj}-\zeta-t_*w_2)}.\] Cauchy–Schwarz and (75) 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 (67). Its density at zero is exactly \[c_r^{-1}\rho_{\mathrm{fc},t_*}(0)=\frac1\pi.\] Since \(c_r\to1\), the enlarged regularity domain proves (65) for \(c_rV_*\) with the original \(g,G\). The strict power margins preserve (68) 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 (2001, 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.

Proof of Theorem 1. Start with any sequence of admissible graph sizes tending to infinity. Proposition 13 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 (83) 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. ◻