QMA-hardness of continuum Coulomb energy with unit nuclear charges

Introduction

An electronic ground-energy problem specifies both an external potential and the states over which the energy is minimized. These choices matter for computational hardness. A finite collection of orbitals can encode a hard spin system, yet an arbitrary continuum state may have lower energy than every state in that collection. Likewise, a freely designed potential can create a useful array of wells without being the potential of positive point nuclei. We prove hardness when the only external data are rational positions of nuclei of charge one, and the minimization is over the full spinful fermionic continuum space.

Schuch and Verstraete proved QMA-completeness for a Hubbard model with local magnetic fields and developed a continuum realization using designed scalar and spin-dependent potentials [10]. Their reduction to Heisenberg exchange and their comparison with the complementary one-particle space are predecessors of the finite and continuum estimates used here. Requiring the external potential to come solely from positive point nuclei imposes a further restriction.

O’Gorman, Irani, Whitfield, and Fefferman proved QMA-completeness for electronic structure in a supplied finite orbital basis [6 Theorem 1]. Their journal version includes a construction with positive unit point charges [6 Appendix D, Corollary 1]. The supplied basis restricts the admissible many-electron states, as the authors emphasize in their discussion of the infinite-dimensional problem in Section V. A projected energy supplies a variational upper bound on the continuum infimum; transferring a NO instance also requires a lower bound on every omitted continuum state. The present theorem retains both unit nuclei and unrestricted spinful continuum minimization.

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

Its Hilbert space and form domain are

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 Theorems 2.13 and 2.19] 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 (H) defines the electronic-energy convention used throughout.

Input or convention	Specification
Nuclear positions	Pairwise distinct rational triples $R_\alpha$
Nuclear charges	$Z_\alpha= 1$ for every $\alpha$
Electron number	$N \ge1$, encoded in unary
Thresholds	Rational numbers $a < b$
Rational encoding	Signed binary numerator and positive binary denominator in lowest terms; fixed self-delimiting encoding
Input length	$L$: total length of the binary and unary data
Promise	$b-a \ge L^{-1}$; YES if $E_0 \le a$, NO if $E_0 \ge b$
Energy convention	Electronic energy; clamped nuclear repulsion omitted

Theorem 1 (Unit-charge Coulomb hardness). The electronic promise problem for (H) on $\mathcal{H}_N$, with distinct rational nuclear positions, all nuclear charges equal to one, unary electron number, and rational thresholds separated by at least $L^{-1}$, is QMA-hard under deterministic classical polynomial-time many-one reductions. The reduction produces polynomially many nuclei and electrons, polynomial rational bit lengths, and threshold separation at least one. The only output data are the nuclear positions, electron number, and thresholds.

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 Theorem 1.1] 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 Theorem 1.1 and Section 2] 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 Theorem 48, Lemma 47, and Section 10.1]. 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 Section 2.4 and Lemma 5], 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.

The lattice source and rational normalization

We first fix the finite input and its error budget. The output of this section is a signed rational Heisenberg Hamiltonian with polynomial geometric extent and coefficients, together with a positive inverse-polynomial decision gap.

Our starting problem is the field-free square-lattice Heisenberg energy problem, which is QMA-hard by Theorem 48 and Section 10.1 of Cubitt, Montanaro and Piddock [3], using the polynomial-weight spatially sparse construction of their Lemma 47. The Hamiltonian on a finite, nonwrapping square-lattice subgraph is

on spins $1/2$, where the components are Pauli matrices and there are no one-site terms. YES instances have ground energy at most the lower threshold, and NO instances have ground energy at least the upper threshold. The cited construction supplies polynomial bounds on lattice size and coefficient magnitudes and an inverse-polynomial promise gap. The following deterministic normalizations give signed rational coefficients, rational thresholds $a_s,b_s$, and polynomial geometric extent.

For completeness, let $\delta> 0$ be a known rational inverse-polynomial lower bound for the incoming promise gap, and let $q$ be the number of weighted edges. Set $\epsilon_0 = \delta/100$. Rounding each weight within $\epsilon_0/\max(1,q)$ changes the Hamiltonian norm by at most $3\epsilon_0$. Write the incoming scalar as $c_0 I$, with $c_0 = 0$ if absent, and the incoming thresholds as $a_0,b_0$. Approximate these three numbers within $\epsilon_0$, obtaining $\widetilde{c}_0,\widetilde{a}_0,\widetilde{b}_0$. The rational Hamiltonian without the scalar then has valid thresholds

The offsets allow $3\epsilon_0$ for the Hamiltonian, $\epsilon_0$ for the scalar, and $\epsilon_0$ for the threshold. Thus both promise implications survive, and $b_s-a_s \ge\delta- 12\epsilon_0 > 0$. Polynomial magnitudes and inverse-polynomial accuracy require only polynomial rational bit lengths. A connected square-grid component with $V$ vertices has coordinate diameter at most $V-1$, by following lattice paths. Translating components into disjoint strips therefore gives polynomial extent without changing their weighted graphs.

Put $n = \max(2,\text{source input size})$, $\gamma= \min(1,b_s-a_s)$. All polynomial choices below can use fixed size, magnitude and inverse-gap bounds for this source. We can translate coordinates to keep them polynomially bounded, omit zero weights, and assume a nonempty vertex set (add a free spin if needed).

The quantity $\gamma$ is the error budget used in the finite reductions below. It concerns the normalized spin instance; the final electronic thresholds will be obtained from an explicitly computed affine energy map in Section 7.3.

Positive spin model

We convert the signed lattice Hamiltonian into a positive-exchange model whose interaction strengths can be encoded by distances. Singlet-pair mediators implement the sign conversion; fixed planar templates then make the resulting link lengths independently adjustable. The same-member versus opposite-member sign mechanism is the isotropic specialization of the mediators in Piddock and Montanaro [9 Section 2.4 and Lemma 5]. We keep the quantitative proof here, including the scalar shifts and the terms unchanged during subdivision.

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 Section 3.2]; we prove precisely the bottom-energy statement needed here. In particular, the first-order block $PBP$ need not vanish.

Lemma 2 (Second-order bottom-energy estimate). Let $A, B$ be self-adjoint operators on a finite-dimensional Hilbert space. Suppose $A$ has kernel projection $P \ne0$, $A \ge gQ$, where $Q = I - P$ and $g > 0$, and $\lVert B\rVert\le g/4$. Then

where the inverse vanishes on the kernel.

Proof. Indeed $|E(A+B)| \le\lVert B\rVert$. At that energy the $Q$-equation can be solved in terms of the nonzero $P$-component of a ground vector. The inverse on $Q$ used there is $(A|_{Q}+QBQ|_{Q}-E(A+B))^{-1}$; it differs from $A^{-1}|_{Q}$ by $O(\lVert B\rVert/g^{2})$ by the resolvent identity. The $P$-equation gives the lower estimate. For the upper estimate take a unit ground vector $p$ for the displayed compression and use $p-A^{-1}Bp$ as trial vector. Its energy numerator is the compressed energy plus $O(\lVert B\rVert^{3}/g^{2})$, and the normalization changes this result by at most the same order. The statement includes $Q=0$. □

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 $\lVert h^{s}_{ij}\rVert=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

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

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 $Jh^{s}_{LR}-(3/2)v_e^2I$ 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_eh^{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

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 (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 $\sqrt{2}$ 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,

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, $\lVert DR\rVert<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.

Localized modes and the finite fermion model

The positive spin model now fixes the target exchange strengths $K_e$. This section constructs smooth localized modes and a finite fermion Hamiltonian whose low-energy spin interaction reproduces $H_+$ up to scale and a scalar shift. Its occupancy penalty retains the intersite direct Coulomb coefficients. Computing the cost of virtual double occupancy determines the required hopping, which we then obtain by choosing the mode separations.

We work in coordinates $y=\lambda x=(r,z)$ with $r\in\mathbb{R}^2$. The common $\lambda$ will be polynomially large. After dividing energies by $\lambda^2$, kinetic energy is the usual $-\Delta/2$, each unit nucleus has coupling $1/\lambda$, and electron repulsion has coefficient $1/\lambda$.

Local modes and Coulomb coefficients

Take $f(r)=(1-|r|^2)_+^{16}$, and fix

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,

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\varphi^2|\nabla(u/\varphi)|^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

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

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 $\lVert R\rVert\le a_*/2$ at our scales, and

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

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$,

Section 5 will produce this one-body term by adjusting the well depths. The finite Hamiltonian to be realized is therefore

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 Eqs. (12)–(13), (18)–(21)]; 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

Indeed the triplets do not hop into same-site pairs, while the singlet hops into two orthogonal doubly occupied states with amplitudes of magnitude $\sqrt{2}t_e$. Apply the creation and annihilation operators to $(c^\dagger_{i\uparrow}c^\dagger_{j\downarrow}-c^\dagger_{i\downarrow}c^\dagger_{j\uparrow})|0\rangle/\sqrt{2}$ 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

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

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(\lvert u_i-u_j\rvert)$ 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

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(\lvert u_i-u_j\rvert)$. Choose the per-link calibration accuracy so that replacing $t_e^{\mathrm{tar}}$ by $t_e$ changes the total Hamiltonian norm by at most $\gamma\tau^2/100$. Together with the second-order remainder, this gives

Polynomial scale convention

Equation (5) 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

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.

Synthesis by a nonnegative charge density

We now construct an attractive Coulomb potential whose low-energy matrix realizes the modes and hopping of Section 4. A broad slab supplies the transverse oscillator near the modes and a positive background density. Adding the planar wells then changes the density by a bounded, compactly supported term, so the total density remains nonnegative. We first construct this potential, then establish a fixed energy cost outside the mode space and compute its matrix on that space.

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

The occupancy identity (C) 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 $\lvert z\rvert\le1/2$; it may be piecewise polynomial with rational data. Define

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 $\lvert\Delta V_\ell/(4\pi)\rvert\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

by compact support and the fundamental solution identity. Thus the density in (4) lies between $\rho/2$ and $3\rho/2$ on $\Omega$. Its total mass is $\rho\lvert\Omega\rvert$, 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$:

Indeed $F_H$ at each $z$ is minimized at $r=0$, and is also nondecreasing in $\lvert z\rvert$ 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 $\lvert z\rvert\le S$ exactly $2\pi\rho\int_{-S}^{S}(\lvert z-w\rvert-\lvert w\rvert)\,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 $\lvert r\rvert\le H/2$, the change relative to the horizontal center is bounded by $CS\lvert r\rvert^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 $\lvert C_H\rvert\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

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 $\lVert\Gamma-I\rVert\le p(n,D)\lambda^{-0.8}$: for example, each off-diagonal overlap is $O(e^{-0.9\lvert u_i-u_j\rvert})$, by leaving part of the decay as an integrable factor. We use the orthonormal modes

In particular, unit vectors in $\operatorname{Ran} P$ have polynomially bounded one-particle kinetic energy.

Proposition 3 (Continuous one-particle control). For the density and separated centers constructed above, with the polynomial scale chosen sufficiently large, there are fixed constants $c>0,C<\infty$ such that

The inequalities are quadratic-form inequalities on the one-particle $H^1$ domain, and the last bound is an operator norm from the range of $P$ to $L^2$. All estimates hold with either spin included.

Proof. For the first inequality compare below, using (5), with the separable planar multiwell $-\Delta_r/2+\sum W_i$ plus the vertical operator with capped potential $2\pi\rho\min(z^2,S^2)$, shifted by $-e_*$. Errors in the lower bound are $O(S^3/H)+O(p(n,D)/\lambda)$, since $\zeta W_i\ge W_i$ and the $b_i/\lambda$ terms are small and bounded everywhere.

For clarity we detail the gap argument. In the plane use a quadratic partition of unity consisting of cutoffs about each site (equal to one inside radius $D/8$, supported at radius $D/3$) and a remaining exterior function. They can be built by smooth angle interpolation from $0$ to $\pi/2$ on the transition annuli using sine and cosine, with squared gradients summing to $C/D^2$ or less. The localization form identity (IMS localization; compare [11 Lemma 3.1]), obtained here by the product rule with kinetic coefficient $1/2$, then costs $O(1/D^2)$. In each site term use the isolated ground and gap bound, and in the exterior use the free kinetic lower bound. This gives

where $P_r$ projects on the planar translates of $\varphi$, $c_1>0$ can be taken below the isolated gap and $1/2$, and the small bound can be made uniformly small. In detail, with site cutoffs $\theta_i$ and $\varphi_i=\varphi(\cdot-u_i)$, the site term evaluated on $\theta_i u$ is at least

the exterior term also bounds its share of the $-1/2+c_1$ term as there is no well on its support. The sum of rank-one operators from the cutoff ground-state modes differs from $P_r$ by at most $O(me^{-c'D})+O(p(n,D)\lambda^{-0.8})$, by the tails and Gram bound.

Similarly localize the vertical problem with an inner cutoff equal to one on $\lvert z\rvert\le S/2$, supported on $\lvert z\rvert\le S$, and an exterior function. Use the harmonic gap and the large $z$ potential lower bound on the respective terms, at gradient cost $O(S^{-2})$, with Gaussian loss on the rank-one mode. This gives a vertical lower bound $\omega/2+c_2(I-P_g)-o(1)I$, where $P_g=\lvert g\rangle\langle g\rvert$. Since $P=P_r\otimes P_g$ on the spatial space and

the tensor sum has a fixed gap on $I-P$ once the errors are small.

The slowly small errors such as $O(D^{-2})$ here are needed only for a fixed gap estimate, not as the energy accuracy per particle.

The kinetic lower bound in (8) follows from (7) or background minimality and bounded $V_\ell$. For the residual estimate, (7), the crude global background bound, and the mode tails give

Indeed moments with $\lvert r\rvert\le H/2$, $\lvert z\rvert\le S$ suffice there, and outside this region the Gaussian or planar exponential beats the global polynomial magnitudes, using $S\ge\log\lambda$. The cutoff defect in $V_\ell$ versus $\sum_i(1+b_i/\lambda)W_i$ acting on any mode is likewise negligible. Since $(-\Delta_y/2+2\pi\rho z^2+W_j-e_*)\psi_j=0$, the remaining residual is a sum of other wells acting on the small planar tails and the term with $b_j/\lambda$. This proves (6). $\square$

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

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(\lvert u_i-u_j\rvert)$. 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

by (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

The preceding estimates give

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.

Replacing the density by unit nuclei

Our input is the nonnegative density in (6). We must replace its Coulomb field by finitely many equal masses $1/\lambda$, so that physical nuclear charges are all one. The construction has two distinct error requirements: a form bound on arbitrary states and a stronger compression bound on the smooth localized modes.

The comparison with $H_d$ uses energies multiplied by $\lambda$. A compressed error must therefore be smaller than $\lambda^{-1}$, up to the required polynomial margin. The error coupling the mode space to its complement can be larger: the fixed complementary gap makes its energy contribution quadratic, as proved in Section 7.1. With $h=(8/(\rho\lambda))^{1/3}$, a cubature rule exact through degree three will give an $O(h^4)$ compressed error and an $O(h^2)$ form error. Both resulting energy errors are $O(\lambda^{-1/3})$ in the units of $H_d$, up to polynomial factors.

Transport and charge placement

The density interpolation uses Moser’s deformation method [5 Section 4], also developed by Dacorogna and Moser [4 Section II, Lemma 3]. 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\geq\rho/2$, the time-dependent vector field on full space

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_tq_t+\nabla\cdot(q_t\mathcal{V}_t)=0$ gives

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\sqrt{3})$. 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

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_*$.

Proposition 4 (Unit-nucleus replacement). For the transported degree-three cubature and rational rounding just specified, with $\delta_R \le(1+H)^{-3}(1+S)^{-2}\lambda^{-3}$, the shifted continuous and point-nucleus operators $L_c,L_p$ satisfy

The first estimate holds for every one-particle $H^1$ function $u$. The second is the norm of the compressed form on the localized-mode space $\operatorname{Ran} P$. Both bounds are in dilated energy units.

Proof. We first work at the exact transported nodes. On an original grid cube $Q$, write $a_{Q,\nu}$, $1\le\nu\le8$, for its Gauss nodes. Taylor’s theorem and exactness through degree three give the common error estimate

whenever $F\circ G$ has four bounded derivatives on $Q$. We apply this first to the singular Coulomb kernel away from its singularity, then to the smooth fields obtained after compression.

General form estimate. For fixed evaluation point $y$, the charge integral is of the kernel $|y-G(b)|^{-1}$ against uniform density in $b$. For cubes whose images touch a ball of radius a fixed large constant times $h$ about $y$, there are only $O(1)$ such cubes. Their continuous contribution is $O(h^2)$ by bounded transformed density. The discrete one is bounded by a sum of $Ch^3/|y-b_\alpha|$ with $|y-b_\alpha|\le C'h$.

On farther cubes, the pure fourth derivative of the kernel is bounded by $C|y-G(b)|^{-5}$. Distances to $y$ are comparable within each such cube, and the bi-Lipschitz change of variables has bounded Jacobian. Summing the pure derivative contributions therefore costs at most

The other chain-rule terms contain a derivative of $G$ of order at least two and a kernel derivative bounded by $C|y-G(b)|^{-j}$, $2\le j\le4$. Such terms occur only on cells meeting the well cylinders, of combined volume $O(m(S+1))$. Below distance one they are bounded by the same inverse fifth power; beyond distance one their contribution is at most $Ch^4m(S+1)$.

This proves a pointwise absolute error bound by $C(h^2+h^4m(S+1))$ plus the close-node sum above. In particular there is no factor the size of the large horizontal area in this bound.

The close-node balls have bounded overlap also upon any fixed dilation. Around each node one has

This follows by a cutoff and the three-dimensional Hardy inequality $\int|u|^2/|y|^2\le4\int|\nabla u|^2$, which follows also by the divergence formula for $y/|y|^2$ and Cauchy–Schwarz. Summing gives the general form bound, since $h^2=O(\lambda^{-2/3})$.

Compressed estimate. For compression use the smooth Coulomb fields

Their derivatives through order four are uniformly bounded (differentiate the translated densities). Fourth derivatives satisfy

Indeed derivatives of the product density decay exponentially about $(u_i,0)$, uniformly by bounded derivatives on the other factor. At a large distance $d=\lvert b-(u_i,0)\rvert$, use a smooth density cutoff supported within distance $d/2$ of that center and equal to one inside $d/3$, with the cutoff radius held fixed when differentiating. There one uses pure kernel derivative decay, and on the remaining tail one differentiates the cut-off density instead, with exponentially small bounds. This proves the assertion. Applying the same degree-three cubature to $K_{ij}\circ G$ now costs $O(h^4(1+m(S+1)))$ for each entry: pure fourth derivatives have summable volume bounds as just proved, and other chain-rule terms live on the exceptional cells. Gram correction and taking matrix norms cost at most polynomial factors. This proves the second estimate for exact nodes.

Rounding. A center shift of size $\delta_R$ changes the inverse distance by at most $\delta_R/(\lvert y-b_\alpha\rvert\lvert y-\widetilde b_\alpha\rvert)$. Cauchy–Schwarz and translated Hardy inequalities bound its form by $4\delta_R\lVert\nabla u\rVert^2$. Summing the masses gives the total bound $4\rho\lvert\Omega\rvert\delta_R\lVert\nabla u\rVert^2$, where

This is smaller than both required errors, also after compression because the modes have polynomial kinetic bounds. It completes (7).

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.

Many-electron continuum estimate

The point nuclei have now been fixed. We must compare their electronic energy with the finite Hamiltonian $H_d$ on all continuum states. First we eliminate the space in which some electron lies outside the localized modes, using the fixed gap and the two replacement estimates. Then we identify the compressed Hamiltonian, keeping direct Coulomb repulsion exactly at its leading order.

In the dilated coordinates the full operator after dividing by $\lambda^2$ and shifting by $-m(C_H+e_*)$ is

We work on the antisymmetric space with spins and its $H^1$ form domain. The Coulomb forms used here have the variational meaning specified in Section 1, whether or not their energy infima are attained. Indeed, fiberwise Hardy and Cauchy–Schwarz give, for every $\eta>0$ and every fixed center $b$,

The same estimate applies to a pair singularity by fixing the other coordinates. For each finite nuclear configuration, summing with sufficiently small $\eta$ gives a closed lower-bounded form whose shifted form norm is equivalent to the $H^1$ norm. Smooth compactly supported antisymmetric functions form a core. The spectral infimum is therefore the infimum of the form over normalized $H^1$ states, and the form estimates below extend from this core by continuity.

Let $\Pi$ project onto the space in which all $m$ particles lie in $P$, and write

The smooth decaying modes lie in the form domain, so their compression is well defined. The one-particle kinetic bound gives $\lVert w\rVert_{1} \le p(n,D)\lVert w\rVert$ for $w \in\operatorname{Ran}\Pi$.

Let

where the compression is defined by the form. This is a finite-dimensional minimum because $\Pi$ projects onto the $m$-electron space built from $2m$ one-particle spin modes.

Control outside the orbital space

Proposition 5 (Eliminating the continuum complement). For the unit-nucleus construction above with $N=m$, and sufficiently large polynomial scale $\lambda$, the full antisymmetric $m$-electron operator in dilated and shifted units satisfies

The lower bound holds on every normalized antisymmetric $H^{1}$ state, including every spin sector.

Proof. The proof needs a lower bound on the complementary space, a bound on its coupling to $\operatorname{Ran}\Pi$, and then a square completion. Write

By (6), as one-particle forms,

for some $c'>0$: the diagonal and cross terms involving $P$ are controlled by $\lVert L_{c}P\rVert$. Since $\sum_{\ell}(I-P)_{\ell}\ge I-\Pi$, summing gives a fixed positive lower bound for $\mathcal{H}^{0}_{c}$ on $\operatorname{Ran}(I-\Pi)$ once $mp(n,D)\lambda^{-0.9}$ is small. We also have the kinetic lower bound $\mathcal{H}^{0}_{c}\ge-\sum_{\ell}\Delta_{\ell}/2-Cm$. A convex combination, with a sufficiently small weight of order $1/m$ on the kinetic bound, retains a positive $L^{2}$ term and gives

for $q\in\operatorname{Ran}(I-\Pi)$. By (7) on each coordinate and positivity of $\mathcal{R}$, this estimate persists with half that coefficient for $\mathcal{H}$.

For $w\in\operatorname{Ran}\Pi$, equation (6) gives $\lVert\mathcal{H}^{0}_{c}w\rVert\le p(n,D)\lambda^{-0.9}\lVert w\rVert$. Hardy applied with one particle coordinate at a time also gives $\lVert\mathcal{R}w\rVert\le p(n,D)\lVert w\rVert$. Finally, polarizing the form bound in (7) and using the $H^{1}$ bound on $w$ controls the point-nucleus replacement in the mixed block. Together these estimates give

The same estimates show that the compressed form has norm less than $c''/(4m)$ at sufficiently large scale. Set $\alpha=c''/(4m)$ and let $\beta=p(n,D)\lambda^{-2/3}$ bound the mixed block. Decompose any normalized form-domain vector as $u=w+q$, where $w=\Pi u$ and $q=(I-\Pi)u$. Subtracting $E_{\Pi}$ makes the $w$-diagonal form nonnegative and leaves at least $\alpha\lVert q\rVert_{1}^{2}$ on the complement. Hence, with all pairings interpreted as forms,

This proves the lower bound in (8); the upper bound follows by testing a minimizing vector in $\operatorname{Ran}\Pi$. The lower-bound argument applies to every normalized antisymmetric state, in any spin sector, without any binding or ground-eigenvector assumption.

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 (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: $\lVert\psi_i\psi_k\rVert_{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

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 (C) combines these terms into the first line of (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.

Proposition 6 (Continuum-to-fermion energy comparison). For the constructed unit nuclei with $N=m$, let $\mathcal{H}$ be the dilated and shifted continuum operator above and let $H_d$ be the finite fermion Hamiltonian (2). At sufficiently large polynomial scale,

Proof. The one- and two-body matrix estimates just obtained compare $\lambda E_{\Pi}$ with $E(H_d)$ within $p(n,D)\lambda^{-0.1}$. Multiplying (8) by $\lambda$ contributes only $p(n,D)\lambda^{-1/3}$, which fits within that bound. □

This analysis treats long-range diagonal effects and possible charge transfer at the virtual-excitation scale in (2)–(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$.

Estimate	Error scale	Role in the comparison
General form replacement, (7)	$\lambda^{-2/3}$	Controls mixed blocks in the $H^1$ norm
Compressed replacement, (7)	$\lambda^{-4/3}$	Direct error on the mode space
Complement elimination, (8)	$\lambda^{-4/3}$	Square of the mixed-block scale
Last two errors after multiplication by $\lambda$	$\lambda^{-1/3}$	Error in the units of $H_d$
Total comparison, (9)	$\lambda^{-0.1}$	Includes the other mode-matrix errors

Physical energies and thresholds

Choose the polynomial scale in $\lambda$ large enough that (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

Combining (1), (3), and (9) now gives a total error at most $7\lambda\tau^2\gamma/100$, and hence at most $\lambda\tau^2\gamma/8$, from

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,

For a NO source instance, the reversed estimate gives

Since $b_s-a_s\ge\gamma$, the threshold separation satisfies

Thus both promises survive rational rounding, and the reduction meets the required output separation $b-a\ge1$.

Implementability and encoding

The preceding construction specifies the nuclear configuration, and Section 7.3 gives threshold accuracies that preserve the source promise. It remains to compute those outputs in deterministic polynomial time. We first fix the order of the scales, then explain how to evaluate the orbital integrals and transported nodes to the required accuracy.

Order and uniformity of the scale choices

Fix the bounded well data and the integer $\rho$ using the positivity bound in (6), 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

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 (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 (10): since $C_H$ is multiplied by $\lambda^2 m$, evaluating it within

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 C_S H$ 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

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 $\epsilon$ by $\epsilon^{-1}$, with $\epsilon$ an appropriately small inverse polynomial. Integrating the local singularity against bounded densities gives an error of order $\epsilon^2$ times the relevant polynomial volume bounds. The capped kernel has Lipschitz constant at most $\epsilon^{-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

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 (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$.

References