The electronic promise problem

Let \(M\ge1\) clamped nuclei have pairwise distinct positions \(R_1,\ldots,R_M\in\mathbb Q^3\), and let \(N\ge1\) be the number of electrons. Each nucleus has charge \(Z_\alpha=1\). In atomic units the electronic Hamiltonian is \[ \tag{H} H=-\frac12\sum_{i=1}^N\Delta_{x_i} -\sum_{i=1}^N\sum_{\alpha=1}^M\frac1{|x_i-R_\alpha|} +\sum_{1\leq i<j\leq N}\frac1{|x_i-x_j|}. \] Its Hilbert space and form domain are \[\mathcal H_N=\bigwedge^N\bigl(L^2(\mathbb R^3)\otimes\mathbb C^2\bigr), \qquad \mathcal Q_N=H^1(\mathbb R^{3N};(\mathbb C^2)^{\otimes N}) \cap\mathcal H_N.\] Antisymmetry exchanges position and spin together. No spin sector is selected. We use the lower-bounded self-adjoint operator associated with the closed Coulomb quadratic form on \(\mathcal Q_N\), and write \(E_0=\inf\operatorname{spec}H\). Equivalently, \(E_0\) is the infimum of this form over normalized vectors in \(\mathcal Q_N\), whether or not it is attained; see [12] for the closed-form representation and variational characterization. Section 7 gives the form bounds and core argument used in the continuum comparison.

The instance and its encoding are specified in Table 1. The promise excludes energies strictly between the thresholds. Nuclear repulsion is omitted: for clamped positions it is a fixed geometry-dependent scalar, whereas [eq:electronic-H] defines the electronic-energy convention used throughout.

Thus the bounded-charge problem holds with the absolute charge bound \(Z_{\max}=1\) and precision exponent \(c=1\). The theorem imposes no neutrality, binding, existence of a ground eigenvector, molecular spectral-gap, or lower nuclear-separation promise. It supplies no orbital basis, electronic state, magnetic field, or additional external potential. Membership in QMA is not asserted. The constructed inequality \(b-a\ge1\) is stronger than the input requirement \(b-a\ge L^{-1}\); neither inequality specifies an excitation gap of the molecule.

The strict-binding excess-charge bound in [8] bounds how many electrons fixed nuclei can bind strictly. Such a bound supplies neither a finite-dimensional truncation nor the continuum lower bound required here.

The companion article [7] develops the finite-gadget framework in greater detail and proves full-continuum hardness with binary-encoded nuclear charges. Its large-charge analytic construction and the polynomial-scale realization here are distinct. Our proof starts directly from the same external spin source and includes the finite estimates it uses; it does not invoke the companion’s hardness theorem.

Proof strategy

The starting point is the QMA-hard energy problem for signed Heisenberg exchange on finite square-lattice subgraphs. We use the spatially sparse construction of Cubitt, Montanaro, and Piddock, whose polynomial interaction scales are essential here [3]. Section 2 puts this source in a rational encoding with an explicit remaining promise gap.

The finite part of the reduction has two stages. Section 3 replaces signed exchange by positive exchange using singlet-pair mediators [9], and gives a planar layout with independently adjustable link lengths. Section 4 associates one localized spatial mode, with two spin states, to each of the resulting \(m\) sites. With \(m\) electrons, a positive occupancy penalty favors one electron per site. Virtual hopping through double occupancy then gives the prescribed positive spin exchange, by the superexchange mechanism of Anderson [1]. The direct Coulomb coefficients enter that penalty and its virtual-excitation denominators explicitly.

The remaining construction realizes these modes using unit nuclei. We first work in dilated coordinates \(y=\lambda x\), with \(\lambda\) a polynomially large scale. After dividing energies by \(\lambda^2\), a unit nucleus has coupling \(1/\lambda\). Section 5 constructs a nonnegative continuous nuclear density producing the planar wells and transverse confinement. A broad, constant-density slab supplies the positive curvature needed to keep the modified density nonnegative. Small changes in the well depths supply a one-site counterterm for the long-range direct electron repulsion.

Section 6 replaces the density by distinct point nuclei. A smooth flow transports uniform density to the constructed density. Applying this flow to a tensor-product two-point Gauss rule preserves equal node masses. The mass is chosen to be exactly \(1/\lambda\), so returning to physical coordinates gives charge one at every nucleus. Rational rounding of the positions is included in the estimates.

The key analytical distinction is between the replacement error on arbitrary form-domain vectors and its compression to the localized modes. The latter is smaller because the modes have smooth Coulomb fields. In Section 7, a uniform one-particle complementary gap excludes low energies from the omitted continuum space, and square completion makes coupling to that space enter the energy quadratically. This gives a lower energy bound on every antisymmetric form-domain state, uniformly in spin, and a matching variational upper bound from the mode space. Together they compare the full continuum infimum with the finite fermion ground energy. Retaining the direct Coulomb interaction and its counterterm is necessary at the polynomial scales used here. Section 8 completes the deterministic computation of the nuclear positions and thresholds, including all scalar shifts.

Energy reduction through singlet pairs

Write \(E\) for the smallest energy in finite dimensions. The following elementary second-order estimate is a form of the usual effective Hamiltonian expansion [2]; we prove precisely the bottom-energy statement needed here. In particular, the first-order block \(PBP\) need not vanish.

First delete all extremely small source weights, below a fixed inverse polynomial in absolute value chosen so the deletion error is at most \(\gamma/100\), by \(\|h^s_{ij}\|=3\). Here and below such cutoffs can be rational. The following operation can be made on a set of edges simultaneously, keeping the other terms unchanged. For each replaced edge \(e=LR\), of weight \(J\ne0\), insert a fresh central pair \(A_e,B_e\), with bond \(\Delta h^s_{A_eB_e}\), and two spokes \[\sqrt\Delta\,v_e(h^s_{LA_e}+h^s_{C_eR}),\qquad v_e=\sqrt{2|J|},\qquad C_e=\begin{cases}B_e& J>0,\\ A_e&J<0.\end{cases}\] Here \(\Delta\) is a common large positive integer. Recall that \(h^s_{ij}=2\operatorname{SWAP}_{ij}-I\) with eigenvalues \(-3,1\) on the singlet and triplets, respectively. The new energy is the previous one minus \[\sum_e (3\Delta+3v_e^2/2)\] up to error that can be made at most \(\gamma/100\). To verify this shift, penalize excitations of the fresh pairs by \(A=\Delta\sum_e(h^s_{A_eB_e}+3I)\), with singlet kernel for each pair and gap \(4\Delta\); the rest \(B\) has norm at most \(\operatorname{poly}(n)(1+\sqrt\Delta)\) if the entering model has polynomial size and weights. The first-order compression is the unreplaced part.

Only two actions on the same fresh pair contribute to the inverse term, with denominator \(4\Delta\). Pauli products on one member of that singlet have expectation \(\delta_{\mu\nu}\) (components indexed by \(\mu,\nu\)), and products across the pair have \(-\delta_{\mu\nu}\), by the singlet Pauli relations. Thus the negative inverse term gives \(J h^s_{LR}-(3/2)v_e^2 I\) per edge. The estimate just proved permits a polynomially large \(\Delta\).

Apply this conversion first to all surviving signed edges. For each of its spokes make a path of length nine by two subdivision stages, each using the same operation for positive weights (a path of length three per replaced edge). In the first subdivision replace both spokes, and in the second replace all edges of their resulting paths, not the original central-pair bond. Each stage can use its own polynomial scale as there are only three stages total. Thus in the final positive model \(H_+=\sum K_e h^s_{ij}\) there are, per input edge still present, one central bond and two nine-link paths, from \(L\) to \(A_e\) and from \(C_e\) to \(R\). Let \(m\) be the total number of sites (including original vertices); every actual link weight \(K_e\) and its reciprocal are polynomially bounded. There is a sum of scalar shifts \(C_s\) from the operations, such that \[|E(H_+)-E(H_{\rm in})+C_s|\le 4\gamma/100. \tag{1}\] All these intermediate weights and shifts are polynomial-time approximable to inverse polynomial error (the square roots need not be exact rational outputs).

Independently adjustable contact lengths

Equation [eq:1] completes the energy conversion. To encode its positive weights by orbital hopping, we now need a geometry in which changing one link length does not constrain the others. Here is a planar contact layout in which all the links have independently adjustable lengths; see Figure 1. Start with original grid spacing 17. Orient the edges in the two positive axial directions, and within each edge gadget use planar coordinates along that direction and in a perpendicular direction, origin at \(L\). In the positive input sign case set \(A_e=(8,0), B_e=(9,0)\); in the negative sign case set \(A_e=(8.5,0), B_e=(8.5,1)\). Recall \(C_e=B_e\) and \(C_e=A_e\), respectively. Each nine-link path joins two on-axis points \(8\) or \(8.5\) apart. Orient it along the edge in order of increasing first coordinate. Make the first two, middle (fifth), and last two steps unit forward steps on that direction, and use two unit forward-sloping steps (steps 3 and 4) with first component \(3/4\) for span \(8\), \(7/8\) for span \(8.5\), balanced by two steps (6 and 7) with transverse components negated.

The link lengths now equal one, and nonlinks have distance at least \(1.4\). Indeed consecutive separations projected along a main path/backbone are at least \(3/4\); nonconsecutive sites on it are separated by at least \(1.5\), and the extra side site for a negative sign is at distance at least \(\sqrt2\) from sites other than its central neighbor (the steps immediately there are straight). These bounds also give separation along each entire axial lattice line across successive gadgets. Parallel lines cannot interfere. Transverse offsets are at most \(2\sqrt{1-(3/4)^2}<1.324\); off-axis internal sites have along-axis projection at least \(2+3/4\) away from bounding original grid coordinates. Against sites of a perpendicular orientation this already gives separation greater than \(1.4\). If both internal sites of perpendicular orientations are on-axis, their coordinate differences in the two directions are at least one each. Original vertices other than a gadget’s two endpoints are far from its internal sites by grid spacing. This verifies the separation.

The unit-length layout already separates links from nonlinks. We now show that its link lengths can vary independently without losing that separation.

There is a fixed \(\varepsilon\in(0,1/100)\) such that any prescribed link lengths in \([1-\varepsilon,1+\varepsilon]\) can be implemented with original vertices fixed and nonlinks still at distance \(>1.2\). To prove this, consider one original-edge gadget. We can first change the central length by moving \(B_e\), then vary each path’s internal sites while its endpoint velocities are prescribed.

For a path with unit step vectors \(e_1,\ldots,e_9\), a linear dependence with coefficients \(a_1,\ldots,a_9\) among the rows of its fixed-end length differential would give, at each internal vertex, \[a_j e_j-a_{j+1}e_{j+1}=0 \qquad (1\le j\le8).\] All vectors \(a_j e_j\) would therefore be equal. Since the path has nonparallel steps, they must all vanish. The differential is thus surjective, as is the full gadget’s link-length differential on its internal-site variables.

Choose a fixed right inverse for each of the two templates and restrict the internal-site displacements to its image. In these coordinates the length increment map is \(q\mapsto q+R(q)\), with \(R(0)=0\) and \(DR(0)=0\). On a small fixed ball, \(\|DR\|<1/2\); for a sufficiently small prescribed increment \(d\), the iteration \(q\mapsto d-R(q)\) is a contraction of that ball. Choose the ball small enough to preserve the nonlink separation, and then fix \(\varepsilon\) accordingly. These solves have fixed dimension and algebraic initial data. Approximating the iterations with guard precision computes positions to any inverse-polynomial error in polynomial time. Because the original vertices stay fixed, separate gadgets can be solved independently; their number does not shrink the allowed length interval.

Local modes and Coulomb coefficients

Take \(f(r)=(1-|r|^2)_+^{16}\), and fix \[\phi=(-\Delta_r+1)^{-1}f,\qquad \varphi=\phi/\|\phi\|_2,\qquad W(r)=-\frac{f(r)}{2\phi(r)} .\] These functions have more bounded continuous derivatives than needed below (in particular through order six). \(W\) is nonpositive and compactly supported, and \(\phi\) is strictly positive, radial, exponentially decaying. For details, \[\phi(r)=\int_0^\infty e^{-t}\int_{\mathbb R^2} (4\pi t)^{-1}e^{-|r-b|^2/(4t)} f(b)\,db\,dt\] solves the equation by the Gaussian heat kernel identity. For \(|r|>2\), it and the derivatives needed below are bounded in absolute value by \(C(1+|r|)^C e^{-|r|}\), by differentiating on \(f\) and using \(t+(|r|-1)^2/(4t)\ge |r|-1\); for example integrate \(dt/t\) with this bound from 1 to \(|r|^2\), bounding the omitted tails directly by Gaussian or exponential decay. Integrating over \(b\) in the inner half disk and \(|r|/2<t<|r|/2+1\) also gives a lower bound on \(\phi\) of the form \(c(1+|r|)^{-C}e^{-|r|}\).

The isolated planar operator \(-\Delta_r/2+W\) has ground state \(\varphi\), energy \(-1/2\), and a fixed positive gap above that mode. Indeed the shifted form on test functions \(u\) is \(\int \phi^2 |\nabla(u/\phi)|^2/2\); by closure and local positivity this gives nonnegativity and a simple kernel. A unit sequence perpendicular to the kernel with form tending to zero would be bounded in \(H^1\) and, upon extraction, converge weakly and locally in \(L^2\) (strongly on compact sets). By compact support of \(W\) its limit would be in the kernel, hence zero. But then the \(1/2\) shift of the form precludes convergence of the form values to zero. For the transverse direction \(z\) we use the normalized Gaussian \(g(z)=(\omega/\pi)^{1/4}e^{-\omega z^2/2}\) of \(-\partial_z^2/2+2\pi\rho z^2\) with energy \(\omega/2\), \(\omega=\sqrt{4\pi\rho}\), where \(\rho>0\) is a large fixed integer to be chosen. This mode too has a fixed gap (same positive ground-state identity, with compactness here from confinement).

We will place the planar centers \(u_i\) at scale \(D=\log\lambda\) using the contact layout. For all sufficiently accurate implementations, we require minimum spacing \(\ge0.95D\), nonlink spacing \(\ge1.2D\), and extent \(\operatorname{poly}(n) O(D)\). Define \[\begin{gathered} \psi_i(y)=\varphi(r-u_i) g(z),\quad e_*=-1/2+\omega/2,\\ v_{ij}=\iint\frac{\psi_i(y)^2\psi_j(y')^2}{|y-y'|}\,dy\,dy' =v(|u_i-u_j|),\qquad U=v(0). \end{gathered}\] These constants measure direct electron repulsion: \(U\) is the same-site coefficient and \(v_{ij}\) the coefficient between sites \(i\) and \(j\). They are bounded uniformly. The function \(v\) is radial as indicated and uniformly Lipschitz, as follows by translating the smooth density with decaying derivatives.

To control departures from one electron per site, we need a uniform lower bound on the Coulomb quadratic form. We prove \((v_{ij})\ge cI\), with fixed \(c>0\), uniformly over these separated geometries once the polynomial scale is sufficiently large. For real coefficients \(a_i\), put \(F=\sum_i a_i\psi_i^2\). Choose a fixed smooth nonnegative, nonzero test bump near the origin, and let \(\chi_i\) be its disjoint translates to the centers \((u_i,0)\). With \(w=\sum_i a_i\chi_i\), Poisson’s equation and integration by parts give \[\left|\int Fw\right|^2\le \frac1{4\pi}\left(\iint\frac{F(y)F(y')}{|y-y'|}\,dy\,dy'\right) \int|\nabla w|^2\] for these smooth decaying densities, or first after truncation at infinity. To bound the left side below even for signed \(a_i\), let \(A_{ij}=\int\psi_i^2\chi_j\). Its diagonal entries equal a fixed \(a_*>0\), while the row and column sums of its off-diagonal entries are exponentially small in \(D\), up to a polynomial factor in \(m\). Thus \(A=a_*I+R\), with \(\|R\|\le a_*/2\) at our scales, and \[\int Fw=a^{\mathsf T}Aa\ge(a_*/2)\sum_i a_i^2, \qquad \int|\nabla w|^2=O\Bigl(\sum_i a_i^2\Bigr).\] The displayed Coulomb inequality proves the required lower bound. Here and below we take \(\lambda\) so large a power of \(n\) that exponential-in-\(-D\) bounds beat the polynomial factors. In particular, applying the bound to two sites and coefficients \((1,-1)\) shows that \(U-v(d)\) is bounded away from zero for every distance \(d\ge0.95D\) in use; it is also uniformly bounded above.

The occupancy penalty and spin exchange

For the moment keep any geometry satisfying the separation bounds, and let the hopping strengths \(t_e\) be free parameters. On \(2m\) canonical fermion modes, two spins at each site, restrict to \(m\) electrons and write \(n_i\) for total occupancy at site \(i\). For this finite model we multiply dilated energies by \(\lambda\), so its Coulomb coefficients have no \(1/\lambda\) factor. The direct Coulomb interaction is \[\frac U2\sum_i n_i(n_i-1)+\sum_{i<j}v_{ij}n_i n_j.\] Set \(s_i=\sum_{j\ne i}v_{ij}\). A diagonal one-body term \(-\sum_i s_i n_i\) cancels the terms linear in the occupancy defects: using \(v_{ii}=U\) and \(\sum_i n_i=m\), \[ \tag{C} \begin{aligned} &\frac U2\sum_i n_i(n_i-1)+\sum_{i<j}v_{ij}n_i n_j-\sum_i s_i n_i\\ &\hspace{2em}=\frac12\sum_{i,j}v_{ij}(n_i-1)(n_j-1)-\sum_{i<j}v_{ij}. \end{aligned} \] Section 5 will produce this one-body term by adjusting the well depths. The finite Hamiltonian to be realized is therefore \[\begin{split} H_d(t)={}&\frac12\sum_{i,j}v_{ij}(n_i-1)(n_j-1)-\sum_{i<j}v_{ij}\\ &-\sum_{e=\{i,j\}}t_e \sum_{\sigma=\uparrow,\downarrow} (c_{i\sigma}^\dagger c_{j\sigma}+c_{j\sigma}^\dagger c_{i\sigma}). \end{split} \tag{2}\] The nonnegative quadratic penalty has kernel exactly the singly occupied spin space. It has a fixed gap: a nonzero integer defect vector \((n_i-1)_i\) has sum zero, hence squared length at least two, so its penalty is at least \(c\) by Coulomb coercivity.

Hopping vanishes on first-order compression to single occupancy. Its second-order effect is the double-occupancy superexchange familiar from Anderson [1]; his transfer-energy discussion also retains intersite Coulomb contributions. The normalization and error bound below are proved for the precise matrix used here.

One allowed hop on \(e=\{i,j\}\) creates a doubly occupied site and an empty site. Its defect vector is \((1,-1)\) on that pair and zero elsewhere, so its excitation energy is exactly \(U-v_{ij}\). Only a hop on the same link can restore single occupancy in one further step. For hopping strength \(t_e\), the negative second-order operator is \[\frac{t_e^2}{U-v_{ij}}(h^s_{ij}-I).\] Indeed the triplets do not hop into same-site pairs, while the singlet hops into two orthogonal doubly occupied states with amplitudes of magnitude \(\sqrt2\,t_e\). Apply the creation and annihilation operators to \((c_{i\uparrow}^\dagger c_{j\downarrow}^\dagger- c_{i\downarrow}^\dagger c_{j\uparrow}^\dagger)|0\rangle/\sqrt2\) to obtain these amplitudes. Other sites can be placed after these two in the ordering: changing the ordering of singly occupied sites gives only an overall sign on the spin space.

With a common positive parameter \(\tau\) to be fixed below, choose target hoppings \[t_e^{\rm tar}=\tau\sqrt{K_e(U-v_{ij})}.\] They give the effective spin term \(\tau^2(H_+-\sum_e K_e)\). Their hopping operator has norm at most \(\tau\operatorname{poly}(n)\), uniformly over admissible separated geometries. By Lemma 2, the energy remainder is at most \(C\tau^3\operatorname{poly}(n)\), since the penalty gap is uniform. Fix a rational inverse-polynomial \(\tau>0\) small enough that the lemma applies and this remainder is at most \(\gamma\tau^2/100\). This choice is independent of \(\lambda\) and of the final centers.

Calibrating the separations

The modes supply the hopping profile \[T(d)=-\int_{\mathbb R^2}\varphi(r)W(r)\varphi(r-d e_1)\,dr =\frac{1}{2\|\phi\|_2^2}\int f(r)\phi(r-d e_1)\,dr .\] It is positive, uniformly Lipschitz, and lies between \(c(1+d)^{-C}e^{-d}\) and \(C(1+d)^Ce^{-d}\) by the tail estimates. The continuum compression will supply hopping \(t_e=\lambda T(|u_i-u_j|)\) in the units of \(H_d\). Although the occupancy penalty involves every pair, the target hopping on link \(e\) depends only on its own distance through \(v_{ij}=v(d)\). Thus each link requires a scalar solve for \[F_e(d)=\lambda T(d)-\tau\sqrt{K_e(U-v(d))}.\] At \(d=(1-\varepsilon/2)D\) and \(d=(1+\varepsilon/2)D\), respectively, \(F_e\) is positive and negative if the fixed polynomial exponent in \(\lambda\) is sufficiently large. Indeed the first term has respective sizes \(\lambda^{\varepsilon/2}\) and \(\lambda^{-\varepsilon/2}\), up to powers of \(D\), whereas the target has polynomial upper and inverse-polynomial lower bounds in \(n\).

Bisection with approximate sign tests finds an inverse-polynomially small residual: stop if the sign is ambiguous within the evaluation error, or continue until the bracket width times a Lipschitz bound is small enough. No monotonicity or derivative lower bound is required. The link-length construction in Section 3, rescaled by \(D\), then implements all these distances at once. Its contraction solves and a final rational approximation give rational centers with the required separations and hopping accuracy. Section 8 details the numerical evaluation.

From now on, \(H_d\) means \(H_d(t)\) at these final centers with \(t_e=\lambda T(|u_i-u_j|)\). Choose the per-link calibration accuracy so that replacing \(t_e^{\rm tar}\) by \(t_e\) changes the total Hamiltonian norm by at most \(\gamma\tau^2/100\). Together with the second-order remainder, this gives \[|E(H_d)+\sum_{i<j}v_{ij} -\tau^2(E(H_+)-\sum_e K_e)|\le 2\gamma\tau^2/100. \tag{3}\]

Polynomial scale convention

Equation [eq:3] fixes the finite-model accuracy required of the continuum construction. All subsequent errors must fit inside an inverse-polynomial fraction of \(\gamma\tau^2\) in the energy units of \(H_d\).

In the estimates below \(p(n,D)\) denotes a polynomial upper bound with uniform coefficients and exponents, which may be enlarged between uses. All requirements of the form "make \(p(n,D)\lambda^{-\alpha}\) small" for fixed \(\alpha>0\), with required smallness at worst inverse polynomial in \(n\), can be imposed simultaneously. Specifically we take an integer \(k\) growing as a sufficiently large fixed power of \(n\) (allowing a fixed factor) and put \[h=1/k,\qquad \lambda=8k^3/\rho .\] Here \(h\) will be the side length of the cubature cells in Section 6; the relation \(\rho h^3/8=1/\lambda\) gives each of their eight nodes the mass needed for a unit nucleus. Other magnitudes such as the large box used below will have stated polynomial bounds in \(\lambda,n\); its volume will enter \(p(n,D)\) only when explicitly accounted for. The fixed \(\rho\) only has to dominate bounded potential derivatives as described next and can be chosen before these scale choices.

The slab and the well-depth counterterm

Use integers \(S\) of size between \(D\) and \(10D\), and \(H=\lceil\lambda^5\rceil\). (For \(S\) one can use a binary bit-length scale.) Let \[\Omega=[-H,H]^2\times[-S,S],\qquad F_H(y)=-\rho\int_\Omega\frac{db}{|y-b|},\qquad C_H=F_H(0).\] The occupancy identity [eq:occupancy-completion] requires the one-body diagonal \(-s_i/\lambda\), where \(s_i=\sum_{j\ne i}v_{ij}\). Put \(W_i=W(\,\cdot-u_i)\). Since \(\langle\psi_i,W_i\psi_i\rangle=\int\varphi^2W<0\), a small increase in each well depth has the required sign and size. Take a fixed sufficiently smooth cutoff \(\zeta\) between 0 and 1, supported inside \((-1,1)\), equal to 1 for \(|z|\le 1/2\); it may be piecewise polynomial with rational data. Define \[b_i=\frac{s_i}{-\int\varphi^2 W},\qquad V_\ell(r,z)=\zeta(z/S)\sum_i(1+b_i/\lambda)W_i(r).\] These \(b_i\ge0\) are polynomially bounded independent of \(D\); hence \(0\le b_i/\lambda<1\) at our scales. For all \(S\ge1\) the derivatives we need of \(V_\ell\) are bounded independently, by disjoint planar supports. Choose the fixed integer \(\rho\) so large, using \(1+b_i/\lambda\le2\), that \(|\Delta V_\ell/(4\pi)|\le\rho/2\). There is no circularity, since the bound based on \(2,W,\zeta\) does not depend on the vertical Gaussian or \(b_i\). The centers fit well inside the horizontal box for our scale. Thus \[F_H+V_\ell =-\int_\Omega \frac{\rho+\Delta V_\ell(b)/(4\pi)}{|y-b|}\,db \tag{4}\] by compact support and the fundamental solution identity.

Thus the density in [eq:4] lies between \(\rho/2\) and \(3\rho/2\) on \(\Omega\). Its total mass is \(\rho|\Omega|\), because \(\int\Delta V_\ell=0\). The next estimates show that the truncated slab approximates the harmonic potential near the modes and provides a global lower bound outside that region.

Here are useful background estimates with constants independent of \(n,\lambda\): \[\begin{split} F_H(y)-C_H &\ge 2\pi\rho \min(z^2,S^2)- C S^3/H,\\ |F_H(y)-C_H-2\pi\rho z^2| &\le C(S^3+S|r|^2)/H \qquad (|r|\le H/2,\ |z|\le S). \end{split} \tag{5}\] Indeed \(F_H\) at each \(z\) is minimized at \(r=0\), and is also nondecreasing in \(|z|\) for fixed \(r\). These statements follow one coordinate at a time by integration of even decreasing kernels on centered intervals. At \(r=0\), extending the horizontal integration to the plane in the difference between \(z\) and 0 gives for \(|z|\le S\) exactly \(2\pi\rho\int_{-S}^S(|z-w|-|w|)\,dw=2\pi\rho z^2\), by polar integration. The omitted difference costs \(CS^3/H\), by inverse cube decay outside horizontal radius \(H\). For \(|r|\le H/2\), the change relative to the horizontal center is bounded by \(CS|r|^2/H\). For example translate the integration square instead of translating the kernel. The first derivatives then come from sides at distance at least \(H/2\), and second derivatives are bounded by \(CS/H\) by integrating first derivatives of the kernel on the sides. Evenness at zero gives the bound. Finally \(0\ge F_H(y)\ge C_H\) and \(|C_H|\le CSH\).

The one-particle gap and residual

The potential has now been expressed as the Coulomb field of a nonnegative density. Before replacing it by point nuclei we need two properties: a fixed energy gap outside the prescribed modes, and a small residual on those modes. The gap will control arbitrary continuum states, while the residual will control mixing with the mode space.

Let the shifted continuous one-particle operator, still before point-charge replacement, be \[L_c=-\Delta_y/2+F_H+V_\ell-C_H-e_* .\] We suppress the tensor identity on spin. Let \(P\) project onto the span of the \(\psi_i\), including both spins when present, and write \(\Gamma_{ij}=\langle\psi_i,\psi_j\rangle\) for their spatial Gram matrix. Exponential tails and the \(0.95D\) spacing give \(\|\Gamma-I\|\le p(n,D)\lambda^{-0.8}\): for example, each off-diagonal overlap is \(O(e^{-0.9|u_i-u_j|})\), by leaving part of the decay as an integrable factor. We use the orthonormal modes \[\widehat\psi_j=\sum_i\psi_i(\Gamma^{-1/2})_{ij},\qquad \|\Gamma^{-1/2}-I\|\le p(n,D)\lambda^{-0.8}.\] In particular, unit vectors in \(\operatorname{Ran}P\) have polynomially bounded one-particle kinetic energy.

The compressed one-particle matrix

The fixed complementary gap alone does not give the required energy accuracy. We therefore compute the matrix on the prescribed modes, including its diagonal counterterm, to a sharper scale.

The continuous one-body compression can now be evaluated more precisely. On the nonorthogonal functions before Gram correction, the column on mode \(j\) comes, up to background and cutoff errors as above, from \[\Big(\sum_{k\ne j} W_k+\lambda^{-1}\sum_k b_k W_k\Big)\psi_j .\]

On row \(j\), the leading contribution is \(-\lambda^{-1}\sum_{i\ne j}v_{ij}\) by construction; the other-well integrals there use two tail factors. On row \(i\ne j\) the first sum at \(k=i\) gives \(-T(|u_i-u_j|)\). The cases \(k\ne i,j\) of the first sum are products of two tail factors at a third site; each such integral is \(O(\lambda^{-1.8})\) by the spacing and compact well support. Terms with the explicit \(\lambda^{-1}\) give only negligible off-diagonal errors (at least one planar tail on each well support, hence bounded by \(p(n,D)\lambda^{-1.9}\) altogether on such an entry). Nonlink hoppings cost at most \(p(n,D)\lambda^{-1.2}\) per term.

To pass to the orthonormal modes, let \(A_{ij}=\langle\psi_i,L_c\psi_j\rangle\). The matrix becomes \(\Gamma^{-1/2}A\Gamma^{-1/2}\), and \[\|\Gamma^{-1/2}A\Gamma^{-1/2}-A\| \le C\|\Gamma-I\|\,\|A\| \le p(n,D)\lambda^{-1.7}\] by [eq:6]. Here \(A\) is the shifted matrix: subtracting \(C_H+e_*\) before making this estimate removes the large background scalar. The correction therefore remains negligible after multiplication by \(\lambda\).

Define the spatial matrix \[\mathsf M_{ij}=\begin{cases} -s_i,&i=j,\\ -\lambda T(|u_i-u_j|),&\{i,j\}\text{ is a link},\\ 0,&\text{otherwise}. \end{cases}\] The preceding estimates give \[\bigl\|\lambda\Gamma^{-1/2}A\Gamma^{-1/2}-\mathsf M\bigr\| \le p(n,D)\lambda^{-0.1}.\] This is the required one-particle compression, with the identity on spin. The diagonal supplies the depth counterterm, and the link entries supply the calibrated hopping.

Transport and charge placement

The density interpolation uses Moser’s deformation method [5], also developed by Dacorogna and Moser [4]. Here the vector field is explicit and compactly supported inside the box; the following continuity-equation calculation proves the required transport directly. Write \(q_t=\rho+t\Delta V_\ell/(4\pi)\), so that \(q_0=\rho\) and \(q_1\) is the desired density on \(\Omega\). Since \(q_t\ge\rho/2\), the time-dependent vector field on full space \[\mathcal V_t(b)=-\frac{\nabla V_\ell(b)}{4\pi q_t(b)}, \qquad 0\le t\le1,\] has uniformly bounded derivatives through order four. It vanishes outside the well cylinders strictly inside the box. Let \(G_t\) be its flow from time zero and put \(G=G_1\). The ODE and its differentiated equations give bounded derivatives of \(G\) through order four and a bounded Lipschitz inverse. The flow is the identity off the cylinders and preserves \(\Omega\).

The identity \(\partial_t q_t+\nabla\cdot(q_t\mathcal V_t)=0\) gives \[q_t(G_t(b))\det DG_t(b)=\rho,\qquad \int_\Omega F(b)q_1(b)\,db =\rho\int_\Omega F(G(b))\,db.\] Thus the flow transports the uniform measure to the target density. The last formula also shows how to retain equal cubature masses: move the nodes by \(G\) and keep their weights unchanged.

Since \(H,S,k\) are integers, cubes of side \(h=1/k\) tile \(\Omega\). On each cube use the eight tensor products of the two equal-weight Gauss points per coordinate, at the midpoint plus or minus \(h/(2\sqrt3)\). The rule integrates polynomials of total degree at most three exactly, by the corresponding moments on each interval. Map these nodes by \(G\). Their common mass remains \(\rho h^3/8=1/\lambda\); polynomial exactness is used in the original cube coordinates, on integrands composed with \(G\).

Each exact mapped point \(b_\alpha\) will be approximated within Euclidean error at most \[\delta_R\le (1+H)^{-3}(1+S)^{-2}\lambda^{-3}\] by rational points \(\tilde b_\alpha\), and the physical nuclear positions are \(\tilde b_\alpha/\lambda\). They remain distinct by node separation and bi-Lipschitz bounds.

Form and compression estimates

Let \(L_p\) be the shifted one-particle operator with these actual point nuclei in dilated units, i.e. with potential \(-\sum_\alpha \lambda^{-1}/|y-\tilde b_\alpha|-C_H-e_*\).

We have obtained the two replacement estimates with constants independent of the horizontal area of the slab. The full mass enters only the rounding bound, where the chosen coordinate accuracy absorbs it. Section 7 now combines these estimates with the fixed complementary gap and the compressed matrix from Section 5.

Control outside the orbital space

Compressed energy

It remains to evaluate \(\lambda E_\Pi\). Section 5 computed the one-particle matrix \(\mathsf M\) in the orthonormal modes \(\widehat\psi_i\), to accuracy \(p(n,D)\lambda^{-0.1}\). Multiplying the compressed replacement error in [eq:7] by \(\lambda\) costs only \(p(n,D)\lambda^{-1/3}\), so the same matrix accuracy holds for the point nuclei. For the repulsion \(\mathcal R\), the pair matrix on the original modes before Gram correction has spatial entries from products of two transition densities. Keeping only \(\psi_i^2,\psi_j^2\) as the respective densities gives the direct coefficients \(v_{ij}\). Any unequal indices within a density give an exponentially small error: \(\|\psi_i\psi_k\|_{L^1}=O(\lambda^{-0.8})\) for \(i\ne k\), and the Coulomb field of the other density product in absolute value is uniformly bounded by orbital boundedness and decay. The matrix of direct interactions in the product basis is bounded, and Gram correction changes its entries negligibly.

These entrywise estimates give a many-electron operator-norm error of at most \(p(n,D)\lambda^{-0.1}\). To see the dimension dependence, write the one- and two-body matrices in the \(2m\) orthonormal spin modes. Each coefficient multiplies a monomial \(c_a^\dagger c_b\) or \(c_a^\dagger c_b^\dagger c_d c_c\), of norm at most one. There are only \(O(m^2)\) and \(O(m^4)\) such monomials, respectively, so summing their coefficient errors costs a polynomial factor on the entire fermion space. In this representation, the retained two-body term is \[\frac{U}{2}\sum_i n_i(n_i-1)+\sum_{i<j}v_{ij} n_i n_j .\] The retained one-body diagonal is \(-\sum_i(\sum_{j\ne i}v_{ij})n_i\). At the filling \(\sum_i n_i=m\), the identity [eq:occupancy-completion] combines these terms into the first line of [eq:2], including its scalar. The hopping matrix gives the remaining line of the calibrated \(H_d\). We have thus identified the compressed Hamiltonian, and the comparison with the full continuum follows.

This analysis treats long-range diagonal effects and possible charge transfer at the virtual-excitation scale in [eq:2]-[eq:3], rather than requiring all intersite Coulomb matrix elements to be negligible compared to the input promise gap.

Table 2 records the different roles of the errors. In particular the arbitrary-state form error cannot be used as a direct energy error after multiplication by \(\lambda\).

Physical energies and thresholds

Choose the polynomial scale in \(\lambda\) large enough that [eq:9] costs at most \(\tau^2\gamma/100\) and that \(\lambda\tau^2\gamma\ge8\). All earlier smallness, separation, and calibration estimates hold simultaneously by the stated growth choices. Undoing the dilation and scalar shift gives the physical energy exactly as \[E_0=\lambda^2\left[m(C_H+e_*)+ \inf\operatorname{Spec}\mathcal H\right].\] Combining [eq:1], [eq:3], and [eq:9] now gives a total error at most \(7\lambda\tau^2\gamma/100\), and hence at most \(\lambda\tau^2\gamma/8\), from \[\mathcal F(E(H_{\rm in})),\qquad \mathcal F(\xi)= \lambda^2 m(C_H+e_*)+ \lambda\left[-\sum_{i<j}v_{ij} +\tau^2\left(\xi-C_s-\sum_e K_e\right)\right]. \tag{10}\] Take output thresholds rationally approximating \(\mathcal F(a_s+\gamma/4)\) and \(\mathcal F(b_s-\gamma/4)\) each within \(\lambda\tau^2\gamma/32\). For a YES source instance, \[E_0\le\mathcal F(a_s)+\lambda\tau^2\gamma/8 <\mathcal F(a_s+\gamma/4)-\lambda\tau^2\gamma/32\le a.\] For a NO source instance, the reversed estimate gives \[E_0\ge\mathcal F(b_s)-\lambda\tau^2\gamma/8 >\mathcal F(b_s-\gamma/4)+\lambda\tau^2\gamma/32\ge b.\] Since \(b_s-a_s\ge\gamma\), the threshold separation satisfies \[b-a\ge\lambda\tau^2\bigl(b_s-a_s-\gamma/2\bigr) -\lambda\tau^2\gamma/16 \ge\frac7{16}\lambda\tau^2\gamma\ge\frac72>1.\] Thus both promises survive rational rounding, and the reduction meets the required output separation \(b-a\ge1\).

Order and uniformity of the scale choices

Fix the bounded well data and the integer \(\rho\) using the positivity bound in [eq:4], under the conditions that the well supports are disjoint and \(b_i/\lambda\le1\). As explained there, that derivative bound uses only \(2,W,\zeta\), so it does not depend on the later Gaussian or Coulomb coefficients. Next fix the source cutoff and three spin gadget penalty scales from the source size and gap bounds. Fix the inverse-polynomial rational \(\tau\) using the uniform Coulomb coercivity constant for sufficiently large \(D\). Finally choose \[k=\lceil Cn^A\rceil,\qquad h=k^{-1},\qquad \lambda=8k^3/\rho,\] with sufficiently large absolute positive integers \(A,C\).

At this last stage, \(K_e,1/K_e\), the layout extent, and \(b_i\) already have fixed polynomial size bounds. Their actual values, including the centers and \(b_i\), may be computed after \(k\) is fixed. Each remaining smallness condition has the form \(p(n,\log\lambda)\lambda^{-\alpha}\) bounded by a specified inverse polynomial in \(n\), with fixed \(\alpha>0\). The endpoint brackets for hopping calibration have the same fixed-power character. A sufficiently large fixed \(A\) therefore meets every asymptotic requirement; increasing \(C\) covers all \(n\ge2\) and the fixed separation and gap conditions. The constants whose existence was proved in the preceding estimates can be dominated by absolute constants in these choices.

The uniformity of this argument depends on the bounds already proved: the polynomials \(p(n,D)\) come from fixed-order derivatives and moments, fixed-dimensional wells and gadgets, and polynomially many matrix entries. The large box is treated separately. Its effects were bounded in [eq:5], by orbital tails, and by the summable fourth derivatives in the cubature argument. Its total mass enters only the much smaller node-rounding error. Thus no exponent needs to grow with the input, and \(H\), the box volume, and \(1/\delta_R\) are all polynomial in \(n\) after the fixed choice of \(A\).

Numerical evaluation and hopping calibration

All required real quantities can be evaluated to inverse-polynomial absolute accuracy. This includes the large scalar shift in [eq:10]: since \(C_H\) is multiplied by \(\lambda^2m\), evaluating it within \[\frac{\tau^2\gamma}{128\lambda m}\] uses at most one quarter of the allowed threshold error \(\lambda\tau^2\gamma/32\). The other scalar terms can be allocated the remaining budget in the same way, because their amplification factors and their number are polynomially bounded. The magnitude \(|C_H|\le CSH\) therefore causes no demand for exponentially fine absolute precision.

We give elementary integration procedures for the quantities used in the construction. In the Gaussian heat formula for \(\phi\) and its derivatives, substitute \(t=s^2\) and rescale the Gaussian integration variable. For a derivative \(D^\alpha\) of any required order, the result is \[D^\alpha\phi(r)=\frac1\pi \int_0^\infty\!\int_{\mathbb R^2} 2s e^{-s^2}e^{-|a|^2}D^\alpha f(r+2sa)\,da\,ds.\] This integrand is regular at \(s=0\). The needed derivatives of \(f\) are bounded piecewise polynomials with bounded first derivatives. Polynomial integration cutoffs make the Gaussian tails smaller than any prescribed inverse-polynomial error; on the remaining domain the integrand and its Lipschitz bound are polynomially bounded. Elementary grid sums in this fixed dimension therefore suffice. Division by \(\phi\) in \(W\) and its derivatives occurs only on the unit disk, where \(\phi\) has a fixed positive lower bound. The other decaying orbital integrals, including normalization and hopping integrals, can likewise be truncated about their centers with polynomial cutoffs.

For Coulomb integrations, including \(v_{ij}\) and \(C_H\), replace \(|x|^{-1}\) inside radius \(\varepsilon\) by \(\varepsilon^{-1}\), with \(\varepsilon\) an appropriately small inverse polynomial. Integrating the local singularity against bounded densities gives an error of order \(\varepsilon^2\) times the relevant polynomial volume bounds. The capped kernel has Lipschitz constant at most \(\varepsilon^{-2}\), so grid sums again give the required accuracy in polynomial time. All integrations have fixed dimension and fixed nesting depth. Exponentials, logarithms, square roots, and the other elementary constants can also be evaluated by rational approximations to the required accuracy; the arguments and reciprocal absolute tolerances are polynomially bounded. Square roots at zero can be handled by bisection. Choosing guard precision for the nested evaluations keeps their accumulated errors within the allotted budgets.

The link hopping equations have strict endpoint brackets, from exponential decay at the two prescribed distance scales and polynomial upper and inverse lower bounds on their targets. In approximate bisection, evaluate the residual to a small fraction of its per-link tolerance. If the sign is ambiguous, the residual is already small enough. Otherwise retain the sign-change bracket until its width times the Lipschitz bound is within tolerance. No derivative lower bound is needed. The fixed layout contractions and subsequent rational center approximations can be computed more accurately than these tolerances require.

Flow evaluation and rational coordinates

Once the rational centers \(u_i\) are fixed, their exact real coefficients \(b_i\) define a reference density and a reference flow \(G\). We approximate the images under that flow; numerical evaluation of \(b_i\) is part of the velocity evaluation error. Thus all numerical and final rational rounding errors can be charged to the single displacement budget \(\delta_R\), without changing the reference density or any cell mass.

The flow lasts for time one and has uniformly bounded velocities and first derivatives in space and time. Euler stepping with an inverse-polynomial time mesh therefore suffices. More explicitly, for mesh \(\nu\) (the reciprocal of the number of steps), initial error \(a_0\), per-step arithmetic error \(a_1\), and velocity evaluation error \(\epsilon_v\), the final position error is \[O(a_0+\nu+\epsilon_v+a_1/\nu).\] Indeed a step multiplies the preceding error by at most \(1+C\nu\) and adds \(O(\nu^2+\nu\epsilon_v+a_1)\). Each parameter can be chosen inverse polynomially small so that the final error is at most \(\delta_R\). Rounding intermediate positions on sufficiently fine rational grids keeps their bit lengths polynomial. Near support boundaries, the specified smooth vanishing and derivative bounds give the same error control. None of these operations uses a many-electron wavefunction.

Output size and encoding

There are \(8|\Omega|h^{-3}\) point nuclei, polynomially many in \(n\), and every nucleus has charge one. Their rational physical positions, the unary electron number, and the rational thresholds from Section 7.3 can all be computed with polynomial encoding lengths. Reducing rational outputs to lowest terms gives the specified input format. To define the map on every input, one may first check the source syntax and the fixed polynomial bounds used in the construction, and output a fixed instance when they fail; this does not affect either promised implication.

The deterministic construction thus supplies the outputs whose energy comparison and threshold margins were proved in [eq:10] and the following inequalities. With the electronic-energy convention that omits nuclear repulsion, this completes the proof of Theorem 1, with \(Z_{\max}=1\) and \(c=1\).