Proof. Let \(K=PBA^{-1}BP\geq0\) and \(\kappa=\lVert K\rVert\). Since \(Bp\in Q\) for \(p\in P\), \[ 0\leq\kappa\leq\frac{\lVert B\rVert^2}{g} \leq g\epsilon^2. \tag{3}\] For \(p\in P\) and \(q\in Q\), the inequality \(QBQ\geq-\epsilon gQ\geq-\epsilon A|_Q\) gives \[\begin{align*} \langle p+q,(A+B)(p+q)\rangle &\geq (1-\epsilon)\langle q,Aq\rangle +2\operatorname{Re}\langle q,Bp\rangle\\ &= (1-\epsilon) \left\|A^{1/2}q+\frac{A^{-1/2}Bp}{1-\epsilon}\right\|^2 -\frac{\langle p,Kp\rangle}{1-\epsilon}. \end{align*}\] On a unit vector this is at least \(-\kappa/(1-\epsilon)\). Consequently \[E(A+B)+\kappa \geq-\frac{\epsilon\kappa}{1-\epsilon} \geq-2g\epsilon^3.\]

For the other direction choose a unit \(p\in P\) with \(Kp=\kappa p\) and set \(q=-A^{-1}Bp\). Then \(\lVert q\rVert\leq\epsilon\) and \[\langle p+q,(A+B)(p+q)\rangle =-\kappa+\langle q,Bq\rangle.\] Writing \(r=\lVert q\rVert^2\), the variational principle therefore gives \[E(A+B)+\kappa \leq \frac{\langle q,Bq\rangle+\kappa r}{1+r} \leq g\epsilon^3+g\epsilon^4 \leq\tfrac32g\epsilon^3.\] The smallest eigenvalue of \(-K\) on \(P\) is \(-\kappa\), proving the claim. If \(Q=0\), the hypothesis \(PBP=0\) forces \(B=0\), and the same formulas hold with all inverse terms zero. ◻

Proof. Take \[A=\Delta\sum_e(h^s_{A_eB_e}+3I),\qquad B=\sqrt\Delta\sum_e v_e(h^s_{LA_e}+h^s_{RC_e}).\] The kernel of \(A\) consists of arbitrary original spins tensored with a singlet on each ancillary pair. Denote its projection by \(P_s\). The gap is \(g=4\Delta\), and \[P_sBP_s=0,\qquad \lVert B\rVert\leq6\sqrt\Delta V.\] A weak interaction from pair \(e\) takes that pair’s singlet into its triplet space and leaves every other ancillary pair in its singlet. Thus every state reached from \(P_s\) in one hop has excitation energy \(4\Delta\). Moreover, terms associated with distinct pairs give zero mixed second-order contributions: a term acting on pair \(f\neq e\) cannot remove the triplet on pair \(e\). This remains true when the two source edges share an original vertex.

For a singlet \(|s\rangle\) on \(A,B\) and Pauli components \(\mu,\nu\), \[\langle s|\sigma_A^\mu\sigma_A^\nu|s\rangle =\delta_{\mu\nu},\qquad \langle s|\sigma_A^\mu\sigma_B^\nu|s\rangle =-\delta_{\mu\nu}.\] Consequently the squares of the two weak terms each compress to \(3I\). Their two ordered cross terms sum to \(2h^s_{LR}\) when \(C_e=A_e\), and to \(-2h^s_{LR}\) when \(C_e=B_e\). The contribution of pair \(e\) to \(-P_sBA^{-1}BP_s\), identified with an operator on the original spins, is therefore \[ -\frac{v_e^2}{4} \begin{cases} 6I+2h^s_{LR},&J_e<0,\\ 6I-2h^s_{LR},&J_e>0 \end{cases} =J_e h^s_{LR}-\frac32v_e^2 I. \tag{7}\] The factor \(4\) in the denominator is the singlet–triplet excitation energy in units of \(\Delta\).

Apply Lemma 3 with \(\epsilon=3V/(2\sqrt\Delta)<1/2\). Its error bound is \(2(4\Delta)\epsilon^3=27V^3/\sqrt\Delta\). Since \(A+B=H_++3k\Delta I\), summing (7) gives the result. If \(k=0\), both perturbation and shift vanish and the assertion is exact. ◻

Proof. Within one gadget the assertion follows directly from (9). For a negative gadget the two distances from \(B_e\) to original endpoints are \(\sqrt2\); for a positive gadget the two nonadjacent endpoint–ancilla distances are \(\sqrt3\). Original vertices are mutually separated by at least \(2\). An original vertex other than the endpoints of an ancilla’s gadget, on that same axis line, has along-line distance at least \(3/2\) to a positive ancilla and at least \(3\) to a negative midpoint. An original vertex off that line differs in a transverse coordinate by at least \(2\).

For ancillas belonging to distinct gadgets on the same axis line, the only along-line separation below \(3/2\) is the separation \(1\) between \(B_e\) and \(A_f\) for consecutive positive gadgets. They have the same height and hence are a contact. All remaining such pairs are at least \(3/2>\sqrt2\) apart. Distinct parallel axis lines have transverse separation at least \(2\).

It remains to compare horizontal and vertical ancillas. Call an ancilla midpoint type if it belongs to a negative gadget, and positive type otherwise. Two midpoint-type projections differ by at least \(1\) in each planar coordinate, so their distance is at least \(\sqrt2\). A midpoint-type and a positive-type projection differ by at least \(1\) and \(1/2\) in the two planar coordinates. Their heights differ by at least \(\sqrt3/2\), since horizontal and vertical heights have opposite signs. Their squared distance is therefore at least \(1+1/4+3/4=2\). Finally two positive-type ancillas from perpendicular directions have height difference \(\sqrt3\), already larger than \(\sqrt2\). These cases exhaust all pairs and also establish distinctness. ◻

Proof. Set the height displacement of every original vertex to zero. The horizontal coordinate of each original vertex will be determined using only its original row, and its vertical planar coordinate using only its original column. These are independent variables even when both directions meet at the same vertex.

Fix one such axis line. Write \(a\) for displacement along it, \(z\) for physical height displacement, and put \(\sigma=1\) for a horizontal line and \(\sigma=-1\) for a vertical line. On every connected run of negative bonds choose the first original \(a\) coordinate to be zero, and accumulate \[ \begin{split} a_R&=a_L+d_{LA_e}+d_{A_eR},\qquad a_{A_e}=a_L+d_{LA_e},\\ a_{B_e}&=a_{A_e},\qquad z_{A_e}=0,\qquad z_{B_e}=\sigma d_{A_eB_e}. \end{split} \tag{11}\] Unconstrained original coordinates may be set to zero. A positive or missing bond places no requirement on the difference of its original endpoints’ \(a\) coordinates at this step, so it separates runs. Each run is a path on a line, and (11) satisfies all three constraints in its negative gadgets.

Now consider the elevated vertices of positive gadgets on this line. Their same-height contacts, including the extra contacts of Lemma 5, form disjoint paths. Set the first \(a\) coordinate on each path to zero and successively set \(a_{\rm next}=a_{\rm current}+d_{\rm edge}\). Every such contact vector is along the axis, so its constraint involves only these \(a\) coordinates. Each elevated vertex has exactly one additional, sloped contact to an original vertex. Writing \(h_0=\sqrt3/2\), set \[ \begin{split} z_{A_e}&=\frac{d_{LA_e}-(a_{A_e}-a_L)/2}{\sigma h_0},\\ z_{B_e}&=\frac{d_{B_eR}-(a_R-a_{B_e})/2}{\sigma h_0}. \end{split} \tag{12}\] Indeed the two unit spoke vectors are \(\tfrac12 e_{\rm axis}+\sigma h_0 e_z\) and \(\tfrac12 e_{\rm axis}-\sigma h_0 e_z\), respectively. Because the original height displacements vanish, these formulas satisfy the remaining equations without imposing any new condition on original coordinates. Set each ancilla’s unused transverse coordinate to zero.

Repeat on all rows and columns. The row procedure changes only original \(x\) coordinates, and the column procedure only original \(y\) coordinates. No ancilla belongs to both procedures. Thus the constructions do not overwrite each other. In particular a square cycle is not a cycle in either scalar accumulation graph; its vector closure is automatic because actual vertex displacements have been assigned.

For the norm bound put \(d_{\max}=\max_{\mathcal C}|d_{ij}|\). Every original along-line coordinate has magnitude at most \(2m d_{\max}\), and every accumulated elevated coordinate has magnitude at most \(m d_{\max}\). Negative ancilla coordinates have magnitude at most \((2m+1)d_{\max}\) along the line and at most \(d_{\max}\) in height. From (12), positive ancilla heights are at most \((1+3m/2)d_{\max}/h_0\). Combining at most two nonzero coordinates at any site gives the stated bound \(6m d_{\max}\). All assignments are additions and multiplication by fixed rational numbers or by \(1/h_0\), proving the algorithmic assertion. ◻

Proof. First use the ideal hoppings \(t^{(0)}\). Let \(A_U=U\sum_i n_{i\uparrow}n_{i\downarrow}\) and let \(T\) be the hopping part of \(H_{\rm Hub}(t^{(0)})\). At filling \(m\), the kernel of \(A_U\) is exactly the singly occupied subspace: having \(m\) fermions, at most one at each of \(m\) sites, forces precisely one per site. Identify this subspace with \(m\) spins, and denote its projection by \(P_{\rm so}\). The complementary excitation energies are positive integer multiples of \(U\). A single hop from \(P_{\rm so}\) creates one vacant site and one doubly occupied site and hence has energy exactly \(U\). In particular \(P_{\rm so}TP_{\rm so}=0\).

A second hop can return to single occupancy only by moving an electron from that doubly occupied site to that vacant site. Thus the two hops must use the same unordered bond; products from different bonds have zero compression, even when the bonds share a vertex. For completeness, consider one bond \(i,j\) of strength \(t\), with local mode order \(i\uparrow,i\downarrow,j\uparrow,j\downarrow\). Write \(|D_i\rangle=c_{i\uparrow}^\dagger c_{i\downarrow}^\dagger|0\rangle\) and similarly for \(|D_j\rangle\). The bond hopping operator \(T_{ij}\) obeys \[T_{ij}|\uparrow,\downarrow\rangle =-t(|D_i\rangle+|D_j\rangle),\qquad T_{ij}|\downarrow,\uparrow\rangle = t(|D_i\rangle+|D_j\rangle).\] It consequently annihilates all triplets and maps the singlet to \(-\sqrt2t(|D_i\rangle+|D_j\rangle)\), of squared norm \(4t^2\). In an ordering containing other sites between \(i\) and \(j\), those sites are singly occupied during this two-hop process; their fixed fermion parity changes transition phases but not the compressed product. Equivalently, permuting the intervening mode pairs gives the same calculation under a fermionic change of basis. It follows that \[ -P_{\rm so}T_{ij}A_U^{-1}T_{ij}P_{\rm so} =-\frac{4t^2}{U}P_{{\rm singlet},ij} =\frac{t^2}{U}(h^s_{ij}-I). \tag{21}\] Here \(P_{{\rm singlet},ij}=(I-h^s_{ij})/4\).

For each spin, \(c_{i\sigma}^\dagger c_{j\sigma} +c_{j\sigma}^\dagger c_{i\sigma}\) has norm one, so \[\lVert T\rVert\leq2\sum_{\mathcal C}t^{(0)}_{ij} =2\eta\sum_{\mathcal C}\sqrt{K_{ij}} \leq2\eta\overline R.\] The last inequality uses \(K_{ij}\leq\Delta\) and \(\Delta\geq1\). Lemma 3, with \(g=U\) and \(\epsilon=2\eta\overline R/U<1/2\), and (21) give \[ \left|E\bigl(H_{\rm Hub}(t^{(0)})\bigr) -\alpha\left(E(H_+)-\sum_{\mathcal C}K_{ij}\right)\right| \leq\frac{16\eta^3\overline R^3}{U^2} \leq\frac{\alpha\gamma}{64}. \tag{22}\] If the complementary subspace is empty, as can happen for \(m=1\), there is no hopping and this assertion is exact.

For the rational hoppings, nonnegativity of \(t^{(0)}_{ij}\) and the accuracy of \(r_{ij}\) imply \(|\max\{f,r_{ij}\}-t^{(0)}_{ij}|\leq f\). Therefore \[ \left\|H_{\rm Hub}(t)-H_{\rm Hub}(t^{(0)})\right\| \leq2q_{\mathcal C}f\leq\frac{\alpha\gamma}{64}. \tag{23}\] Also \(t^{(0)}_{ij}\leq\eta\overline R\leq U/1024\) and \(f\leq U/(128\cdot1024^2)\), which verify \(f\leq t_{ij}<1\). The floor \(f\) is an inverse polynomial. This argument covers an ideal hopping much smaller than \(f\) and an ideal hopping equal to zero; no coefficient lower bound has been used. If \(\mathcal C\) is empty, (23) is simply zero.

Finally, Proposition 4 with (8) contributes at most another \(\alpha\gamma/64\). Adding the three errors proves (20). All sizes, magnitudes, and reciprocal precision budgets in the construction are polynomial in \(n\). Each \(K_{ij}\) is either \(\Delta\), a square root of the rational number \(2\Delta|J_e|\), or zero. The displayed scalar shift is thus efficiently approximable to the asserted precision, completing the proof. ◻

Proof. On the nonzero support of \(\chi_k\), \(r_k\ge D/16\). Since \(\rho_k\le4\), \[|y-X_k-\rho_k e_z|\ge r_k-4\ge r_k/2\] for large \(D\). The secondary singularity is inside the region where the cap is identically zero. Differentiation through order two, including derivatives of the cap, gives \(\|D^j(\chi_kV_k)\|_\infty\le CD^{-j-1}\). Summing proves the derivative bounds. For example, \(q\le1\), \(q'\le1\), and sufficiently large \(D\) give a bound below \(100/D\) for each capped potential, proving the stated function bound.

Put \(H_i^0=-\Delta/2-r_i^{-1}\) and \(B_i=h_i-H_i^0\). Hardy’s inequality and (31) give \[ \|B_i u\|_2\le C\frac mD\|u\|_{H^1}. \tag{32}\] For any potential in these single-well operators, Hardy and Fourier interpolation give, for each \(a>0\), a uniform bound \(\|Vu\|_2\le a\|\Delta u\|_2+C_a\|u\|_2\). The Kato–Rellich theorem for relatively bounded perturbations (Teschl 2009, Theorem 6.4, p. 135) therefore gives domain \(H^2\). The same estimates, with the Laplacian term absorbed, prove uniform graph-norm equivalence: \[ C^{-1}\|u\|_{H^2} \le \|h_i u\|_2+C\|u\|_2 \le C\|u\|_{H^2}. \tag{33}\] The corresponding statement holds for \(H_i^0\). At form level a fixed positive shift of either operator is bounded above and below by constant multiples of \(\|u\|_{H^1}^2\), giving the asserted form domains.

The hydrogenic ground state and its isolation follow without an excited-state formula. Direct integration by parts gives \[ \langle u,(H_i^0+1/2)u\rangle =\frac12\left\|\nabla u+\frac{y-X_i}{r_i}u\right\|_2^2. \tag{34}\] Its zero vectors are exactly the multiples of \(\phi_i\). If normalized vectors orthogonal to \(\phi_i\) had energies tending to \(-1/2\), they would be bounded in \(H^1\). The Coulomb form is compact under weak convergence of such a sequence: within radius \(a\) of the pole its absolute value is bounded by \(Ca\|u\|_{H^1}^2\), by Sobolev and Hölder; beyond radius \(R\) it is bounded by \(R^{-1}\|u\|_2^2\); on the intervening annulus local compactness applies. A weak limit would therefore obey \[\langle u,H_i^0u\rangle\le-1/2,\qquad \langle u,H_i^0u\rangle\ge-\tfrac12\|u\|_2^2\ge-1/2.\] Equality forces norm one and a multiple of \(\phi_i\), contrary to orthogonality. This proves a fixed gap \(g_0>0\).

By (32), after a fixed positive shift the perturbed form lies between \(1-Cm/D\) and \(1+Cm/D\) times the unperturbed shifted form. The variational principle puts its first level within \(Cm/D\) of \(-1/2\), and its second variational level at least \(-1/2+g_0-Cm/D\). The negative first level is attained. The preceding truncation argument gives form compactness of both own Coulomb potentials and of the bounded decaying \(W_i\). A normalized minimizing sequence is bounded in \(H^1\). Passing to a weakly convergent subsequence with limit \(u\), compactness and lower semicontinuity give \[\lambda_i\|u\|_2^2\le\langle u,h_i u\rangle\le\lambda_i<0.\] Since \(\|u\|_2\le1\), these inequalities force \(\|u\|_2=1\) and attainment. For small enough \(m/D\) the first level is simple and separated from the rest by a fixed \(g>0\). Taking the absolute value of a real minimizer cannot increase its energy. The Euler equation and the strong maximum principle on the connected complement of the poles make it positive there.

The form bounds give \(\|\psi_i\|_{H^1}\le C\). Its equation implies \[(H_i^0+1/2)\psi_i =(\lambda_i+1/2)\psi_i-B_i\psi_i,\] with right side of norm at most \(Cm/D\). On \(\phi_i^\perp\), the inverse of \(H_i^0+1/2\) is bounded from \(L^2\) to \(H^2\): use the gap for the \(L^2\) bound, then (33). The orthogonal component of \(\psi_i\) consequently has \(H^2\) norm \(O(m/D)\). Normalization and the positive overlap with \(\phi_i\) put its coefficient along \(\phi_i\) within \(O((m/D)^2)\) of one. This proves (30) and the uniform sup bound. Positivity also holds at the poles by continuity and this sup-norm approximation, since both poles lie within distance four of \(X_i\). ◻

Proof. Taking the inner product of the eigen-equation with \(\phi_i\) gives \[ \lambda_i=-\frac12+(1-q)-q'f(\rho_i) -\langle\phi_i,W_i\phi_i\rangle+O((m/D)^2), \tag{37}\] where \[f(\rho)=\int_{\mathbb R^3}\frac{\pi^{-1}e^{-2|y|}}{|y-\rho e_z|}\,\mathrm dy.\] Here \(\langle\phi_i,\psi_i\rangle=1+O((m/D)^2)\), and replacing \(\psi_i\) by \(\phi_i\) in a perturbing expectation costs \(O((m/D)^2)\), by (32) and (30). Also \(\langle\phi_i,r_i^{-1}\phi_i\rangle=1\). The spherical average of \(|r\omega-\rho e_z|^{-1}\) is \(1/\max(r,\rho)\), as direct integration in \(\omega\cdot e_z\) shows. Consequently, \[\begin{align*} f(\rho)&=\rho^{-1}\int_0^\rho4r^2e^{-2r}\,\mathrm dr +\int_\rho^\infty4r e^{-2r}\,\mathrm dr, \tag{38}\\ f'(\rho)&=-\rho^{-2}\int_0^\rho4r^2e^{-2r}\,\mathrm dr. \tag{39}\end{align*}\] In particular, \[f(1)=1-2e^{-2}>0.7,\quad f(4)\le\tfrac14,\quad \frac{1-5e^{-2}}{16}\le-f'(\rho)\le1\quad(1\le\rho\le4).\] At the first endpoint the leading shift in (37) is at most \((500-700)m/D\); at the second it is at least \((500-250-100)m/D\). The quadratic and rounding errors fit inside these strict margins, proving bracketing.

For the secant bound, fix a ground function \(u\) at either endpoint, and put \(A(t)=X_i+t e_z\). The kernel derivative is dominated by \(|y-A(t)|^{-2}\), and Hardy gives \[\begin{align*} \int\frac{|u^2-\phi_i^2|}{|y-A(t)|^2}\,\mathrm dy &\le\left\|\frac{u-\phi_i}{|\,\cdot-A(t)\,|}\right\|_2 \left\|\frac{u+\phi_i}{|\,\cdot-A(t)\,|}\right\|_2 \\ &\le4\|\nabla(u-\phi_i)\|_2 \bigl(\|\nabla u\|_2+\|\nabla\phi_i\|_2\bigr) \le C\frac mD. \tag{40}\end{align*}\] The same bound with \(u^2\) justifies absolute integration in \(y\) and \(t\). Thus replacing the density by \(\phi_i^2\) when integrating the moving-pole kernel derivative incurs at most \(Cq'(m/D)(b-a)\) error. The variational principle gives \[\langle\psi_i(b),(h_i(b)-h_i(a))\psi_i(b)\rangle \le\lambda_i(b)-\lambda_i(a) \le\langle\psi_i(a),(h_i(b)-h_i(a))\psi_i(a)\rangle.\] Both bounding quantities differ from \(-q'(f(b)-f(a))\) by the error just estimated, proving (35).

For \(k\ne i\), the moving secondary pole in \(\chi_kV_k\) stays at least \(D/32\) from the region where its cap is nonzero. Its bounded potential changes in sup norm by at most \(Cq'D^{-2}|\rho_k-\widetilde\rho_k|\). Summing and applying the variational principle for a bounded additive perturbation proves (36). ◻

Proof. On \(r_i\ge .01D\), both own poles are far away and \(V_i+W_i=O(m/D)\). For \(a=.99\) and \(b(y)=C\exp[-a(r_i-.01D)]\), direct calculation gives \[\frac{(h_i-\lambda_i)b}{b} =-\frac{a^2}{2}+\frac a{r_i}-V_i-W_i-\lambda_i.\] Its constant part is \((1-a^2)/2+O(m/D)>0\) at the chosen scales. The zero-order coefficient of \(h_i-\lambda_i\) is also positive throughout this exterior region. The uniform sup bound fixes \(C\) so that \(b\ge|\psi_i|\) at the inner boundary. An \(H^2(\mathbb R^3)\) function is uniformly continuous and vanishes at infinity; otherwise uniform continuity on a sequence of disjoint balls would contradict its square integrability. On a large finite annulus add a constant \(\varepsilon>0\) to \(b\) to dominate both signs of \(\psi_i\) at the outer boundary. The supersolution inequality is preserved. Testing the negative part of \(b+\varepsilon\mp\psi_i\) proves comparison. Send the outer radius to infinity and then \(\varepsilon\) to zero. Inside \(.01D\), the same sup bound proves (41).

For the gradient estimate choose a smooth cutoff \(\zeta\) that vanishes on \(r_i\le .645D\), equals one on \(r_i\ge .65D\), and has gradient \(O(D^{-1})\). Testing the eigen-equation with \(\zeta^2\psi_i\) gives \[\frac12\|\nabla(\zeta\psi_i)\|_2^2 =\int(V_i+W_i+\lambda_i)|\zeta\psi_i|^2\,\mathrm dy +\frac12\int|\nabla\zeta|^2|\psi_i|^2\,\mathrm dy.\] All coefficients on this support are bounded. Thus (41) bounds the exterior \(H^1\) norm by a polynomial in \(D\) times \(\exp[-.99(.645-.01)D]\). Since \(.99(.645-.01)=.62865>.61\), the polynomial factor is absorbed for large \(D\). ◻

Proof. For two concentric balls of radii between \(1/2\) and \(1\), separated by \(h\), the local estimates needed below are \[\begin{align*} \|D^2v\|_{2,\mathrm{inner}} &\le C\bigl(\|\Delta v\|_{2,\mathrm{outer}} +h^{-2}\|v\|_{2,\mathrm{outer}}\bigr), \tag{45}\\ \|\nabla v\|_{2,\mathrm{inner}} &\le C\bigl(h\|\Delta v\|_{2,\mathrm{outer}} +h^{-1}\|v\|_{2,\mathrm{outer}}\bigr). \tag{46}\end{align*}\] To obtain them, insert smooth cutoffs with first and second derivatives \(O(h^{-1})\) and \(O(h^{-2})\). Integration by parts against a squared cutoff and Cauchy–Schwarz, with \(ab\le (h a)^2/2+(b/h)^2/2\), bound first derivatives on an intermediate ball as in (46). Apply the whole-space Fourier inequality \(\|D^2w\|_2\le C\|\Delta w\|_2\) to the cutoff function and substitute that first-derivative bound in its commutator terms to obtain (45). Approximation proves these estimates for weak solutions. Successive localized Laplacian regularity makes \(u\) smooth when \(F\) is smooth, justifying the following differentiations.

Fix an integer \(\ell\ge2\), set \(r_j=1-j/(2\ell)\), and write \[T_j=\max_{|\alpha|=j}\|\partial^\alpha u\|_{L^2(B_{r_j})}, \qquad M=\|u\|_{L^2(B_1)}.\] The first-derivative bound gives \(T_1\le C(\ell+A/\ell)M\). For \(j\le\ell-2\), apply (45) to an order-\(j\) derivative between \(B_{r_j}\) and \(B_{r_{j+2}}\). Leibniz’ rule yields \[ T_{j+2}\le C\left( \sum_{b=0}^j\binom jb A^{b+1}b!\,T_{j-b} +\ell^2T_j\right). \tag{47}\] Lower-order derivatives are available on the larger balls \(B_{r_{j-b}}\). The multi-index coefficient sum of order \(b\) is \(\binom jb\), by the multinomial identity. If preceding orders satisfy \(T_s\le M(K\ell)^s\), the sum is at most \[AM(K\ell)^j\sum_{b=0}^j(A/K)^b \le2AM(K\ell)^j\] for \(K\ge2A\), using \(\binom jb b!\le\ell^b\). Increasing the fixed \(K\) absorbs the remaining \(C\ell^2\) factor into \((K\ell)^2\). Both parity chains therefore satisfy \(T_j\le M(K\ell)^j\).

Set \(\ell=k+2\). Sobolev embedding on fixed smaller balls bounds the sup norm of an order-\(k\) derivative on \(B_{1/4}\) by \(CM[K(k+2)]^{k+2}\). The inequality \(k!\ge(k/e)^k\), with the small cases absorbed in the constant, converts this into (44). Taylor’s formula, summed over multi-indices, bounds the degree-\(L\) remainder on \(|z|\le a\) by a constant times \((C_1a)^{L+1}\). The first-derivative remainder has an additional factor polynomial in \(L\) and the same geometric decay. Fixing \(a\) small enough that \(C_1a<1\) proves the last assertion. ◻

Proof. Work in local coordinates \(z=y-X_i\) on \([-D,D]^3\), with poles at \(0\) and \(\rho_i e_z\). Put \(\delta=2^{-8D}\) and fix a sufficiently small rational \(\xi>0\), as specified below. On each coordinate axis let \(B\) be the set of projections of the two poles. Starting from \([-D,D]\), bisect dyadically until every interval \(I\) obeys \[ |I|\le\xi\min\{1,\operatorname{dist}(I,B)+\delta\}. \tag{50}\] Every leaf length is at least \(\xi\delta/2\), since its parent failed the stopping rule. Above a fixed spacing there are \(O(D)\) intervals. At a finer level of length \(h\), only intervals within \(h/\xi\) of one of the at most two projections can still split, so a fixed number of intervals per level are split there. There are \(O(\log(D/\delta))=O(D+\log D)\) levels, hence \(O(D)\) intervals per axis. This remains true for poles arbitrarily close to dyadic boundaries. The tensor product has \(O(D^3)\) nodes and cells.

Let \(b_w\) be the tensor-linear hat at node \(w\). The full hats form a partition of unity in the cube. Retain nodes satisfying \[ |w|<.73D,\qquad |w|>\delta,\qquad |w-\rho_i e_z|>\delta, \tag{51}\] and take \(b_w(z)(z-w)^\alpha\) for all \(|\alpha|\le L_0\), where \(L_0\) is a sufficiently large fixed multiple of \(D\). These functions span \(\mathcal V_i\), whose dimension before removing dependencies is \(O(D^6)\).

Here is the support estimate that avoids any aspect-ratio assumption. For either pole \(a\), every coordinate interval adjacent to \(w\) obeys \[|I|\le\xi\min\{1,|w_j-a_j|+\delta\}.\] The support radius \(R_w^{\mathrm{hat}}\) of a retained hat thus satisfies \[ \begin{split} R_w^{\mathrm{hat}} &\le\xi\bigl(|w-a|+\sqrt3\delta\bigr) \le(1+\sqrt3)\xi|w-a|,\\ R_w^{\mathrm{hat}}&\le\sqrt3\xi. \end{split} \tag{52}\] For \(d_w=\min\{1,|w|,|w-\rho_i e_z|\}\) this gives \(R_w^{\mathrm{hat}}\le C\xi d_w\). Each support avoids both poles and, for large \(D\), lies in \(|z|<.8D\).

Choose a fixed \(c<1/10\) and set \(R=c d_w\). The whole ball \(B(w,R)\) avoids both poles. Throughout it every remote primary has distance at least \(D-2S-.73D-c>D/8\), so it also avoids every remote cap transition. On this enlarged ball the eigen-equation contains only pure Coulomb potentials. Rescaling by \(R\) gives \(\Delta u=Fu\) on the unit ball. Its coefficient satisfies (43) with an absolute \(A\): local inverse-distance expansions give factorial derivative bounds because \(R/\operatorname{dist}(w,a)\le c\); the extra \(R^2\) in the coefficient makes \(R^2/\operatorname{dist}(w,a)\) bounded; and remote contributions sum to at most \(Cm/D\). The eigenvalue term is bounded too. The sup norm of \(u\) is uniform by 8.

Choose \(\xi\) sufficiently small compared to \(c\) and the fixed analytic radius in 11, putting the hat supports inside that smaller ball after rescaling. With \(L_0\) a sufficiently large fixed multiple of \(D\), its Taylor polynomial \(p_w\) satisfies \[ |\psi_i-p_w|+R|\nabla\psi_i-\nabla p_w|\le e^{-20D} \quad\text{on }\operatorname{supp}b_w. \tag{53}\] Increasing the fixed multiple absorbs the polynomial factor in the differentiated remainder. The Taylor coefficients establish existence of an approximant; the algorithm will not compute them.

To remove tiny neighborhoods of the poles and the exterior, let \(a_{\mathrm{out}}\) equal one up to radius \(.66D\), decrease linearly to zero by \(.72D\), and vanish thereafter. Let \(a_{\mathrm{in}}\) vanish up to \(\delta\), increase linearly to one by \(2\delta\), and equal one thereafter. Sample \[a_{\mathrm{out}}(|z|)\, a_{\mathrm{in}}(|z|)\, a_{\mathrm{in}}(|z-\rho_i e_z|)\] at all grid nodes to get weights \(c_w\in[0,1]\), and put \(b=\sum_w c_wb_w\). Nonzero weights occur only at retained nodes. The function \(b\) equals one away from \(O(\delta)\)-balls about the poles and away from \(|z|>.65D\): a node within \(2\delta\) of a pole has adjacent lengths \(O(\xi\delta)\), whereas all intervals have length at most \(\xi\). These facts localize the inner and outer transitions, respectively.

Each first derivative of a tensor-linear interpolant is a convex combination of edge secants in that coordinate. Consequently \(|\nabla b|\le C/\delta\) on the inner transitions and \(|\nabla b|\le C/D\) on the outer transition, irrespective of aspect ratios. The bounded sup norm of \(\psi_i\) makes the inner cutoff-gradient cost \(C\delta^{1/2}\). Also, by its \(H^2\) bound, \[\|\nabla\psi_i\|_{L^2(B_{C\delta}(a))} \le C\delta\|\nabla\psi_i\|_6\le C\delta.\] The exterior cost follows from (42). Thus \[ \|\psi_i-b\psi_i\|_{H^1} \le C\bigl(\sqrt\delta+e^{-.61D}\bigr). \tag{54}\] The cutoff weights, like the Taylor coefficients, are only existence coefficients in the explicitly generated space.

Set \(v=\sum_w c_wb_wp_w\in\mathcal V_i\). At almost every point at most eight hats are nonzero. Their derivatives cost at most \(C/\delta\), by the minimum interval length, and \(R^{-1}\le C/\delta\). On a region of volume \(O(D^3)\), (53) gives \[\|b\psi_i-v\|_{H^1}\le C D^{3/2}\delta^{-1}e^{-20D}.\] Since \(\sqrt\delta=e^{-4(\log2)D}\) and \(\delta^{-1}e^{-20D}=e^{-(20-8\log2)D}\), this and (54) prove (48) for large \(D\).

Finally, the form of \(h_i-\lambda_i\) is bounded in absolute value by a constant times the \(H^1\) norm squared. For \(e=v-\psi_i\), the eigen-equation cancels the cross terms, giving \[\frac{\langle v,h_iv\rangle}{\|v\|_2^2}-\lambda_i =\frac{\langle e,(h_i-\lambda_i)e\rangle}{\|v\|_2^2} \le C\|e\|_{H^1}^2.\] Here \(\|v\|_2\) is bounded away from zero. The variational principle supplies the opposite inequality for the lowest Ritz value, proving (49). ◻

Proof. We divide the construction into exact polynomial arithmetic, certified Coulomb integration, and rational eigenvalue enclosure.

Exact matrices and a mass lower bound. Construct the basis in 12. All grid decisions are rational comparisons: one-dimensional distances in (50) are rational, and node membership in (51) is checked by comparing squared distances. Grid depth is \(O(D+\log D)\). Endpoints, hat coefficients, and coefficients of the degree-\(O(D)\) monomials have polynomial bit length. For instance, raising a rational coordinate of \(b\) bits to degree \(L\) costs at most \(O(Lb)\) bits, before polynomial sums; here \(b=O(B+D+\log D)\) and \(L=O(D)\). Only polynomially many coefficients and rational operations occur in any polynomial expansion or box integral. Summing rational terms by a product of their denominators gives a conservative polynomial bit bound even without exploiting common denominators.

Let \(f_1,\ldots,f_{r_0}\) be this possibly redundant basis, \(r_0=O(D^6)\). The mass and kinetic matrices \[ M_{ab}=\int f_af_b\,\mathrm dy,\qquad T_{ab}=\frac12\int\nabla f_a\cdot\nabla f_b\,\mathrm dy \tag{56}\] are exact rationals obtained by box integration. A product of two hats has only a fixed number of box pieces; alternatively, summing over the entire polynomial-size grid also suffices. Remove redundant functions by exact rational Gram-matrix elimination, retaining an independent subset. This gives a positive definite matrix \(M\) of dimension \(r\le r_0\), with corresponding kinetic matrix \(T\). The subset is nonempty by the approximation property.

Let \(B_0\) be the largest bit length of a numerator or positive denominator of the entries of this reduced \(M\), enlarged to at least one. It is computed from the exact entries and is polynomial in \(D,m,B\). If \(Q\) is the product of all entry denominators, then \(Q\le2^{r^2B_0}\) and \(QM\) is an integer positive definite matrix. Therefore \(\det(QM)\ge1\), while \(\lambda_{\max}(M)\le r2^{B_0}\). The explicitly computable number \[ \mu=2^{-K_M},\qquad K_M=r^3B_0+(r-1)\bigl(B_0+\lceil\log_2r\rceil\bigr) \tag{57}\] satisfies \(\lambda_{\min}(M)\ge\mu\). This intentionally coarse determinant bound is enough: \(K_M\) is polynomial. Exact rank selection is itself polynomial-bit arithmetic. After clearing denominators, fraction-free elimination has entries given by minors of an integer matrix. Hadamard’s determinant bound bounds their bit lengths by a polynomial in its dimension and entry bit length. Equivalently, ordinary rational elimination with fractions reduced at each step has the same polynomial bound through the bordered-minor formulas.

Certified Coulomb entries. The preceding determinant bound shows how many extra bits compensate for a poorly conditioned basis. It remains to compute each potential entry to that absolute accuracy. The graded supports avoid every pole, which permits a uniformly convergent kernel expansion on each box. On the supports of the basis functions every other cap equals one, so every potential entry is a finite sum of pure Coulomb integrals against polynomials on rational boxes. For each such box \(Q_0\), each pole \(a\), and its midpoint \(c_0\), put \[s=|c_0-a|^2,\qquad z(y)=\frac{2(c_0-a)\cdot(y-c_0)+|y-c_0|^2}{s}.\] Every contributing box is in the support of a retained hat. By (52), its diameter is a fixed small fraction of its distance from either own pole. Remote poles are farther still. Decreasing the fixed \(\xi\) if necessary gives \(|z(y)|\le1/2\) throughout every such box. Thus \[ \frac1{|y-a|} =s^{-1/2}\sum_{j=0}^J\binom{-1/2}{j}z(y)^j +\mathcal R_J(y),\qquad |\mathcal R_J(y)|\le s^{-1/2}2^{-J}. \tag{58}\] The coefficient bound \(\left|\binom{-1/2}{j}\right| =4^{-j}\binom{2j}{j}\le1\) proves this geometric remainder bound by summing the tail.

Here \(s\) and all coefficients of \(z\) are exact rationals. Both \(\log_2^+(s^{-1/2})\) and the bit lengths of these rationals are polynomial in \(D,m,B\); indeed the own-pole distances on supports are bounded below by a fixed multiple of \(\delta\). The positive scalar \(s^{-1/2}\) can be approximated to any requested polynomial number of bits by rational bisection and exact squaring.

We spell out an absolute-error certificate, so cancellation among polynomial terms causes no difficulty. If \(P\) is the polynomial multiplying a kernel on \(Q_0\), a rational coefficient-and-volume bound gives \[A_{Q_0}\ge\int_{Q_0}|P(y)|\,\mathrm dy.\] One may take the volume times the sum of absolute coefficients, each multiplied by the appropriate coordinate bound to its degree. Its bit length is polynomial. Include the charge in this bound. If \(s^{-1/2}\le2^b\), the truncation error is at most \(A_{Q_0}2^{b-J}\). An approximation to \(s^{-1/2}\) of absolute error \(2^{-q_0}\) adds at most \(2A_{Q_0}2^{-q_0}\), since the absolute sum in (58) is at most two. For \(T_0\) box–pole summands, taking \[\begin{align*} J&\ge p+b+\lceil\log_2^+ A_{Q_0}\rceil +\lceil\log_2 T_0\rceil+3,\\ q_0&\ge p+\lceil\log_2^+ A_{Q_0}\rceil +\lceil\log_2 T_0\rceil+4 \end{align*}\] makes their total error smaller than \(2^{-p}\). Use maxima over the summands to choose common orders. Here \(\log_2^+t=\max\{0,\log_2t\}\); integer upper bounds can be read directly from rational numerator and denominator bit lengths.

The polynomial \(z(y)^j\) has degree \(2j\) in three variables. Even after multiplication by the basis polynomials it has only \(O((D+J+1)^3)\) possible monomials. Polynomial multiplication, integration, and coefficient bit lengths are polynomial in \(D,m,B,J,q_0\). These bounds establish a deterministic entry algorithm with any prescribed polynomial number \(p\) of absolute precision bits. No matrix element is computed by integrating across a Coulomb singularity or a cap transition.

Matrix error and exact eigenvalue decisions. Let \(A\) be the exact form matrix of \(h_i\) in the independent basis. Compute a symmetric rational matrix \(\widetilde A\) by evaluating each potential entry, once per unordered index pair, to error at most \[\tau=\frac{\varepsilon_D\mu}{16r}.\] The kinetic matrix is exact. Since a symmetric matrix’s operator norm is at most \(r\) times its largest absolute entry, \[ \|M^{-1/2}(\widetilde A-A)M^{-1/2}\| \le\varepsilon_D/16. \tag{59}\] The required entry precision has \(O(D+K_M+\log r)\) bits and is therefore polynomial. Let \(\widetilde\lambda\) be the minimum of the generalized Rayleigh quotient for \((\widetilde A,M)\). The variational characterization and (59) give \[|\widetilde\lambda-\lambda_i^{\mathrm R}| \le\varepsilon_D/16.\]

An initial rational bracket is \([-R_0,R_0]\), where \[R_0=1+\frac{r\max_{a,b}|\widetilde A_{ab}|}{\mu}.\] Its bit length is polynomial. To bisect it at a rational \(z\), test the inertia of \(\widetilde A-zM\) exactly. For clarity, an elementary exact inertia algorithm repeatedly uses congruence and Schur complements. A nonzero diagonal entry is a one-dimensional pivot and contributes its sign. If every remaining diagonal vanishes but an off-diagonal entry \(a\ne0\) remains, pivot on the block \(\left(\begin{smallmatrix}0&a\\a&0\end{smallmatrix}\right)\), which contributes one positive and one negative direction. A zero remaining matrix contributes only zero directions. These operations give the full inertia. After denominators are cleared, each Schur entry is a ratio of bordered minors; determinant bounds again give polynomial bit growth. Thus this test is a polynomial-time exact rational calculation, including the case of zero eigenvalues.

Because \(M\) is positive definite, \(\widetilde A-zM\) is positive semidefinite exactly when \(z\le\widetilde\lambda\). Update the lower or upper endpoint accordingly. Polynomially many bisections give a rational enclosure \([\ell,u]\) of \(\widetilde\lambda\) of width at most \(\varepsilon_D/8\). Arbitrary proximity to a bisection point incurs no further precision requirement and no eigenvalue-gap assumption.

The fixed lower bound on \(D\) can be enlarged so that \(Ce^{-1.16D}\le\varepsilon_D/8\); this is possible since \((8/5)\log2<1.16\). Combining the one-sided Ritz bias and the two-sided matrix error gives the certified interval \[ \lambda_i\in \left[\ell-\frac{\varepsilon_D}{16} -\frac{\varepsilon_D}{8}, u+\frac{\varepsilon_D}{16}\right]. \tag{60}\] Its width is at most \(3\varepsilon_D/8\). Notice that the Ritz bias is subtracted from the lower endpoint; a Ritz upper bound alone would not suffice for tuning. This proves the proposition.

All dimensions, operation counts, and working bit lengths are polynomial in \(D,m,B\). The mesh and the independent subset may change at different rational samples, but every such sample obeys the same estimates. Constants controlling the polynomial degree, the cutoff size, and the minimum admissible \(D\) are fixed once from the uniform inequalities above. They are not input-dependent advice or a search for an unknown spectral gap. ◻

Proof. For frozen parameters other than \(\rho_i\), the endpoint brackets, continuity, and strict secant bound in 9 give a unique root of \(\lambda_i=-1/2\) in \([1,4]\). Let \(T_i(\rho)\) be that root and let \(T=(T_1,\ldots,T_m)\), using the old other coordinates in every update. Combining (35) and (36) gives \[ \|T(\rho)-T(\widetilde\rho)\|_\infty \le\theta\|\rho-\widetilde\rho\|_\infty,\qquad \theta=Cm/D^2\le1/2. \tag{62}\] Thus \(T\) has a unique fixed point \(\rho^*\) in the cube. This also follows directly by iterating the map: its successive differences are geometrically summable, and the limit is fixed by continuity.

We implement each coordinate update with a certified scalar bisection. Start with \([1,4]\). At its midpoint invoke 13. If the whole resulting eigenvalue interval lies below \(-1/2\), retain the upper half of the parameter interval; if it lies above \(-1/2\), retain the lower half. If it contains \(-1/2\), stop and return that sample. In this last case the true energy differs from \(-1/2\) by at most \(\varepsilon_D\), so (35) puts the sample within \(C\varepsilon_D/q'\) of the exact root. Otherwise, after \(J=\lceil4D\rceil\) bisections return the midpoint of the remaining parameter interval. All sign decisions were certified, so that interval still contains the root and has length \(3\,2^{-J}\). Since \(q'\le1\) and \(2^{-4D}\ll\varepsilon_D\), the same update-error bound holds: \[ \|\widehat T(\rho)-T(\rho)\|_\infty \le C\varepsilon_D/q'. \tag{63}\] This comparison-or-stop rule never requests an exact sign of an unknown irrational eigenvalue at a putative zero.

Start with all coordinates equal to two and perform \(J\) approximate simultaneous updates. If \(e_j\) is the sup-norm distance to \(\rho^*\), then \[e_{j+1}\le\tfrac12 e_j+C\varepsilon_D/q', \qquad e_J\le3\,2^{-J}+2C\varepsilon_D/q'.\] The own-coordinate upper slope and the other-coordinate variation bound together give \[|\lambda_i(\rho)-\lambda_i(\rho^*)| \le Cq'\|\rho-\rho^*\|_\infty.\] Consequently the output has \[|\lambda_i+1/2| \le Cq'2^{-J}+C\varepsilon_D \le e^{-1.01D}\] for large enough \(D\), since \((8/5)\log2>1.1>1.01\).

Each scalar bisection restarts from the rational interval \([1,4]\) and makes at most \(O(D)\) midpoint operations. Its samples, including early-return samples, therefore have \(O(D)\) fractional bits independently of the number of simultaneous iterations. There are \(O(mD^2)\) oracle calls. Their data bit lengths consist of the fixed-center and rounded-charge bit lengths together with these \(O(D)\) parameter bits, so 13 proves the total polynomial-time bound. In particular the charge-rounding error is part of every evaluated operator and is compensated by tuning; it is not introduced as an uncontrolled perturbation afterward. ◻

Proof. Fix a site, suppress its subscript, and write \(r=|y-X|\) and \(\omega=(y-X)/r\) when \(r>0\). The capped-field estimates in Proposition 8 give \[ \|W\|_\infty\le \frac{Cm}{D},\qquad \|\nabla W\|_\infty\le \frac{Cm}{D^2},\qquad \|D^2W\|_\infty\le \frac{Cm}{D^3}. \tag{67}\] Set \(A(y)=\int_0^1W(X+t(y-X))\,\,\mathrm dt\), so \(F=rA\). For \(0<r\le3D\), differentiation gives \[\begin{align*} \nabla F &=\omega A+r\int_0^1t\nabla W(X+t(y-X))\,\,\mathrm dt,\\ \Delta F &=\frac{2A}{r} +2\omega\cdot\int_0^1t\nabla W(X+t(y-X))\,\,\mathrm dt\\ &\quad+r\int_0^1t^2\Delta W(X+t(y-X))\,\,\mathrm dt. \end{align*}\] Consequently, \[ |\nabla F|\le \frac{Cm}{D},\qquad |\Delta F|\le C\left(\frac{m}{Dr}+\frac{m}{D^2}\right). \tag{68}\] Also \(F(r,\omega)=\int_0^rW(X+t\omega)\,\,\mathrm dt\), and hence \(\partial_rF=W(y)\) exactly.

Let \(f=\pi^{-1/2}e^{-r+F}\). A direct calculation on the annulus \(10\le r\le3D\) gives \[\begin{align*} \frac{(h_i-\lambda_i)f}{f} &=-(\lambda_i+1/2)+\frac1r-V_i -\frac12\bigl(\Delta F+|\nabla F|^2\bigr), \tag{69}\\ \left|\frac{(h_i-\lambda_i)f}{f}\right| &\le C\left(\frac{m}{Dr}+\frac{m^2}{D^2}\right). \tag{70}\end{align*}\] In the first equality the term \(\partial_rF\) cancels \(W_i\). For the second, note that \(\rho_i\le4\) and \(r\ge10\), so \[\left|\frac1r-V_i(y)\right| \le \frac{|1-q|}{r}+\frac{q'}{r-4} \le \frac{Cm}{Dr}.\] The charge-rounding errors and \(e^{-1.01D}\) are absorbed by the right-hand side of (70).

Choose a sufficiently large fixed constant \(A_0\) and put \[ g(r)=A_0\left[ \frac{m}{D}\bigl(1+\log(r/10)\bigr) +\frac{m^2r}{D^2}\right]. \tag{71}\] Then \[g'=A_0\left(\frac{m}{Dr}+\frac{m^2}{D^2}\right),\qquad \Delta g=A_0\left(\frac{m}{Dr^2}+\frac{2m^2}{D^2r}\right), \qquad 0\le\Delta g\le\frac{2g'}r.\] Replacing \(f\) by \(fe^{\pm g}\) adds to its operator quotient the term \[ \pm(1-\partial_rF)g' \mp\frac12\Delta g-\frac12(g')^2. \tag{72}\] Since \(r\ge10\), \(\partial_rF=O(m/D)\), and \(g'\) is uniformly small for sufficiently large \(D\), the added term is at least \(c g'\) for the plus sign and at most \(-c g'\) for the minus sign. Increasing \(A_0\) dominates (70). Thus \(fe^g\) is a supersolution and \(fe^{-g}\) a subsolution.

At \(r=10\), Proposition 8 and Sobolev embedding give \(\psi_i=\pi^{-1/2}e^{-10}+O(m/D)\), while \(F=O(m/D)\). The reference value \(e^{-10}\) is a fixed positive number. A further fixed increase of \(A_0\) therefore ensures \(fe^{-g}\le\psi_i\le fe^g\) on the inner boundary. On the outer boundary \(r=3D\), Lemma 10 gives \[|\psi_i|\le C e^{-0.99(3-0.01)D}=C e^{-2.9601D}.\] Moreover \(0\le F\le Cm\) there. With \(\kappa=e^{-2.5D}\), the functions \[f_+=fe^g+\kappa,\qquad f_-=fe^{-g}-\kappa\] satisfy \(f_-\le\psi_i\le f_+\) on both boundary spheres for large \(D\). On the annulus the zero-order coefficient obeys \[-V_i-W_i-\lambda_i \ge \frac12-\frac1{10}-O(m/D)\ge\frac14.\] Adding the positive constant \(\kappa\) preserves the supersolution inequality, and subtracting it preserves the subsolution inequality. The maximum principle therefore yields the same inequalities throughout the annulus. For completeness, its weak form here follows by testing a supersolution-minus-solution with its nonnegative negative part: the resulting identity is the negative of a sum of a gradient norm and a strictly positive weighted \(L^2\) norm, so the negative part vanishes. The same argument applies to a solution-minus-subsolution.

On \(10\le r\le2D\), one has \(g(r)=O(\varepsilon_{\mathrm{rel}})\) and, since \(F\ge0\), \[\frac{\kappa}{f}\le \sqrt\pi\,e^{-D/2}.\] The comparison inequalities imply (66) on this region. On \(r\le10\), the reference exponential is bounded below by a fixed positive constant, \(F=O(m/D)\), and the uniform absolute orbital estimate gives the same relative conclusion directly. ◻

Proof. For large \(D\), both poles are strictly inside \(|z|<D/16\), where the kernels vanish. On the nonzero part of either kernel, \(|z|\ge D/16\) and \(|z-\rho e_z|\ge |z|/2\). Rescale by \(z=D\zeta\). Both kernels are \(D^{-1}\) times a function whose derivatives through order two are uniformly bounded: the shifted pole in the second rescaled kernel has position \((\rho/D)e_z\), a fixed positive distance from the support of \(\chi\). The \(C^2\) matching of \(\chi\) at its transition endpoints makes these global bounds, including at the hole boundaries. This proves (75).

Since \(|q-1|+q'\le Cm/D\) after rounding, \[W_i-W_i^*= \sum_{k\ne i}\bigl[(q-1)G_D(y-X_k) +q'J_{D,\rho_k}(y-X_k)\bigr]\] has sup norm \(O(m^2/D^2)\). Integration over a ray of length at most \(3D\) proves the second bound in (76). Exponentiating this bound in Proposition 15 proves the final assertion. None of the kernel estimates require tuning; they hold uniformly on the whole parameter cube. ◻

Proof. On the support of \(V_j^{\mathrm{hole}}\), write \(y=X_j+v\). Then \(|v|\le D/8\) and \(|y-X_i|\le D+2S+D/8<2D\). Proposition 15 and Lemma 16 apply to both wave functions. Their product is \(\pi^{-1}e^{L_X(v)}\) times a uniform \(1+O(\varepsilon_{\mathrm{rel}})\) factor. The primary kernel has coefficient \(q=1+O(m/D)\), so it remains to compare the secondary kernel with the primary one in the resulting positive integral.

For \(|v|>10\), \[\frac{|v|}{|v-\rho_j e_z|} \le\frac{|v|}{|v|-4}\le\frac53.\] For \(|v|\le10\), the hole weight is one when \(D\) is large. The ray derivative estimates, also valid for \(F^*\), show that \(L_X\) is Lipschitz on this ball with a fixed constant. This statement requires no differentiability of \(|v|\) at zero; it follows either by integrating the bounded derivatives away from zero or directly from the ray formula. Thus \(e^{L_X(v)}\) is within fixed factors of \(e^{L_X(0)}\) throughout the ball. Moreover, \[\int_{|v|\le10}\frac{\,\mathrm d^3v}{|v-\rho_j e_z|} \le\int_{|w|\le14}\frac{\,\mathrm d^3w}{|w|}<\infty, \qquad \int_{|v|\le1}\frac{\,\mathrm d^3v}{|v|}>0.\] The secondary-kernel integral on the inner ball is therefore at most a fixed multiple of the primary-kernel integral. Combining the two regions and multiplying by \(q'=O(m/D)\) gives a relative \(O(m/D)\) correction. Positivity allows the uniform relative wave-function errors to be taken outside both integrals. These errors are included in (79). ◻

Proof. Put \(u=u_{ij}=p_j-p_i\), so \(|u|=1\). The triangle inequality gives \(|Du+v|+|v|\ge D\). Each ray correction in the integral is nonnegative and at most \(Cm\), by (75) and its ray length. The integral of \(|v|^{-1}\) over the ball of radius \(D/8\) is \(O(D^2)\), proving the upper bound. On \(|v|\le1\) the hole weight is one, \(|Du+v|+|v|\le D+2\), and the corrections are nonnegative. Integration over this fixed ball gives the lower bound; for example \(2e^{-D-2}\) is a valid lower bound for large \(D\). ◻

Proof. Compare the two integrals at the same \(v\). In the ray defining \(F_i^*(X_j+v)\), the reference path point \(Dp_i+t(Dp_j+v-Dp_i)\) moves by \((1-t)s_i+ts_j\), of norm at most \(S\). Each field center moves by at most \(S\) as well. The globally Lipschitz bound in (75) therefore changes the field along the path by at most \(CmS/D^2\). The ray has length \(O(D)\) and its length changes by at most \(2S\). Since the field itself is \(O(m/D)\), \[ |F_i^*(X_j+v;X)-F_i^*(Dp_j+v;X^0)| \le\frac{CmS}{D}. \tag{83}\] The same bound holds for the \(F_j^*\) term; its ray length \(|v|\) does not change. The global nature of the capped-kernel bound is essential here: it also applies when a moving ray crosses a hole or a transition interface.

Write \(u=u_{ij}\), \(\delta=s_j-s_i\), and \(v=zu+w\) with \(w\perp u\). Since \(|v|\le D/8\), \(|Du+v|\ge7D/8\). The Hessian of the norm on the segment from \(Du+v\) to \(Du+v+\delta\) has norm \(O(1/D)\). Its gradient at \(Du+v\) differs from \(u\) by \(O(|w|/D)\). Taylor’s formula gives \[ |Du+v+\delta|-|Du+v| =u\cdot\delta+O\left(\frac{S|w|}{D}+\frac{S^2}{D}\right). \tag{84}\] This estimate includes the whole backward segment \(v=-tu\), \(0\le t\le D/8\): there \(Du+v=(D-t)u\) and the gradient is exactly \(u\). Thus the segment’s non-negligible longitudinal extent creates neither a different coefficient nor an \(O(S)\) transverse error.

On the region \(|w|\le D^{3/4}\), Equations (83) and (84) show that \[ L_X(v)-L_{X^0}(v)+u\cdot\delta =O\left(\frac{mS}{D}+\frac{S^2}{D}+S D^{-1/4}\right). \tag{85}\] The nonnegative integration weight is identical in the two integrals, so exponentiation controls their contributions on this region by the same relative factor.

On the remaining region \(|w|>D^{3/4}\), one has \[ \begin{split} |Du+v|+|v|-D &\ge z+|v| =\frac{|w|^2}{|v|-z}\\ &\ge\frac{4|w|^2}{D} \ge4\sqrt D. \end{split} \tag{86}\] The identity is used only for \(w\ne0\), and the denominator is at most \(2|v|\le D/4\). The field corrections are bounded by \(Cm\), without any angular constancy assumption. Displacement changes a distance by at most \(2S\), and \(|u\cdot\delta|\le2S\). It follows that the omitted reference and displaced contributions, after division by \(e^{-D}\) and \(e^{-D-u\cdot\delta}\) respectively, are at most \[ C D^2\exp\bigl[-c\sqrt D+C(m+S)\bigr]. \tag{87}\] By (80), these are also bounds on their contributions relative to the reference integral, with the appropriate displacement factor. The assumption \(m+S\le D^{1/4}\) and sufficiently large \(D\) absorb both the possibly large anisotropic factor \(e^{Cm}\) and the polynomial volume factor into \(e^{-c_1\sqrt D}\). Combining this with (85), and using smallness of its right-hand side, proves the lemma. ◻

Proof. Let \(L=L_{X^0}\) in (78). On the integration ball, the two distance terms have Lipschitz constants one, and the ray formulas and (75) give a Lipschitz constant \(Cm/D\) for each correction. Hence \[ |L(v)-L(v')|\le C|v-v'|, \qquad |L(v)|\le C(D+m). \tag{89}\] Only the exponential factor has been included in this logarithm. No lower bound on \(1-\chi\) or differentiability of \(\log(1-\chi)\) is needed.

We describe a mesh with exact positive weights. Parameterize the unit sphere, up to a set of area zero, by \[\omega(\mu,\varphi) =\bigl(\sqrt{1-\mu^2}\cos\varphi, \sqrt{1-\mu^2}\sin\varphi,\mu\bigr).\] Here \(-1\le\mu\le1\), \(0\le\varphi\le2\pi\), and \(\,\mathrm d\omega=\,\mathrm d\mu\,\,\mathrm d\varphi\). Use \(M\) equal intervals in \(\mu\), \(N\) equal intervals in \(\varphi\), and a rational radial partition with mesh at most \(h\). Insert all the finitely many breakpoints of \(\chi(r/D)\) into that partition. On a radial interval \(I_a\), the weight \[w_a=\int_{I_a}(1-\chi(r/D))r\,\,\mathrm dr\] is a nonnegative exact rational, computable by polynomial integration. Discard intervals with zero weight. The normalized angular weight of every angular rectangle is exactly \[ \frac1\pi\,\frac2M\,\frac{2\pi}{N}=\frac4{MN}. \tag{90}\] Choose one sample \(v_{abc}\) in each resulting spatial cell, for example at the parameter midpoints.

The square-root parameterization satisfies \[|\omega(\mu,\varphi)-\omega(\mu',\varphi')| \le C|\mu-\mu'|^{1/2}+C|\varphi-\varphi'|.\] Therefore \(M=O((D/h)^2)\) and \(N=O(D/h)\), with sufficiently large fixed constants, make every spatial cell have diameter \(O(h)\), including cells adjacent to the poles. There are \(O((D/h)^4)\) cells, a polynomial number. Their exact weighted sample sum is \[Q=\sum_{a,b,c}\frac{4w_a}{MN}\,e^{L(v_{abc})}.\] Positivity and (89) imply \[ e^{-Ch}Q\le T_{ij}(X^0)\le e^{Ch}Q. \tag{91}\] Thus a fixed multiple of \(\epsilon\) is sufficient for \(h\). This relative estimate is independent of the size \(e^{-D}\) of the integral.

It remains to specify evaluation and bit accuracy. For a ray of length \(r\le2D\), the integrand \(t\mapsto W_i^*(X_i+t(y-X_i))\) has Lipschitz constant at most \(Cmr/D^2\). A midpoint rule with \(K\) equal panels therefore approximates its ray primitive with error at most \[ C\frac{mr^2}{D^2K}\le\frac{Cm}{K}. \tag{92}\] Taking \(K=O(m/\epsilon)\) and evaluating each field sample to absolute error \(O(\epsilon/D)\) gives an \(O(\epsilon)\) error in each primitive. Each field sample is a sum of \(m-1\) globally Lipschitz capped kernels. At small radius the kernel is evaluated as zero, and elsewhere its denominator is bounded below by a fixed multiple of \(D\); equivalently one evaluates the fixed continuous function \(\chi(a)/a\), extended by zero on \([0,1/16]\), at \(a=|z|/D\). Its global Lipschitz bound controls errors even if the approximate radius lies on the other side of a transition breakpoint. No crossing location needs to be found.

The coordinates \(Dp_i\), sample coordinates, square roots, and trigonometric values can be computed to inverse-polynomial absolute accuracy with polynomially many bits. The global kernel and ray bounds just proved control the resulting errors. Thus all sample logarithms can be approximated to an arbitrarily specified \(O(\epsilon)\) absolute error with polynomial work. Such logarithm errors change the positive sample sum by a factor \(e^{O(\epsilon)}\), not by an additive error multiplied by an uncontrolled exponential.

For the remaining arithmetic, compute the scaled sum \(e^D Q\), using sample exponentials \(e^{D+L(v_{abc})}\). Its exact value is bounded below by a positive constant up to the factor in (91), because (80) gives \(e^D T_{ij}(X^0)\ge2e^{-2}\). The sum of its non-exponential weights is \(O(D^2)\), and all exponent arguments have magnitude \(O(D+m)\). After the logarithms have been approximated as above, each elementary exponential can be rounded to sufficiently small absolute error \(O(\epsilon/D^2)\), clipping a negative rounded value to zero; this adds \(O(\epsilon)\) to the scaled sum. These last approximations require only \(O(D+m+\log(D/\epsilon))\) precision bits, with polynomial working precision. This step does not require exponentially accurate ray quadrature: the earlier ray errors have already been bounded multiplicatively. Standard range reduction and Taylor series compute the square roots, trigonometric values, exponentials, and logarithms to these precisions in polynomial time. The exact rational weights have polynomial bit lengths, and the sum contains polynomially many terms.

Allocate a sufficiently small fixed fraction of \(\epsilon\) to each of the spatial, ray, coordinate, and arithmetic errors. The resulting positive rational sum is \(\widehat I_{ij}\) and satisfies (88). ◻

Proof. The single-well operators have common domain \(H^2\), and \[h_i-h_j=V_j^{\mathrm{hole}}-V_i^{\mathrm{hole}}.\] Pair this identity between the real eigenfunctions \(\psi_j\) and \(\psi_i\), using self-adjointness of \(h_j\), to obtain (93). All the pairings are well-defined by the Coulomb domain bounds of Proposition 8. Equivalently, in the two expansions of the full matrix element \(\langle\psi_j,(h+1/2)\psi_i\rangle\), contributions from every hole other than \(i\) and \(j\) cancel identically, because they integrate the same real product against the same multiplication potential.

To record an explicit overlap bound, Lemma 10 gives \[|\psi_i(y)\psi_j(y)| \le C e^{0.0198D}e^{-0.99(r_i+r_j)}.\] Use \(r_i+r_j\ge D-2S\) on a factor with exponent \(0.90\), and retain \(e^{-0.09r_i}\) for integrability. It follows that \[|\langle\psi_j,\psi_i\rangle| \le C e^{-0.8802D+1.8S} \int_{\mathbb R^3}e^{-0.09|y-X_i|}\,\,\mathrm d^3y \le e^{-0.85D}\] for sufficiently large \(D\). Finally \(|\lambda_i-\lambda_j|\le2e^{-1.01D}\) proves the second bound. ◻

Proof. Compute, to additive error at most \(\tau/4\), the quantities \[ d_{ij}=\log\left(\frac{Z_*T_{ij}(X^0)}{t_{ij}}\right), \qquad X^0=(Dp_i)_i. \tag{98}\] Lemma 20 supplies the requisite relative integral accuracy. One can evaluate the logarithm as \[\log Z_*-D+\log\bigl(e^D T_{ij}(X^0)\bigr)-\log t_{ij},\] so all elementary-function calculations have polynomial bit cost. The two-sided bound (80) yields \[ \max_{(i,j)\in\mathcal C}|d_{ij}| \le C\bigl(m+\log D+\log(1/t_{\min})+1\bigr). \tag{99}\]

Apply the constructive linear right inverse in Lemma 6 to the approximated data. Its operator bound is \(Cm\) from the maximum edge datum to the maximum site displacement. Its formulas use rational operations and fixed elements of \(\mathbb Q(\sqrt3)\), so they can be evaluated with polynomially many bits. Round the resulting positions \(Dp_i+s_i\) to rational positions with Euclidean site error at most \(\tau/8\), allocating another such allowance if approximate arithmetic is used in the right inverse. Define the final displacements from these rational positions by \(s_i=X_i-Dp_i\). This definition preserves the exact contact vector decomposition even though \(s_i\) need not be rational. Since every \(u_{ij}\) is a unit vector, the resulting equations satisfy \[ |u_{ij}\cdot(s_j-s_i)-d_{ij}|\le\tau. \tag{100}\] The right-inverse bound and (99) imply (96). For example, the \(\sqrt3\) arithmetic may be done to absolute accuracy \(O(\tau/[D+\operatorname{poly}(m)\max(1,\max|d_{ij}|)])\); the required number of bits is polynomial. Noncontacts retain separation at least \(\sqrt2D-2S\) by the triangle inequality.

There is no requirement to realize prescribed distances exactly. For a contact with \(\delta=s_j-s_i\) the actual distance is \[|Du_{ij}+\delta| =D+u_{ij}\cdot\delta+O(S^2/D).\] Actual coordinates automatically satisfy all vector-cycle identities; the displayed per-edge correction is already included in Lemma 19 and is not accumulated through a further nonlinear solve.

Put \(e_{ij}=u_{ij}\cdot(s_j-s_i)-d_{ij}\). The key cancellation is the exact identity \[ Z_*T_{ij}(X^0)e^{-u_{ij}\cdot(s_j-s_i)} =t_{ij}e^{-e_{ij}},\qquad |e_{ij}|\le\tau. \tag{101}\] Lemma 19 therefore gives \(Z_*T_{ij}(X)=t_{ij}(1+O(\varepsilon_{\mathrm{disp}}+\tau))\). After tuning, Lemma 17 adds only \(O(\varepsilon_{\mathrm{rel}})\) relative error. This proves (97). In particular no factor \(e^{Cm}\) or \(e^{|d_{ij}|}\) multiplies its absolute error, even for the smallest floored hopping. Since \(t_{ij}\le1\), the absolute error per contact is at most \(C\mathcal E_{\mathrm{cal}}\); summing hopping operators incurs only polynomial factors in the number of sites.

All computations of \(X\) used the surrogate fields \(W_i^*\). They depend on the primary centers and on \(D\), but not on any eigenfunction or on the secondary parameters \(\rho_i\). Thus the algorithm first fixes the rounded charges, computes and rounds these centers, and only then invokes Proposition 14 at the final rational centers. The surrogate-field replacement bounds hold uniformly over the tuning cube, and the relative orbital comparisons hold uniformly over all tuned outputs. Thus tuning at these fixed centers requires no recalibration and introduces no subsequent position-rounding error.

Finally Lemma 21 and \(Z_*\le2e^D\) give \(Z_*|B_{ij}-B_{ji}|\le4e^{-0.86D}\). The common target energy \(-1/2\) is sufficient; exact equality of the \(\lambda_i\) is not required. The full one-body operator is Hermitian, so one may use either orientation for its off-diagonal entry after the negligible other-hole terms are accounted for in Section 5. ◻

Proof. Since \(h=h_i-\sum_{k\ne i}V_k^{\mathrm{hole}}\), \[ A\psi_i=\epsilon_i\psi_i- \sum_{k\ne i}V_k^{\mathrm{hole}}\psi_i. \tag{110}\] The \(k\)-th hole is supported in \(r_k\le D/8\). On it, \[r_i\ge 7D/8-2S\quad(k\ne i).\] Both Coulomb poles of \(V_k\) lie in this hole for large \(D\). The charges \(q,q'\) are bounded, and direct radial integration, enlarging the integration ball by at most four for the secondary pole, gives \[ \|V_k^{\mathrm{hole}}\|_2\le C D^{1/2}, \qquad \|V_k^{\mathrm{hole}}\|_1\le C D^2. \tag{111}\] Consequently the sum in (110) has norm at most \[CmD^{1/2}\exp(-.85635D+1.98S).\] Together with the tuning error, this is at most \(e^{-.75D}\) for the stated regime and sufficiently large \(D\).

For \(i\ne j\), put \(\ell_{ij}=|X_i-X_j|\). The product of the two tail bounds is at most \[C\exp\bigl(-.99(r_i+r_j)+.0198D\bigr).\] Use \(r_i+r_j\ge \ell_{ij}\) on a factor with exponent \(.90\); the remaining factor is integrable uniformly in the centers, since \(\int e^{-.09(r_i+r_j)}\,\,\mathrm dy\le \int e^{-.09r_i}\,\,\mathrm dy\le C\). Thus \[|G_{ij}|\le C e^{-.90\ell_{ij}+.0198D} \le C e^{-.8802D+1.8S}.\] The diagonal entries are one. The row-sum estimate and \(m+S\le D^{1/4}\) prove the second inequality in (107).

Let \(M=\mathcal U^*A\mathcal U\), which is real symmetric. We next estimate its entries before orthonormalization. If \(k\ne i,j\), both \(r_i\) and \(r_j\) are at least \(7D/8-2S\) on the \(k\)-th hole. Equations (103) and (111) therefore imply \[\sum_{k\ne i,j} |\langle\psi_j,V_k^{\mathrm{hole}}\psi_i\rangle| \le CmD^2e^{-1.7127D+3.96S}.\] The same estimate bounds \(\sum_{k\ne i}\langle\psi_i,V_k^{\mathrm{hole}}\psi_i\rangle\). For a noncontact pair, the remaining hole at \(j\) is bounded using \(r_i+r_j\ge\sqrt2D-2S\), giving \[|\langle\psi_j,V_j^{\mathrm{hole}}\psi_i\rangle| \le CD^2 e^{-(.99\sqrt2-.0198)D+1.98S}.\] Here \(.99\sqrt2-.0198>1.38\). All these bounds are smaller than \(e^{-1.1D}\), after absorbing their polynomial factors and fixed constants. For \(i\ne j\), the term \(\epsilon_iG_{ji}\) is at most \(e^{-1.86D}\). It follows, adjusting a fixed constant, that \[\begin{align*} |M_{ii}-\epsilon_i|&\le Ce^{-1.1D},\\ |M_{ji}+B_{ij}|&\le Ce^{-1.1D} &&(i\to j\text{ a chosen contact}),\tag{112}\\ |M_{ij}|&\le Ce^{-1.1D} &&(i\ne j,\ \{i,j\}\notin\mathcal C). \end{align*}\] The exactly symmetric matrix \(M\) requires only one of the two orientations of each contact. Its diagonal tuning contribution has norm at most \(e^{-1.01D}\), not \(m e^{-1.01D}\). Entrywise comparison with \(\mathsf T\) and a row-sum bound now give \[\|Z_*M-\mathsf T\| \le m\epsilon_{\mathrm{hop}} +Z_*e^{-1.01D}+CmZ_*e^{-1.1D}.\]

For completeness, the residual estimate controls the entire matrix, without using a dimension-dependent bound on a many-body space: \[\|A\mathcal U\|\le\sqrt m\,e^{-.75D},\qquad \|\mathcal U\|=\|G\|^{1/2}\le2,\qquad \|M\|\le C\sqrt m\,e^{-.75D}.\] Functional calculus on the positive matrix \(G\) gives \(\|G^{-1/2}-I\|\le Ce^{-.85D}\) and \(\|G^{-1/2}\|\le2\). Consequently \[ \|G^{-1/2}MG^{-1/2}-M\| \le C\sqrt m\,e^{-1.6D} \le Cm e^{-1.6D}. \tag{113}\] This orthonormalizes the matrix of the shifted operator \(h+1/2\). In the original nonorthogonal modes, \(M=\mathcal U^*h\mathcal U+\tfrac12G\). The shift is essential in this estimate: one must add \(\tfrac12G\), rather than \(\tfrac12I\), before orthonormalization. Equation (108) follows.

The same residual estimate applied to \(A\mathcal J_0=A\mathcal U G^{-1/2}\) proves the first part of (109). Finally, as a map into \(H^1\), \(\mathcal U\) has norm at most \(C\sqrt m\), by (102). Thus \[\|\omega_i-\psi_i\|_{H^1} \le C\sqrt m\,e^{-.85D}.\] Combining this with \(\|\psi_i-\phi_i\|_{H^1}\le Cm/D\) gives the remaining assertion. ◻

Proof. Choose smooth real functions \(\theta_0,\theta_1,\ldots,\theta_m\) with \[\sum_{a=0}^m\theta_a^2=1,\qquad \theta_i=1\ \text{on }r_i\le D/5,\qquad \operatorname{supp}\theta_i\subseteq\{r_i\le D/4\}.\] The balls of radius \(D/4\) are disjoint by (105). In the annulus \(D/5<r_i<D/4\), take \(\theta_i=\cos f(r_i/D)\) and \(\theta_0=\sin f(r_i/D)\), where \(f\) is a fixed smooth transition from zero to \(\pi/2\), constant near both endpoints. Set \(\theta_0=0\) in the inner balls and one outside the outer balls. At most one transition annulus is present at any point, so \[ \sum_{a=0}^m|\nabla\theta_a|^2\le CD^{-2}. \tag{114}\] This construction also verifies smoothness at the boundaries.

Let \(u\in H^1(\mathbb R^3;\mathbb C^2)\), \(Pu=0\), and \(\|u\|=1\). In spin component \(\sigma\), global orthogonality implies \[b_{i\sigma}:=\langle\psi_i,\theta_i u_\sigma\rangle =\langle(\theta_i-1)\psi_i,u_\sigma\rangle.\] The tail estimate, integrated over \(r_i\ge D/5\), gives \[\|(\theta_i-1)\psi_i\| \le C\operatorname{poly}(D)e^{-.99(.2-.01)D} \le e^{-.18D}.\] Here the polynomial comes only from radial integration and has fixed degree. Consequently \[ \sum_{i,\sigma}|b_{i\sigma}|^2\le me^{-.36D}. \tag{115}\] There is no assumption that a localized vector is exactly orthogonal to the local ground state.

On \(\operatorname{supp}\theta_i\), every \(k\ne i\) satisfies \[r_k\ge D-2S-D/4=3D/4-2S>D/8.\] Hence \(\chi_k=1\) throughout that support and \(h=h_i\) there. Choose a fixed \(0<g_{\mathrm{loc}}\le\min(g_{\mathrm w}/2,1/4)\). Equations (104) and (103), for large \(D\), imply \[\begin{align*} \langle\theta_i u,A\theta_i u\rangle &\ge g_{\mathrm{loc}}\|\theta_i u\|^2 -(g_{\mathrm{loc}}+e^{-1.01D})\sum_\sigma|b_{i\sigma}|^2 . \end{align*}\]

On \(\operatorname{supp}\theta_0\), each primary distance is at least \(D/5\). Each secondary distance is at least \(D/5-\rho_i\ge D/5-4\ge D/6\). Since \(q+q'\le2\) for large \(D\), the magnitude of the full attraction on this support is at most \(Cm/D\). Thus \[\langle\theta_0u,A\theta_0u\rangle \ge(1/2-Cm/D)\|\theta_0u\|^2 \ge g_{\mathrm{loc}}\|\theta_0u\|^2.\]

The kinetic product rule yields the IMS identity \[\langle u,Au\rangle =\sum_{a=0}^m\langle\theta_au,A\theta_au\rangle -\frac12\int\sum_{a=0}^m|\nabla\theta_a|^2|u|^2.\] It holds first for smooth functions and then on \(H^1\) by continuity of the Coulomb forms. Summing the preceding bounds and using (114) and (115) gives \[\langle u,Au\rangle \ge g_{\mathrm{loc}}-Cme^{-.36D}-CD^{-2}\ge g_{\mathrm{loc}}/2.\] Taking \(g=g_{\mathrm{loc}}/2\) proves the proposition. ◻

Proof. Let \((K_{\mathrm{raw}})_{ij} =\langle\nabla\psi_i,\nabla\psi_j\rangle\). This is positive semidefinite, with trace at most \(Cm\). In the orthonormal spatial modes the kinetic matrix is \[K_0=G^{-1/2}K_{\mathrm{raw}}G^{-1/2}.\] Hence \[\operatorname{Tr}K_0 =\operatorname{Tr}(K_{\mathrm{raw}}G^{-1}) \le\|G^{-1}\|\operatorname{Tr}K_{\mathrm{raw}}\le Cm.\] The spin matrix \(K=K_0\otimes I_2\) has the same bound after adjusting the constant.

For any normalized fermionic \(u\), its one-particle density matrix \(\gamma\) in these modes satisfies \(0\le\gamma\le I\). Indeed, for a coefficient vector \(v\), the quadratic form of \(\gamma\) is \(\|c(v)u\|^2\); the anticommutation relation \(c(v)c(v)^\dagger+c(v)^\dagger c(v)=\|v\|^2I\) places it between zero and \(\|v\|^2\). Thus, for an arbitrary superposition in the finite-mode space, \[\sum_\ell\|\nabla_{y_\ell}u\|_2^2 =\operatorname{Tr}(\gamma K)\le\operatorname{Tr}K\le Cm.\]

For a pair \(\ell<k\), applying the three-dimensional Hardy inequality in \(y_\ell\), with \(y_k\) and all other variables held fixed, gives \[\|\,|y_\ell-y_k|^{-1}u\|_2 \le2\|\nabla_{y_\ell}u\|_2.\] Integration over the remaining variables and spin components justifies this inequality for the full wavefunction. Put \(b_m=m(m-1)/2\). Cauchy–Schwarz in the sum of pair multipliers gives, when \(m\ge2\), \[\begin{align*} \|W_{\mathrm{ee}}u\|_2^2 &\le b_m\sum_{\ell<k} \|\,|y_\ell-y_k|^{-1}u\|_2^2\\ &\le4b_m\sum_{\ell<k}\|\nabla_{y_\ell}u\|_2^2\\ &\le4b_m(m-1)\sum_\ell\|\nabla_{y_\ell}u\|_2^2 \le Cm^4. \end{align*}\] For \(m=1\) the repulsion is zero. This proves the operator bound on the whole finite-mode space, rather than only on its individual Slater determinants. ◻

Proof. Let \(\varepsilon=\|AP\|\le C\sqrt m\,e^{-.75D}\). The one-particle off-block operator \[R_{\mathrm{off}}=PA(I-P)+(I-P)AP\] extends to a bounded self-adjoint operator, of norm \(\|(I-P)AP\|\le\varepsilon\). The adjoint identifies the two blocks even though \(A\) itself is unbounded. Its block-diagonal remainder satisfies \[A-R_{\mathrm{off}} =PAP+(I-P)A(I-P)\ge-\varepsilon P+g(I-P)\] by Proposition 24. Consequently, on the full many-particle form domain, \[ d\Gamma(A)\ge gN_\perp-2m\varepsilon I, \qquad N_\perp=d\Gamma(I-P). \tag{120}\] Here the low diagonal can cost at most \(m\varepsilon\), and \(\|d\Gamma(R_{\mathrm{off}})\|\le m\varepsilon\). The commuting one-particle projections \((I-P)^{(\ell)}\) show that \(N_\perp\) has nonnegative integer eigenvalues, with kernel exactly \(\operatorname{Ran}\Pi\). Therefore \(N_\perp\ge Q_m\), including on superpositions of states with different numbers of high-mode particles.

The repulsion form is nonnegative on every wavefunction, so its compression to \(Q_m\) is nonnegative even though \(W_{\mathrm{ee}}\) need not commute with \(\Pi\). Equation (120) gives \[ Q_m\mathcal H Q_m\ge(g-2m\varepsilon)Q_m \ge cQ_m,\qquad c=g/2, \tag{121}\] for large \(D\).

By Lemma 25, the coupling from the finite subspace has operator norm at most \[ \beta:=\|Q_m\mathcal H\Pi\| \le m\varepsilon+\frac{Cm^2}{Z_*} \le Cm^{3/2}e^{-.75D}+\frac{Cm^2}{Z_*}. \tag{122}\] The first term follows because the block-diagonal one-particle sum preserves \(\operatorname{Ran}\Pi\). The compressed operator also obeys \[\|\Pi\mathcal H\Pi\|\le m\varepsilon+Cm^2/Z_*=o(1).\]

For any form-domain vector \(u=p+v\), where \(p=\Pi u\) and \(v=Q_m u\), the finite-mode operator-domain inclusion and the bounded coupling justify \[\begin{align*} \langle u,\mathcal H u\rangle &\ge\langle p,\mathcal H p\rangle +c\|v\|^2-2\beta\|p\|\|v\|\\ &\ge\langle p,\mathcal H p\rangle -\frac{2\beta^2}{c}\|p\|^2+\frac c2\|v\|^2. \end{align*}\] Equivalently, this is the form inequality \[\mathcal H\ge \Pi\mathcal H\Pi-\frac{2\beta^2}{c}\Pi+\frac c2Q_m.\] Since the finite block has norm \(o(1)\), the bottom of the right side is \(E_\Pi-2\beta^2/c\) for large \(D\). The opposite variational inequality follows by restricting test vectors to \(\operatorname{Ran}\Pi\). This proves the explicit error bound. In the present regime, with \(Z_*\ge e^D/2\), it is at most \(e^{-1.2D}\). All lower bounds apply to quadratic forms and hence to spectral infima; existence of an eigenvector at the bottom of the full continuum spectrum is not used. ◻

Proof. For four spatial functions define \[\mathcal V(a,b;c,d)= \iint\frac{\overline{a(x)b(y)}\,c(x)d(y)} {|x-y|}\,\,\mathrm dx\,\,\mathrm dy.\] Hardy in the \(x\) coordinate and Cauchy–Schwarz show that \[ |\mathcal V(a,b;c,d)| \le2\|a\|_2\|b\|_2\|\nabla c\|_2\|d\|_2. \tag{125}\] Telescoping a change of the four arguments therefore makes this coefficient Lipschitz in their \(H^1\) norms on uniformly bounded \(H^1\) sets. The same statement holds for fixed spin orbitals, with the usual spin Kronecker factors.

The coefficients of \(W_{\mathrm{eff}}\) are these integrals in the modes \(\omega_i\) and their two spin copies. Equation (109) permits their replacement by \(\phi_i\) with error at most \(Cm/D\) per coefficient. Indeed all four \(H^1\) norms stay uniformly bounded. In the standard second-quantized expression each coefficient multiplies \(\tfrac12c_a^\dagger c_b^\dagger c_d c_c\), a monomial of norm at most one. There are at most \((2m)^4\) coefficients. Thus this replacement changes the operator by at most \(Cm^5/D\). We are comparing coefficient arrays in the fixed abstract orthonormal site space; no orthogonality of the replacement functions \(\phi_i\) is asserted or needed.

Choose a fixed smooth radial cutoff equal to one on \(r_i\le D/5\), zero on \(r_i\ge D/4\), and with gradient bounded by \(C/D\), and let \(\zeta_i\) be its product with \(\phi_i\). Direct integration of the hydrogenic tail gives \[\|\zeta_i-\phi_i\|_{H^1} \le C\operatorname{poly}(D)e^{-D/5}\le e^{-.18D}.\] The supports of the \(\zeta_i\) are disjoint. By (125), replacing all \(\phi_i\) by \(\zeta_i\) changes the coefficient operator by at most \(Cm^4e^{-.18D}\).

For the truncated functions a nonzero coefficient must have the same site label for both factors at \(x\), and likewise at \(y\). Therefore the only surviving terms are onsite and intersite density terms. If \(i\ne j\), their spatial coefficient satisfies \[0\le V_{ij}^{\mathrm{dir}} =\iint\frac{|\zeta_i(x)|^2|\zeta_j(y)|^2}{|x-y|} \,\,\mathrm dx\,\,\mathrm dy \le\frac{1}{D/2-2S}\le\frac CD,\] using \(\|\zeta_i\|_2\le1\). The intersite operator is \(\sum_{i<j}V_{ij}^{\mathrm{dir}}n_i n_j\), where \(n_i=n_{i\uparrow}+n_{i\downarrow}\) has norm two. It consequently has norm at most \(Cm^2/D\).

The onsite spatial integral differs from \(\mathcal V(\phi_i,\phi_i;\phi_i,\phi_i)=U\) by at most \(Ce^{-.18D}\). The two equal-spin monomials vanish by anticommutation. The two opposite-spin monomials, together with the factor \(1/2\) in the interaction, give exactly \(U n_{i\uparrow}n_{i\downarrow}\). Combining these estimates yields \[\left\|W_{\mathrm{eff}} -U\sum_i n_{i\uparrow}n_{i\downarrow}\right\| \le C\left(\frac{m^5}{D} +m^4e^{-.18D}+\frac{m^2}{D}\right) \le\frac{Cm^5}{D}.\] ◻

Proof. The exact compression, in the abstract site space, is \[\mathcal J_m^*(Z_*\mathcal H)\mathcal J_m =d\Gamma\bigl(Z_*\mathcal J^*A\mathcal J\bigr) +W_{\mathrm{eff}}.\] The Hubbard hopping term is \(d\Gamma(\mathsf T\otimes I_2)\). On \(m\) particles, second quantization of a bounded one-particle difference increases its norm by at most a factor \(m\). Lemma 23 and Lemma 27 therefore give \[\left\|\mathcal J_m^*(Z_*\mathcal H)\mathcal J_m -H_{\mathrm{Hub}}(t)\right\| \le m\delta_1+Cm^5/D.\] For finite self-adjoint matrices, the difference of their spectral minima is bounded by the operator norm difference, by the variational principle. Lemma 26 adds at most \(Z_*e^{-1.2D}\) when the restriction to \(\operatorname{Ran}\Pi\) is removed. Finally, (117) gives the exact energy identity \[Z_*\inf\operatorname{Spec}\mathcal H =\frac{E_0}{Z_*}+\frac{mZ_*}{2}.\] This proves (126). Substitute (108) and \(Z_*\le2e^D\) to obtain (127).

The right side involves only fixed polynomial factors in \(m\), the calibrated hopping error, inverse powers of \(D\), and the displayed decaying exponentials. Equation (106) and Proposition 22 allow the hopping error to be made an arbitrarily small fixed inverse polynomial. The simultaneous choices, their order, and polynomial bit cost are specified in Section 6. No selected-orbital assumption has been imposed on the spectral infimum in this comparison. ◻