==== TABLE OF CONTENTS ====
Continuum-Coulomb-hardness-with-binary-nuclear-charges//00-introduction.tex
Continuum-Coulomb-hardness-with-binary-nuclear-charges//01-finite-models.tex
Continuum-Coulomb-hardness-with-binary-nuclear-charges//02-wells.tex
Continuum-Coulomb-hardness-with-binary-nuclear-charges//03-hopping.tex
Continuum-Coulomb-hardness-with-binary-nuclear-charges//04-continuum.tex
Continuum-Coulomb-hardness-with-binary-nuclear-charges//05-completion.tex
Continuum-Coulomb-hardness-with-binary-nuclear-charges//paper.tex
Continuum-Coulomb-hardness-with-binary-nuclear-charges//preamble.tex
Continuum-Coulomb-hardness-with-binary-nuclear-charges//references.bib
===========================
==== BEGIN Continuum-Coulomb-hardness-with-binary-nuclear-charges//00-introduction.tex ====
\section{The problem and the result}\label{sec:introduction}
The exact molecular electronic Hamiltonian has unusually rigid input data: its external potential is generated by a list of positive nuclei. A hardness proof for a matrix obtained by projecting onto supplied orbitals does not by itself establish hardness for this operator. Projection provides an upper bound on the ground energy, while a reduction also needs a lower bound that excludes all omitted continuum states. We give such a reduction for the integer-charge, rational-position encoding below.
\subsection{The precise promise problem}
Fix a positive integer $c$. An input consists of distinct positions $R_1,\ldots,R_M\in\Q^3$, positive integers $Z_1,\ldots,Z_M$, an electron number $N\ge1$, and rational thresholds $a0$ are independent of the large scale parameters and may change between estimates. A polynomial bound always means a fixed polynomial with uniform coefficients and degree. Exponentially small bounds are stated explicitly when they are later amplified. All operator inequalities involving unbounded Hamiltonians are quadratic-form inequalities unless their bounded compressions are specified.
==== END Continuum-Coulomb-hardness-with-binary-nuclear-charges//00-introduction.tex ====
==== BEGIN Continuum-Coulomb-hardness-with-binary-nuclear-charges//01-finite-models.tex ====
\section{Finite models and a controllable contact geometry}
\label{sec:finite}
We first reduce a signed spin Hamiltonian to a half-filled Hubbard model.
The geometric construction is part of this reduction: every Hubbard
hopping will later be set by a small displacement of a nuclear site.
For a finite-dimensional self-adjoint operator $A$, write $E(A)$ for its
smallest eigenvalue. On two spin-$1/2$ sites set
\[
h^s_{ij}=X_iX_j+Y_iY_j+Z_iZ_j
=\boldsymbol\sigma_i\cdot\boldsymbol\sigma_j.
\]
Its eigenvalues are $-3$ on the singlet and $1$ on the triplet, so
$\lVert h^s_{ij}\rVert=3$.
\Needspace{14\baselineskip}
\subsection{The source problem and a quantitative perturbation bound}
\begin{theorem}[Source problem]\label{thm:source}
There is a QMA-complete promise problem whose Hamiltonians are
\[
H_{\rm in}=\sum_{e=LR}J_e h^s_{LR}
\]
on finite, nonwrapping subgraphs of the square lattice, with no
one-qubit terms. The signed weights $J_e$ are rational. The number of
vertices and edges, the lattice extent, the weight magnitudes, and the
threshold magnitudes are bounded by fixed polynomials in the source
size $n\geq2$. The rational thresholds $a_s0$ be a known rational inverse-polynomial lower bound on the
incoming gap, let $q_s$ be the number of weighted edges, and put
$\varepsilon_0=\delta/100$. Round each edge coefficient within
$\varepsilon_0/\max\{1,q_s\}$. Since $\|h^s_{ij}\|=3$, the resulting
Hamiltonian changes by norm at most $3\varepsilon_0$. If the source also
has an identity coefficient $c_0$ and thresholds $a_0,b_0$, approximate
each of those three numbers within $\varepsilon_0$, obtaining
$\widetilde c_0,\widetilde a_0,\widetilde b_0$. For the rounded Hamiltonian
with its identity term removed, set
\[
a_s=\widetilde a_0-\widetilde c_0+5\varepsilon_0,
\qquad
b_s=\widetilde b_0-\widetilde c_0-5\varepsilon_0.
\]
The two threshold approximation errors and the Hamiltonian error sum to
at most $5\varepsilon_0$, so both promises survive. Moreover,
$b_s-a_s\ge\delta-12\varepsilon_0>0$. All rational descriptions have
polynomial length. A connected grid component with $v$ vertices has
coordinate diameter at most $v-1$, by following its lattice paths;
translating components into disjoint strips gives polynomial extent.
These are normalization steps applied to the cited source, not additional
conclusions imported from its statement.
This imported source theorem is the only hardness theorem we need.
In what follows zero-weight edges are omitted. No lower bound
on the magnitudes of the remaining weights is assumed, and no total-spin
sector is selected.
By appending an isolated original vertex if necessary, at a lattice
position outside the existing finite drawing, we may assume that there
is at least one original vertex. Tensoring with its uncoupled spin
preserves the minimum energy and all polynomial bounds.
We use the following finite-dimensional estimate twice. It controls the
bottom energy and does not require a gap within the encoded subspace.
The effective second-order operator is the familiar one from
Schrieffer--Wolff perturbation theory; see
Bravyi, DiVincenzo, and Loss~\cite[Section~3.2]{BravyiDiVincenzoLoss2011}.
We prove the particular bound needed here directly, so its constants and
the empty-complement case are explicit.
\begin{lemma}[Second-order bottom-energy estimate]
\label{lem:second-order}
Let $A\geq0$ be a finite-dimensional self-adjoint operator, let
$P\neq0$ project onto its kernel, and put $Q=I-P$. Suppose
$A\geq gQ$ for some $g>0$. Let $B=B^*$ satisfy $PBP=0$ and
$\lVert B\rVert\leq\epsilon g$, where $0\leq\epsilon<1/2$.
Define $A^{-1}$ to vanish on $P$. Then
\begin{equation}\label{eq:second-order-bound}
\left|E(A+B)-E\bigl(-PBA^{-1}BP\big|_{P}\bigr)\right|
\leq 2g\epsilon^3.
\end{equation}
The conclusion also holds when $Q=0$.
\end{lemma}
\begin{proof}
Let $K=PBA^{-1}BP\geq0$ and $\kappa=\lVert K\rVert$.
Since $Bp\in Q$ for $p\in P$,
\begin{equation}\label{eq:second-order-size}
0\leq\kappa\leq\frac{\lVert B\rVert^2}{g}
\leq g\epsilon^2.
\end{equation}
For $p\in P$ and $q\in Q$, the inequality
$QBQ\geq-\epsilon gQ\geq-\epsilon A|_Q$ gives
\begin{align*}
\langle p+q,(A+B)(p+q)\rangle
&\geq (1-\epsilon)\langle q,Aq\rangle
+2\operatorname{Re}\langle q,Bp\rangle\\
&= (1-\epsilon)
\left\|A^{1/2}q+\frac{A^{-1/2}Bp}{1-\epsilon}\right\|^2
-\frac{\langle p,Kp\rangle}{1-\epsilon}.
\end{align*}
On a unit vector this is at least $-\kappa/(1-\epsilon)$.
Consequently
\[
E(A+B)+\kappa
\geq-\frac{\epsilon\kappa}{1-\epsilon}
\geq-2g\epsilon^3.
\]
For the other direction choose a unit $p\in P$ with $Kp=\kappa p$
and set $q=-A^{-1}Bp$. Then $\lVert q\rVert\leq\epsilon$ and
\[
\langle p+q,(A+B)(p+q)\rangle
=-\kappa+\langle q,Bq\rangle.
\]
Writing $r=\lVert q\rVert^2$, the variational principle therefore gives
\[
E(A+B)+\kappa
\leq \frac{\langle q,Bq\rangle+\kappa r}{1+r}
\leq g\epsilon^3+g\epsilon^4
\leq\tfrac32g\epsilon^3.
\]
The smallest eigenvalue of $-K$ on $P$ is $-\kappa$, proving the
claim. If $Q=0$, the hypothesis $PBP=0$ forces $B=0$, and the same
formulas hold with all inverse terms zero.
\end{proof}
\subsection{Replacing signed bonds by positive Heisenberg bonds}
Piddock and Montanaro's singlet mediator construction changes the sign of
an effective interaction by coupling to the same or to opposite members
of a strongly bound pair~\cite[Section~2.4 and Lemma~5]{PiddockMontanaro2017}.
We use its isotropic specialization below, keeping the scalar shift and
the simultaneous-gadget error explicit. This stage produces only positive
spin bonds; the next geometric stage assigns each of them a controllable
unit contact.
For every source edge $e=LR$ introduce two ancillary spins $A_e,B_e$.
There are $k$ such pairs. Define
\begin{equation}\label{eq:positive-gadget}
\begin{split}
v_e&=\sqrt{2|J_e|},\qquad
C_e=\begin{cases}A_e,&J_e<0,\\ B_e,&J_e>0,\end{cases}\\
H_+&=\Delta\sum_e h^s_{A_eB_e}
+\sqrt\Delta\sum_e v_e
\bigl(h^s_{LA_e}+h^s_{RC_e}\bigr).
\end{split}
\end{equation}
Every nonzero coefficient in $H_+$ is positive. Put
\begin{equation}\label{eq:spin-shift}
V=\sum_e v_e,\qquad
C_0=3k\Delta+\frac32\sum_e v_e^2.
\end{equation}
\begin{proposition}[Positive-spin reduction]\label{prop:positive-spin}
If $\sqrt\Delta>3V$, then
\begin{equation}\label{eq:positive-spin-error}
\bigl|E(H_+)-E(H_{\rm in})+C_0\bigr|
\leq\frac{27V^3}{\sqrt\Delta}.
\end{equation}
In particular, a polynomially large integer $\Delta$ makes this error
an arbitrarily small prescribed inverse polynomial.
\end{proposition}
\begin{proof}
Take
\[
A=\Delta\sum_e(h^s_{A_eB_e}+3I),\qquad
B=\sqrt\Delta\sum_e v_e(h^s_{LA_e}+h^s_{RC_e}).
\]
The kernel of $A$ consists of arbitrary original spins tensored with a
singlet on each ancillary pair. Denote its projection by $P_s$.
The gap is $g=4\Delta$, and
\[
P_sBP_s=0,\qquad
\lVert B\rVert\leq6\sqrt\Delta V.
\]
A weak interaction from pair $e$ takes that pair's singlet into its
triplet space and leaves every other ancillary pair in its singlet.
Thus every state reached from $P_s$ in one hop has excitation energy
$4\Delta$. Moreover, terms associated with distinct pairs give zero
mixed second-order contributions: a term acting on pair $f\neq e$
cannot remove the triplet on pair $e$. This remains true when the two
source edges share an original vertex.
For a singlet $|s\rangle$ on $A,B$ and Pauli components $\mu,\nu$,
\[
\langle s|\sigma_A^\mu\sigma_A^\nu|s\rangle
=\delta_{\mu\nu},\qquad
\langle s|\sigma_A^\mu\sigma_B^\nu|s\rangle
=-\delta_{\mu\nu}.
\]
Consequently the squares of the two weak terms each compress to $3I$.
Their two ordered cross terms sum to $2h^s_{LR}$ when $C_e=A_e$,
and to $-2h^s_{LR}$ when $C_e=B_e$. The contribution of pair $e$ to
$-P_sBA^{-1}BP_s$, identified with an operator on the original spins,
is therefore
\begin{equation}\label{eq:spin-gadget-effective}
-\frac{v_e^2}{4}
\begin{cases}
6I+2h^s_{LR},&J_e<0,\\
6I-2h^s_{LR},&J_e>0
\end{cases}
=J_e h^s_{LR}-\frac32v_e^2 I.
\end{equation}
The factor $4$ in the denominator is the singlet--triplet excitation
energy in units of $\Delta$.
Apply Lemma~\ref{lem:second-order} with
$\epsilon=3V/(2\sqrt\Delta)<1/2$. Its error bound is
$2(4\Delta)\epsilon^3=27V^3/\sqrt\Delta$.
Since $A+B=H_++3k\Delta I$, summing
\eqref{eq:spin-gadget-effective} gives the result. If $k=0$, both
perturbation and shift vanish and the assertion is exact.
\end{proof}
Here is one explicit scale choice, which we retain for later bookkeeping.
Set
\begin{equation}\label{eq:spin-scale}
\begin{split}
\gamma&=\min\{1,\delta_s\},\qquad
W=\max\bigl(\{1\}\cup\{|J_e|:e\}\bigr),\\
\overline V&=\max\{1,2kW\},\qquad
\Delta=\left\lceil
\left(\frac{1728\overline V^3}{\gamma}\right)^2
\right\rceil.
\end{split}
\end{equation}
Then $V\leq\overline V$, and
\eqref{eq:positive-spin-error} is at most $\gamma/64$.
Furthermore $v_e\leq\sqrt\Delta$, so every coefficient of $H_+$ is
at most $\Delta$. All these scales have polynomial magnitude in $n$.
Intermediate square roots cause no encoding issue: their nonnegative
values admit polynomial-precision rational approximations by bisection.
\subsection{Contact geometry and independent length control}
Place the original lattice vertices at spacing $2$ in the $x,y$ plane.
Orient each source edge in the positive $x$ or positive $y$ direction.
Use local coordinates $(\ell,\zeta)$, where $\ell$ is distance from its
left endpoint along that direction, and $\zeta=z$ for a horizontal edge
but $\zeta=-z$ for a vertical edge. The two ancillas are placed at
\begin{equation}\label{eq:gadget-coordinates}
\begin{array}{c|cc}
& A_e & B_e\\ \hline
J_e<0 &(1,0)&(1,1)\\[2pt]
J_e>0 &(1/2,\sqrt3/2)&(3/2,\sqrt3/2).
\end{array}
\end{equation}
The local endpoints are always $L=(0,0)$ and $R=(2,0)$.
Let $p_1,\ldots,p_m$ denote all original and ancillary sites.
Their coordinates belong to $\mathbb Q(\sqrt3)^3$ and have polynomial
bit descriptions.
\Cref{fig:gadgets} shows these placements and the additional contact
between consecutive positive gadgets.
\begin{figure}[tbp]
\centering
\begin{tikzpicture}[
x=1cm,y=1cm,line cap=round,
original/.style={circle,fill=black,draw=black,inner sep=1.7pt},
ancilla/.style={circle,fill=white,draw=black,inner sep=1.9pt},
weak/.style={line width=0.65pt},
strong/.style={line width=1.5pt},
extra/.style={densely dashed,line width=0.75pt},
every node/.style={font=\small}]
\begin{scope}
\node at (1,2.15) {(a) $J_e<0$};
\draw[weak] (0,0)--(1,0)--(2,0);
\draw[strong] (1,0)--(1,1);
\node[original] at (0,0) {};
\node[original] at (2,0) {};
\node[ancilla] at (1,0) {};
\node[ancilla] at (1,1) {};
\node[anchor=north] at (0,-0.08)
{\shortstack{$L$\\[-1pt]\scriptsize$(0,0)$}};
\node[anchor=north] at (1,-0.08)
{\shortstack{$A_e$\\[-1pt]\scriptsize$(1,0)$}};
\node[anchor=north] at (2,-0.08)
{\shortstack{$R$\\[-1pt]\scriptsize$(2,0)$}};
\node[anchor=south] at (1,1.08)
{\shortstack{$B_e$\\[-1pt]\scriptsize$(1,1)$}};
\node[anchor=west] at (1.10,0.62) {$\Delta$};
\node[font=\scriptsize,anchor=south] at (0.45,0.10) {$w_e$};
\node[font=\scriptsize,anchor=south] at (1.65,0.10) {$w_e$};
\end{scope}
\begin{scope}[shift={(4,0)}]
\node at (1,2.15) {(b) $J_e>0$};
\draw[weak] (0,0)--(0.5,0.8660254);
\draw[weak] (1.5,0.8660254)--(2,0);
\draw[strong] (0.5,0.8660254)--(1.5,0.8660254);
\node[original] at (0,0) {};
\node[original] at (2,0) {};
\node[ancilla] at (0.5,0.8660254) {};
\node[ancilla] at (1.5,0.8660254) {};
\node[anchor=north] at (0,-0.08)
{\shortstack{$L$\\[-1pt]\scriptsize$(0,0)$}};
\node[anchor=north] at (2,-0.08)
{\shortstack{$R$\\[-1pt]\scriptsize$(2,0)$}};
\node[anchor=south east] at (0.52,0.97)
{\shortstack{$A_e$\\[-1pt]\scriptsize$(1/2,h)$}};
\node[anchor=south west] at (1.48,0.97)
{\shortstack{$B_e$\\[-1pt]\scriptsize$(3/2,h)$}};
\node[anchor=north] at (1,0.79) {$\Delta$};
\node[font=\scriptsize,anchor=east] at (0.18,0.45) {$w_e$};
\node[font=\scriptsize,anchor=west] at (1.82,0.45) {$w_e$};
\end{scope}
\begin{scope}[shift={(8,0)}]
\node at (2,2.15) {(c) Successive positive bonds};
\draw[weak] (0,0)--(0.5,0.8660254);
\draw[weak] (1.5,0.8660254)--(2,0)--(2.5,0.8660254);
\draw[weak] (3.5,0.8660254)--(4,0);
\draw[strong] (0.5,0.8660254)--(1.5,0.8660254);
\draw[strong] (2.5,0.8660254)--(3.5,0.8660254);
\draw[extra] (1.5,0.8660254)--(2.5,0.8660254);
\foreach \x in {0,2,4} {\node[original] at (\x,0) {};}
\foreach \x in {0.5,1.5,2.5,3.5}
{\node[ancilla] at (\x,0.8660254) {};}
\node[anchor=north] at (0,-0.08) {$L$};
\node[anchor=north] at (2,-0.08) {$R=L'$};
\node[anchor=north] at (4,-0.08) {$R'$};
\node[anchor=south] at (0.5,0.97) {$A_e$};
\node[anchor=south] at (1.5,0.97) {$B_e$};
\node[anchor=south] at (2.5,0.97) {$A_f$};
\node[anchor=south] at (3.5,0.97) {$B_f$};
\draw[thin] (-0.2,0)--(-0.2,0.8660254);
\draw[thin] (-0.26,0)--(-0.14,0);
\draw[thin] (-0.26,0.8660254)--(-0.14,0.8660254);
\node[anchor=east,font=\scriptsize] at (-0.24,0.433) {$h$};
\draw[thin] (2,1.30)--(2,0.91);
\node[font=\scriptsize,anchor=south] at (2,1.32)
{extra contact};
\end{scope}
\end{tikzpicture}
\caption{The contact construction in the local $(\ell,\zeta)$ plane,
with $h=\sqrt3/2$ and $w_e=\sqrt\Delta\,v_e$.
Filled vertices are original spins and open vertices are ancillas.
Every solid segment has length one; thick segments carry weight
$\Delta$ and thin segments the indicated weak weight. The dashed
segment in (c) is an additional unit contact, with weight zero in
$H_+$; it receives only the small positive hopping floor in the Hubbard
model. Horizontal gadgets use $\zeta=z$, and vertical gadgets use
$\zeta=-z$.}
\label{fig:gadgets}
\end{figure}
\begin{lemma}[Contact classification]\label{lem:contact-layout}
All the sites $p_i$ are distinct. Every pair is either at distance one
or at distance at least $\sqrt2$. The complete set $\mathcal C$ of
unit-distance pairs consists of:
\begin{enumerate}
\item $L A_e$, $A_e R$, and $A_e B_e$ for a negative source edge;
\item $L A_e$, $A_e B_e$, and $B_e R$ for a positive source edge;
\item $B_e A_f$ for two consecutive positive source edges $e,f$ on
the same original lattice line.
\end{enumerate}
Thus every interaction of $H_+$ is a contact. Its weight on a contact
will be denoted $K_{ij}\geq0$, using $K_{ij}=0$ for contacts in the
third class.
\end{lemma}
\begin{proof}
Within one gadget the assertion follows directly from
\eqref{eq:gadget-coordinates}. For a negative gadget the two distances
from $B_e$ to original endpoints are $\sqrt2$; for a positive gadget
the two nonadjacent endpoint--ancilla distances are $\sqrt3$.
Original vertices are mutually separated by at least $2$.
An original vertex other than the endpoints of an ancilla's gadget,
on that same axis line, has along-line distance at least $3/2$ to a
positive ancilla and at least $3$ to a negative midpoint.
An original vertex off that line differs in a transverse coordinate
by at least $2$.
For ancillas belonging to distinct gadgets on the same axis line,
the only along-line separation below $3/2$ is the separation $1$
between $B_e$ and $A_f$ for consecutive positive gadgets. They have
the same height and hence are a contact. All remaining such pairs
are at least $3/2>\sqrt2$ apart. Distinct parallel axis lines have
transverse separation at least $2$.
It remains to compare horizontal and vertical ancillas. Call an
ancilla \emph{midpoint type} if it belongs to a negative gadget,
and \emph{positive type} otherwise. Two midpoint-type projections
differ by at least $1$ in each planar coordinate, so their distance
is at least $\sqrt2$. A midpoint-type and a positive-type projection
differ by at least $1$ and $1/2$ in the two planar coordinates.
Their heights differ by at least $\sqrt3/2$, since horizontal and
vertical heights have opposite signs. Their squared distance is
therefore at least $1+1/4+3/4=2$. Finally two positive-type ancillas
from perpendicular directions have height difference $\sqrt3$,
already larger than $\sqrt2$. These cases exhaust all pairs and
also establish distinctness.
\end{proof}
Orient each contact arbitrarily and write
$u_{ij}=p_j-p_i$, so $|u_{ij}|=1$. Its first-order length change under
site displacements $s_i$ is $u_{ij}\cdot(s_j-s_i)$. Reversing the
orientation reverses both factors and leaves this scalar unchanged.
The following stronger-than-counting fact is what permits independent
hopping calibration.
\begin{lemma}[A bounded right inverse for contact lengths]
\label{lem:displacement}
For every real assignment $(d_{ij})_{\{i,j\}\in\mathcal C}$ there
are displacements $s_i\in\mathbb R^3$ such that
\begin{equation}\label{eq:contact-displacement}
u_{ij}\cdot(s_j-s_i)=d_{ij},\qquad
\max_i|s_i|\leq6m\max_{\mathcal C}|d_{ij}|.
\end{equation}
For an empty contact set take all displacements zero. The construction
is linear in $d$, uses polynomially many arithmetic operations with
coefficients in $\mathbb Q(\sqrt3)$, and has no consistency condition
around lattice cycles.
\end{lemma}
\begin{proof}
Set the height displacement of every original vertex to zero.
The horizontal coordinate of each original vertex will be determined
using only its original row, and its vertical planar coordinate using
only its original column. These are independent variables even when
both directions meet at the same vertex.
Fix one such axis line. Write $a$ for displacement along it, $z$
for physical height displacement, and put $\sigma=1$ for a horizontal
line and $\sigma=-1$ for a vertical line. On every connected run of
negative bonds choose the first original $a$ coordinate to be zero,
and accumulate
\begin{equation}\label{eq:negative-displacement}
\begin{split}
a_R&=a_L+d_{LA_e}+d_{A_eR},\qquad
a_{A_e}=a_L+d_{LA_e},\\
a_{B_e}&=a_{A_e},\qquad
z_{A_e}=0,\qquad z_{B_e}=\sigma d_{A_eB_e}.
\end{split}
\end{equation}
Unconstrained original coordinates may be set to zero. A positive
or missing bond places no requirement on the difference of its
original endpoints' $a$ coordinates at this step, so it separates
runs. Each run is a path on a line, and
\eqref{eq:negative-displacement} satisfies all three constraints in
its negative gadgets.
Now consider the elevated vertices of positive gadgets on this line.
Their same-height contacts, including the extra contacts of
Lemma~\ref{lem:contact-layout}, form disjoint paths. Set the first
$a$ coordinate on each path to zero and successively set
$a_{\rm next}=a_{\rm current}+d_{\rm edge}$.
Every such contact vector is along the axis, so its constraint
involves only these $a$ coordinates. Each elevated vertex has
exactly one additional, sloped contact to an original vertex.
Writing $h_0=\sqrt3/2$, set
\begin{equation}\label{eq:positive-displacement}
\begin{split}
z_{A_e}&=\frac{d_{LA_e}-(a_{A_e}-a_L)/2}{\sigma h_0},\\
z_{B_e}&=\frac{d_{B_eR}-(a_R-a_{B_e})/2}{\sigma h_0}.
\end{split}
\end{equation}
Indeed the two unit spoke vectors are
$\tfrac12 e_{\rm axis}+\sigma h_0 e_z$ and
$\tfrac12 e_{\rm axis}-\sigma h_0 e_z$, respectively.
Because the original height displacements vanish, these formulas
satisfy the remaining equations without imposing any new condition
on original coordinates. Set each ancilla's unused transverse
coordinate to zero.
Repeat on all rows and columns. The row procedure changes only
original $x$ coordinates, and the column procedure only original
$y$ coordinates. No ancilla belongs to both procedures. Thus the
constructions do not overwrite each other. In particular a square
cycle is not a cycle in either scalar accumulation graph; its vector
closure is automatic because actual vertex displacements have been
assigned.
For the norm bound put $d_{\max}=\max_{\mathcal C}|d_{ij}|$.
Every original along-line coordinate has magnitude at most
$2m d_{\max}$, and every accumulated elevated coordinate has
magnitude at most $m d_{\max}$. Negative ancilla coordinates
have magnitude at most $(2m+1)d_{\max}$ along the line and at
most $d_{\max}$ in height. From
\eqref{eq:positive-displacement}, positive ancilla heights are at
most $(1+3m/2)d_{\max}/h_0$. Combining at most two nonzero
coordinates at any site gives the stated bound $6m d_{\max}$.
All assignments are additions and multiplication by fixed rational
numbers or by $1/h_0$, proving the algorithmic assertion.
\end{proof}
Only approximate length equations will be needed. To make their bit
complexity explicit, suppose the data are known within $\tau/4$ and
apply the displayed right inverse to those approximations. Approximate
its output and $Dp_i$, then round to rational positions $X_i$, with
total Euclidean error at most $\tau/4$ per site. Defining
$s_i=X_i-Dp_i$ absorbs this rounding even though the $p_i$ may be
irrational. Every contact equation then has error at most
$\tau/4+2(\tau/4)<\tau$, and
\begin{equation}\label{eq:rounded-displacement}
\max_i|s_i|\leq6m\bigl(\max_{\mathcal C}|d_{ij}|+\tau\bigr)+\tau.
\end{equation}
The number of arithmetic operations is polynomial, and all accumulated
coefficients have polynomial size. Computing $\sqrt3$ and the other
real inputs to sufficiently many additional polynomially bounded bits
therefore gives the stated absolute accuracy in polynomial time.
In particular this does not call for an ill-conditioned general
linear-system solve.
For later use, actual distances satisfy the exact vector identities
and the elementary estimates
\begin{equation}\label{eq:displaced-geometric-bounds}
\begin{aligned}
X_j-X_i&=D u_{ij}+s_j-s_i
&&(\{i,j\}\in\mathcal C),\\
|X_j-X_i|&=D+u_{ij}\cdot(s_j-s_i)+O(S^2/D)
&&(\{i,j\}\in\mathcal C),\\
|X_j-X_i|&\geq\sqrt2D-2S
&&(\{i,j\}\notin\mathcal C),
\end{aligned}
\end{equation}
provided $\max_i|s_i|\leq S$ and $S\leq D/4$. The middle estimate
is the Taylor bound for the Euclidean norm away from zero. The last
is the triangle inequality and Lemma~\ref{lem:contact-layout}.
No exact nonlinear prescribed-length problem is being solved; the
quadratic remainder will be included in the hopping estimates.
\subsection{The half-filled Hubbard model}
We now realize those positive bonds as second-order hopping processes.
Virtual double occupancy gives the antiferromagnetic superexchange
mechanism of Anderson~\cite{Anderson1959}. In our Pauli-matrix
normalization the coefficient is $t^2/U$ multiplying
$h^s_{ij}-I$. We calculate the fermionic signs below and retain the
constant term, because it contributes to the final decision thresholds.
Associate two fermionic modes, labelled $\uparrow,\downarrow$, to each
of the $m$ sites and restrict to exactly $m$ fermions. Let
$c_{i\sigma}^\dagger,c_{i\sigma}$ satisfy the canonical
anticommutation relations and put
$n_{i\sigma}=c_{i\sigma}^\dagger c_{i\sigma}$.
For nonnegative hopping strengths define
\begin{equation}\label{eq:hubbard-definition}
\begin{split}
H_{\rm Hub}(t)
&=-\sum_{\{i,j\}\in\mathcal C}\sum_{\sigma=\uparrow,\downarrow}
t_{ij}\bigl(c_{i\sigma}^\dagger c_{j\sigma}
+c_{j\sigma}^\dagger c_{i\sigma}\bigr)
+U\sum_i n_{i\uparrow}n_{i\downarrow},\\
U&=\iint_{\mathbb R^3\times\mathbb R^3}
\frac{\phi(y)^2\phi(y')^2}{|y-y'|}\,\dd y\,\dd y',
\qquad \phi(y)=\pi^{-1/2}e^{-|y|}.
\end{split}
\end{equation}
The radial probability density of $\phi^2$ is $p(r)=4r^2e^{-2r}$.
Angular integration gives the kernel $1/\max(r,s)$: for $r,s>0$,
\[
\frac12\int_{-1}^{1}\frac{\dd u}{\sqrt{r^2+s^2-2rsu}}
=\frac1{\max(r,s)}.
\]
Since $\int_0^r p(s)\,\dd s=1-e^{-2r}(1+2r+2r^2)$,
\begin{equation}\label{eq:onsite-integral}
\begin{split}
U&=8\int_0^\infty r e^{-2r}
\bigl[1-e^{-2r}(1+2r+2r^2)\bigr]\,\dd r\\
&=8\left(\frac14-\frac1{16}-\frac1{16}-\frac3{64}\right)
=\frac58.
\end{split}
\end{equation}
Thus $U$ is a fixed positive rational constant.
Let $q_{\mathcal C}=|\mathcal C|$ and define the rational scales
\begin{equation}\label{eq:hubbard-scales}
\overline R=\max\{1,q_{\mathcal C}\Delta\},\qquad
\eta=\frac{U\gamma}{1024\overline R^3},\qquad
\alpha=\frac{\eta^2}{U},\qquad
f=\frac{\alpha\gamma}{128\max\{1,q_{\mathcal C}\}}.
\end{equation}
The ideal hoppings are $t^{(0)}_{ij}=\eta\sqrt{K_{ij}}$.
For each contact compute a rational $r_{ij}$ with
$|r_{ij}-t^{(0)}_{ij}|\leq f$, and output
\begin{equation}\label{eq:hopping-floor}
t_{ij}=\max\{f,r_{ij}\}.
\end{equation}
Polynomial-precision bisection on nonnegative square roots (at most
two nested square roots here) suffices to compute $r_{ij}$; no exact
comparison between irrational numbers is required. The maximum in
\eqref{eq:hopping-floor} is a comparison of rationals. All final
hoppings satisfy $1/\poly(n)\leq t_{ij}<1$, including those with
$K_{ij}=0$.
\begin{proposition}[Hubbard reduction with rational hoppings]
\label{prop:hubbard}
The construction above is deterministic and polynomial time. With
\begin{equation}\label{eq:finite-scalar-shift}
C_s=C_0+\sum_{\{i,j\}\in\mathcal C}K_{ij}
=3k\Delta+\frac32\sum_e v_e^2
+\sum_{\{i,j\}\in\mathcal C}K_{ij},
\end{equation}
it satisfies
\begin{equation}\label{eq:finite-combined-error}
\left|E\bigl(H_{\rm Hub}(t)\bigr)
-\alpha\bigl(E(H_{\rm in})-C_s\bigr)\right|
\leq\frac{3\alpha\gamma}{64}.
\end{equation}
The scalar $C_s$ has polynomial magnitude and can be approximated to
any polynomial-bit absolute precision in polynomial time.
\end{proposition}
\begin{proof}
First use the ideal hoppings $t^{(0)}$. Let
$A_U=U\sum_i n_{i\uparrow}n_{i\downarrow}$ and let $T$ be the hopping
part of $H_{\rm Hub}(t^{(0)})$. At filling $m$, the kernel of $A_U$
is exactly the singly occupied subspace: having $m$ fermions, at most
one at each of $m$ sites, forces precisely one per site. Identify this
subspace with $m$ spins, and denote its projection by $P_{\rm so}$.
The complementary excitation energies are positive integer multiples
of $U$. A single hop from $P_{\rm so}$ creates one vacant site
and one doubly occupied site and hence has energy exactly $U$.
In particular $P_{\rm so}TP_{\rm so}=0$.
A second hop can return to single occupancy only by moving an electron
from that doubly occupied site to that vacant site. Thus the two hops
must use the same unordered bond; products from different bonds have
zero compression, even when the bonds share a vertex.
For completeness, consider one bond $i,j$ of strength $t$, with local
mode order $i\uparrow,i\downarrow,j\uparrow,j\downarrow$.
Write $|D_i\rangle=c_{i\uparrow}^\dagger
c_{i\downarrow}^\dagger|0\rangle$ and similarly for $|D_j\rangle$.
The bond hopping operator $T_{ij}$ obeys
\[
T_{ij}|\uparrow,\downarrow\rangle
=-t(|D_i\rangle+|D_j\rangle),\qquad
T_{ij}|\downarrow,\uparrow\rangle
= t(|D_i\rangle+|D_j\rangle).
\]
It consequently annihilates all triplets and maps the singlet to
$-\sqrt2t(|D_i\rangle+|D_j\rangle)$, of squared norm $4t^2$.
In an ordering containing other sites between $i$ and $j$, those sites
are singly occupied during this two-hop process; their fixed fermion
parity changes transition phases but not the compressed product.
Equivalently, permuting the intervening mode pairs gives the same
calculation under a fermionic change of basis. It follows that
\begin{equation}\label{eq:hubbard-exchange}
-P_{\rm so}T_{ij}A_U^{-1}T_{ij}P_{\rm so}
=-\frac{4t^2}{U}P_{{\rm singlet},ij}
=\frac{t^2}{U}(h^s_{ij}-I).
\end{equation}
Here $P_{{\rm singlet},ij}=(I-h^s_{ij})/4$.
For each spin, $c_{i\sigma}^\dagger c_{j\sigma}
+c_{j\sigma}^\dagger c_{i\sigma}$ has norm one, so
\[
\lVert T\rVert\leq2\sum_{\mathcal C}t^{(0)}_{ij}
=2\eta\sum_{\mathcal C}\sqrt{K_{ij}}
\leq2\eta\overline R.
\]
The last inequality uses $K_{ij}\leq\Delta$ and $\Delta\geq1$.
Lemma~\ref{lem:second-order}, with $g=U$ and
$\epsilon=2\eta\overline R/U<1/2$, and
\eqref{eq:hubbard-exchange} give
\begin{equation}\label{eq:ideal-hubbard-error}
\left|E\bigl(H_{\rm Hub}(t^{(0)})\bigr)
-\alpha\left(E(H_+)-\sum_{\mathcal C}K_{ij}\right)\right|
\leq\frac{16\eta^3\overline R^3}{U^2}
\leq\frac{\alpha\gamma}{64}.
\end{equation}
If the complementary subspace is empty, as can happen for $m=1$,
there is no hopping and this assertion is exact.
For the rational hoppings, nonnegativity of $t^{(0)}_{ij}$ and the
accuracy of $r_{ij}$ imply
$|\max\{f,r_{ij}\}-t^{(0)}_{ij}|\leq f$.
Therefore
\begin{equation}\label{eq:hopping-rounding-error}
\left\|H_{\rm Hub}(t)-H_{\rm Hub}(t^{(0)})\right\|
\leq2q_{\mathcal C}f\leq\frac{\alpha\gamma}{64}.
\end{equation}
Also $t^{(0)}_{ij}\leq\eta\overline R\leq U/1024$ and
$f\leq U/(128\cdot1024^2)$, which verify $f\leq t_{ij}<1$.
The floor $f$ is an inverse polynomial. This argument covers an ideal
hopping much smaller than $f$ and an ideal hopping equal to zero; no
coefficient lower bound has been used. If $\mathcal C$ is empty,
\eqref{eq:hopping-rounding-error} is simply zero.
Finally, Proposition~\ref{prop:positive-spin} with
\eqref{eq:spin-scale} contributes at most another
$\alpha\gamma/64$. Adding the three errors proves
\eqref{eq:finite-combined-error}. All sizes, magnitudes, and reciprocal
precision budgets in the construction are polynomial in $n$.
Each $K_{ij}$ is either $\Delta$, a square root of the rational number
$2\Delta|J_e|$, or zero. The displayed scalar shift is thus efficiently
approximable to the asserted precision, completing the proof.
\end{proof}
==== END Continuum-Coulomb-hardness-with-binary-nuclear-charges//01-finite-models.tex ====
==== BEGIN Continuum-Coulomb-hardness-with-binary-nuclear-charges//02-wells.tex ====
\section{Isolated wells and deterministic energy tuning}
\label{sec:wells}
This section supplies an exponentially accurate classical procedure for
tuning the individual well energies. The procedure operates on
one-electron operators only; it does not diagonalize any of the spin or
Hubbard Hamiltonians from \Cref{sec:finite}.
There are two separate tasks. First we must show that moving a secondary
positive nucleus changes one isolated energy monotonically, with much
weaker effects on the other sites. Then we must evaluate those energies
accurately enough to implement the resulting contraction. The required
accuracy is exponential in $D$: an inverse-polynomial energy error in
these dilated units would grow after the final amplification. The
analytic and arithmetic estimates below address that precision cost.
Fix centers \(X_1,\ldots,X_m\in\R^3\) satisfying
\begin{equation}
|X_i-X_j|\ge D-2S\quad(i\ne j),\qquad S\le D/100.
\label{eq:wells-separation}
\end{equation}
Here \(D\) is a sufficiently large integer and \(m/D\) is smaller than a
sufficiently small absolute constant. All constants in this section are
uniform in the centers and the parameters \(\rho_i\in[1,4]\).
The stronger scale requirements used later imply these hypotheses.
Computational statements additionally assume rational centers. Their
bit lengths, rather than their numerical magnitudes, enter the
complexity bounds.
Choose an integer \(Z_*\) with \(e^D/2\le Z_*\le2e^D\). Round the following
two numbers to nearest integers, once and for all, before tuning:
\begin{equation}
\begin{split}
Z_*q&=\operatorname{round}\!\left(
Z_*\left(1-\frac{500m}{D}\right)\right),\\
Z_*q'&=\operatorname{round}\!\left(Z_*\frac{1000m}{D}\right).
\end{split}
\label{eq:wells-charges}
\end{equation}
Both integers are positive for the scales considered here. In particular,
\begin{equation}
\begin{split}
q&=1-\frac{500m}{D}+O(Z_*^{-1}),\\
q'&=\frac{1000m}{D}+O(Z_*^{-1}),\qquad
c\frac mD\le q'\le C\frac mD.
\end{split}
\label{eq:wells-charge-bounds}
\end{equation}
The primary and secondary nuclei of site \(i\) are at \(X_i\) and
\(X_i+\rho_i e_z\), respectively, in scaled coordinates. Write
\begin{equation}
\begin{split}
r_i&=|y-X_i|,\qquad
V_i(y)=\frac q{r_i}+\frac{q'}{|y-X_i-\rho_i e_z|},\\
h&=-\frac12\Delta-\sum_{i=1}^m V_i.
\end{split}
\label{eq:wells-full-one-body}
\end{equation}
These are the one-electron operators obtained from physical coordinates
by \(y=Z_*x\) and division of energies by \(Z_*^2\).
Fix a nondecreasing \(C^2\) function \(\chi:[0,\infty)\to[0,1]\) that
vanishes on \([0,1/16]\) and equals one on \([1/8,\infty)\). It can have
rational piecewise-polynomial data: integrate a normalized positive
multiple of \((t-1/16)^2(1/8-t)^2\) on the transition interval and extend
constantly outside. Set
\begin{equation}
\chi_i(y)=\chi(r_i/D),\qquad
W_i=\sum_{k\ne i}\chi_k V_k,\qquad
h_i=-\frac12\Delta-V_i-W_i.
\label{eq:wells-capped}
\end{equation}
The caps remove the other nuclear singularities from \(h_i\) while
retaining their fields near site \(i\).
\subsection{Uniform ground-state estimates}
We use the three-dimensional Sobolev bounds
\(\|u\|_6\le C\|u\|_{H^1}\) and
\(\|u\|_\infty\le C\|u\|_{H^2}\), together with their localized versions.
We also use the translated Hardy inequality
\begin{equation}
\left\|\frac{u}{|\,\cdot-a\,|}\right\|_2
\le2\|\nabla u\|_2\qquad(a\in\R^3).
\label{eq:wells-hardy}
\end{equation}
Indeed, integrate
\(\operatorname{div}((y-a)/|y-a|^2)=|y-a|^{-2}\)
against \(|u|^2\) for a smooth compactly supported \(u\). Integration by
parts and Cauchy--Schwarz give
\[
\|u/|\,\cdot-a\,|\|_2^2
\le2\|u/|\,\cdot-a\,|\|_2\|\nabla u\|_2.
\]
Approximation proves \eqref{eq:wells-hardy} on \(H^1\).
Thus Coulomb multiplication maps \(H^1\) continuously into \(L^2\),
uniformly in the location of the pole.
\begin{proposition}[Uniform isolated-well estimates]
\label{prop:isolated-well}
The operators \(h_i\) are self-adjoint on \(H^2(\R^3)\), with form domain
\(H^1(\R^3)\). Each has a simple ground eigenvalue \(\lambda_i\) and a
real, positive, normalized eigenfunction \(\psi_i\). There is a constant
\(g>0\) such that the remaining spectrum is at least \(\lambda_i+g\).
For \(\phi_i(y)=\pi^{-1/2}e^{-r_i}\),
\begin{equation}
|\lambda_i+1/2|\le C\frac mD,\qquad
\|\psi_i-\phi_i\|_{H^2}\le C\frac mD.
\label{eq:wells-h2}
\end{equation}
Moreover, for derivative orders \(j=0,1,2\),
\begin{equation}
\|D^jW_i\|_\infty\le C\frac m{D^{j+1}},
\qquad 0\le W_i\le\frac{100m}{D}.
\label{eq:wells-field-bounds}
\end{equation}
These assertions hold before tuning and uniformly throughout
\([1,4]^m\).
\end{proposition}
\begin{proof}
On the nonzero support of \(\chi_k\), \(r_k\ge D/16\). Since \(\rho_k\le4\),
\[
|y-X_k-\rho_k e_z|\ge r_k-4\ge r_k/2
\]
for large \(D\). The secondary singularity is inside the region where
the cap is identically zero. Differentiation through order two,
including derivatives of the cap, gives
\(\|D^j(\chi_kV_k)\|_\infty\le CD^{-j-1}\).
Summing proves the derivative bounds. For example, \(q\le1\), \(q'\le1\),
and sufficiently large \(D\) give a bound below \(100/D\) for each capped
potential, proving the stated function bound.
Put \(H_i^0=-\Delta/2-r_i^{-1}\) and \(B_i=h_i-H_i^0\).
Hardy's inequality and \eqref{eq:wells-field-bounds} give
\begin{equation}
\|B_i u\|_2\le C\frac mD\|u\|_{H^1}.
\label{eq:wells-perturbation}
\end{equation}
For any potential in these single-well operators, Hardy and Fourier
interpolation give, for each \(a>0\), a uniform bound
\(\|Vu\|_2\le a\|\Delta u\|_2+C_a\|u\|_2\).
The Kato--Rellich theorem for relatively bounded perturbations
\cite[Theorem~6.4, p.~135]{Teschl2009} therefore gives domain \(H^2\). The same estimates, with the Laplacian term
absorbed, prove uniform graph-norm equivalence:
\begin{equation}
C^{-1}\|u\|_{H^2}
\le \|h_i u\|_2+C\|u\|_2
\le C\|u\|_{H^2}.
\label{eq:wells-graph-norm}
\end{equation}
The corresponding statement holds for \(H_i^0\).
At form level a fixed positive shift of either operator is bounded
above and below by constant multiples of \(\|u\|_{H^1}^2\), giving
the asserted form domains.
The hydrogenic ground state and its isolation follow without an
excited-state formula. Direct integration by parts gives
\begin{equation}
\langle u,(H_i^0+1/2)u\rangle
=\frac12\left\|\nabla u+\frac{y-X_i}{r_i}u\right\|_2^2.
\label{eq:wells-hydrogen-square}
\end{equation}
Its zero vectors are exactly the multiples of \(\phi_i\).
If normalized vectors orthogonal to \(\phi_i\) had energies tending to
\(-1/2\), they would be bounded in \(H^1\). The Coulomb form is compact
under weak convergence of such a sequence: within radius \(a\) of the
pole its absolute value is bounded by \(Ca\|u\|_{H^1}^2\), by Sobolev and
H\"older; beyond radius \(R\) it is bounded by \(R^{-1}\|u\|_2^2\);
on the intervening annulus local compactness applies. A weak limit
would therefore obey
\[
\langle u,H_i^0u\rangle\le-1/2,\qquad
\langle u,H_i^0u\rangle\ge-\tfrac12\|u\|_2^2\ge-1/2.
\]
Equality forces norm one and a multiple of \(\phi_i\), contrary to
orthogonality. This proves a fixed gap \(g_0>0\).
By \eqref{eq:wells-perturbation}, after a fixed positive shift the
perturbed form lies between \(1-Cm/D\) and \(1+Cm/D\) times the
unperturbed shifted form. The variational principle puts its first
level within \(Cm/D\) of \(-1/2\), and its second variational level at
least \(-1/2+g_0-Cm/D\). The negative first level is attained.
The preceding truncation argument gives form compactness of both own
Coulomb potentials and of the bounded decaying \(W_i\). A normalized
minimizing sequence is bounded in \(H^1\). Passing to a weakly convergent
subsequence with limit \(u\), compactness and lower semicontinuity give
\[
\lambda_i\|u\|_2^2\le\langle u,h_i u\rangle\le\lambda_i<0.
\]
Since \(\|u\|_2\le1\), these inequalities force \(\|u\|_2=1\)
and attainment. For small enough \(m/D\) the
first level is simple and separated from the rest by a fixed \(g>0\).
Taking the absolute value of a real minimizer cannot increase its
energy. The Euler equation and the strong maximum principle on the
connected complement of the poles make it positive there.
The form bounds give \(\|\psi_i\|_{H^1}\le C\). Its equation implies
\[
(H_i^0+1/2)\psi_i
=(\lambda_i+1/2)\psi_i-B_i\psi_i,
\]
with right side of norm at most \(Cm/D\). On \(\phi_i^\perp\), the inverse
of \(H_i^0+1/2\) is bounded from \(L^2\) to \(H^2\): use the gap for the
\(L^2\) bound, then \eqref{eq:wells-graph-norm}.
The orthogonal component of \(\psi_i\) consequently has \(H^2\) norm
\(O(m/D)\). Normalization and the positive overlap with \(\phi_i\) put
its coefficient along \(\phi_i\) within \(O((m/D)^2)\) of one.
This proves \eqref{eq:wells-h2} and the uniform sup bound. Positivity
also holds at the poles by continuity and this sup-norm approximation,
since both poles lie within distance four of \(X_i\).
\end{proof}
\subsection{Moving the secondary nucleus}
The derivative of a moving Coulomb potential is not an \(L^2\)-bounded
multiplier. We instead control its expectations by Hardy's inequality.
\begin{lemma}[Bracketing and uniform secant bounds]
\label{lem:wells-secants}
With centers and rounded charges fixed, for every choice of the other
parameters,
\(\lambda_i(\rho_i=1)<-1/2<\lambda_i(\rho_i=4)\).
For \(1\le a0.7,\quad f(4)\le\tfrac14,\quad
\frac{1-5e^{-2}}{16}\le-f'(\rho)\le1\quad(1\le\rho\le4).
\]
At the first endpoint the leading shift in
\eqref{eq:wells-energy-expansion} is at most \((500-700)m/D\);
at the second it is at least \((500-250-100)m/D\).
The quadratic and rounding errors fit inside these strict margins,
proving bracketing.
For the secant bound, fix a ground function \(u\) at either endpoint,
and put \(A(t)=X_i+t e_z\). The kernel derivative is dominated by
\(|y-A(t)|^{-2}\), and Hardy gives
\begin{align}
\int\frac{|u^2-\phi_i^2|}{|y-A(t)|^2}\dd y
&\le\left\|\frac{u-\phi_i}{|\,\cdot-A(t)\,|}\right\|_2
\left\|\frac{u+\phi_i}{|\,\cdot-A(t)\,|}\right\|_2 \notag\\
&\le4\|\nabla(u-\phi_i)\|_2
\bigl(\|\nabla u\|_2+\|\nabla\phi_i\|_2\bigr)
\le C\frac mD.
\label{eq:wells-density-hardy}
\end{align}
The same bound with \(u^2\) justifies absolute integration in \(y\)
and \(t\). Thus replacing the density by \(\phi_i^2\) when integrating
the moving-pole kernel derivative incurs at most
\(Cq'(m/D)(b-a)\) error.
The variational principle gives
\[
\langle\psi_i(b),(h_i(b)-h_i(a))\psi_i(b)\rangle
\le\lambda_i(b)-\lambda_i(a)
\le\langle\psi_i(a),(h_i(b)-h_i(a))\psi_i(a)\rangle.
\]
Both bounding quantities differ from \(-q'(f(b)-f(a))\) by the error
just estimated, proving \eqref{eq:wells-own-slope}.
For \(k\ne i\), the moving secondary pole in \(\chi_kV_k\) stays at
least \(D/32\) from the region where its cap is nonzero.
Its bounded potential changes in sup norm by at most
\(Cq'D^{-2}|\rho_k-\widetilde\rho_k|\).
Summing and applying the variational principle for a bounded additive
perturbation proves \eqref{eq:wells-other-slope}.
\end{proof}
\Needspace{10\baselineskip}
\subsection{Uniform tails}
\begin{lemma}[Crude tails before tuning]
\label{lem:well-tail}
Uniformly for all \(\rho\in[1,4]^m\),
\begin{equation}
|\psi_i(y)|\le C\exp[-.99(r_i-.01D)],
\label{eq:wells-tail}
\end{equation}
and, for sufficiently large \(D\),
\begin{equation}
\|\psi_i\|_{H^1(\{r_i>.65D\})}\le e^{-.61D}.
\label{eq:wells-h1-tail}
\end{equation}
\end{lemma}
\begin{proof}
On \(r_i\ge .01D\), both own poles are far away and
\(V_i+W_i=O(m/D)\). For \(a=.99\) and
\(b(y)=C\exp[-a(r_i-.01D)]\), direct calculation gives
\[
\frac{(h_i-\lambda_i)b}{b}
=-\frac{a^2}{2}+\frac a{r_i}-V_i-W_i-\lambda_i.
\]
Its constant part is \((1-a^2)/2+O(m/D)>0\) at the chosen scales.
The zero-order coefficient of \(h_i-\lambda_i\) is also positive
throughout this exterior region. The uniform sup bound fixes \(C\)
so that \(b\ge|\psi_i|\) at the inner boundary. An \(H^2(\R^3)\) function
is uniformly continuous and vanishes at infinity; otherwise uniform
continuity on a sequence of disjoint balls would contradict its
square integrability. On a large finite annulus add a constant
\(\varepsilon>0\) to \(b\) to dominate both signs of \(\psi_i\) at the
outer boundary. The supersolution inequality is preserved.
Testing the negative part of \(b+\varepsilon\mp\psi_i\) proves
comparison. Send the outer radius to infinity and then
\(\varepsilon\) to zero. Inside \(.01D\), the same sup bound proves
\eqref{eq:wells-tail}.
For the gradient estimate choose a smooth cutoff \(\zeta\) that vanishes
on \(r_i\le .645D\), equals one on \(r_i\ge .65D\), and has gradient
\(O(D^{-1})\). Testing the eigen-equation with \(\zeta^2\psi_i\) gives
\[
\frac12\|\nabla(\zeta\psi_i)\|_2^2
=\int(V_i+W_i+\lambda_i)|\zeta\psi_i|^2\dd y
+\frac12\int|\nabla\zeta|^2|\psi_i|^2\dd y.
\]
All coefficients on this support are bounded. Thus
\eqref{eq:wells-tail} bounds the exterior \(H^1\) norm by a polynomial in
\(D\) times \(\exp[-.99(.645-.01)D]\).
Since \(.99(.645-.01)=.62865>.61\), the polynomial factor is absorbed
for large \(D\).
\end{proof}
\subsection{An analytic estimate with uniform constants}
We give the local regularity argument used in the numerical procedure.
It explains why the shrinking length scales near a moving Coulomb pole
do not increase the required polynomial degree beyond \(O(D)\).
\begin{lemma}[Uniform analytic patches]
\label{lem:wells-analytic}
Suppose \(u\) is bounded on the unit ball in \(\R^3\) and satisfies
\(\Delta u=Fu\) there, where
\begin{equation}
\|\partial^\beta F\|_\infty\le A^{|\beta|+1}|\beta|!
\quad\text{for every multi-index }\beta.
\label{eq:wells-analytic-coefficient}
\end{equation}
There is a constant \(A_1\), depending only on \(A\), such that
\begin{equation}
\max_{|\alpha|=k}\|\partial^\alpha u\|_{L^\infty(B_{1/4})}
\le \|u\|_{L^2(B_1)}A_1^{k+1}k!.
\label{eq:wells-analytic-derivatives}
\end{equation}
On a sufficiently small fixed ball about the origin, Taylor
polynomials consequently have geometrically decreasing errors in
both value and first derivative, with constants depending only on \(A\).
\end{lemma}
\begin{proof}
For two concentric balls of radii between \(1/2\) and \(1\), separated
by \(h\), the local estimates needed below are
\begin{align}
\|D^2v\|_{2,\mathrm{inner}}
&\le C\bigl(\|\Delta v\|_{2,\mathrm{outer}}
+h^{-2}\|v\|_{2,\mathrm{outer}}\bigr),
\label{eq:wells-local-h2}\\
\|\nabla v\|_{2,\mathrm{inner}}
&\le C\bigl(h\|\Delta v\|_{2,\mathrm{outer}}
+h^{-1}\|v\|_{2,\mathrm{outer}}\bigr).
\label{eq:wells-local-h1}
\end{align}
To obtain them, insert smooth cutoffs with first and second derivatives
\(O(h^{-1})\) and \(O(h^{-2})\). Integration by parts against a squared
cutoff and Cauchy--Schwarz, with \(ab\le (h a)^2/2+(b/h)^2/2\),
bound first derivatives on an intermediate ball as in
\eqref{eq:wells-local-h1}. Apply the whole-space Fourier inequality
\(\|D^2w\|_2\le C\|\Delta w\|_2\) to the cutoff function and substitute
that first-derivative bound in its commutator terms to obtain
\eqref{eq:wells-local-h2}. Approximation proves these estimates for
weak solutions. Successive localized Laplacian regularity makes \(u\)
smooth when \(F\) is smooth, justifying the following differentiations.
Fix an integer \(\ell\ge2\), set \(r_j=1-j/(2\ell)\), and write
\[
T_j=\max_{|\alpha|=j}\|\partial^\alpha u\|_{L^2(B_{r_j})},
\qquad M=\|u\|_{L^2(B_1)}.
\]
The first-derivative bound gives \(T_1\le C(\ell+A/\ell)M\).
For \(j\le\ell-2\), apply \eqref{eq:wells-local-h2} to an order-\(j\)
derivative between \(B_{r_j}\) and \(B_{r_{j+2}}\).
Leibniz' rule yields
\begin{equation}
T_{j+2}\le C\left(
\sum_{b=0}^j\binom jb A^{b+1}b!\,T_{j-b}
+\ell^2T_j\right).
\label{eq:wells-analytic-recurrence}
\end{equation}
Lower-order derivatives are available on the larger balls
\(B_{r_{j-b}}\). The multi-index coefficient sum of order \(b\)
is \(\binom jb\), by the multinomial identity.
If preceding orders satisfy \(T_s\le M(K\ell)^s\), the sum is at most
\[
AM(K\ell)^j\sum_{b=0}^j(A/K)^b
\le2AM(K\ell)^j
\]
for \(K\ge2A\), using \(\binom jb b!\le\ell^b\).
Increasing the fixed \(K\) absorbs the remaining \(C\ell^2\) factor
into \((K\ell)^2\). Both parity chains therefore satisfy
\(T_j\le M(K\ell)^j\).
Set \(\ell=k+2\). Sobolev embedding on fixed smaller balls bounds the
sup norm of an order-\(k\) derivative on \(B_{1/4}\) by
\(CM[K(k+2)]^{k+2}\).
The inequality \(k!\ge(k/e)^k\), with the small cases absorbed in
the constant, converts this into \eqref{eq:wells-analytic-derivatives}.
Taylor's formula, summed over multi-indices, bounds the degree-\(L\)
remainder on \(|z|\le a\) by a constant times \((C_1a)^{L+1}\).
The first-derivative remainder has an additional factor polynomial
in \(L\) and the same geometric decay. Fixing \(a\) small enough that
\(C_1a<1\) proves the last assertion.
\end{proof}
\subsection{A polynomial-size graded Ritz space}
Graded meshes with increasing polynomial degree are a standard way to
combine singular local behavior with exponential approximation away from
the singularities. For related eigenvalue approximation results with
singular potentials, see Maday and Marcati~\cite{MadayMarcati2019}.
Their bounded-domain discontinuous-Galerkin estimates are not an oracle
for the present moving attractive Coulomb wells. We construct the needed
conforming space and then certify its matrix arithmetic separately.
More specifically, the products of tensor-linear hats with local
polynomials below form a partition-of-unity approximation space, as in
Melenk and Babu\v{s}ka~\cite[Section~2]{MelenkBabuska1996}.
We give the cutoff and derivative estimates explicitly because the tensor
cells can have arbitrarily large aspect ratios.
This is a form-domain approximation: the piecewise polynomial
functions need only belong to \(H^1\). Their kinetic matrix uses first
weak derivatives, not second derivatives across element boundaries.
\begin{lemma}[Exponentially accurate Ritz approximation]
\label{lem:wells-ritz-approximation}
For each \(i\), each rational choice of centers and secondary parameters,
and each sufficiently large \(D\), one can explicitly construct a space
\(\mathcal V_i\subset H^1(\R^3)\) of rational piecewise polynomials,
of dimension polynomial in \(D\), supported in \(r_i<.8D\), such that
\begin{equation}
\inf_{v\in\mathcal V_i}\|v-\psi_i\|_{H^1}\le e^{-.58D}.
\label{eq:wells-ritz-h1}
\end{equation}
Its lowest Ritz value \(\lambda_i^{\mathrm R}\) satisfies
\begin{equation}
0\le\lambda_i^{\mathrm R}-\lambda_i\le C e^{-1.16D}.
\label{eq:wells-ritz-bias}
\end{equation}
The construction and constants are uniform over rational parameters
throughout \([1,4]^m\).
\end{lemma}
\begin{proof}
Work in local coordinates \(z=y-X_i\) on \([-D,D]^3\), with poles
at \(0\) and \(\rho_i e_z\). Put \(\delta=2^{-8D}\) and fix a sufficiently
small rational \(\xi>0\), as specified below. On each coordinate axis
let \(B\) be the set of projections of the two poles. Starting from
\([-D,D]\), bisect dyadically until every interval \(I\) obeys
\begin{equation}
|I|\le\xi\min\{1,\dist(I,B)+\delta\}.
\label{eq:wells-grid-rule}
\end{equation}
Every leaf length is at least \(\xi\delta/2\), since its parent failed
the stopping rule. Above a fixed spacing there are \(O(D)\) intervals.
At a finer level of length \(h\), only intervals within \(h/\xi\) of
one of the at most two projections can still split, so a fixed number
of intervals per level are split there.
There are \(O(\log(D/\delta))=O(D+\log D)\) levels, hence \(O(D)\)
intervals per axis. This remains true for poles arbitrarily close to
dyadic boundaries. The tensor product has \(O(D^3)\) nodes and cells.
Let \(b_w\) be the tensor-linear hat at node \(w\). The full hats form
a partition of unity in the cube. Retain nodes satisfying
\begin{equation}
|w|<.73D,\qquad |w|>\delta,\qquad
|w-\rho_i e_z|>\delta,
\label{eq:wells-retained-nodes}
\end{equation}
and take \(b_w(z)(z-w)^\alpha\) for all \(|\alpha|\le L_0\), where \(L_0\)
is a sufficiently large fixed multiple of \(D\). These functions span
\(\mathcal V_i\), whose dimension before removing dependencies is
\(O(D^6)\).
Here is the support estimate that avoids any aspect-ratio assumption.
For either pole \(a\), every coordinate interval adjacent to \(w\) obeys
\[
|I|\le\xi\min\{1,|w_j-a_j|+\delta\}.
\]
The support radius \(R_w^{\mathrm{hat}}\) of a retained hat thus satisfies
\begin{equation}
\begin{split}
R_w^{\mathrm{hat}}
&\le\xi\bigl(|w-a|+\sqrt3\delta\bigr)
\le(1+\sqrt3)\xi|w-a|,\\
R_w^{\mathrm{hat}}&\le\sqrt3\xi.
\end{split}
\label{eq:wells-patch-radius}
\end{equation}
For \(d_w=\min\{1,|w|,|w-\rho_i e_z|\}\) this gives
\(R_w^{\mathrm{hat}}\le C\xi d_w\). Each support avoids both poles and,
for large \(D\), lies in \(|z|<.8D\).
Choose a fixed \(c<1/10\) and set \(R=c d_w\).
The whole ball \(B(w,R)\) avoids both poles. Throughout it every remote
primary has distance at least \(D-2S-.73D-c>D/8\), so it also avoids
every remote cap transition. On this enlarged ball the eigen-equation
contains only pure Coulomb potentials. Rescaling by \(R\) gives
\(\Delta u=Fu\) on the unit ball. Its coefficient satisfies
\eqref{eq:wells-analytic-coefficient} with an absolute \(A\):
local inverse-distance expansions give factorial derivative bounds
because \(R/\dist(w,a)\le c\); the extra \(R^2\) in the coefficient
makes \(R^2/\dist(w,a)\) bounded; and remote contributions sum to at
most \(Cm/D\). The eigenvalue term is bounded too.
The sup norm of \(u\) is uniform by \Cref{prop:isolated-well}.
Choose \(\xi\) sufficiently small compared to \(c\) and the fixed
analytic radius in \Cref{lem:wells-analytic}, putting the hat supports
inside that smaller ball after rescaling. With \(L_0\) a sufficiently
large fixed multiple of \(D\), its Taylor polynomial \(p_w\) satisfies
\begin{equation}
|\psi_i-p_w|+R|\nabla\psi_i-\nabla p_w|\le e^{-20D}
\quad\text{on }\supp b_w.
\label{eq:wells-taylor-error}
\end{equation}
Increasing the fixed multiple absorbs the polynomial factor in the
differentiated remainder. The Taylor coefficients establish existence
of an approximant; the algorithm will not compute them.
To remove tiny neighborhoods of the poles and the exterior, let
\(a_{\mathrm{out}}\) equal one up to radius \(.66D\), decrease linearly
to zero by \(.72D\), and vanish thereafter. Let \(a_{\mathrm{in}}\) vanish
up to \(\delta\), increase linearly to one by \(2\delta\), and equal one
thereafter. Sample
\[
a_{\mathrm{out}}(|z|)\,
a_{\mathrm{in}}(|z|)\,
a_{\mathrm{in}}(|z-\rho_i e_z|)
\]
at all grid nodes to get weights \(c_w\in[0,1]\), and put
\(b=\sum_w c_wb_w\). Nonzero weights occur only at retained nodes.
The function \(b\) equals one away from \(O(\delta)\)-balls about the
poles and away from \(|z|>.65D\): a node within \(2\delta\) of a pole
has adjacent lengths \(O(\xi\delta)\), whereas all intervals have
length at most \(\xi\). These facts localize the inner and outer
transitions, respectively.
Each first derivative of a tensor-linear interpolant is a convex
combination of edge secants in that coordinate. Consequently
\(|\nabla b|\le C/\delta\) on the inner transitions and
\(|\nabla b|\le C/D\) on the outer transition, irrespective of aspect
ratios. The bounded sup norm of \(\psi_i\) makes the inner
cutoff-gradient cost \(C\delta^{1/2}\). Also, by its \(H^2\) bound,
\[
\|\nabla\psi_i\|_{L^2(B_{C\delta}(a))}
\le C\delta\|\nabla\psi_i\|_6\le C\delta.
\]
The exterior cost follows from \eqref{eq:wells-h1-tail}. Thus
\begin{equation}
\|\psi_i-b\psi_i\|_{H^1}
\le C\bigl(\sqrt\delta+e^{-.61D}\bigr).
\label{eq:wells-cutoff-error}
\end{equation}
The cutoff weights, like the Taylor coefficients, are only existence
coefficients in the explicitly generated space.
Set \(v=\sum_w c_wb_wp_w\in\mathcal V_i\).
At almost every point at most eight hats are nonzero. Their derivatives
cost at most \(C/\delta\), by the minimum interval length, and
\(R^{-1}\le C/\delta\).
On a region of volume \(O(D^3)\), \eqref{eq:wells-taylor-error} gives
\[
\|b\psi_i-v\|_{H^1}\le C D^{3/2}\delta^{-1}e^{-20D}.
\]
Since \(\sqrt\delta=e^{-4(\log2)D}\) and
\(\delta^{-1}e^{-20D}=e^{-(20-8\log2)D}\),
this and \eqref{eq:wells-cutoff-error} prove
\eqref{eq:wells-ritz-h1} for large \(D\).
Finally, the form of \(h_i-\lambda_i\) is bounded in absolute value
by a constant times the \(H^1\) norm squared. For \(e=v-\psi_i\),
the eigen-equation cancels the cross terms, giving
\[
\frac{\langle v,h_iv\rangle}{\|v\|_2^2}-\lambda_i
=\frac{\langle e,(h_i-\lambda_i)e\rangle}{\|v\|_2^2}
\le C\|e\|_{H^1}^2.
\]
Here \(\|v\|_2\) is bounded away from zero.
The variational principle supplies the opposite inequality for the
lowest Ritz value, proving \eqref{eq:wells-ritz-bias}.
\end{proof}
\subsection{Certified polynomial-bit eigenvalue evaluation}
A small Ritz space by itself does not establish bit complexity.
We next bound the conditioning loss and give certified matrix
entries and eigenvalue enclosures. All arithmetic below is classical.
\begin{proposition}[Uniform single-well eigenvalue oracle]
\label{prop:ritz-oracle}
Suppose the centers and the parameters in \([1,4]^m\) are rational,
and all their numerators and denominators, as well as those of
\(q,q'\), have at most \(B\) bits. Under the hypotheses of this section,
a deterministic algorithm returns rational numbers
\(\ell_i\le\lambda_i\le u_i\) with
\begin{equation}
u_i-\ell_i\le\varepsilon_D,\qquad
\varepsilon_D=2^{-\lceil8D/5\rceil}1.1>1.01\).
Each scalar bisection restarts from the rational interval \([1,4]\)
and makes at most \(O(D)\) midpoint operations. Its samples,
including early-return samples, therefore have \(O(D)\) fractional
bits independently of the number of simultaneous iterations.
There are \(O(mD^2)\) oracle calls. Their data bit lengths consist
of the fixed-center and rounded-charge bit lengths together with
these \(O(D)\) parameter bits, so \Cref{prop:ritz-oracle} proves
the total polynomial-time bound.
In particular the charge-rounding error is part of every evaluated
operator and is compensated by tuning; it is not introduced as an
uncontrolled perturbation afterward.
\end{proof}
==== END Continuum-Coulomb-hardness-with-binary-nuclear-charges//02-wells.tex ====
==== BEGIN Continuum-Coulomb-hardness-with-binary-nuclear-charges//03-hopping.tex ====
\section{Relative tails and computable hopping calibration}
\label{sec:hopping}
The absolute orbital estimate of Proposition~\ref{prop:isolated-well}
does not determine an overlap of order $e^{-D}$. This section proves a
relative tail estimate, identifies positive integrals that determine the
hoppings, and shows how to calibrate those integrals before tuning the
single-well energies. All constants below are independent of the site
layout, $D$, the displacements, and the parameters $\rho_i\in[1,4]$.
The cap profile $\chi$ is fixed throughout.
We use the positions and operators from Section~\ref{sec:wells}. In
particular,
\[
X_i=Dp_i+s_i,\qquad \max_i|s_i|\le S,
\qquad |X_i-X_j|\ge D-2S\quad(i\ne j).
\]
We take $D$ sufficiently large and, for the estimates below, may assume
$m+S\le D^{1/4}$. Where a relative error is expanded, its displayed
bound is assumed smaller than a fixed positive constant. The final
parameter choice in Section~\ref{sec:completion} satisfies these
conditions and makes every such error as small as required.
\subsection{Relative comparison of single-well tails}
For $r_i=|y-X_i|$ define the ray primitive
\begin{equation}
F_i(y)=r_i\int_0^1 W_i\bigl(X_i+t(y-X_i)\bigr)\,\dd t,
\qquad F_i(X_i)=0.
\label{eq:ray-primitive}
\end{equation}
The integral follows the entire segment from $X_i$ to $y$; it is not
a local approximation to the field at either endpoint. Write
\begin{equation}
\varepsilon_{\mathrm{rel}}
=\frac{m\log D+m^2}{D}.
\label{eq:relative-tail-error}
\end{equation}
\begin{proposition}[Uniform relative tail]
\label{prop:relative-tail}
Suppose the actual rounded charges and the parameters $\rho_i$ obey
the conclusions of Proposition~\ref{prop:tuning}, in particular
$|\lambda_i+1/2|\le e^{-1.01D}$. Uniformly for $r_i\le2D$,
\begin{equation}
\psi_i(y)=\pi^{-1/2}e^{-r_i+F_i(y)}
\bigl(1+O(\varepsilon_{\mathrm{rel}})\bigr).
\label{eq:relative-tail}
\end{equation}
\end{proposition}
\begin{proof}
Fix a site, suppress its subscript, and write $r=|y-X|$ and
$\omega=(y-X)/r$ when $r>0$. The capped-field estimates in
Proposition~\ref{prop:isolated-well} give
\begin{equation}
\|W\|_\infty\le \frac{Cm}{D},\qquad
\|\nabla W\|_\infty\le \frac{Cm}{D^2},\qquad
\|D^2W\|_\infty\le \frac{Cm}{D^3}.
\label{eq:hopping-field-derivatives}
\end{equation}
Set $A(y)=\int_0^1W(X+t(y-X))\,\dd t$, so $F=rA$.
For $010$,
\[
\frac{|v|}{|v-\rho_j e_z|}
\le\frac{|v|}{|v|-4}\le\frac53.
\]
For $|v|\le10$, the hole weight is one when $D$ is large. The
ray derivative estimates, also valid for $F^*$, show that $L_X$
is Lipschitz on this ball with a fixed constant. This statement
requires no differentiability of $|v|$ at zero; it follows either
by integrating the bounded derivatives away from zero or directly
from the ray formula. Thus $e^{L_X(v)}$ is within fixed factors
of $e^{L_X(0)}$ throughout the ball. Moreover,
\[
\int_{|v|\le10}\frac{\dd^3v}{|v-\rho_j e_z|}
\le\int_{|w|\le14}\frac{\dd^3w}{|w|}<\infty,
\qquad
\int_{|v|\le1}\frac{\dd^3v}{|v|}>0.
\]
The secondary-kernel integral on the inner ball is therefore at
most a fixed multiple of the primary-kernel integral. Combining
the two regions and multiplying by $q'=O(m/D)$ gives a relative
$O(m/D)$ correction. Positivity allows the uniform relative
wave-function errors to be taken outside both integrals. These
errors are included in \eqref{eq:hole-surrogate-comparison}.
\end{proof}
\begin{lemma}[Size of the reference integrals]
\label{lem:reference-hopping}
Let $X^0=(Dp_i)_i$. For every contact,
\begin{equation}
c e^{-D}\le T_{ij}(X^0)\le C D^2 e^{-D+Cm}.
\label{eq:reference-hopping-bounds}
\end{equation}
\end{lemma}
\begin{proof}
Put $u=u_{ij}=p_j-p_i$, so $|u|=1$. The triangle inequality
gives $|Du+v|+|v|\ge D$. Each ray correction in the integral is
nonnegative and at most $Cm$, by \eqref{eq:capped-kernel-bounds}
and its ray length. The integral of $|v|^{-1}$ over the ball of
radius $D/8$ is $O(D^2)$, proving the upper bound. On $|v|\le1$
the hole weight is one, $|Du+v|+|v|\le D+2$, and the corrections
are nonnegative. Integration over this fixed ball gives the lower
bound; for example $2e^{-D-2}$ is a valid lower bound for large $D$.
\end{proof}
\Needspace{10\baselineskip}
\subsection{Sensitivity to displacement}
\begin{lemma}[Uniform displacement law]
\label{lem:hopping-displacement}
Let $X_i=Dp_i+s_i$ with $|s_i|\le S$. There is a constant
$c_1>0$ such that, uniformly on all contacts,
\begin{align}
\frac{T_{ij}(X)}{T_{ij}(X^0)}
&=e^{-u_{ij}\cdot(s_j-s_i)}
\bigl(1+O(\varepsilon_{\mathrm{disp}})\bigr),
\label{eq:hopping-displacement-law}\\
\varepsilon_{\mathrm{disp}}
&=\frac{mS}{D}+\frac{S^2}{D}
+S D^{-1/4}+e^{-c_1\sqrt D}.
\label{eq:hopping-displacement-error}
\end{align}
In particular the coefficient of the projected displacement in the
leading exponent is one.
\end{lemma}
\begin{proof}
Compare the two integrals at the same $v$. In the ray defining
$F_i^*(X_j+v)$, the reference path point
$Dp_i+t(Dp_j+v-Dp_i)$ moves by $(1-t)s_i+ts_j$, of norm at
most $S$. Each field center moves by at most $S$ as well.
The globally Lipschitz bound in \eqref{eq:capped-kernel-bounds}
therefore changes the field along the path by at most $CmS/D^2$.
The ray has length $O(D)$ and its length changes by at most $2S$.
Since the field itself is $O(m/D)$,
\begin{equation}
|F_i^*(X_j+v;X)-F_i^*(Dp_j+v;X^0)|
\le\frac{CmS}{D}.
\label{eq:ray-displacement}
\end{equation}
The same bound holds for the $F_j^*$ term; its ray length $|v|$
does not change. The global nature of the capped-kernel bound is
essential here: it also applies when a moving ray crosses a hole
or a transition interface.
Write $u=u_{ij}$, $\delta=s_j-s_i$, and $v=zu+w$ with
$w\perp u$. Since $|v|\le D/8$,
$|Du+v|\ge7D/8$. The Hessian of the norm on the segment from
$Du+v$ to $Du+v+\delta$ has norm $O(1/D)$. Its gradient at
$Du+v$ differs from $u$ by $O(|w|/D)$. Taylor's formula gives
\begin{equation}
|Du+v+\delta|-|Du+v|
=u\cdot\delta+O\left(\frac{S|w|}{D}+\frac{S^2}{D}\right).
\label{eq:distance-displacement}
\end{equation}
This estimate includes the whole backward segment $v=-tu$,
$0\le t\le D/8$: there $Du+v=(D-t)u$ and the gradient is
exactly $u$. Thus the segment's non-negligible longitudinal
extent creates neither a different coefficient nor an $O(S)$
transverse error.
On the region $|w|\le D^{3/4}$, Equations
\eqref{eq:ray-displacement} and \eqref{eq:distance-displacement}
show that
\begin{equation}
L_X(v)-L_{X^0}(v)+u\cdot\delta
=O\left(\frac{mS}{D}+\frac{S^2}{D}+S D^{-1/4}\right).
\label{eq:main-region-log-error}
\end{equation}
The nonnegative integration weight is identical in the two
integrals, so exponentiation controls their contributions on
this region by the same relative factor.
On the remaining region $|w|>D^{3/4}$, one has
\begin{equation}
\begin{split}
|Du+v|+|v|-D
&\ge z+|v|
=\frac{|w|^2}{|v|-z}\\
&\ge\frac{4|w|^2}{D}
\ge4\sqrt D.
\end{split}
\label{eq:transverse-excess}
\end{equation}
The identity is used only for $w\ne0$, and the denominator is
at most $2|v|\le D/4$. The field corrections are bounded by
$Cm$, without any angular constancy assumption. Displacement
changes a distance by at most $2S$, and $|u\cdot\delta|\le2S$.
It follows that the omitted reference and displaced contributions,
after division by $e^{-D}$ and $e^{-D-u\cdot\delta}$ respectively,
are at most
\begin{equation}
C D^2\exp\bigl[-c\sqrt D+C(m+S)\bigr].
\label{eq:transverse-integral-error}
\end{equation}
By \eqref{eq:reference-hopping-bounds}, these are also bounds on
their contributions relative to the reference integral, with
the appropriate displacement factor. The assumption
$m+S\le D^{1/4}$ and sufficiently large $D$ absorb both the
possibly large anisotropic factor $e^{Cm}$ and the polynomial
volume factor into $e^{-c_1\sqrt D}$. Combining this with
\eqref{eq:main-region-log-error}, and using smallness of its
right-hand side, proves the lemma.
\end{proof}
\subsection{Deterministic relative quadrature}
We have reduced hopping calibration to positive integrals whose size is
exponentially small. Direct absolute quadrature at that scale would
obscure the complexity bound. Instead we control the oscillation of the
logarithm of the exponential factor, while integrating its nonnegative
radial weight exactly. Positivity then converts local errors into a
relative error for the whole integral.
\begin{lemma}[Computing the positive reference integrals]
\label{lem:hopping-quadrature}
For $0<\epsilon<1/10$, one can compute a positive rational
$\widehat I_{ij}$ satisfying
\begin{equation}
\left|\log\frac{\widehat I_{ij}}{e^D T_{ij}(X^0)}\right|
\le\epsilon
\label{eq:reference-quadrature-accuracy}
\end{equation}
by a deterministic algorithm polynomial in $D$, $m$, the bit
lengths of the layout data, and $\epsilon^{-1}$. In particular
the algorithm is polynomial in the source size when $D$ is a
fixed polynomial and $\epsilon$ an inverse polynomial.
\end{lemma}
\begin{proof}
Let $L=L_{X^0}$ in \eqref{eq:hopping-log-integrand}. On the
integration ball, the two distance terms have Lipschitz constants
one, and the ray formulas and \eqref{eq:capped-kernel-bounds} give
a Lipschitz constant $Cm/D$ for each correction. Hence
\begin{equation}
|L(v)-L(v')|\le C|v-v'|,
\qquad |L(v)|\le C(D+m).
\label{eq:hopping-log-regularity}
\end{equation}
Only the exponential factor has been included in this logarithm.
No lower bound on $1-\chi$ or differentiability of $\log(1-\chi)$
is needed.
We describe a mesh with exact positive weights. Parameterize
the unit sphere, up to a set of area zero, by
\[
\omega(\mu,\varphi)
=\bigl(\sqrt{1-\mu^2}\cos\varphi,
\sqrt{1-\mu^2}\sin\varphi,\mu\bigr).
\]
Here $-1\le\mu\le1$, $0\le\varphi\le2\pi$, and
$\dd\omega=\dd\mu\,\dd\varphi$. Use $M$ equal intervals
in $\mu$, $N$ equal intervals in $\varphi$, and a rational radial
partition with mesh at most $h$. Insert all the finitely many
breakpoints of $\chi(r/D)$ into that partition. On a radial
interval $I_a$, the weight
\[
w_a=\int_{I_a}(1-\chi(r/D))r\,\dd r
\]
is a nonnegative exact rational, computable by polynomial
integration. Discard intervals with zero weight. The normalized
angular weight of every angular rectangle is exactly
\begin{equation}
\frac1\pi\,\frac2M\,\frac{2\pi}{N}=\frac4{MN}.
\label{eq:exact-angular-weight}
\end{equation}
Choose one sample $v_{abc}$ in each resulting spatial cell, for
example at the parameter midpoints.
The square-root parameterization satisfies
\[
|\omega(\mu,\varphi)-\omega(\mu',\varphi')|
\le C|\mu-\mu'|^{1/2}+C|\varphi-\varphi'|.
\]
Therefore $M=O((D/h)^2)$ and $N=O(D/h)$, with sufficiently
large fixed constants, make every spatial cell have diameter
$O(h)$, including cells adjacent to the poles. There are
$O((D/h)^4)$ cells, a polynomial number. Their exact weighted
sample sum is
\[
Q=\sum_{a,b,c}\frac{4w_a}{MN}\,e^{L(v_{abc})}.
\]
Positivity and \eqref{eq:hopping-log-regularity} imply
\begin{equation}
e^{-Ch}Q\le T_{ij}(X^0)\le e^{Ch}Q.
\label{eq:positive-quadrature-bracket}
\end{equation}
Thus a fixed multiple of $\epsilon$ is sufficient for $h$.
This relative estimate is independent of the size $e^{-D}$ of
the integral.
It remains to specify evaluation and bit accuracy. For a ray of
length $r\le2D$, the integrand
$t\mapsto W_i^*(X_i+t(y-X_i))$ has Lipschitz constant
at most $Cmr/D^2$. A midpoint rule with $K$ equal panels
therefore approximates its ray primitive with error at most
\begin{equation}
C\frac{mr^2}{D^2K}\le\frac{Cm}{K}.
\label{eq:ray-quadrature-error}
\end{equation}
Taking $K=O(m/\epsilon)$ and evaluating each field sample to
absolute error $O(\epsilon/D)$ gives an $O(\epsilon)$ error
in each primitive. Each field sample is a sum of $m-1$
globally Lipschitz capped kernels. At small radius the kernel
is evaluated as zero, and elsewhere its denominator is bounded
below by a fixed multiple of $D$; equivalently one evaluates the
fixed continuous function $\chi(a)/a$, extended by zero on
$[0,1/16]$, at $a=|z|/D$. Its global Lipschitz bound controls
errors even if the approximate radius lies on the other side of
a transition breakpoint. No crossing location needs to be found.
The coordinates $Dp_i$, sample coordinates, square roots, and
trigonometric values can be computed to inverse-polynomial
absolute accuracy with polynomially many bits. The global
kernel and ray bounds just proved control the resulting errors.
Thus all sample logarithms can be approximated to an arbitrarily
specified $O(\epsilon)$ absolute error with polynomial work.
Such logarithm errors change the positive sample sum by a factor
$e^{O(\epsilon)}$, not by an additive error multiplied by
an uncontrolled exponential.
For the remaining arithmetic, compute the scaled sum $e^D Q$,
using sample exponentials $e^{D+L(v_{abc})}$. Its exact value
is bounded below by a positive constant up to the factor in
\eqref{eq:positive-quadrature-bracket}, because
\eqref{eq:reference-hopping-bounds} gives
$e^D T_{ij}(X^0)\ge2e^{-2}$. The sum of its non-exponential
weights is $O(D^2)$, and all exponent arguments have magnitude
$O(D+m)$. After the logarithms have been approximated as above,
each elementary exponential can be rounded to sufficiently small
absolute error $O(\epsilon/D^2)$, clipping a negative rounded
value to zero; this adds $O(\epsilon)$ to the scaled sum. These
last approximations require only
$O(D+m+\log(D/\epsilon))$ precision bits, with polynomial
working precision. This step does not require exponentially
accurate ray quadrature: the earlier ray errors have already
been bounded multiplicatively. Standard range reduction and
Taylor series compute the square roots, trigonometric values,
exponentials, and logarithms to these precisions in polynomial
time. The exact rational weights have polynomial bit lengths,
and the sum contains polynomially many terms.
Allocate a sufficiently small fixed fraction of $\epsilon$ to
each of the spatial, ray, coordinate, and arithmetic errors.
The resulting positive rational sum is $\widehat I_{ij}$ and
satisfies \eqref{eq:reference-quadrature-accuracy}.
\end{proof}
\subsection{Calibration, rounding, and orientation}
The oriented attraction integrals are close to a symmetric matrix
for an exact reason that is useful independently of their asymptotic
formula.
\begin{lemma}[Orientation identity]
\label{lem:hopping-orientation}
For distinct sites,
\begin{equation}
B_{ij}-B_{ji}
=(\lambda_i-\lambda_j)\langle\psi_j,\psi_i\rangle.
\label{eq:hopping-orientation-identity}
\end{equation}
For tuned wells and large $D$ under the stated separation bounds,
\begin{equation}
|\langle\psi_j,\psi_i\rangle|\le e^{-0.85D},\qquad
|B_{ij}-B_{ji}|\le2e^{-1.86D}.
\label{eq:hopping-orientation-error}
\end{equation}
\end{lemma}
\begin{proof}
The single-well operators have common domain $H^2$, and
\[
h_i-h_j=V_j^{\mathrm{hole}}-V_i^{\mathrm{hole}}.
\]
Pair this identity between the real eigenfunctions $\psi_j$ and
$\psi_i$, using self-adjointness of $h_j$, to obtain
\eqref{eq:hopping-orientation-identity}. All the pairings are
well-defined by the Coulomb domain bounds of
Proposition~\ref{prop:isolated-well}. Equivalently, in the two
expansions of the full matrix element
$\langle\psi_j,(h+1/2)\psi_i\rangle$, contributions from every
hole other than $i$ and $j$ cancel identically, because they
integrate the same real product against the same multiplication
potential.
To record an explicit overlap bound, Lemma~\ref{lem:well-tail}
gives
\[
|\psi_i(y)\psi_j(y)|
\le C e^{0.0198D}e^{-0.99(r_i+r_j)}.
\]
Use $r_i+r_j\ge D-2S$ on a factor with exponent $0.90$, and
retain $e^{-0.09r_i}$ for integrability. It follows that
\[
|\langle\psi_j,\psi_i\rangle|
\le C e^{-0.8802D+1.8S}
\int_{\mathbb R^3}e^{-0.09|y-X_i|}\,\dd^3y
\le e^{-0.85D}
\]
for sufficiently large $D$. Finally
$|\lambda_i-\lambda_j|\le2e^{-1.01D}$ proves the second bound.
\end{proof}
For a tolerance $\tau>0$, define the total calibration error
\begin{equation}
\begin{split}
\mathcal E_{\mathrm{cal}}
={}&\frac{m\log D+m^2}{D}
+\frac{mS}{D}+\frac{S^2}{D}\\
&+S D^{-1/4}+\tau+e^{-c_1\sqrt D}.
\end{split}
\label{eq:total-calibration-error}
\end{equation}
\begin{proposition}[Deterministic hopping calibration]
\label{prop:calibration}
Let the contact layout be as in Lemma~\ref{lem:contact-layout},
choose one orientation per contact, and let the prescribed rational
hoppings satisfy $00\), independent of the
construction parameters, such that
\begin{equation}
h_i-\lambda_i\ge
g_{\mathrm w}\bigl(I-|\psi_i\rangle\langle\psi_i|\bigr).
\label{eq:continuum-well-gap}
\end{equation}
All operator inequalities involving an unbounded operator are
understood as quadratic-form inequalities. By
Lemma~\ref{lem:contact-layout} and the displacement bound,
\begin{equation}
|X_i-X_j|\ge D-2S\quad(i\ne j),\qquad
|X_i-X_j|\ge\sqrt2D-2S\quad(\{i,j\}\notin\mathcal C).
\label{eq:continuum-separations}
\end{equation}
The second inequality concerns distinct noncontact sites.
\subsection{The one-particle matrix}
Recall the hole potentials from \cref{sec:hopping}, and write
\[
A=h+\tfrac12,\qquad
V_i^{\mathrm{hole}}=(1-\chi_i)V_i,\qquad
\epsilon_i=\lambda_i+\tfrac12.
\]
For the chosen orientation \(i\to j\) of a contact, define
\[
B_{ij}=\langle\psi_j,V_j^{\mathrm{hole}}\psi_i\rangle,\qquad
\epsilon_{\mathrm{hop}}
=\max_{i\to j}\bigl|Z_*B_{ij}-t_{ij}\bigr|.
\]
The maximum is zero if there are no contacts.
Proposition~\ref{prop:calibration} gives
\begin{equation}
\epsilon_{\mathrm{hop}}\le C\mathcal E_{\mathrm{cal}},
\label{eq:continuum-calibration-interface}
\end{equation}
because all the target hoppings are at most one.
Let \(\mathsf T\) be the real symmetric \(m\)-by-\(m\) matrix with
entry \(-t_{ij}\) on contacts and zero elsewhere, including its
diagonal.
Define \(\mathcal U:\mathbb C^m\to L^2(\mathbb R^3)\) by
\(\mathcal U e_i=\psi_i\), and let \(G=\mathcal U^*\mathcal U\)
be its Gram matrix. The notation \(\mathcal U\) denotes this map;
the scalar onsite repulsion remains \(U=5/8\).
\begin{lemma}[Residuals and the orthonormalized matrix]
\label{lem:one-body-matrix}
For the preceding parameters and sufficiently large \(D\),
\begin{equation}
\|A\psi_i\|_2\le e^{-.75D},\qquad
\|G-I\|\le e^{-.85D}.
\label{eq:continuum-residual-gram}
\end{equation}
In particular \(G\) is positive definite.
The isometry
\(\mathcal J_0=\mathcal U G^{-1/2}\) has orthonormal columns
\(\omega_i=\mathcal J_0e_i\), spanning the same spatial orbital
space. Its matrix satisfies
\begin{align}
\|Z_*\mathcal J_0^*A\mathcal J_0-\mathsf T\|
&\le \delta_1, \label{eq:continuum-one-body-error}\\
\delta_1
&:=m\epsilon_{\mathrm{hop}}+Z_*e^{-1.01D}
+CmZ_*e^{-1.1D}+CmZ_*e^{-1.6D}. \notag
\end{align}
If \(P_0=\mathcal J_0\mathcal J_0^*\), then
\begin{equation}
\|AP_0\|\le C\sqrt m\,e^{-.75D},
\qquad
\|\omega_i-\phi_i\|_{H^1}\le Cm/D.
\label{eq:continuum-projected-residual}
\end{equation}
\end{lemma}
\begin{proof}
Since \(h=h_i-\sum_{k\ne i}V_k^{\mathrm{hole}}\),
\begin{equation}
A\psi_i=\epsilon_i\psi_i-
\sum_{k\ne i}V_k^{\mathrm{hole}}\psi_i.
\label{eq:continuum-residual-identity}
\end{equation}
The \(k\)-th hole is supported in \(r_k\le D/8\). On it,
\[
r_i\ge 7D/8-2S\quad(k\ne i).
\]
Both Coulomb poles of \(V_k\) lie in this hole for large \(D\).
The charges \(q,q'\) are bounded, and direct radial integration,
enlarging the integration ball by at most four for the secondary
pole, gives
\begin{equation}
\|V_k^{\mathrm{hole}}\|_2\le C D^{1/2},
\qquad
\|V_k^{\mathrm{hole}}\|_1\le C D^2.
\label{eq:continuum-hole-norms}
\end{equation}
Consequently the sum in \eqref{eq:continuum-residual-identity}
has norm at most
\[
CmD^{1/2}\exp(-.85635D+1.98S).
\]
Together with the tuning error, this is at most \(e^{-.75D}\)
for the stated regime and sufficiently large \(D\).
For \(i\ne j\), put \(\ell_{ij}=|X_i-X_j|\).
The product of the two tail bounds is at most
\[
C\exp\bigl(-.99(r_i+r_j)+.0198D\bigr).
\]
Use \(r_i+r_j\ge \ell_{ij}\) on a factor with exponent \(.90\);
the remaining factor is integrable uniformly in the centers,
since
\(\int e^{-.09(r_i+r_j)}\,\dd y\le
\int e^{-.09r_i}\,\dd y\le C\).
Thus
\[
|G_{ij}|\le C e^{-.90\ell_{ij}+.0198D}
\le C e^{-.8802D+1.8S}.
\]
The diagonal entries are one. The row-sum estimate and
\(m+S\le D^{1/4}\) prove the second inequality in
\eqref{eq:continuum-residual-gram}.
Let \(M=\mathcal U^*A\mathcal U\), which is real symmetric.
We next estimate its entries before orthonormalization.
If \(k\ne i,j\), both \(r_i\) and \(r_j\) are at least
\(7D/8-2S\) on the \(k\)-th hole. Equations
\eqref{eq:continuum-tail-input} and
\eqref{eq:continuum-hole-norms} therefore imply
\[
\sum_{k\ne i,j}
|\langle\psi_j,V_k^{\mathrm{hole}}\psi_i\rangle|
\le CmD^2e^{-1.7127D+3.96S}.
\]
The same estimate bounds
\(\sum_{k\ne i}\langle\psi_i,V_k^{\mathrm{hole}}\psi_i\rangle\).
For a noncontact pair, the remaining hole at \(j\) is bounded using
\(r_i+r_j\ge\sqrt2D-2S\), giving
\[
|\langle\psi_j,V_j^{\mathrm{hole}}\psi_i\rangle|
\le CD^2
e^{-(.99\sqrt2-.0198)D+1.98S}.
\]
Here \(.99\sqrt2-.0198>1.38\).
All these bounds are smaller than \(e^{-1.1D}\), after
absorbing their polynomial factors and fixed constants.
For \(i\ne j\), the term \(\epsilon_iG_{ji}\) is at most
\(e^{-1.86D}\). It follows, adjusting a fixed constant, that
\begin{align}
|M_{ii}-\epsilon_i|&\le Ce^{-1.1D},\notag\\
|M_{ji}+B_{ij}|&\le Ce^{-1.1D}
&&(i\to j\text{ a chosen contact}),\label{eq:continuum-raw-entries}\\
|M_{ij}|&\le Ce^{-1.1D}
&&(i\ne j,\ \{i,j\}\notin\mathcal C).\notag
\end{align}
The exactly symmetric matrix \(M\) requires only one of the two
orientations of each contact. Its diagonal tuning contribution
has norm at most \(e^{-1.01D}\), not \(m e^{-1.01D}\).
Entrywise comparison with \(\mathsf T\) and a row-sum bound now give
\[
\|Z_*M-\mathsf T\|
\le m\epsilon_{\mathrm{hop}}
+Z_*e^{-1.01D}+CmZ_*e^{-1.1D}.
\]
For completeness, the residual estimate controls the entire matrix,
without using a dimension-dependent bound on a many-body space:
\[
\|A\mathcal U\|\le\sqrt m\,e^{-.75D},\qquad
\|\mathcal U\|=\|G\|^{1/2}\le2,\qquad
\|M\|\le C\sqrt m\,e^{-.75D}.
\]
Functional calculus on the positive matrix \(G\) gives
\(\|G^{-1/2}-I\|\le Ce^{-.85D}\) and
\(\|G^{-1/2}\|\le2\). Consequently
\begin{equation}
\|G^{-1/2}MG^{-1/2}-M\|
\le C\sqrt m\,e^{-1.6D}
\le Cm e^{-1.6D}.
\label{eq:continuum-gram-correction}
\end{equation}
This orthonormalizes the matrix of the shifted operator
\(h+1/2\). In the original nonorthogonal modes,
\(M=\mathcal U^*h\mathcal U+\tfrac12G\).
The shift is essential in this estimate: one must add
\(\tfrac12G\), rather than \(\tfrac12I\), before
orthonormalization.
Equation~\eqref{eq:continuum-one-body-error} follows.
The same residual estimate applied to
\(A\mathcal J_0=A\mathcal U G^{-1/2}\) proves the first part of
\eqref{eq:continuum-projected-residual}.
Finally, as a map into \(H^1\), \(\mathcal U\) has norm at most
\(C\sqrt m\), by \eqref{eq:continuum-orbital-input}. Thus
\[
\|\omega_i-\psi_i\|_{H^1}
\le C\sqrt m\,e^{-.85D}.
\]
Combining this with \(\|\psi_i-\phi_i\|_{H^1}\le Cm/D\)
gives the remaining assertion.
\end{proof}
We now include spin. Set
\[
\mathfrak h=L^2(\mathbb R^3;\mathbb C^2),\qquad
P=P_0\otimes I_2,\qquad
\mathcal J=\mathcal J_0\otimes I_2.
\]
Thus \(P_0\) has rank \(m\), whereas \(P\) has rank \(2m\).
The operator \(A\) acts trivially on spin; we use the same symbol
for this extension. Its residual and matrix bounds are unchanged.
\subsection{A gap on the entire one-particle complement}
Spatial localization separates the nuclei from the exterior, where the
attraction is uniformly weak. We use the IMS product-rule identity
(compare Simon~\cite[Lemma~3.1]{Simon1983}); its kinetic-error coefficient
is $1/2$ for our normalization $-\Delta/2$. The identity is derived below.
The additional point to verify is that a globally orthogonal function
has only exponentially small components along the local ground states
after localization.
\begin{proposition}[Continuum complement gap]
\label{prop:complement-gap}
There is a constant \(g>0\) such that, for all sufficiently large
\(D\),
\[
(I-P)A(I-P)\ge g(I-P).
\]
The inequality holds on the full one-particle form domain,
including functions supported between the nuclei or at infinity.
\end{proposition}
\begin{proof}
Choose smooth real functions \(\theta_0,\theta_1,\ldots,\theta_m\)
with
\[
\sum_{a=0}^m\theta_a^2=1,\qquad
\theta_i=1\ \text{on }r_i\le D/5,\qquad
\supp\theta_i\subseteq\{r_i\le D/4\}.
\]
The balls of radius \(D/4\) are disjoint by
\eqref{eq:continuum-separations}. In the annulus
\(D/5D/8.
\]
Hence \(\chi_k=1\) throughout that support and \(h=h_i\) there.
Choose a fixed
\(04(b+h+r+2),
\end{equation}
and put
\begin{equation}\label{eq:choice-D}
D=A(n+2)^K,
\qquad
\tau=c_*(n+2)^{-h-r-2},
\end{equation}
where $A$ is a sufficiently large fixed integer and $c_*>0$ a sufficiently small fixed rational. The factors $1+\log D$ in \cref{eq:uniform-displacement-size} do not make this choice circular: $K$ is chosen first, and then $A$ is fixed. Each $D$-dependent expression tends to zero relative to its error budget as $A$ increases, even with its powers of $\log A$. The remaining growth in $n$ is controlled by \cref{eq:choice-K}. In particular,
\[
n^h S D^{-1/4}=o(n^{-r}),\qquad
n^h S^2/D=o(n^{-r}),\qquad
m+S=o(D^{1/4}).
\]
Here the little-oh statements describe the uniform large-scale estimates used to choose the fixed constants; increasing $A$ covers every finite small-$n$ case. The same choices bound all the terms in \cref{eq:calibration-budget} and the Coulomb-compression error by their assigned fractions of \cref{eq:continuum-budget}. No search for a successful $D$ is part of the algorithm.
It is important to keep the energy units explicit. After dilation and subtraction of the isolated-well ground energy,
\begin{equation}\label{eq:units-final}
\mathcal H=H/Z_*^2+m/2,
\qquad Z_*\mathcal H=H/Z_*+mZ_*/2.
\end{equation}
\Needspace{10\baselineskip}
The exponential margins that are needed after multiplication by $Z_*\asymp e^D$ are:
\begin{center}
\begin{tabular}{@{}lll@{}}
\toprule
Source of error & In $\mathcal H$ units & In $Z_*\mathcal H$ units\\
\midrule
Tuned well energies & $m e^{-1.01D}$ & $O(m e^{-.01D})$\\
Discarded one-body entries & $O(m^2 e^{-1.1D})$ & $O(m^2 e^{-.1D})$\\
Orthonormalization & $\poly(m)e^{-1.6D}$ & $\poly(m)e^{-.6D}$\\
Continuum complement & $e^{-1.2D}$ & $O(e^{-.2D})$\\
\bottomrule
\end{tabular}
\end{center}
The polynomial multiplicities are retained in \cref{eq:continuum-simulation-final}. The positive exponential margins allow all these terms to fit their assigned error budgets by increasing the fixed $A$. The one-body residual $e^{-.75D}$ is used quadratically against a fixed complementary gap; multiplying that residual linearly by $Z_*$ would not be a valid estimate. \Cref{prop:continuum-simulation} supplies the required quadratic bound.
With these choices, that proposition gives
\begin{equation}\label{eq:full-final-error}
\left|E_0/Z_*+mZ_*/2-\alpha(E_{\rm in}-C_s)\right|
\le \frac{\alpha\gamma}{16}
\le \frac{\alpha\delta_s}{16}.
\end{equation}
\subsection{The deterministic construction and its bit cost}
First check the source encoding and the fixed polynomial bounds in \cref{thm:source}, including the inverse-gap bound. If a check fails, output the fixed instance with one electron, one unit-charge nucleus at the origin, and thresholds $-1/2$ and $1/2$. Otherwise carry out the following steps. The guard is polynomial-time and makes the reduction a total polynomial-time map; it never rejects a promised source instance.
\begin{enumerate}[label=\arabic*.,leftmargin=*]
\item \emph{Finite models.} Read the signed rational Heisenberg instance. Form the ancillary graph, the algebraic reference positions $p_i$, the complete contact set $\mathcal C$, the explicit scales $\Delta,\eta$, and the positive rational hoppings $t_{ij}$ of \cref{prop:hubbard}. The finite contact classification avoids deciding arbitrary geometric equalities numerically. All square roots here are approximated with certified rational intervals to the stated inverse-polynomial errors.
\item \emph{Charge scale.} Choose $D$ as in \cref{eq:choice-D} and compute a positive integer $Z_*$ with
\[
e^D/2\le Z_*\le2e^D.
\]
This does not require an exact comparison with a transcendental number. For example, compute the rational sum $s_D=\sum_{j=0}^{4D}D^j/j!$ and take $Z_*=\lfloor s_D\rfloor$. The ratio of successive omitted terms is at most $1/4$, and $j!\ge(j/e)^j$, so $0b.
\]
These inequalities concern the exact continuum spectral infimum by \cref{prop:continuum-simulation}, rather than only a projected eigenvalue.
The threshold separation is exactly
\begin{equation}\label{eq:output-gap}
b-a=\tfrac12 Z_*\alpha\delta_s.
\end{equation}
The reciprocal of $\alpha\delta_s$ is polynomially bounded in $n$, whereas $Z_*\ge e^D/2$. Increasing the same fixed scale constant if necessary makes \cref{eq:output-gap} at least one for every input. Since the output length $L\ge1$, this implies $b-a\ge L^{-1}$ without any circular estimate involving $L$.
We have constructed a legal nuclear instance in deterministic polynomial time, preserving YES and NO and permitting all spin sectors and all continuum states. This proves \cref{thm:main}.
==== END Continuum-Coulomb-hardness-with-binary-nuclear-charges//05-completion.tex ====
==== BEGIN Continuum-Coulomb-hardness-with-binary-nuclear-charges//paper.tex ====
\documentclass[11pt]{article}
\input{preamble}
\title{Continuum Coulomb hardness with binary nuclear charges}
\author{OpenAI}
\date{September 24, 2026}
\begin{document}
\maketitle
\begin{abstract}
We prove that approximating the electronic Coulomb spectral infimum in three-dimensional space is QMA-hard when positive integer nuclear charges are encoded in binary. The nuclei have distinct rational positions, the electron number is unary, and the energy is minimized over all antisymmetric continuum states and spin sectors. A deterministic classical polynomial-time reduction produces instances with threshold separation at least one. The nuclear charges may be exponentially large, but every output has polynomial bit length.
\end{abstract}
\input{sections/00-introduction}
\input{sections/01-finite-models}
\input{sections/02-wells}
\input{sections/03-hopping}
\input{sections/04-continuum}
\input{sections/05-completion}
\bibliographystyle{plain}
\bibliography{references}
\end{document}
==== END Continuum-Coulomb-hardness-with-binary-nuclear-charges//paper.tex ====
==== BEGIN Continuum-Coulomb-hardness-with-binary-nuclear-charges//preamble.tex ====
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage{lmodern}
\usepackage[margin=1.05in]{geometry}
\usepackage{amsmath,amssymb,amsthm,mathtools,bm}
\usepackage{graphicx,tikz,booktabs}
\usetikzlibrary{arrows.meta,calc,patterns,decorations.pathreplacing}
\usepackage{microtype}
\usepackage{enumitem}
\usepackage{xcolor}
\usepackage{etoolbox,needspace}
\usepackage[hyphens]{url}
\usepackage[colorlinks=true,linkcolor=blue!45!black,citecolor=blue!45!black,urlcolor=blue!45!black]{hyperref}
\usepackage[nameinlink,capitalise,noabbrev]{cleveref}
\hypersetup{pdftitle={Continuum Coulomb hardness with binary nuclear charges},pdfauthor={OpenAI},pdfsubject={A reduction for the exact all-spin electronic Coulomb spectral infimum},pdfkeywords={quantum complexity, Coulomb Hamiltonian, electronic structure, QMA-hardness}}
\newcommand{\R}{\mathbb R}
\newcommand{\N}{\mathbb N}
\newcommand{\Q}{\mathbb Q}
\newcommand{\C}{\mathbb C}
\newcommand{\dd}{\,\mathrm d}
\newcommand{\supp}{\operatorname{supp}}
\newcommand{\Tr}{\operatorname{Tr}}
\newcommand{\poly}{\operatorname{poly}}
\newcommand{\dist}{\operatorname{dist}}
\newcommand{\diag}{\operatorname{diag}}
\newcommand{\Spec}{\operatorname{Spec}}
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Corollary}
\theoremstyle{definition}
\newtheorem{definition}[theorem]{Definition}
\theoremstyle{remark}
\newtheorem{remark}[theorem]{Remark}
\crefname{theorem}{Theorem}{Theorems}
\crefname{lemma}{Lemma}{Lemmas}
\crefname{proposition}{Proposition}{Propositions}
\crefname{corollary}{Corollary}{Corollaries}
\crefname{definition}{Definition}{Definitions}
\crefname{remark}{Remark}{Remarks}
\crefname{section}{Section}{Sections}
\crefname{figure}{Figure}{Figures}
\setlist{topsep=5pt,itemsep=3pt,parsep=0pt}
\setlength{\emergencystretch}{1.5em}
\allowdisplaybreaks[1]
\BeforeBeginEnvironment{theorem}{\Needspace{4\baselineskip}}
\BeforeBeginEnvironment{lemma}{\Needspace{4\baselineskip}}
\BeforeBeginEnvironment{proposition}{\Needspace{4\baselineskip}}
\BeforeBeginEnvironment{corollary}{\Needspace{4\baselineskip}}
\BeforeBeginEnvironment{definition}{\Needspace{4\baselineskip}}
\BeforeBeginEnvironment{remark}{\Needspace{4\baselineskip}}
==== END Continuum-Coulomb-hardness-with-binary-nuclear-charges//preamble.tex ====
==== BEGIN Continuum-Coulomb-hardness-with-binary-nuclear-charges//references.bib ====
@misc{CMP2019,
author = {Toby Cubitt and Ashley Montanaro and Stephen Piddock},
title = {Universal quantum {Hamiltonians}},
year = {2019},
howpublished = {arXiv:1701.05182v4 [quant-ph]},
note = {Full manuscript, version 4, 31 May 2019.
\href{https://arxiv.org/pdf/1701.05182v4}{Online version}},
eprint = {1701.05182v4},
archivePrefix = {arXiv},
primaryClass = {quant-ph},
url = {https://arxiv.org/abs/1701.05182v4}
}
@article{CMP2018,
author = {Toby S. Cubitt and Ashley Montanaro and Stephen Piddock},
title = {Universal quantum {Hamiltonians}},
journal = {Proceedings of the National Academy of Sciences},
volume = {115},
number = {38},
pages = {9497--9502},
year = {2018},
doi = {10.1073/pnas.1804949115}
}
@article{PiddockMontanaro2017,
author = {Stephen Piddock and Ashley Montanaro},
title = {The complexity of antiferromagnetic interactions and {2D} lattices},
journal = {Quantum Information and Computation},
volume = {17},
number = {7--8},
pages = {636--672},
year = {2017},
eprint = {1506.04014v2},
archivePrefix = {arXiv},
url = {https://www.rintonpress.com/xxqic17/qic-17-78/0636-0672.pdf}
}
@article{LiuChristandlVerstraete2007,
author = {Yi-Kai Liu and Matthias Christandl and Frank Verstraete},
title = {Quantum computational complexity of the {$N$}-representability problem: {QMA} complete},
journal = {Physical Review Letters},
volume = {98},
pages = {110503},
year = {2007},
doi = {10.1103/PhysRevLett.98.110503}
}
@article{WhitfieldLoveAspuruGuzik2013,
author = {James Daniel Whitfield and Peter John Love and Al{\'a}n Aspuru-Guzik},
title = {Computational complexity in electronic structure},
journal = {Physical Chemistry Chemical Physics},
volume = {15},
pages = {397--411},
year = {2013},
doi = {10.1039/C2CP42695A}
}
@article{SchuchVerstraete2009,
author = {Norbert Schuch and Frank Verstraete},
title = {Computational complexity of interacting electrons and fundamental limitations of density functional theory},
journal = {Nature Physics},
volume = {5},
pages = {732--735},
year = {2009},
doi = {10.1038/nphys1370},
eprint = {0712.0483v2},
archivePrefix = {arXiv}
}
@article{Anderson1959,
author = {P. W. Anderson},
title = {New approach to the theory of superexchange interactions},
journal = {Physical Review},
volume = {115},
pages = {2--13},
year = {1959},
doi = {10.1103/PhysRev.115.2}
}
@article{BravyiDiVincenzoLoss2011,
author = {Sergey Bravyi and David P. DiVincenzo and Daniel Loss},
title = {{Schrieffer--Wolff} transformation for quantum many-body systems},
journal = {Annals of Physics},
volume = {326},
number = {10},
pages = {2793--2826},
year = {2011},
doi = {10.1016/j.aop.2011.06.004}
}
@article{MelenkBabuska1996,
author = {J. M. Melenk and I. Babu{\v{s}}ka},
title = {The partition of unity finite element method: Basic theory and applications},
journal = {Computer Methods in Applied Mechanics and Engineering},
volume = {139},
number = {1--4},
pages = {289--314},
year = {1996},
doi = {10.1016/S0045-7825(96)01087-0}
}
@article{MadayMarcati2019,
author = {Yvon Maday and Carlo Marcati},
title = {Regularity and {$hp$} discontinuous {Galerkin} finite element approximation of linear elliptic eigenvalue problems with singular potentials},
journal = {Mathematical Models and Methods in Applied Sciences},
volume = {29},
number = {8},
pages = {1585--1617},
year = {2019},
eprint = {1810.09010},
archivePrefix = {arXiv},
url = {https://arxiv.org/abs/1810.09010}
}
@article{Simon1983,
author = {Barry Simon},
title = {Semiclassical analysis of low lying eigenvalues. {I}. Non-degenerate minima: asymptotic expansions},
journal = {Annales de l'Institut Henri Poincar{\'e}, Section A, Physique Th{\'e}orique},
volume = {38},
number = {3},
pages = {295--308},
year = {1983},
url = {https://www.numdam.org/item/AIHPA_1983__38_3_295_0/}
}
@misc{UnitCoulomb2026,
author = {{OpenAI}},
title = {{QMA-hardness of continuum Coulomb energy with unit nuclear charges}},
howpublished = {OpenAI Math Release preprint
\href{https://github.com/openai/math/blob/main/preprints/QMA-hardness-of-continuum-Coulomb-energy-with-unit-nuclear-charges-September-24-2026/QMA-hardness-of-continuum-Coulomb-energy-with-unit-nuclear-charges-September-24-2026.pdf}{OAI:QMA-hardness-of-continuum-Coulomb-energy-with-unit-nuclear-charges-September-24-2026}},
year = {2026}
}
@misc{UniformExcessCharge2026,
author = {{OpenAI}},
title = {{Uniform excess charge for Coulomb molecules and the outer radius of neutral atoms}},
howpublished = {OpenAI Math Release preprint
\href{https://github.com/openai/math/blob/main/preprints/Uniform-excess-charge-for-Coulomb-molecules-and-the-outer-radius-of-neutral-atoms-September-24-2026/paper.pdf}{OAI:Uniform-excess-charge-for-Coulomb-molecules-and-the-outer-radius-of-neutral-atoms-September-24-2026}},
year = {2026}
}
@article{OGorman2022,
author = {Bryan O'Gorman and Sandy Irani and James Whitfield
and Bill Fefferman},
title = {Intractability of Electronic Structure in a Fixed Basis},
journal = {PRX Quantum},
volume = {3},
number = {2},
pages = {020322},
year = {2022},
doi = {10.1103/PRXQuantum.3.020322},
note = {\href{https://journals.aps.org/prxquantum/pdf/10.1103/PRXQuantum.3.020322}{Online version}},
url = {https://journals.aps.org/prxquantum/abstract/10.1103/PRXQuantum.3.020322}
}
@article{Schleich2026,
author = {Philipp Schleich and Lasse Bj{\o}rn Kristensen
and Jorge A. Campos-Gonzalez-Angulo
and Abdulrahman Aldossary and Davide Avagliano
and Mohsen Bagherimehrab and Christoph Gorgulla
and Joe Fitzsimons and Al{\'a}n Aspuru-Guzik},
title = {Chemically motivated simulation problems are efficiently
solvable on a quantum computer},
journal = {Digital Discovery},
volume = {5},
number = {1},
pages = {64--87},
year = {2026},
doi = {10.1039/D5DD00377F},
note = {\href{https://pubs.rsc.org/en/content/articlepdf/2025/dd/d5dd00377f}{Online version}},
url = {https://doi.org/10.1039/D5DD00377F}
}
@book{Teschl2009,
author = {Teschl, Gerald},
title = {Mathematical Methods in Quantum Mechanics: With Applications to {Schr\"odinger} Operators},
series = {Graduate Studies in Mathematics},
volume = {99},
publisher = {American Mathematical Society},
year = {2009},
note = {First edition, Theorem 6.4.
\href{https://www.mat.univie.ac.at/~gerald/ftp/book-schroe/schroe.pdf}{Author-hosted text}}
}
==== END Continuum-Coulomb-hardness-with-binary-nuclear-charges//references.bib ====