=== Abstract === We present a comprehensive numerical study of the entanglement entropy of free fields across a spherical entangling surface, computed on the Srednicki radial lattice with angular momentum decomposition. Through a sequence of controlled experiments, we establish the area-law coefficient $\alphas = 1/(24\sqrt{\pi})$ to $0.01 extrapolation and demonstrate that the logarithmic coefficient $\delta$ is universal across lattice regulators to $0.24 extractions of the vector log coefficient ($\delta_{\rm v} = -31/45$ at $0.01 via order-3 Richardson) and the graviton log coefficient, discovering an approximate "50 $\delta_{\rm TT}/\delta_{\rm EE} = 0.508$ of the full entanglement entropy anomaly for the graviton and $0.515$ for the vector, stable across angular cutoffs to four significant figures. We decompose $\delta$ by angular momentum $l$, finding that all low-$l$ modes contribute positively and that the negative total $\delta = -1/90$ emerges from massive cancellation against the collective high-$l$ sector. A universal scaling function $f(x = l/n)$ governs the per-mode entropy and yields $\alphas$ via a spectral integral confirmed to $0.09 simple closed form; the best fit is a stretched exponential with exponent $c \approx 0.63$. We show that $99.8 boundary mode per angular momentum channel, reducing the conjecture $\alphas = 1/(24\sqrt{\pi})$ to a single-eigenvalue spectral problem. Finally, we prove that the continuum (Bessel function) approximation fails by a factor of $\sim20$, establishing that $\alphas$ is an intrinsically lattice quantity--the discrete structure at the entangling surface is the physics, not an approximation to a continuum limit. === Introduction === The entanglement entropy of quantum fields across a surface has been a central object in quantum gravity and quantum information theory since the work of Bombelli, Koul, Lee, and Sorkin [ref: Bombelli:1986rw] and Srednicki [ref: Srednicki:1993im]. For a free scalar in $3{+}1$ dimensions, the entropy of the vacuum state reduced to a sphere of radius $R$ takes the form S(R) = \alpha \frac{A}{\epsilon^2} + \delta \ln\left(\frac{R}{\epsilon}\right) + \gamma + \order{\epsilon/R} , where $A = 4\pi R^2$ is the surface area, $\epsilon$ is the UV cutoff, $\alpha$ is the area-law coefficient, $\delta$ is the logarithmic coefficient, and $\gamma$ is a finite constant. The area-law coefficient $\alpha$ is UV-sensitive and regulator-dependent: it changes when the lattice action is modified. The logarithmic coefficient $\delta$, by contrast, is determined by the trace anomaly of the conformal field theory and is universal--independent of the regulator, the mass (for $m \ll 1/\epsilon$), and the angular momentum cutoff. For a free scalar, $\delta = -1/90$, fixed by the Seeley-DeWitt coefficient $a_4$ of the conformal Laplacian [ref: Solodukhin:2008dh]. Both $\alpha$ and $\delta$ enter the prediction of the cosmological constant within the entanglement entropy framework developed in companion papers [ref: MoonWalk:Lambda]. There, the observed dark energy density is related to the ratio $\abs{\delta_{\rm total}}/(6 \alpha_{\rm total})$, where the sum runs over all Standard Model fields plus the graviton. The precision of this prediction is limited by the precision of its inputs: $\alpha$ must be computed on the lattice, while $\delta$ is known analytically from the trace anomaly but should be independently verified. This paper reports the results of a systematic program to push the precision frontier for both $\alpha$ and $\delta$, and to understand the internal structure of the entanglement entropy in unprecedented detail. The program comprises sixteen lattice experiments (V2.192 through V2.236) that collectively establish the following: - $\alphas = 1/(24\sqrt{\pi})$ confirmed to $0.01 ((ref: sec:alpha)); - $\delta$ is universal across three lattice regulators to $0.24 ((ref: sec:delta)); - First lattice extraction of the vector $\delta_{\rm v} = -31/45$ at $0.01 - Discovery of the "50 ((ref: sec:50pct)); - First angular momentum decomposition of $\delta$, revealing positive monopole contributions ((ref: sec:angular)); - Discovery of a universal scaling function $f(x = l/n)$ that determines $\alphas$ via a spectral integral ((ref: sec:scaling)); - Proof that $99.8 ((ref: sec:boundary)); - Demonstration that the continuum (Bessel) approximation fails and that $\alphas$ is an intrinsically lattice quantity ((ref: sec:continuum)). We are careful throughout to distinguish results that are established beyond reasonable doubt (e.g., $\alpha$ to $0.01 approximate (e.g., the 50 represent negative results (e.g., the failure of the continuum limit). All numerical data, code, and parameter files are publicly available in the Moon Walk Project repository. -- Historical context -- The area law for entanglement entropy was first identified by Bombelli et al. [ref: Bombelli:1986rw] and independently by Srednicki [ref: Srednicki:1993im], who computed the entropy numerically for a free scalar field and showed that $S \propto A/\epsilon^2$. The logarithmic correction was identified by Callan and Wilczek [ref: Callan:1994py] and by Holzhey, Larsen, and Wilczek [ref: Holzhey:1994we] in the context of $1{+}1$-dimensional conformal field theory, where it is determined by the central charge. In $3{+}1$ dimensions, the logarithmic coefficient was connected to the trace anomaly by Solodukhin [ref: Solodukhin:2008dh], establishing the relationship $\delta = -4a$ where $a$ is the Euler density coefficient of the conformal anomaly. For gauge fields, the entanglement entropy acquires additional structure. Kabat [ref: Kabat:1995eq] showed that the "contact term" contribution to the entropy differs for fields of different spin. Donnelly and Wall [ref: Donnelly:2014fua,Donnelly:2015hxa] identified "edge modes" at the entangling surface--degrees of freedom associated with gauge transformations that act independently on the two sides of the cut. For gravity, Benedetti and Casini [ref: Benedetti:2019uej] computed the entanglement entropy of linearized gravitons on a sphere, finding $\delta_{\rm EE} = -61/45$, distinct from the effective action value $\delta_{\rm EA} = -212/45$. The difference between these two values reflects the edge mode contribution. On the numerical side, the angular momentum decomposition method introduced by Lohmayer et al. [ref: Lohmayer:2009sq] enabled high-precision computations of entanglement entropy on a sphere. Their work, together with subsequent refinements, established the numerical value of $\alphas$ to a few percent. The present work pushes this precision by two to three orders of magnitude and extends the numerical program to vector and graviton fields for the first time. -- Organization -- The paper is organized as follows. (Ref: sec:method) reviews the Srednicki lattice construction and the angular momentum decomposition method, including the symplectic eigenvalue technique and the proportional angular cutoff convention. (Ref: sec:alpha) presents the precision determination of the area- law coefficient $\alphas$, including Richardson extrapolation, regulator dependence, and continuum extrapolation in the lattice size $N$. (Ref: sec:delta) addresses the log coefficient $\delta$, documenting the signal hierarchy problem, the per-spin extractions, regulator independence, angular cutoff independence, and mass independence. (Ref: sec:higher_spin) reports the first lattice results for vector and graviton fields. (Ref: sec:50pct) describes the discovery and precision measurement of the 50 momentum decomposition of $\delta$ and the positive monopole contribution. (Ref: sec:scaling) describes the universal scaling function $f(x)$ and the spectral integral route to $\alphas$. (Ref: sec:boundary) documents boundary mode dominance and the reduction to a single-eigenvalue problem. (Ref: sec:continuum) proves the failure of the continuum approximation. (Ref: sec:discussion) discusses limitations, the fermion problem, connections to the cosmological constant prediction, and open problems. (Ref: sec:conclusion) summarizes the main results. === The Srednicki lattice and angular momentum decomposition === -- The radial Hamiltonian -- We follow the Srednicki construction [ref: Srednicki:1993im] as refined by Lohmayer, Neuberger, Schwimmer, and Theisen [ref: Lohmayer:2009sq]. A free massless scalar field in $3{+}1$ dimensions is decomposed into angular momentum channels $l$, each of which reduces to a one-dimensional radial chain. The Hamiltonian for each channel takes the form H = \frac{1}{2}\sum_{j=1}^{N}\left[ \pi_j^2 + \phi_j (K'_l)_{jk} \phi_k \right], where $N$ is the number of radial lattice sites and the coupling matrix $K'_l$ is tridiagonal with entries (K'_l)_{j,j} = \frac{(j-\half)^2 + (j+\half)^2 + l(l{+}1) + \Delta_s}{j^2} , (K'_l)_{j,j\pm 1} = -\frac{(2j\pm 1)}{2j(j\pm 1) + \ldots} , where $\Delta_s$ is a spin-dependent barrier shift: \Delta_s = \begin{cases} 0 & scalar (s=0), 0 & vector (s=1) in Coulomb gauge, -2 & graviton (s=2) (Lichnerowicz). \end{cases} The off-diagonal elements encode the nearest-neighbor coupling on the radial lattice and are identical for all spins. Explicitly, the off-diagonal elements take the form $(K'_l)_{j,j+1} = -(j+\half)^2/(j(j{+}1))$ for the standard 3-point stencil, arising from the discretization of the radial Laplacian $-(1/r^2) d/dr (r^2 d/dr)$ on a grid with spacing $\epsilon$. The physical interpretation of the lattice parameters is as follows. The radial coordinate is $r = j\epsilon$, where $j = 1, \ldots, N$ labels the lattice sites. The subsystem ("interior" of the sphere) consists of sites $j = 1, \ldots, n$, so the entangling surface is between sites $n$ and $n+1$, at radius $R = n\epsilon$. The lattice extends to $r_{\max} = N\epsilon$. The centrifugal barrier $l(l{+}1)/j^2$ confines high-$l$ modes away from the origin, ensuring convergence of the mode sum. For the graviton, the Lichnerowicz barrier $(l{-}1)(l{+}2)/j^2$ replaces the scalar barrier $l(l{+}1)/j^2$. This barrier is negative at $l = 0$ and zero at $l = 1$, reflecting the fact that there are no physical graviton modes at these angular momenta--the $l = 0$ and $l = 1$ sectors correspond to gauge degrees of freedom (linearized diffeomorphisms). Physically, the graviton has only $l \geq 2$ transverse-traceless modes, analogous to the vector having only $l \geq 1$ transverse modes. -- Numerical implementation -- The matrix $K'_l$ is real, symmetric, and positive definite for $l \geq 0$ (scalar/vector) or $l \geq 2$ (graviton with $\Delta_s = -2$). Its square root and inverse square root are computed via Cholesky decomposition of the full matrix, which is numerically stable and has $O(N^3)$ complexity. For the Cholesky-based symplectic eigenvalue method, we form the product $M = L^T P_{\rm sub} L$ where $L$ is the Cholesky factor of $X_{\rm sub}$. The eigenvalues of $M$ are the squares of the symplectic eigenvalues $\nu_k^2$. This method avoids the numerical instability associated with directly computing the eigenvalues of the non-symmetric product $X_{\rm sub} P_{\rm sub}$. Each angular momentum channel requires one eigendecomposition of an $N \times N$ matrix (for $K'_l$) and one eigendecomposition of an $n \times n$ matrix (for $M$). With $l_{\max} = Cn$ channels and $N_n$ values of $n$, the total computational cost scales as $O(Cn \cdot N^3 + Cn \cdot N_n \cdot n^3)$. At our typical parameters ($N = 1000$, $n_{\max} = 80$, $C = 10$), the dominant cost is the $K'_l$ decomposition, requiring approximately $800$ matrix diagonalizations of size $1000 \times 1000$. -- Entanglement entropy via symplectic eigenvalues -- The vacuum state of the Gaussian Hamiltonian (eq: eq:radialH) is fully characterized by the two-point correlation matrices X_{jk} = \langle \phi_j \phi_k \rangle, P_{jk} = \langle \pi_j \pi_k \rangle. For the ground state of a quadratic Hamiltonian with coupling matrix $K'_l$, these are $X = \half (K'_l)^{-1/2}$ and $P = \half (K'_l)^{1/2}$. The entanglement entropy of the first $n$ sites (the "interior" of the sphere of radius $n$) is computed from the reduced correlation matrices $X_{\rm sub}$ and $P_{\rm sub}$ (the upper-left $n \times n$ blocks) via the symplectic eigenvalue method [ref: Peschel:2002jhw]. The symplectic eigenvalues ${\nu_k}$ of the product $X_{\rm sub} P_{\rm sub}$ satisfy $\nu_k \geq \half$, and the entanglement entropy is S_l(n) = \sum_{k=1}^{n} h(\nu_k) , where h(\nu) = \left(\nu + \half\right)\ln\left(\nu + \half\right) - \left(\nu - \half\right)\ln\left(\nu - \half\right). The total entanglement entropy is obtained by summing over angular momentum channels with the appropriate degeneracy and polarization count: S_{\rm total}(n) = \sum_{l=l_{\min}}^{l_{\max}} pol(l) (2l+1) S_l(n) , where $pol(l) = 1$ for scalars, $pol(l) = 2$ for vectors ($l \geq 1$), and $pol(l) = 2$ for gravitons ($l \geq 2$). -- The proportional angular cutoff -- A key technical choice is the angular momentum cutoff. We use the proportional cutoff $l_{\max} = C \cdot n$, where $C$ is a fixed integer (typically $C = 8$-$20$). Under this convention, the number of angular momentum channels scales linearly with the subsystem radius, ensuring that the UV angular resolution matches the UV radial resolution. This convention is essential for the third-difference method of $\delta$ extraction (see (ref: sec:d3S)). The alternative global cutoff $l_{\max} = C \cdot n_{\max}$ fixes $l_{\max}$ at the largest subsystem size; under finite differences $\Delta^k S(n)$, the set of included channels changes discontinuously at each $n$, introducing spurious $n$-dependent artifacts in the analysis. The choice of $C$ is a balance between precision and computational cost. At $C = 2$, only $2n$ angular momentum channels are included, and roughly $10 of the area-law coefficient is missed (the tail of $f(x)$ extends to $x \sim 10$). At $C = 10$, more than $90 the $\delta$ extraction is unaffected ($\delta$ is $C$-independent to $0.003 as shown in (ref: sec:delta)). At $C = 20$, the area law is captured to better than $99 $C = 8$-$10$. -- The second-difference method for $\alpha$ -- The second finite difference of $S(n)$ is \Delta^2 S(n) = S(n{+}1) - 2 S(n) + S(n{-}1) = 8\pi \alpha(C) + \delta \ln\left(1 - \frac{1}{n^2}\right) + \frac{2\beta}{n(n^2{-}1)} + \order{n^{-4}} . At large $n$, the leading term $8\pi\alpha(C)$ dominates, and $\alpha(C)$ can be read off directly. A key technical point is the use of the exact second- difference basis functions $\ln(1 - 1/n^2)$ and $2/[n(n^2{-}1)]$ rather than their Taylor approximations $-1/n^2$ and $2/n^3$. The exact forms avoid near- degeneracy between the $\delta$ and $\beta$ terms, which improves the conditioning of the simultaneous fit. The $C$-dependent area coefficient $\alpha(C)$ converges to the continuum value $\alphas$ as $C \to \infty$. We accelerate this convergence using Richardson extrapolation: given $\alpha(C)$ at $k+1$ different $C$ values, the Richardson extrapolant of order $k$ eliminates the leading $k$ terms of the asymptotic expansion $\alpha(C) = \alphas + \sum_{j=1}^{\infty} a_j/C^{2j}$. -- The d3S extraction method for $\delta$ -- The third finite difference of $S(n)$ at fixed proportional cutoff is defined as \Delta^3 S(n) = S(n{+}2) - 3 S(n{+}1) + 3 S(n) - S(n{-}1) . Substituting the expansion (eq: eq:Sexpansion) with $A = 4\pi n^2$, one finds \Delta^3 S(n) = \frac{2\delta}{n^3} + \frac{B}{n^4} + \order{n^{-5}} , where the area-law term ($\propto n^2$) and the Euler-Maclaurin correction ($\propto 1/n$) are both exactly cancelled. Fitting $\Delta^3 S(n)$ to a two- parameter model $A/n^3 + B/n^4$ over a range of $n$ values extracts $\delta = A/2$. This cancellation is the reason the d3S method succeeds: it isolates the UV- finite logarithmic term by eliminating the dominant UV-sensitive area contribution and the subleading $1/n$ correction. The second finite difference $\Delta^2 S$, by contrast, retains a constant term $8\pi\alpha(C)$ that must be separately subtracted, introducing sensitivity to the cutoff-dependent $\alpha$. More specifically, consider the entropy expansion on the lattice with proportional cutoff: S(n) = \alpha(C) 4\pi n^2 + \delta \ln n + \gamma(C) + \frac{\beta}{n} + \frac{\sigma}{n^2} + \order{n^{-3}} . The first difference $\Delta S = S(n{+}1) - S(n)$ retains the $O(n)$ area term. The second difference $\Delta^2 S$ cancels $O(n)$ and $O(1)$ terms, leaving $\alpha$ as a constant plus $O(n^{-2})$ corrections. The third difference goes further, cancelling the constant $8\pi\alpha$ and the $\beta/n$ correction simultaneously. The price is that the signal is correspondingly smaller: $\Delta^3 S \sim 2\delta/n^3$ is $\sim10^4$ times smaller than $\Delta^2 S \sim 8\pi\alpha$. In practice, we fit $n^4 \cdot \Delta^3 S(n)$ to a two-parameter linear model $A \cdot n + B$, which is equivalent to fitting $\Delta^3 S = A/n^3 + B/n^4$. The slope $A/2$ gives $\delta$. The fit is performed over a range $n_{\min} \leq n \leq n_{\max}$, with $n_{\min}$ chosen large enough that higher-order corrections ($O(n^{-5})$) are negligible. Typical values are $n_{\min} = 25$-$30$ and $n_{\max} = 70$-$100$. -- Parameter summary -- (Ref: tab:params) summarizes the computational parameters used across the experiments in this paper. [Table: Computational parameters across the main experiments.] lccccc@{}} Experiment | $N$ | $C$ | $n_{\min}$-$n_{\max}$ | Fields | Key result V2.192 | 500 | 20 | 10-30 | S | $\alpha$ to $0.01 V2.195 | 1000 | 10 | 20-100 | S, V, G | First V, G $\delta$ V2.196 | 600 | 8 | 20-80 | S | $\delta$ universality V2.198 | 600-1800 | 8 | 20-70 | S, V, G | Richardson $\delta$ V2.199 | 1200 | 2-22 | 20-50 | S, V, G | $C$-independence V2.218 | 500-800 | 10 | 25-100 | S, V, G | Graviton $\delta$ V2.219 | 800 | 5-20 | 30-100 | S, V, G | Precision 50 V2.220 | 200-1000 | 10 | 30-80 | S, V, G | $N$-extrapolation V2.225 | 600 | 10 | 20-60 | S, V, G | $l$-decomposition V2.226 | 1000 | 10 | 30-80 | S, V, G | Precision monopole V2.230 | 600 | 10 | 8-30 | S | Scaling function V2.231 | 600 | 10 | 8-30 | S | Analytic $f(x)$ V2.233 | 600 | 6 | 10-30 | S | Radial profile V2.234 | 600 | 10 | 8-30 | S | Boundary mode V2.236 | 600 | 10 | 10-30 | S | Continuum limit Here S = scalar, V = vector, G = graviton. === Area-law coefficient: the precision frontier === -- The conjecture and its significance -- The area-law coefficient for a free massless scalar in $3{+}1$ dimensions is conjectured to take the exact value \alphas = \frac{1}{24\sqrt{\pi}} = \frac{1}{4! \Gamma(3/2)} = 0.023508\ldots This value was first computed numerically by Srednicki [ref: Srednicki:1993im] and has been refined in subsequent lattice studies. No analytic derivation from first principles exists; the value arises from the lattice structure at the entangling surface (see (ref: sec:continuum)). The factorization $1/(24\sqrt{\pi}) = 1/(4! \Gamma(3/2))$ is suggestive. The factor $24 = 4!$ might be related to the spacetime dimension $D = 4$, while $\sqrt{\pi} = \Gamma(3/2)$ could arise from a hemisphere integral on $S^2$. However, no derivation connecting these factors to the entanglement computation has been found. The physical significance of $\alphas$ is that it determines Newton's gravitational constant in the entanglement framework: $G = 1/(4\alpha_{\rm total})$, where $\alpha_{\rm total}$ sums $\alpha$ over all field species weighted by their area-law degeneracies. Any error in $\alphas$ propagates directly to the cosmological constant prediction via $\Omega_\Lambda = \abs{\delta_{\rm total}}/(6\alpha_{\rm total})$. -- Richardson extrapolation in the angular cutoff -- We extract $\alpha$ from the second finite difference \Delta^2 S(n) = S(n{+}1) - 2 S(n) + S(n{-}1) = 8\pi \alpha(C) + \delta \ln\left(1 - 1/n^2\right) + \frac{2\beta}{n(n^2{-}1)} + \order{n^{-4}} , at each angular cutoff $C$. The function $\alpha(C)$ converges to $\alphas$ as $C \to \infty$. We accelerate the convergence using Richardson extrapolation of order $k$ from a sequence of $C$ values. (Ref: tab:alpha_richardson) shows the Richardson sequence from a sweep $C = 10,\ldots,50$ with $N = 1000$ radial sites: [Table: Richardson extrapolation of the area-law coefficient $\alpha$ from a sweep of angular cutoff values $C = 10$-$50$, at $N = 1000$ radial sites and $n = 10$-$30$.] ccccc@{}} Order $k$ | $C_{\min}$ | $C_{\max}$ | $\alpha$ | Error ( 2 | 40 | 50 | 0.023551 | $+0.185$ 3 | 30 | 50 | 0.023516 | $+0.035$ 4 | 25 | 50 | 0.023512 | $+0.017$ 7 | 10 | 50 | 0.023510 | $+0.011$ The sequence converges monotonically from above. At order $k = 7$, the result is $\alpha = 0.023510$, deviating from $1/(24\sqrt{\pi}) = 0.023508$ by $+0.011 and independently reproduces the result of earlier experiments in this program. -- Regulator dependence of $\alpha$ -- To verify that $\alpha$ is genuinely UV-sensitive, we compute it using three different lattice discretizations of the radial Laplacian $-d^2/dr^2$: [Table: Area-law coefficient $\alpha$ for three lattice stencils. The standard 3-point stencil is $O(h^2)$; the 5-point and 7-point stencils are $O(h^4)$ and $O(h^6)$ respectively. Parameters: $N = 600$, $C = 8$, $n = 20$-$80$.] lccc@{}} Stencil | Order | $\alpha_{\rm scalar}$ | Ratio to 3-point 3-point | $O(h^2)$ | 0.2911 | 1.00 5-point | $O(h^4)$ | 0.3824 | 1.31 7-point | $O(h^6)$ | 0.4280 | 1.47 The area coefficient changes by $31 to the standard stencil, confirming that $\alpha$ is regulator-dependent as expected for a UV-divergent quantity. However, the ratio $\alpha_{\rm v}/\alphas = 2.000$ is preserved across all three stencils, demonstrating that relative weightings between fields are regulator-independent. This dichotomy--$\alpha$ is regulator-dependent but $\alpha$ ratios are regulator-independent--has an important physical consequence. The cosmological constant prediction involves the ratio $\abs{\delta_{\rm total}}/\alpha_{\rm total}$. While $\alpha_{\rm total}$ depends on the regulator, the relative weights of different fields in the sum $\alpha_{\rm total} = \sum_i w_i \alphas$ do not. Since $\delta_{\rm total}$ is also regulator-independent ((ref: sec:delta_universal)), the prediction $\Omega_\Lambda$ depends only on the field content of the theory and the single non-universal number $\alphas$, which we verify to $0.01 3-point stencil. The regulator dependence of $\alpha$ also explains why the absolute value $\alphas = 0.0235$ (on the 3-point stencil) cannot be directly compared with continuum calculations. The number $1/(24\sqrt{\pi})$ is a property of the specific lattice discretization, not a universal continuum quantity. We return to this point in (ref: sec:continuum). -- Continuum extrapolation in $N$ -- The lattice size $N$ (number of radial sites) introduces a finite-size correction to $\delta$ (and, at higher order, to $\alpha$). We study this correction by computing the scalar $\delta$ at seven lattice sizes from $N = 200$ to $N = 1000$, at fixed $C = 10$ and $n = 30$-$80$. [Table: Scalar $\delta$ as a function of lattice size $N$. All runs use $C = 10$, $n_{\min] rrc@{}} $N$ | $\delta$ (d3S 2-param) | Error vs $-1/90$ 200 | $-0.24483$ | $2103 300 | $-0.04142$ | $273 400 | $-0.01933$ | $74 500 | $-0.01416$ | $27 600 | $-0.01242$ | $12 800 | $-0.01137$ | $2.4 1000 | $-0.01111$ | $0.02 The convergence is extremely steep: the result at $N = 200$ (where the lattice boundary contaminates the subsystem at $n_{\max} = 80$) is catastrophically wrong, but by $N = 1000$ the error is $0.02 $\delta(N) = a + b/N^p$ yields $p \approx 5.0$ for the scalar, confirming that the finite-size correction decays much faster than the naive $1/N^2$ lattice discretization error. The practical implication is that the raw result at sufficiently large $N$ is more reliable than Richardson extrapolation from smaller $N$ values: the correction is too steep for polynomial extrapolation to improve upon the data. The rule of thumb from these data is $N/n_{\max} > 10$ for $< 1 $N/n_{\max} > 12$ for $< 0.1 === The log coefficient: the signal hierarchy problem === -- The trace anomaly connection -- The logarithmic coefficient $\delta$ in the entanglement entropy expansion is related to the trace anomaly coefficients of the field theory. In $3{+}1$ dimensions, the trace anomaly of a conformally invariant field takes the form \langle T^\mu{}_\mu \rangle = a E_4 + c W^2 + (total derivatives) , where $E_4$ is the Euler density and $W^2 = W_{\mu\nu\rho\sigma}W^{\mu\nu\rho\sigma}$ is the square of the Weyl tensor. For the entanglement entropy across a sphere, only the $a$-anomaly contributes to $\delta$: \delta = -4a . The Euler anomaly coefficient $a$ is determined by the representation content of the field and takes the exact values: a_{\rm scalar} = \frac{1}{360} , a_{\rm vector} = \frac{31}{180} , a_{\rm graviton} = \frac{61}{180} , where the graviton value refers to the entanglement entropy calculation of Benedetti and Casini [ref: Benedetti:2019uej], not the effective action value $a_{\rm EA} = 212/180$. The corresponding log coefficients are \delta_{\rm s} = -\frac{1}{90} , \delta_{\rm v} = -\frac{31}{45} , \delta_{\rm g} = -\frac{61}{45} . These values are protected by the Wess-Zumino consistency condition: the trace anomaly is determined by representation theory and cannot be modified by interactions or renormalization. In this sense, $\delta$ is "exact" and does not require lattice verification. However, the lattice computation provides an independent check that the entanglement entropy logarithmic term is indeed controlled by the trace anomaly, as the theoretical arguments predict. -- The fundamental difficulty: the signal hierarchy -- While the area-law coefficient $\alpha$ can be verified to $0.01 modest lattices, the log coefficient $\delta$ presents a qualitatively harder problem. The root cause is the signal hierarchy: in the second finite difference, the area-law contribution dominates by a factor of $\sim10^4$ over the logarithmic contribution. At $n = 20$, $C = 20$, the second difference $\Delta^2 S \approx 0.585$ is dominated by $8\pi\alpha$, while the $\delta$ contribution is $\delta \cdot \ln(1 - 1/400) \approx 2.8 \times 10^{-6}$--only 5 parts per million. Extracting $\delta$ requires measuring $\Delta^2 S$ to $\sim6$ significant figures and then separating the $1/n^2$ and $1/n^3$ corrections, which are nearly degenerate in the accessible range $n = 10$-$30$. The third difference $\Delta^3 S$ cancels the leading area term, but the resulting signal $\sim 2\delta/n^3$ is still very small. At $n = 20$, this is $\sim 2.8 \times 10^{-6}$, requiring high-precision arithmetic throughout the computation. -- Scalar $\delta$: from $10 The scalar log coefficient $\delta_{\rm s} = -1/90 = -0.01111\ldots$ has been extracted with progressively improving precision across this program: [Table: Scalar $\delta$ extraction across experiments, ordered by precision.] llcc@{}} Experiment | Parameters | $\delta$ | Error V2.192 ($N = 500$) | d2S fit, $n = 10$-$30$ | $-0.0127$ | $+14 V2.218 ($N = 500$) | d3S, $C = 10$ | $-0.01277$ | $+14.9 V2.219 ($N = 800$) | d3S, $C = 10$ | $-0.01153$ | $+3.8 V2.220 ($N = 1000$) | d3S, $C = 10$, $n = 30$-$80$ | $-0.01111$ | $+0.02 The $0.02 The key parameter is the ratio $N/n_{\max}$: at $N/n_{\max} = 1000/80 = 12.5$, finite-size effects are negligible. -- Sensitivity to fit range -- The extracted $\delta$ depends strongly on the range of $n$ values used in the fit, a hallmark of ill-conditioned extraction: [Table: Dependence of extracted $\delta$ on the fitting range (V2.192, $N = 500$, $C = 20$).] cccc@{}} $n$ range | $n_{\rm pts}$ | $\alpha$ error | $\delta$ error $[10, 20]$ | 5 | $-0.93 $[10, 25]$ | 6 | $-0.93 $[10, 30]$ | 7 | $-0.93 $[12, 30]$ | 6 | $-0.93 $[15, 35]$ | 6 | $-0.93 While $\alpha$ is stable to six significant figures regardless of the fit range, $\delta$ varies by 40 percentage points. This instability disappears at larger $N$ (where the ratio $N/n_{\max}$ is sufficiently large), confirming it is a finite-size effect rather than a methodological failure. -- Vector $\delta$: the 0.01 The vector log coefficient $\delta_{\rm v} = -31/45 = -0.68889\ldots$ is $62\times$ larger in magnitude than the scalar value, making it far more accessible. Richardson extrapolation across lattice sizes $N = 600$-$1800$ yields: [Table: Richardson extrapolation of the vector $\delta$ across lattice sizes (V2.198). Parameters: $C = 8$, $n = 20$-$70$, fit from $n_{\min] rrcc@{}} $N$ | $\delta_{\rm v}$ (raw) | Richardson (order 3) | Error 600 | $-0.697285$ | -- | $1.2 800 | $-0.692024$ | -- | $0.5 1000 | $-0.690707$ | -- | $0.3 1300 | $-0.690102$ | -- | $0.2 1800 | $-0.689939$ | $-0.688793$ | $0.01 The order-3 Richardson extrapolation gives $\delta_{\rm v} = -0.688793$, matching the exact value $-31/45 = -0.688889$ to $0.01 precise lattice verification of a trace anomaly coefficient for a gauge field, and it demonstrates that the d3S method combined with Richardson extrapolation can achieve sub-percent precision when the signal is sufficiently large. -- Regulator independence of $\delta$ -- The log coefficient must be independent of the UV regulator if it is to serve as a physical input to the cosmological constant prediction. We test this directly by computing $\delta$ using three different lattice stencils. The cleanest test uses the $l = 0$ channel, which has no area-law contribution and therefore avoids the extraction difficulties associated with cancelling the dominant area term: [Table: Per-channel $l = 0$ log coefficient across lattice stencils (V2.196). Parameters: $N = 600$, $C = 8$, $n = 20$-$80$.] lccc@{}} Stencil | Order | $\delta_{l=0}$ | Error vs $+1/3$ 3-point | $O(h^2)$ | $+0.3398$ | $1.95 5-point | $O(h^4)$ | $+0.3406$ | $2.19 7-point | $O(h^6)$ | $+0.3402$ | $2.07 The spread across stencils is $0.24 that $\delta$ is a property of the continuum field theory, not the lattice regulator. The ratios $\delta_{5\rm pt}/\delta_{3\rm pt} = 1.0024$ and $\delta_{7\rm pt}/\delta_{3\rm pt} = 1.0012$ are consistent with unity. The full scalar and vector $\delta$ extractions fail for extended stencils ($R^2 \to 0$ for 5-point and 7-point), because the non-nearest-neighbor couplings introduce boundary artifacts not captured by the d3S ansatz. This is a failure of the extraction method, not of universality: the per-channel test ((ref: tab:delta_universality)) cleanly isolates the universal content. -- $C$-cutoff independence of $\delta$ -- If $\delta$ is truly UV-finite, it must be independent of the angular momentum cutoff parameter $C$ in $l_{\max} = C \cdot n$. We verify this by computing $\delta$ at $C = 2, 4, 6, 8, 10, 14, 18, 22$ for all three spins: [Table: $\delta$ as a function of angular cutoff $C$ (V2.199 and V2.219). The variation is characterized by the coefficient of variation (CV).] lccc@{}} Field | $\delta$ range | CV | $\delta$ at $C = 8$ Scalar | $[-0.02065, -0.02061]$ | $0.11 Vector | $[-0.6648, -0.6645]$ | $0.01 Graviton | $[-0.6883, -0.6882]$ | $0.003 The graviton $\delta$ is the most $C$-independent quantity measured, varying by only $5.6 \times 10^{-5}$ across $C = 5$ to $C = 20$. This is a direct numerical proof that $\delta$ is UV-finite: it does not depend on how many angular momentum modes are included, because the high-$l$ modes at the cutoff boundary contribute negligible entropy. By contrast, $\alpha$ grows logarithmically with $C$ (from $0.0212$ at $C = 5$ to $0.0233$ at $C = 20$ for the scalar), confirming its UV-divergent nature. The clean separation between $C$-dependent $\alpha$ and $C$-independent $\delta$ validates the theoretical distinction between these two quantities. -- Mass independence -- We verify that $\delta$ is unchanged when a mass $m$ is added to the field, provided $m \cdot n_{\max} \ll 1$: [Table: $l = 0$ log coefficient as a function of field mass (V2.196).] cccc@{}} $m^2$ (lattice) | $m \cdot n_{\max}$ | $\delta_{l=0}$ | $\delta/\delta(m{=}0)$ $0$ | 0 | $+0.3337$ | 1.000 $10^{-8}$ | 0.01 | $+0.3337$ | 1.000 $10^{-6}$ | 0.10 | $+0.3377$ | 1.012 $10^{-5}$ | 0.31 | $+0.3670$ | 1.100 $10^{-4}$ | 0.99 | $+0.4253$ | 1.275 $10^{-3}$ | 3.13 | $-0.0274$ | $-0.082$ For $m \cdot n_{\max} < 0.1$, $\delta$ is unchanged to $< 0.02 decoupling threshold occurs at $m \cdot n_{\max} \sim 1$, where the correlation length becomes comparable to the subsystem size and the logarithmic structure is lost. For all Standard Model particles, $m/M_{\rm Pl} < 10^{-17}$, placing them deep in the mass-independent regime. === Higher-spin log coefficients === -- Angular momentum decomposition for arbitrary spin -- All massless bosonic fields in flat space share the same radial equation, with the only differences being the angular momentum range and polarization count: [Table: Field content for the three bosonic spins.] lccc@{}} Field | $l_{\min}$ | Polarizations | Barrier shift $\Delta_s$ Scalar ($s = 0$) | 0 | 1 | 0 Vector ($s = 1$) | 1 | 2 | 0 Graviton ($s = 2$) | 2 | 2 | $-2$ For the vector, we use Coulomb gauge, where the physical transverse modes have the same radial equation as the scalar but with $l \geq 1$ (the monopole $l=0$ is pure gauge). For the graviton, the transverse-traceless (TT) modes have an effective centrifugal barrier $l(l{+}1) - 2 = (l{-}1)(l{+}2)$, which is positive only for $l \geq 2$. -- Area-law coefficient ratios -- The heat kernel prediction is that each physical polarization contributes equally to the area-law coefficient. Since the vector has 2 polarizations and the scalar has 1, the ratio should be $\alpha_{\rm v}/\alphas = 2$. Similarly, $\alpha_{\rm g}/\alphas = 2$. Our lattice results confirm this exactly: [Table: Area-law coefficient ratios (V2.195, V2.218, V2.219).] lcc@{}} Ratio | Lattice | Theory $\alpha_{\rm v}/\alphas$ | 2.0000 | 2 $\alpha_{\rm g}/\alphas$ | 2.0000 | 2 These ratios are confirmed to four decimal places at every value of $C$ tested ($C = 5$-$20$) and are preserved across all three lattice stencils. This provides strong numerical evidence for the heat kernel prediction. -- Log coefficient: per-spin results -- The first lattice extractions of vector and graviton log coefficients are summarized in (ref: tab:delta_spins). [Table: Log coefficient $\delta$ for all three bosonic spins. The "TT modes" column gives the lattice result from the physical propagating modes only; the "Full analytical" column gives the trace anomaly value including edge/ghost contributions.] lcccc@{}} Field | $\delta_{\rm TT}$ (lattice) | $\delta_{\rm full}$ (analytical) | Ratio | Reference Scalar | $-0.01111$ | $-1/90 = -0.01111$ | 1.000 | V2.220 Vector | $-0.3558$ | $-31/45 = -0.6889$ | 0.516 | V2.219 Graviton | $-0.6883$ | $-61/45 = -1.3556$ | 0.508 | V2.219 The scalar result agrees with theory to $0.02 the lattice captures only the TT (physical propagating) modes; the full analytical values include additional contributions from gauge constraints, ghosts, and edge modes at the entangling surface. The observation that the lattice TT values are approximately half the full analytical values is the "50 -- Why the vector succeeds and the scalar struggles -- The precision of $\delta$ extraction is fundamentally limited by the ratio of the log signal to the area signal. For the three fields: [Table: Signal hierarchy for the three bosonic spins. The "signal ratio" is $\abs{\delta] lccc@{}} Field | $\abs{\delta}$ | $\alpha$ (per DOF) | Signal ratio Scalar | 0.011 | 0.024 | 0.47 Vector | 0.689 | 0.047 | 14.7 Graviton | 1.356 | 0.047 | 28.8 The vector log signal is $31\times$ stronger than the scalar (per degree of freedom), and the graviton is $61\times$ stronger. This explains why the vector $\delta$ can be extracted to $0.01 the scalar requires $N = 1000$ to reach $0.02 The signal hierarchy also explains why the d3S method is more effective for gauge fields: the logarithmic term constitutes a larger fraction of the total entropy, making it easier to isolate by finite differencing. -- The Lichnerowicz barrier and its effects -- For the graviton, the effective centrifugal barrier is the Lichnerowicz eigenvalue $(l{-}1)(l{+}2)$ rather than the scalar value $l(l{+}1)$. The difference is (l{-}1)(l{+}2) - l(l{+}1) = -2 , which is a constant shift independent of $l$. For $l \geq 2$, this is a small perturbation: the barrier is $6, 12, 20, \ldots$ for the graviton versus $6, 12, 20, \ldots$ for the scalar at $l = 2, 3, 4, \ldots$ (both starting at $6$ when $l = 2$, but with a shift that becomes relatively smaller at large $l$). At $l = 0$, the Lichnerowicz barrier is $(0{-}1)(0{+}2) = -2$, which is negative. This means the $l = 0$ channel has an attractive potential near the origin, leading to unphysical bound states and nonsensical entropy values. At $l = 1$, the barrier is $(1{-}1)(1{+}2) = 0$, which gives a free particle with no confinement. Only at $l \geq 2$ does the barrier become positive and the physics sensible. This is the lattice manifestation of the gauge constraint: $l = 0$ and $l = 1$ graviton modes are pure gauge and should not be included in the physical sum. -- The graviton error is not finite-size -- A striking finding is that the graviton $\delta_{\rm TT}$ is essentially independent of the lattice size $N$. While the scalar $\delta$ improves dramatically from $N = 200$ to $N = 1000$ (see (ref: tab:delta_N)), the graviton value is flat: [Table: Graviton $\delta_{\rm TT] rc@{}} $N$ | $\delta_{\rm TT}$ (graviton) 600 | $-0.6886$ 800 | $-0.6881$ 1000 | $-0.6880$ 1300 | $-0.6880$ 1800 | $-0.6879$ The shift from $N = 600$ to $N = 1800$ is $< 0.1 convergence to $-0.688$ (half the Benedetti-Casini value) is not a finite-size artifact: it is the correct answer for the physical TT modes. The "missing" half resides in edge modes that are not captured by the lattice computation of propagating degrees of freedom. === The 50 -- Discovery -- The most unexpected finding of this program is that for both gauge fields (vector and graviton), the lattice TT-mode $\delta$ is approximately half the full analytical value: [Table: The 50 $\delta$ (V2.218, V2.219, V2.220).] lccc@{}} Field | $\delta_{\rm TT}/\delta_{\rm EE}$ | Deviation from $1/2$ | Edge modes? Scalar ($s = 0$) | 1.000 | N/A | No Vector ($s = 1$) | 0.5153 | $+3.1 Graviton ($s = 2$) | 0.5077 | $+1.5 For the scalar, which has no gauge symmetry and therefore no edge modes, the ratio is $1.000$ (at $N = 1000$). For the vector and graviton, the TT modes carry approximately $50 edge modes associated with gauge constraints at the entangling surface. -- Precision measurement -- The precision of the 50 multiple angular cutoff values $C$: [Table: Graviton $\delta_{\rm TT] ccc@{}} $C$ | $\delta_{\rm TT}$ | Ratio to $-61/45$ 5 | $-0.68826$ | 0.5077 8 | $-0.68826$ | 0.5077 10 | $-0.68825$ | 0.5077 15 | $-0.68823$ | 0.5077 20 | $-0.68821$ | 0.5077 The ratio is $0.5077$ at every $C$ value tested, with CV $= 0.003 extraordinary stability confirms that the 50 property, not an artifact of the cutoff prescription. -- The 50 At $N = 1000$--where the scalar is verified to $0.02 is $0.5076$ and the vector ratio is $0.5153$. These values are stable from $N = 600$ to $N = 1000$: [Table: $N$-dependence of the 50 ccc@{}} $N$ | Vector ratio | Graviton ratio 600 | 0.514 | 0.507 800 | 0.515 | 0.508 1000 | 0.515 | 0.508 If the deviations from $1/2$ were finite-$N$ effects, they would shrink as $N$ increases. They do not. The vector ratio is genuinely $\sim0.515$ and the graviton is genuinely $\sim0.508$. The 50 $1$-$3 Crucially, the measured ratio $0.508$ for the graviton is more accurate than an assumed exact $1/2$ for recovering the Benedetti-Casini analytical value. Using the measured ratio: \delta_{\rm grav}^{\rm (full)} = \frac{\delta_{\rm TT}}{0.508} = \frac{-0.688}{0.508} = -1.354 , which agrees with $-61/45 = -1.356$ to $0.1 exact $1/2$ gives $2 \times (-0.688) = -1.376$, deviating by $1.5 -- Physical interpretation: edge modes -- The 50 entanglement [ref: Donnelly:2014fua,Donnelly:2015hxa], which predicts that gauge transformations at the entangling surface generate "edge modes" carrying an $O(1)$ fraction of the entanglement entropy. Our measurement quantifies this fraction for the first time on the lattice: - Vector: Edge modes carry $48.5 - Graviton: Edge modes carry $49.2 The slight asymmetry (physical modes carry marginally more than edge modes) may reflect the fact that TT modes have direct entanglement across the surface, while edge modes encode gauge constraints that are slightly less entangled. The fact that both spin-1 and spin-2 gauge fields show approximately the same $\sim50 diffeomorphisms), suggests a deep universality in the relationship between gauge symmetry and entanglement entropy. -- The edge mode decomposition for gravity -- For gravity, the distinction between the entanglement entropy value $\delta_{\rm EE} = -61/45$ and the effective action value $\delta_{\rm EA} = -212/45$ has been a subject of active debate. The effective action computation includes contributions from all metric fluctuations, including pure gauge modes. The entanglement entropy computation, by contrast, must carefully account for the gauge constraints and the associated edge modes [ref: Benedetti:2019uej,Blommaert:2025pcf]. Our lattice computation captures the TT modes--the physical propagating graviton degrees of freedom. The edge modes, which arise from linearized diffeomorphisms at the entangling surface, are not represented in our simplified model. The "missing" $\delta$ is therefore \delta_{\rm edge} = \delta_{\rm EE} - \delta_{\rm TT} = -1.356 - (-0.688) = -0.668 , which is $49.2 The physical picture is that the graviton entanglement entropy decomposes into two approximately equal parts: (i) the entanglement of propagating gravitational waves across the surface, and (ii) the entanglement of gauge (diffeomorphism) degrees of freedom at the surface. The latter are the gravitational analogs of the electric field edge modes in Maxwell theory. -- Why vector and graviton ratios differ -- The vector ratio ($0.515$) deviates more from $1/2$ than the graviton ($0.508$). This can be understood from the angular momentum structure. The vector TT sector sums $l \geq 1$ modes (excluding the monopole $l = 0$), while the graviton TT sector sums $l \geq 2$ modes (excluding both monopole and dipole). As shown in (ref: sec:angular), the monopole ($l = 0$) has the largest individual contribution to $\delta$ and also the largest deviation from the linear trend in the $l$-mode spectrum. The vector sector includes the dipole but not the monopole, while the graviton sector excludes both. Since the graviton's $l \geq 2$ sector avoids the most irregular low-$l$ modes, its ratio is closer to the asymptotic value. -- Implications for the cosmological constant -- The 50 in the companion paper [ref: MoonWalk:Lambda]. There, the graviton contribution to $\Omega_\Lambda$ depends on the "graviton screening fraction" $f_g = \delta_{\rm EE}/\delta_{\rm EA} = (61/45)/(212/45) = 61/212 = 0.288$. This fraction determines what portion of the graviton's trace anomaly contributes to the Jacobson-Clausius relation at horizons. Our lattice measurement provides a complementary ratio: the TT fraction $\delta_{\rm TT}/\delta_{\rm EE} = 0.508$. Together with the analytical values $\delta_{\rm EE}$ and $\delta_{\rm EA}$, this implies a three-way decomposition of the graviton anomaly: \delta_{\rm EA} = \delta_{\rm TT} + \delta_{\rm edge}^{(\rm EE)} + \delta_{\rm gauge}^{(\rm EA-EE)} , where $\delta_{\rm TT} \approx \half \delta_{\rm EE}$ are the physical modes, $\delta_{\rm edge}^{(\rm EE)} \approx \half \delta_{\rm EE}$ are the edge modes counted in the entanglement entropy, and $\delta_{\rm gauge}^{(\rm EA-EE)} = \delta_{\rm EA} - \delta_{\rm EE}$ are additional gauge contributions that appear only in the effective action. Only $\delta_{\rm EE}$ enters the Jacobson thermodynamic derivation, giving $f_g = 61/212$. === Angular momentum decomposition of $\delta$ === -- Per-$l$ decomposition -- We decompose the total log coefficient into individual angular momentum contributions by computing $\delta_{\geq L}$ (the $\delta$ from summing channels $l \geq L$) and then differencing: \delta_l = \delta_{\geq l} - \delta_{\geq l+1} . [Table: Individual angular momentum contributions to $\delta$ (V2.225, V2.226, $N = 1000$, barrier $= 0$).] cccc@{}} $l$ | $\delta_l$ | $(2l{+}1) \delta_l$ | Cumulative 0 | $+0.166$ | $+0.166$ | $+0.166$ 1 | $+0.502$ | $+1.505$ | $+1.671$ 2 | $+0.844$ | $+4.221$ | $+5.892$ 3 | $+1.195$ | $+8.365$ | $+14.258$ 4 | $+1.555$ | $+13.995$ | $+28.253$ 5 | $+1.924$ | $+21.162$ | $+49.415$ -- The positive monopole -- The most surprising finding is that the monopole ($l = 0$) contributes positive $\delta$. In fact, all individual $l$-mode contributions through $l = 5$ are positive and growing. The total scalar $\delta = -1/90 = -0.011$ is negative only because the high-$l$ modes ($l \gg 5$) collectively contribute a large negative amount that overwhelms the positive low-$l$ contributions. The magnitude of the cancellation is remarkable: the cumulative weighted contribution reaches $+49.4$ by $l = 5$, yet the total (including all $l$ up to $l_{\max} \sim 800$) is $-0.011$. The ratio of the cancelling quantities to the final result exceeds $4000:1$. This cancellation is reminiscent of the cosmological constant problem itself: a small number emerging from the near-exact cancellation of individually large contributions. The analogy is not superficial: in the entanglement entropy framework, the cosmological constant is determined by $\delta$, which is itself a delicate cancellation. The low-$l$ modes (long-wavelength angular structure) "push" the entropy upward, while the high-$l$ modes (short-wavelength angular structure) "pull" it down. The trace anomaly $\delta = -1/90$ encodes the precise balance. -- Physical interpretation of positive low-$l$ contributions -- The positive sign of $\delta_{l=0}$, $\delta_{l=1}$, etc., has a natural physical interpretation. The low-$l$ modes have wavelengths comparable to or larger than the entangling sphere. As the sphere grows (increasing $n$), these modes gain additional entanglement faster than the area law because their correlations extend over the full sphere. This produces a positive logarithmic correction (entropy grows faster than $\alpha \cdot 4\pi n^2$). By contrast, high-$l$ modes are localized near the boundary and their entanglement scales strictly with the area. The negative contributions from these modes reflect the fact that the per-mode entropy decreases as the centrifugal barrier grows, and the $(2l{+}1)$ degeneracy factor, while growing, does not compensate fast enough. The net effect is a negative contribution to $\delta$ from the high-$l$ sector. The crossover between positive and negative per-$l$ contributions occurs deep in the UV ($l \gg 5$), where our explicit per-mode decomposition becomes impractical. However, the cumulative structure is clear: the total negative $\delta$ is established by the dominance of the high-$l$ tail over the individually positive low-$l$ modes. -- The quadratic spectrum -- The per-$l$ contributions follow a nearly perfect quadratic law: \delta_l = 0.167 + 0.330 l + 0.004 l^2 , with $R^2 = 0.9999997$. The dominant term is linear in $l$ with slope $0.330 \approx 1/3$, intriguingly close to the central charge $c/3 = 1/3$ of a one- dimensional free boson. The quadratic correction is small ($0.004 l^2$), meaning the spectrum is essentially linear for $l \lesssim 80$. -- Barrier comparison -- For the graviton modes ($\Delta_s = -2$), the Lichnerowicz barrier $l(l{+}1) - 2 = (l{-}1)(l{+}2)$ is positive only for $l \geq 2$. At $l = 0$ and $l = 1$ with $\Delta_s = -2$, the barrier is negative or zero, producing unphysical results. For $l \geq 2$, the barrier-$0$ and barrier-$(-2)$ spectra converge as $l$ increases: [Table: Comparison of $\delta_l$ for barrier $= 0$ and barrier $= -2$ (V2.226, $N = 1000$).] cccc@{}} $l$ | $\delta_l$ (barrier $= 0$) | $\delta_l$ (barrier $= -2$) | Difference 2 | 0.844 | 0.840 | $0.44 3 | 1.195 | 1.191 | $0.34 4 | 1.555 | 1.551 | $0.28 5 | 1.924 | 1.920 | $0.22 The convergence rate scales as $\sim 2/[l(l{+}1)]$, consistent with the barrier shift being a first-order perturbation to the centrifugal potential. -- Cross-checks -- The angular decomposition provides exact internal consistency checks: [Table: Cross-check identities (V2.225, $N = 600$).] lccc@{}} Relation | Expected | Measured | Match $\delta_{\rm vec}^{\rm TT} = 2 \times \delta_{\rm scalar}(l \geq 1)$ | $-0.354$ | $-0.354$ | Exact $\delta_{\rm grav}^{\rm TT} = 2 \times \delta_{\rm scalar}(l \geq 2, \Delta_s = -2)$ | $-0.687$ | $-0.687$ | Exact These identities hold to six or more significant figures, confirming that the lattice correctly decomposes the field content by angular momentum and that the code is internally consistent. -- The monopole prediction from the 50 If the 50 specific value for the monopole contribution. The vector TT sector sums $l \geq 1$ modes, so: \delta_{\rm vec}^{\rm TT} = \frac{\delta_{\rm vec}^{\rm full}}{2} = -\frac{31}{90} , which implies \delta_{\rm scalar}(l \geq 1) = -\frac{31}{180} , and therefore \delta_{l=0} = -\frac{1}{90} - \left(-\frac{31}{180}\right) = \frac{29}{180} = 0.16111\ldots Our measurement at $N = 1000$ gives $\delta_{l=0} = +0.166$, which is $3.3 above the prediction $29/180$. Since this $3.3 $N = 600$ to $N = 1000$ (it was $3.1 50 exactly equal to, $29/180$. === The universal scaling function $f(x = l/n)$ === -- Discovery of scaling collapse -- The per-mode entanglement entropy $S_l(n)$ depends on $l$ and $n$ individually. However, when plotted as a function of the ratio $x = l/n$, the data collapse onto a single universal curve: S_l(n) \to f(x) as n \to \infty , x = l/n . The quality of the collapse is excellent. At $x = 1.0$, the spread across $n = 8$-$30$ is only $0.2 [Table: Scaling collapse of $S_l$ at fixed $x = l/n$ (V2.230).] ccccccc@{}} $x$ | $n = 8$ | $n = 12$ | $n = 16$ | $n = 20$ | $n = 30$ | Spread 0.5 | 0.1273 | 0.1300 | 0.1313 | 0.1321 | 0.1332 | $8 1.0 | 0.0558 | 0.0559 | 0.0559 | 0.0559 | 0.0559 | $0.2 2.0 | 0.0146 | 0.0143 | 0.0141 | 0.0140 | 0.0139 | $5 3.0 | 0.0051 | 0.0049 | 0.0048 | 0.0047 | 0.0047 | $8 The deviations at $x \neq 1$ are finite-$n$ corrections that decrease as $n$ increases. The physical interpretation is that $x = l/n$ is the ratio of angular wavelength to the entangling sphere radius, and the universal function $f(x)$ encodes how vacuum correlations project onto each angular scale. The scaling collapse is not exact at finite $n$--it is an asymptotic property of the $n \to \infty$ limit. The finite-$n$ corrections are of order $1/n$ for generic $x$, but at $x = 1$ the corrections are anomalously small ($0.2 at $n = 8$). This suggests that $x = 1$ (the mode whose angular wavelength equals the sphere radius) has special significance, perhaps related to the crossover between the "interior" and "exterior" regimes of the mode. The existence of a scaling function means that the per-mode entropy is determined by a single dimensionless parameter $x = l/n$ in the thermodynamic limit. This is reminiscent of scaling functions in critical phenomena, where universal functions emerge at second-order phase transitions. Here, the "critical point" is the entangling surface itself, and the scaling variable $x$ measures the ratio of the mode's angular size to the size of the entangling region. -- Properties of $f(x)$ -- The scaling function at $n = 30$ has the following features: [Table: Properties of the universal scaling function $f(x)$ at $n = 30$ (V2.230, V2.231).] lc@{}} Property | Value $f(0)$ | $0.651$ (grows as $\ln n$; diverges as $n \to \infty$) $f(0.5)$ | 0.133 $f(1.0)$ | 0.056 $f(2.0)$ | 0.014 $f(5.0)$ | 0.001 Half-maximum | $x = 0.13$ The function $f(x)$ is rapidly decreasing: it drops by a factor of 12 from $x = 0$ to $x = 1$ and by a factor of 650 from $x = 0$ to $x = 5$. The half-maximum at $x = 0.13$ means most of the per-mode entropy comes from very low angular momenta. However, the area law is determined by the integral weighted by $(2l{+}1) \approx 2nx$, which peaks at $x \sim 0.5$. -- Logarithmic divergence at $x = 0$ -- The $s$-wave ($l = 0$) entropy $f(0, n) = S_0(n)$ grows logarithmically with $n$: f(0, n) = B \ln n + const , B = 0.161 \approx \frac{1}{2\pi} = 0.159 . This connection between the $s$-wave entropy and the logarithmic correction $\delta$ is significant: the same physics that determines the trace anomaly coefficient appears in the $n$-dependence of the monopole mode. The coefficient $B$ matches $1/(2\pi)$ to $1.2 -- Alpha from the spectral integral -- The area-law coefficient can be independently computed from $f(x)$ via the spectral integral: \alphas = \frac{1}{4\pi}\int_0^{C} f(x) (2x + 1/n) \dd x \xrightarrow{n \to \infty} \frac{1}{2\pi}\int_0^{\infty} f(x) x \dd x . [Table: Area-law coefficient from the spectral integral (V2.230).] cccc@{}} $n$ | $C$ | $\alpha_{\rm integral}$ | Error 10 | 8 | 0.024419 | $+3.9 15 | 8 | 0.023790 | $+1.2 20 | 8 | 0.023466 | $-0.2 25 | 8 | 0.023269 | $-1.0 30 | 10 | 0.023488 | $-0.09 At $n = 30$, $C = 10$, the spectral integral gives $\alphas = 0.02349$, within $0.09 the conjectured value by a completely independent computational route. The cumulative integral reveals which angular momenta contribute to the area law: - $50 - $90 The area law is not dominated by the $s$-wave or low-$l$ modes; it comes from a broad range of angular momenta peaking around $l \sim n/2$ but with a long tail extending to $l \sim 10n$. The spectral integral convergence as a function of $x_{\max}$ is: [Table: Convergence of the spectral integral $\int_0^{x_{\max] ccc@{}} $x_{\max}$ | $\int_0^{x_{\max}} x f(x) \dd x$ | Error from $\sqrt{\pi}/12$ 5 | 0.1340 | $-9.3 10 | 0.1447 | $-2.1 15 | 0.1472 | $-0.3 20 | 0.1482 | $+0.4 The slow convergence (requiring $x_{\max} > 15$ for $< 1 the extended tail of $f(x)$. This is why large $C$ values are needed in the direct lattice computation: modes with $l$ up to $10n$ contribute significantly to the area law. -- The identity to prove -- The spectral integral identity can be stated precisely: \frac{1}{2\pi}\int_0^{\infty} x f(x) \dd x = \frac{1}{24\sqrt{\pi}} = \frac{\sqrt{\pi}}{12 \cdot 2\pi} . Equivalently, \int_0^{\infty} x f(x) \dd x = \frac{\sqrt{\pi}}{12} = 0.14770\ldots This is a well-defined integral identity connecting the scaling function $f(x)$ (a numerically known but analytically uncharacterized function) to the special value $\sqrt{\pi}/12$. A proof would require either (i) finding the analytic form of $f(x)$ and evaluating the integral, or (ii) computing the integral directly from the spectral theory of the tridiagonal matrix $K'_l$ without explicitly constructing $f(x)$. -- No simple closed form -- We systematically fit $f(x)$ to twelve candidate analytic forms. The results are summarized in (ref: tab:fx_fits). [Table: Fits of $f(x)$ to candidate analytic forms (V2.231, $n = 30$). The "integral error" is the discrepancy between the integral of the fitted function (weighted by $x$) and the target $\sqrt{\pi] lccc@{}} Function | $R^2$ | Integral error | Parameters Stretched exponential | 0.99987 | $-5.0 Double exponential | 0.99889 | $-24.6 Rational exponential | 0.99776 | $+9.4 $\exp(-b\sqrt{x})$ | 0.98938 | $+67.2 Simple exponential | 0.97112 | $-52.0 Heat kernel | 0.96148 | $-58.4 Power law | 0.95658 | $+1183 Gaussian | 0.88811 | $-69.6 No candidate function simultaneously achieves $R^2 > 0.999$ and integral error $< 1 f(x) \approx 0.654 \exp\left(-2.47 x^{0.626}\right), which captures the body of the function ($R^2 = 0.99987$) but misses the integral by $5 tail at $x > 10$. -- Tail behavior -- The tail of $f(x)$ decays more slowly than any simple exponential. The log- derivative $d(\ln f)/dx$ decreases monotonically with $x$: [Table: Log-derivative of $f(x)$ (V2.231).] cc@{}} $x$ | $d(\ln f)/dx$ 0.5 | $-1.95$ 1.0 | $-1.58$ 2.0 | $-1.23$ 3.0 | $-0.98$ 5.0 | $-0.67$ 7.0 | $-0.50$ A pure exponential would give a constant log-derivative. The monotonic decrease rules out simple exponential decay and points to a function that transitions from exponential-like at small $x$ to power-law-like at large $x$. -- Convergence of the stretching exponent -- The stretching exponent $c$ in the best-fit function (eq: eq:stretched_exp) is still evolving at $n = 30$: [Table: Stretching exponent $c$ as a function of $n$ (V2.231).] ccc@{}} $n$ | $c$ | $R^2$ 8 | 0.744 | 0.99988 12 | 0.705 | 0.99990 16 | 0.679 | 0.99991 20 | 0.659 | 0.99990 25 | 0.640 | 0.99989 30 | 0.626 | 0.99987 Extrapolating to $n \to \infty$ suggests $c \sim 0.55$-$0.60$, but this remains uncertain without larger-$n$ data. The value is tantalizingly close to simple fractions such as $5/8$ or $3/5$, though this may be coincidental. -- Moment analysis -- The $x$-weighted distribution $g(x) = x f(x)/\int$ that determines $\alphas$ has the following moments: [Table: Moments of $x f(x)$ (V2.231).] lccc@{}} Moment | $x f(x)$ | Gaussian | Exponential Mean | 1.88 | -- | -- Std dev | 1.78 | -- | -- Skewness | 1.85 | 0 | 2 Excess kurtosis | 3.69 | 0 | 6 The skewness ($1.85$) is close to but less than the exponential value ($2$), and the excess kurtosis ($3.69$) is between the Gaussian ($0$) and exponential ($6$) values. The distribution is strongly right-skewed but softer than exponential. === Boundary mode dominance === -- Single-mode fraction -- The modular spectrum of the reduced density matrix for each angular momentum channel $l$ consists of a set of symplectic eigenvalues ${\nu_k}$, $k = 1, \ldots, n$, ordered by decreasing entropy contribution $h(\nu_k)$. We find that the largest eigenvalue--the "boundary mode"--overwhelmingly dominates: [Table: Fraction of entropy carried by the boundary mode, weighted by $(2l{+] cc@{}} $x = l/n$ | Boundary mode fraction 0.0 | 95.1 0.5 | 99.7 1.0 | 100.0 $\geq 2.0$ | 100.0 $(2l{+}1)$-weighted average | 99.8 The boundary mode carries $99.8 The only channel where sub-dominant modes contribute appreciably is the $s$-wave ($l = 0$), where the centrifugal barrier vanishes and deeper radial modes retain some entanglement. -- Spatial localization -- The boundary mode is spatially localized at the entangling surface. Using the Williamson eigenvector center-of-mass to define the "distance from boundary" $d$ of each mode, we find: [Table: Modular spectrum of the $l = 0$ channel at $n = 25$ (V2.233).] ccccc@{}} Mode | $\varepsilon_k$ | $T = 1/\varepsilon$ | $d_{\rm com}$ | Entropy fraction Hottest | 1.69 | 0.59 | 0.23 | 94.4 2nd | 5.23 | 0.19 | 1.05 | 5.4 3rd | 9.09 | 0.11 | 1.98 | 0.2 All else | $> 13$ | $< 0.08$ | $> 2.9$ | $\sim 0 The boundary mode sits at $d \approx 0.23$ lattice spacings from the entangling surface. The correlation between modular energy and distance is nearly perfect: $\rho(\varepsilon, d) = 0.997$. -- Radial entropy profile -- Summing over all angular momentum channels, the spatial distribution of entropy is unambiguous: - $99.86 - $100.00 - The intermediate zone ($d > 3$) and deep interior ($d > n/2$) contribute exactly zero to numerical precision. This confirms that the entanglement entropy is entirely boundary-localized. Both the area-law contribution ($\alpha$) and the log correction ($\delta$) originate from the same ultraviolet boundary layer. Their distinction is not spatial but geometrical: $\alpha$ is the dominant $n^2$ scaling and $\delta$ is the sub- dominant $\ln n$ correction from the curvature of the entangling surface. This finding refutes the intuitive picture in which the area law comes from boundary modes and the log correction comes from modes at intermediate distance. Instead, both $\alpha$ and $\delta$ originate from the same ultraviolet degrees of freedom. The log correction is not a separate "bulk" contribution but a geometrical property of the boundary mode: as the sphere grows, the boundary mode's eigenvalue changes in a way that encodes the curvature of the entangling surface. The result has a direct analog in black hole thermodynamics. There, the Bekenstein-Hawking entropy $S_{\rm BH} = A/(4G)$ is entirely localized at the horizon. The quantum corrections--including the logarithmic correction $\delta \ln A$--arise from the same horizon degrees of freedom, not from fields in the black hole interior. Our lattice computation provides the first explicit numerical demonstration of this localization in a controlled quantum field theory setting. -- Bisognano-Wichmann and Casini-Huerta-Myers tests -- We compare the modular energies to two theoretical predictions: - Bisognano-Wichmann (BW): For a half-space, $\varepsilon = 2\pi d$, predicting a linear relationship. - Casini-Huerta-Myers (CHM): For a sphere in a CFT, $\varepsilon = \pi(n^2 - r^2)/n$, predicting a parabolic relationship. [Table: Comparison of modular spectrum with BW and CHM predictions (V2.233).] ccc@{}} $n$ | $R^2$ (BW) | $R^2$ (CHM) 15 | 0.769 | 0.883 20 | 0.714 | 0.841 25 | 0.691 | 0.815 30 | 0.711 | 0.814 The CHM (sphere) prediction consistently outperforms the BW (half-space) prediction, confirming that the geometry of the entangling surface matters for the modular spectrum. The fact that the log correction $\delta$ is a curvature effect is consistent with this finding: CHM incorporates the sphere curvature that BW neglects. -- Reduction to a single-eigenvalue problem -- The boundary mode dominance reduces the conjecture $\alphas = 1/(24\sqrt{\pi})$ to a spectral theory problem. Define the boundary mode function g(x) = \nu_{\max}(l = xn, n) as n \to \infty . Then \alphas \approx \frac{1}{2\pi}\int_0^{\infty} h\left(g(x)\right) x \dd x , where $h$ is the entropy function (eq: eq:hnu). The function $g(x)$ shows universal scaling for $x \geq 0.5$ (0.2 across $n = 8$-$30$) but diverges logarithmically at $x = 0$: g(0, n) = 0.076 \ln n + 0.482 . Like $f(x)$, the function $g(x)$ has three regimes: - $x \to 0$: logarithmic divergence (from the $s$-wave); - $0 < x < 1$: rapid decay (stretched exponential, $c \approx 0.6$); - $x > 1$: exponential approach to $\nu = 1/2$ (vacuum; centrifugal barrier suppresses entanglement). The best analytic fit for $g(x) - 1/2$ is a stretched exponential with $R^2 = 0.99967$, but its integral gives $\alpha = 0.0203$, $13.6 target. The missing contribution comes from the $x = 0$ singularity and the extended tail. The single-mode integral at $n = 30$ gives $\alpha = 0.02309$, capturing $98.2 modes in the modular spectrum. -- Physical implications -- The boundary mode reduction has three important consequences: - The area law is a single-mode effect. Despite the entropy being a trace over $n$ modes per angular momentum channel, $99.8 single mode per channel. The Srednicki entropy is not a collective phenomenon of many modes; it is dominated by the mode with the strongest coupling across the entangling surface. - Both $\alpha$ and $\delta$ share the same boundary layer. The area-law coefficient and the log correction both originate from the boundary mode. Their distinction is in the $n$-dependence of this mode's eigenvalue: the leading $n^2$ scaling gives $\alpha$ and the subleading $\ln n$ correction gives $\delta$. - The conjecture is a well-posed spectral problem. Proving $\alphas = 1/(24\sqrt{\pi})$ reduces to finding the largest symplectic eigenvalue of a structured tridiagonal matrix with centrifugal potential $l(l{+}1)/j^2$. This connects to the theory of Jacobi matrices and Szeg\H{o}-type limit theorems. === The continuum limit and the lattice as physics === -- The Bessel function approximation -- In the continuum limit, the eigenvectors of the radial coupling matrix $K'_l$ approach Bessel functions $J_{l+1/2}(kr)$. The boundary correlator at the entangling surface ($r = n$) would then be given by a Bessel integral: X[n,n] \sim \frac{1}{\pi}\int_0^{\infty} n J_{l+1/2}^2(kn) \frac{\dd k}{k} . Using Debye (WKB) asymptotics for Bessel functions in the limit of large order, this yields a closed-form expression for the scaling function: \nu_{\rm WKB}(x) = \frac{1}{\pi^2}\sqrt{\frac{\arccos(x/\pi) \sqrt{\pi^2 - x^2}}{x}} . This formula has appealing properties: it diverges as $x \to 0$, vanishes at $x = \pi$ (the UV cutoff), and has the correct qualitative shape. -- The continuum approximation fails by a factor of 20 -- Quantitatively, the Bessel/WKB prediction is catastrophically wrong: [Table: Lattice vs Bessel/WKB boundary correlator (V2.236, $n = 20$).] cccc@{}} $x$ | $\nu_{\rm lattice}$ | $\nu_{\rm Bessel}$ | Ratio 0 | 0.756 | 1.420 | 0.53 0.5 | 0.551 | 0.290 | 1.90 1.0 | 0.518 | 0.191 | 2.71 2.0 | 0.504 | 0.102 | 4.95 At $x = 0$, the Bessel prediction is too large by a factor of 2. At $x \geq 0.5$, it is too small by factors of 2-5. The discrepancy grows with $x$. The area-law coefficient from the WKB formula is $\alpha_{\rm WKB} = 0.0012$, a factor of $\sim20$ below the lattice value $\alphas = 0.0235$: [Table: Area-law coefficient: lattice vs continuum predictions (V2.236).] lcc@{}} Method | $\alpha$ | Error Lattice (2nd difference) | 0.02244 | $-4.5 Bessel integral ($n = 20$) | 0.00128 | $-94.6 WKB numerical | 0.00121 | $-94.9 WKB closed form | 0.00122 | $-94.8 Target: $1/(24\sqrt{\pi})$ | 0.02351 | -- The Bessel and WKB predictions are converged with each other (they agree to $5 continuum approximation fundamentally fails to capture the boundary correlator. -- Why the continuum limit fails -- The area-law coefficient is a UV-sensitive quantity that depends on the physics at the entangling surface--specifically, on the mode structure at the lattice boundary. The Bessel function approximation is valid for smooth modes far from the boundary, but fails precisely for the modes that contribute most to the entanglement: those localized within one or two lattice spacings of the entangling surface. The lattice boundary condition (the entangling surface sits between two discrete lattice sites) is fundamentally different from the continuum boundary condition (smooth matching of wavefunctions). The lattice eigenvector amplitudes at the boundary are 3-8 times smaller than the Bessel prediction, because the discrete structure modifies the mode functions at the UV scale. -- The lattice is the physics -- This result has a profound implication for the interpretation of the area-law coefficient. In the entanglement entropy framework, the entangling surface is at the Planck scale. There is no "continuum limit" at the boundary--the discrete structure is the physical structure. This explains three puzzles: - Why no closed form for $f(x)$ or $g(x)$ exists ((ref: sec:no_closed_form,sec:boundary)): These functions are determined by the discrete lattice eigenvectors near the boundary, not by continuum Bessel functions. They are intrinsically lattice quantities. - Why the scaling function converges rapidly with $n$: The lattice structure at the boundary is self-similar as $n$ increases. Once the system is large enough that the boundary at $r = N$ is far from the entangling surface at $r = n$, the boundary mode is determined entirely by local lattice structure. - Why $\alphas$ is independent of the system size: The area coefficient is a property of the local lattice structure at the entangling surface, not the global geometry. It converges to its final value as soon as $N/n$ is large enough. -- Implications for proving $\alphas = 1/(24\sqrt{\pi -- )$} The failure of the continuum limit redirects the search for an analytic proof. The proof must work within the lattice framework, using methods specific to discrete systems: - Tridiagonal matrix theory: $\alphas$ comes from the eigenvalue structure of the lattice matrix $K'_l$, not its continuum limit. The result might follow from Szeg\H{o}-type theorems applied to the specific tridiagonal structure. - Combinatorial identity: The number $1/(24\sqrt{\pi}) = 1/(4!\cdot\Gamma(3/2))$ might emerge from combinatorial counting of lattice boundary configurations. - Functional equation: The scaling function $g(x)$ satisfies a discrete recursion from the tridiagonal structure that might have an exact solution via finite difference methods. The key lesson is that $\alphas$ encodes Planck-scale physics. A proof must derive this quantity from the discrete entangling surface, not from a continuum approximation to it. -- What the failure means for the framework -- The failure of the continuum limit might seem to undermine the physical significance of $\alphas$. If the number depends on the lattice structure, how can it have physical meaning? The resolution is that in the entanglement entropy framework, the UV cutoff is physical--it is the Planck length. The lattice is not an approximation to an underlying continuum; it is the fundamental description at the Planck scale. The area-law coefficient $\alpha = 1/(4G\epsilon^2)$ relates Newton's constant to the cutoff, and $\alphas$ (the coefficient per scalar degree of freedom) is the lattice quantity that makes this relation precise. The regulator dependence of $\alpha$ ((ref: tab:alpha_regulator)) is then not a bug but a feature: different lattice actions correspond to different microscopic theories of quantum gravity, each with its own value of $G$. The ratios between different fields ($\alpha_{\rm v}/\alphas = 2$, $\alpha_{\rm g}/\alphas = 2$) are universal because they do not depend on the microscopic details of the cutoff, only on the representation content of the fields. Similarly, $\delta$ is regulator-independent because it is determined by the trace anomaly, a topological quantity that depends only on the field content. The ratio $\Omega_\Lambda = \abs{\delta}/({6\alpha})$ is therefore determined by the Standard Model field content and the single lattice number $\alphas$, which is verified to $0.01 The continuum failure thus strengthens, rather than weakens, the framework's physical interpretation. The area-law coefficient is not a "divergent artifact" that should be subtracted away (as in the cosmological constant problem's standard formulation); it is a physical quantity that encodes the relationship between the UV cutoff and the gravitational coupling. The lattice computation directly measures this relationship. === Discussion and open problems === -- Summary of precision achievements -- The program has established the following precision benchmarks: [Table: Summary of precision achievements.] llcc@{}} Quantity | Method | Precision | Experiment $\alphas = 1/(24\sqrt{\pi})$ | Richardson ($C$-extrap.) | $0.01 $\alphas$ (spectral integral) | $\int f(x) x \dd x$ | $0.09 $\delta_{\rm s} = -1/90$ | d3S at $N = 1000$ | $0.02 $\delta_{\rm v} = -31/45$ | Richardson ($N$-extrap.) | $0.01 $\delta$ universality | 3 stencils, $l = 0$ | $0.24 $\delta$ $C$-independence | $C = 2$-$22$ | $0.003 50 $\alpha$ ratios | multi-$C$, multi-stencil | $0.05 -- Comparison with previous numerical work -- Several groups have computed entanglement entropy on the lattice, and it is useful to compare our results with the existing literature. Srednicki's original computation [ref: Srednicki:1993im] established the area law and gave $\alphas \approx 0.03$ with $\sim30 small lattice with $N = 100$, $n_{\max} = 10$). Lohmayer et al. [ref: Lohmayer:2009sq] improved this to $\alphas \approx 0.0236$ with a few-percent precision using the angular momentum decomposition and larger lattices. Our result $\alphas = 0.023510$ (error $0.011 these by two to three orders of magnitude. For the log coefficient, Lohmayer et al.\ extracted $\delta_{\rm s} \approx -0.012$ with $\sim10 $0.02 is using large $N/n_{\max}$ ratios ($> 12$) to eliminate finite-size effects, combined with the d3S extraction method. No previous numerical work has extracted the vector or graviton log coefficients on the lattice. The results reported here ($\delta_{\rm v}$ to $0.01 $\delta_{\rm g}^{\rm TT}$ with the 50 -- Consistency checks across experiments -- An important feature of this program is the internal consistency across sixteen independent experiments. The same quantities are extracted by different methods, at different lattice sizes, and with different parameter choices: - $\alphas$: confirmed by $C$-extrapolation (V2.192), spectral integral (V2.230), and $N$-extrapolation (V2.220). All three methods agree to $0.1 - $\delta_{\rm s} = -1/90$: confirmed at $0.02 (V2.220), and at $0.6 - $\delta_{\rm v} = -31/45$: confirmed at $0.01 $N$-extrapolation (V2.198), at $3.5 (V2.199), and at $4.1 - The 50 in V2.219, and independently measured in V2.220 and V2.226. - Area-law ratios $\alpha_{\rm v}/\alphas = \alpha_{\rm g}/\alphas = 2$: confirmed in V2.195, V2.218, V2.219, across all $C$ values and all three stencils. The consistency of these cross-checks gives confidence that the results are robust and not artifacts of specific parameter choices. -- Limitations and caveats -- We emphasize several limitations of this work: - Free fields only. All computations use free (non-interacting) fields. The trace anomaly coefficients $\delta$ are exact UV quantities that are not modified by interactions [ref: Duff:1993wm], but the area-law coefficient $\alpha$ could in principle receive interaction corrections. - Fermion $\delta$ cannot be verified. The angular momentum decomposition fails for fermions: the per-mode entropy decay ($\sim \kappa^{-1.8}$) is too slow for the sum to converge. The 45 Weyl fermions in the Standard Model contribute $\delta_{\rm ferm} = -2.75$ (25 be independently verified on a bosonic lattice. The value relies on the continuum heat kernel result. - The graviton model is simplified. We model the graviton as 2 polarizations with the Lichnerowicz barrier $l(l{+}1) - 2$, starting at $l = 2$. This captures the transverse-traceless modes but not the full BRST structure of linearized gravity. Computing edge modes directly would require the Donnelly- Wall extended Hilbert space formalism, which is beyond the scope of this work. - The 50 $\delta_{\rm TT}/\delta_{\rm EE} = 0.508$ deviates from $1/2$ by $1.5 While this is consistent with the level of finite-$N$ corrections observed in the scalar calibration, it could also indicate a genuine asymmetry between physical and edge mode contributions. A definitive resolution requires either much larger lattices ($N \geq 2000$) or an independent analytical argument. - The scaling function has no closed form. The universal function $f(x)$ is well-characterized numerically but resists analytic identification. The proof that $\alphas = 1/(24\sqrt{\pi})$ remains an open problem. -- The fermion problem -- The inability to verify fermion entanglement entropy on the lattice is a fundamental limitation. Earlier experiments in this program found that the Wilson fermion parameter $r_W$ controls the ratio $\alpha_{\rm Weyl}/\alphas$, which ranges from 0.76 ($r_W = 3.0$) to 18.3 ($r_W = 0.5$) and passes through $2.0$ at $r_W \approx 1.7$. The root cause is that $\alpha$ is UV-sensitive and any fermion discretization modifies the UV structure. The bosonic area-law ratios ($\alpha_{\rm v}/\alphas = 2$ and $\alpha_{\rm g}/\alphas = 2$) are verified because the radial equation is identical for all bosonic spins; only the angular momentum range differs. Fermions have a qualitatively different radial structure, and the lattice cannot faithfully reproduce it. For the log coefficient $\delta$, the situation is reversed: $\delta$ is UV- finite and in principle verifiable. However, the per-mode entropy decay for fermions is too slow ($\kappa^{-1.8}$) for the angular momentum sum to converge on accessible lattices. The 90 Weyl effective scalar degrees of freedom (45 Weyl fermions $\times$ 2 from the heat kernel relation) must be taken from the continuum result. -- Connection to the cosmological constant prediction -- The precision results of this paper serve as inputs to the cosmological constant prediction developed in companion papers of the Moon Walk Project. The prediction has the form \Omega_\Lambda = \frac{\abs{\delta_{\rm total}}}{6 \alpha_{\rm total}} , where $\delta_{\rm total}$ is the sum of trace anomaly coefficients for all fields and $\alpha_{\rm total}$ is the corresponding sum of area-law coefficients. The inputs verified in this paper are: - $\alphas = 1/(24\sqrt{\pi}) = 0.02351$ (to $0.01 - $\delta_{\rm s} = -1/90$ (to $0.02 - $\delta_{\rm v} = -31/45$ (to $0.01 - $\alpha_{\rm v}/\alphas = \alpha_{\rm g}/\alphas = 2$ (to $0.05 - $\delta$ is regulator-independent (to $0.24 - $\delta$ is mass-independent for $m \ll 1/\epsilon$ (to $0.02 - $\delta$ is $C$-independent (to $0.003 The graviton contribution to the prediction depends on the edge-mode fraction $f_g = \delta_{\rm EE}/\delta_{\rm EA}$, which the companion paper derives analytically as $f_g = 61/212$. Our lattice measurement of $\delta_{\rm TT}/\delta_{\rm EE} = 0.508$ provides independent evidence that the TT modes carry approximately half the entanglement entropy anomaly, which is consistent with--but does not directly determine--the edge-mode fraction. -- Open problems -- Several questions raised by this work remain open. We organize them by estimated difficulty and potential impact. - Near-term (likely solvable with current methods) - - Push $\delta_{\rm s$ to $0.001 ($0.02 $N = 2000$, $n_{\max} = 150$. This would require $\sim10$ hours of computation on a laptop and would verify the scalar trace anomaly to unprecedented precision. - Richardson-extrapolate the graviton TT $\delta$. The graviton value is already stable across $N$, but Richardson extrapolation in $N$ (as done for the vector in V2.198) might sharpen the 50 $0.508$ to determine whether it converges to a simple fraction. - Extend the per-$l$ decomposition to $l = 10$-$20$. The quadratic law (eq: eq:delta_l_fit) is fitted to $l = 0$-$5$. Extending to higher $l$ would determine whether the spectrum remains quadratic or develops new structure. - Medium-term (require new methods or substantial computation) - - Derive the 50 entanglement entropy that predicts the TT-mode fraction to be approximately $1/2$? The Donnelly-Wall framework predicts an $O(1)$ edge mode contribution but does not predict the specific fraction. A derivation would connect our numerical result to the algebraic structure of gauge constraints at the entangling surface. - Determine the analytic form of $f(x)$. The scaling function has a logarithmic singularity at $x = 0$, a stretched-exponential body, and a power- law tail. It likely involves special functions related to the spectral theory of the tridiagonal matrix $K'_l$. Possible approaches include: the inverse spectral theory of Jacobi matrices, the Marchenko-Pastur distribution for structured random matrices, or the Szeg\H{o} limit theorem for Toeplitz-like operators. - Compute edge modes directly. Implementing the Donnelly-Wall extended Hilbert space on the lattice would allow a direct computation of the edge mode contribution, testing the 50 would require implementing the gauge- invariant factorization of the Hilbert space at the entangling surface. - Extend to fermions. Either develop a lattice fermion formulation that preserves the correct UV structure for entanglement entropy, or find an alternative method (e.g., mutual information, which is UV-insensitive) to verify the fermion $\delta$ independently of the heat kernel. - Long-term (fundamental theoretical challenges) - - Prove $\alphas = 1/(24\sqrt{\pi)$ analytically.} The problem is now cleanly formulated as a spectral theory question ((ref: sec:boundary)), but no proof exists. The failure of the continuum limit ((ref: sec:continuum)) means the proof must use lattice-specific methods. The factorization $1/(4! \Gamma(3/2))$ hints at a combinatorial or special-function identity, but the connection to the tridiagonal matrix structure remains unclear. - Understand the per-$l$ spectrum analytically. The quadratic law $\delta_l \approx 0.167 + 0.330 l + 0.004 l^2$ and the fact that all low-$l$ modes contribute positively demand an explanation. The linear slope $0.330 \approx 1/3$ may have a connection to the central charge of a 1D boson. Understanding this spectrum would illuminate the angular momentum structure of quantum vacuum fluctuations. - Extend to curved spacetime. All computations in this paper use flat space. The cosmological constant prediction applies to de Sitter space, where the entangling surface is the cosmological horizon. While the trace anomaly coefficients are unchanged, the area-law coefficient could in principle receive curvature corrections. Earlier experiments in this program (V2.68, V2.69) found that $\alpha$ is unchanged by de Sitter curvature to $< 1 === Conclusion === We have presented a comprehensive lattice study of entanglement entropy for free bosonic fields across a spherical surface, achieving precisions ranging from $0.003 The main results are: - The area-law coefficient $\alphas = 1/(24\sqrt{\pi})$ is confirmed to $0.01 integral), establishing it as one of the most precisely known non-trivial quantities in lattice quantum field theory. - The log coefficient $\delta$ is a universal, regulator-independent, mass-independent, cutoff-independent quantity. These properties, verified numerically at the sub-percent level, are precisely what is required for $\delta$ to serve as a physical input to the cosmological constant prediction. - The first lattice extraction of the vector trace anomaly coefficient achieves $0.01 first graviton lattice measurement reveals the "50 modes carry approximately half the full entanglement entropy anomaly. - The angular momentum decomposition of $\delta$ reveals a rich internal structure: all low-$l$ modes contribute positively, and the negative total emerges from massive cancellation with the collective high-$l$ sector. - A universal scaling function $f(x = l/n)$ governs the per-mode entropy. It has no simple closed form, diverges logarithmically at $x = 0$, and has a stretched-exponential body with a slow power-law tail. - The entanglement entropy is overwhelmingly ($99.8 single boundary mode per angular momentum channel, reducing the conjecture $\alphas = 1/(24\sqrt{\pi})$ to a single-eigenvalue spectral problem. - The continuum (Bessel function) approximation fails by a factor of 20, proving that $\alphas$ is an intrinsically lattice quantity. The discrete structure at the entangling surface is not an approximation to a continuum limit --it is the physics. Taken together, these results establish the Srednicki lattice as a precision tool for entanglement entropy computation and provide the numerical foundations for the cosmological constant prediction. Three results stand out as particularly significant for the broader program: First, the area-law coefficient is verified to $0.01 methods. This is the single non-topological input to the cosmological constant prediction, and its precision directly limits the prediction's accuracy. The $0.01 by the graviton edge-mode fraction, not by the lattice computation. Second, the 50 contributions to the graviton entanglement entropy. While the rule is approximate (not exact $1/2$), the measured ratio $0.508$ recovers the Benedetti-Casini analytical value $-61/45$ to $0.1 the full graviton $\delta$. This provides independent numerical support for the edge-mode decomposition that enters the cosmological constant prediction. Third, the failure of the continuum limit proves that $\alphas$ is an intrinsically lattice quantity. This negative result is physically important: it means that the area-law coefficient encodes information about the discrete structure at the entangling surface, consistent with the framework's identification of this surface with the Planck scale. The proof of $\alphas = 1/(24\sqrt{\pi})$, if it exists, must come from discrete mathematics rather than continuum field theory. The key open challenge remains proving $\alphas = 1/(24\sqrt{\pi})$ analytically. The boundary mode reduction ((ref: sec:boundary)) formulates this as a concrete spectral theory problem: find the largest symplectic eigenvalue of a structured tridiagonal matrix with centrifugal potential. The universal scaling function $f(x)$ and the boundary mode function $g(x)$ provide complementary characterizations of the same problem. Progress on either front would elevate the cosmological constant prediction from a precise numerical result to a mathematical theorem. \noindentData availability. All numerical data, parameter files, and analysis code are available in the Moon Walk Project repository. Each experiment (V2.192 through V2.236) includes a complete results file, unit tests, and a self- contained run script for reproducibility. \noindentReproducibility. Every numerical result reported in this paper can be reproduced from a single Python script per experiment, using only NumPy, SciPy, and standard library dependencies. The computational cost of the full program is approximately 50 CPU-hours on a modern laptop, with the most expensive single run (V2.220 at $N = 1000$) requiring roughly 3 hours. === References === [Bombelli:1986rw] L. Bombelli, R. K. Koul, J. Lee and R. D. Sorkin, "Quantum source of entropy for black holes," Phys.\ Rev.\ D 34 (1986) 373. [Srednicki:1993im] M. Srednicki, "Entropy and area," Phys.\ Rev.\ Lett.\ 71 (1993) 666 [arXiv:hep-th/9303048]. [Lohmayer:2009sq] R. Lohmayer, H. Neuberger, A. Schwimmer and S. Theisen, "Numerical determination of entanglement entropy for a sphere," Phys.\ Lett.\ B 685 (2010) 222 [arXiv:0911.4283]. [Peschel:2002jhw] I. Peschel, "Calculation of reduced density matrices from correlation functions," J.\ Phys.\ A 36 (2003) L205 [arXiv:cond-mat/0212631]. [Solodukhin:2008dh] S. N. Solodukhin, "Entanglement entropy, conformal invariance and extrinsic geometry," Phys.\ Lett.\ B 665 (2008) 305 [arXiv:0802.3117]. [Donnelly:2014fua] W. Donnelly and A. C. Wall, "Entanglement entropy of electromagnetic and gravitational fields," Phys.\ Rev.\ D 91 (2015) 044032. [Donnelly:2015hxa] W. Donnelly and A. C. Wall, "Geometric entropy and edge modes of the electromagnetic field," Phys.\ Rev.\ D 94 (2016) 104053. [Duff:1993wm] M. J. Duff, "Twenty years of the Weyl anomaly," Class.\ Quant.\ Grav.\ 11 (1994) 1387 [arXiv:hep-th/9308075]. [Casini:2011kv] H. Casini, M. Huerta and R. C. Myers, "Towards a derivation of holographic entanglement entropy," JHEP 1105 (2011) 036 [arXiv:1102.0440]. [Benedetti:2019uej] V. Benedetti and H. Casini, "Entanglement entropy of linearized gravitons in a sphere," Phys.\ Rev.\ D 101 (2020) 045004 [arXiv:1908.01800]. [Kabat:1995eq] D. Kabat, "Black hole entropy and entropy of entanglement," Nucl.\ Phys.\ B 453 (1995) 281 [arXiv:hep-th/9503016]. [Blommaert:2025pcf] A. Blommaert and G. Colin-Ellerin, "Edge modes and gravitational entropy," arXiv:2501.xxxxx (2025). [Callan:1994py] C. G. Callan and F. Wilczek, "On geometric entropy," Phys.\ Lett.\ B 333 (1994) 55 [arXiv:hep-th/9401072]. [Holzhey:1994we] C. Holzhey, F. Larsen and F. Wilczek, "Geometric and renormalized entropy in conformal field theory," Nucl.\ Phys.\ B 424 (1994) 443 [arXiv:hep-th/9403108]. [MoonWalk:Lambda] Moon Walk Project, "The cosmological constant from entanglement entropy: a lattice derivation," companion paper (2024-2026). \appendix === Convergence tables === For reference, we collect the detailed convergence data from the main experiments. -- Scalar $\delta$ convergence in $N$ -- (Ref: tab:app_scalar_N) gives the full scalar $\delta$ convergence data, including both the raw d3S extraction and the N-dependence fit parameters. [Table: Full scalar $\delta$ convergence data (V2.220, $C = 10$, $n = 30$-$80$).] rcccc@{}} $N$ | $N/n_{\max}$ | $\delta$ | Error | $\delta$ (corrected) 200 | 2.5 | $-0.24483$ | $2103 300 | 3.75 | $-0.04142$ | $273 400 | 5.0 | $-0.01933$ | $74 500 | 6.25 | $-0.01416$ | $27 600 | 7.5 | $-0.01242$ | $12 800 | 10.0 | $-0.01137$ | $2.4 1000 | 12.5 | $-0.01111$ | $0.02 The "corrected" column applies a free-exponent Richardson extrapolation using only the three largest $N$ values. The correction improves the $N = 600$ result from $12 raw $N = 1000$ result ($0.02 -- Vector $\delta$ convergence in $N$ -- [Table: Vector $\delta$ convergence data (V2.198, $C = 8$, $n = 20$-$70$).] rccc@{}} $N$ | $\delta_{\rm v}$ | Error vs $-31/45$ | Ratio to theory 600 | $-0.69729$ | $1.22 800 | $-0.69202$ | $0.45 1000 | $-0.69071$ | $0.26 1300 | $-0.69010$ | $0.18 1800 | $-0.68994$ | $0.15 \multicolumn{2}{@{}l}{Richardson (order 3)} | $0.01 -- Graviton $\delta_{\rm TT -- $ stability} [Table: Graviton $\delta_{\rm TT] cccc@{}} $N$ | $C$ | $\delta_{\rm TT}$ | Ratio to $\delta_{\rm EE}$ 500 | 10 | $-0.6881$ | 0.508 600 | 10 | $-0.6886$ | 0.508 800 | 5 | $-0.6883$ | 0.508 800 | 10 | $-0.6883$ | 0.508 800 | 20 | $-0.6882$ | 0.508 1000 | 10 | $-0.6880$ | 0.508 === The spectral integral: detailed derivation === The total entanglement entropy at subsystem size $n$ with proportional cutoff $C$ is S(n) = \sum_{l=0}^{Cn} (2l+1) S_l(n) . Introducing the scaling variable $x = l/n$ and converting the sum to an integral: S(n) \approx n \int_0^C (2xn + 1) f(x) \dd x = 2n^2 \int_0^C x f(x) \dd x + n \int_0^C f(x) \dd x . The area-law contribution $\alpha \cdot 4\pi n^2$ comes from the first term: \alpha \cdot 4\pi = 2\int_0^C x f(x) \dd x , giving \alpha = \frac{1}{2\pi}\int_0^{\infty} x f(x) \dd x , in the limit $C \to \infty$. The conjectured identity is \int_0^{\infty} x f(x) \dd x = \frac{\sqrt{\pi}}{12} = 0.14770\ldots The second term $n \int f(x) \dd x$ contributes to the $O(n)$ piece of the entropy, which is subdominant to the area law. The log correction $\delta$ comes from the $n$-dependence of $f(x)$ itself: the scaling function is not exactly $n$-independent but has corrections of order $\ln(n)/n$ at $x = 0$ (from the $s$-wave logarithmic divergence) that produce the $\delta \ln n$ term. === Error budget === We summarize the error budget for the main quantities: [Table: Error budget for the precision frontier results.] lccccc@{}} Quantity | Statistical | Finite-$N$ | Finite-$n$ | $C$-cutoff | Total $\alphas$ | negligible | $< 0.001 $\delta_{\rm s}$ | negligible | $0.02 $\delta_{\rm v}$ (Richardson) | negligible | eliminated | $0.01 50 $\delta$ universality | negligible | $\sim2 All computations are deterministic (no Monte Carlo), so statistical errors are absent. The dominant systematic for $\alphas$ is the $C$-extrapolation (eliminated by Richardson); for $\delta_{\rm s}$, it is the finite-$N$ correction (controlled by using $N = 1000$); for the 50 yet eliminated, as the ratio is stable but may have a residual $\sim1