Conventions. Atomic units throughout. |A| is the volume of A. p(n,D) is some polynomial in the named variables with absolute constants; it can be enlarged from line to line. E(H) is the bottom of the spectrum of H.
A single electron near a proton has a lowest energy; any fixed nuclear configuration gives one real number whose exact computation is hard, and even its threshold decision is extraordinarily hard.
QMA is the class of decision problems whose YES instances admit a polynomial-size quantum proof. It contains NP; no efficient algorithm is known for QMA-hard problems.
Aharonov-Naveh 2002, arXiv:quant-ph/0210077. Kitaev-Shen-Vyalyi 2002, AMS Grad Studies 47.
Feynman (1982) and Lloyd (1996) framed physical simulation as a natural quantum-computational task; QMA-hardness is the case where even this remains intractable.
Feynman, R. P., "Simulating physics with computers", Int. J. Theor. Phys.
21 (1982), 467-488. https://doi.org/10.1007/BF02650179
Lloyd, S., "Universal quantum simulators", Science 273 (1996), 1073-1078.
https://doi.org/10.1126/science.273.5278.1073
The problem decided. Given N electrons interacting with M nuclei of charge
exactly one, positioned at pairwise distinct rational triples, and given
two rational thresholds a < b with b - a >= 1/L where L is the input
length, decide whether the electronic ground energy E_0 lies in
(-inf, a] or in [b, inf). The promise excludes energies strictly between
a and b. Deciding this is QMA-hard under deterministic classical
polynomial-time many-one reductions. The reduction produces polynomially
many nuclei and electrons, rational positions and thresholds whose bit
lengths are polynomial in the source length, and a threshold separation
b - a >= 1.
What makes the statement sharp. The nuclei are genuine point charges of
charge one, not an invented potential with wells and oscillators built
in. The minimization is over every continuum many-electron state of
either spin, not over a supplied finite basis. No magnetic field, extra
external potential, neutrality, binding, existence of a ground
eigenvector, or excitation-gap assumption is imposed. Membership in QMA
is not asserted.
Schuch and Verstraete (2009) proved QMA-completeness for the Hubbard model
with local magnetic fields and outlined a continuum realization using
designed scalar and spin-dependent potentials. O'Gorman, Irani, Whitfield,
and Fefferman (2022) proved QMA-completeness for electronic structure in a
supplied finite orbital basis, which fixes the Hilbert space in advance.
Cubitt, Montanaro, and Piddock (2019) proved QMA-hardness of the signed
exchange problem on finite square-lattice graphs, which provides the
starting point. Each relaxation is removed at once: the free continuum,
with real unit point nuclei, reached from the same signed-exchange source.
Schuch, N. and Verstraete, F., "Computational complexity of interacting
electrons and fundamental limitations of density functional theory",
Nature Physics 5 (2009), 732-735. arXiv:0712.0483.
local: 0712.0483.pdf
O'Gorman, B., Irani, S., Whitfield, J., and Fefferman, B., "Intractability
of electronic structure in a fixed basis", PRX Quantum 3 (2022), 020322.
doi:10.1103/PRXQuantum.3.020322.
Cubitt, T., Montanaro, A., and Piddock, S., "Universal quantum Hamiltonians",
pnas 115 (2018), 9497-9502. arXiv:1701.05182v4.
local: 1701.05182v4.pdf
Electrons in R^3. Each electron also has a spin, which is a direction in a
two-dimensional complex space C^2. So one electron lives in L^2(R^3) tensor
c^2. n electrons live in the N-fold antisymmetric tensor power
H_N = Lambda^N ( L^2(R^3) tensor C^2 ).
"Antisymmetric" means two electrons may not occupy the same (position,spin)
state: the wavefunction changes sign when you swap two electrons. Position
and spin are swapped together, not separately. No spin sector is fixed.
For M nuclei at positions R_1,...,R_M with charges Z_alpha (here all equal
to 1) the electronic Hamiltonian in atomic units is
H = -1/2 sum_{i=1}^N Delta_{x_i}
- sum_{i=1}^N sum_{alpha=1}^M 1/|x_i - R_alpha|
+ sum_{1<=i<j<=N} 1/|x_i - x_j| . (1)
The three pieces are: electron kinetic energy; electron-nucleus attraction;
electron-electron repulsion. The nuclear repulsion sum_{alpha<beta}
Z_alpha Z_beta/|R_alpha-R_beta| is a fixed number once the nuclei are fixed,
so it is omitted. The energy in (1) is the electronic energy. This is
the convention used everywhere below.
The natural space of states is the H^1 form domain
Q_N = H^1(R^{3N};(C^2)^{tensor N}) intersect H_N .
H is unbounded below on all of L^2 if this is careless, but on Q_N it is
bounded below and has a well-defined self-adjoint realization obtained from
the closed Coulomb quadratic form. Write
E_0 = inf spec H = inf{ <psi,H psi> : psi in Q_N, ||psi|| = 1 },
whether or not the infimum is attained by an actual eigenvector. This
"form" viewpoint is exactly what is needed later: energies are compared
by testing the quadratic form on arbitrary trial states, never assuming
that a ground eigenvector exists.
Teschl, G., "Mathematical Methods in Quantum Mechanics with Applications
to Schrodinger Operators", Graduate Studies in Mathematics 99, AMS, 2009.
https://www.mat.univie.ac.at/~gerald/ftp/book-schroe/
Reed, M. and Simon, B., "Methods of Modern Mathematical Physics, Vol. II:
Fourier Analysis, Self-Adjointness", Academic Press, 1975.
ISBN: 978-0125850025.
A "promise problem" is a decision question that is only asked on inputs
obeying some rule (the promise). The inputs are described exactly here.
Input
M >= 1 nuclear positions R_1,...,R_M in Q^3, pairwise distinct.
All charges Z_alpha = 1.
N >= 1 electrons, written in unary (N ones, so a large N costs length N).
Two rational numbers a < b, the thresholds.
encoding
Each rational is a signed binary numerator over a positive binary
denominator, in lowest terms, inside a fixed self-delimiting encoding.
length L
Total bit length of all the above (the unary count contributes N).
Promise
b - a >= 1/L.
Answer
YES if E_0 <= a. NO if E_0 >= b.
(The promise b-a >= 1/L makes the two cases disjoint.)
The promise b - a >= 1/L concerns decision precision, not a molecular spectral gap; the theorem assumes no excitation gap, binding, neutrality, or existence of a ground eigenvector.
Energy convention: electronic energy; nuclear repulsion omitted.
Main theorem. The promise problem is QMA-hard under deterministic classical polynomial-time many-one reductions, with output consisting of polynomially many unit-charge rational nuclei at distinct positions, polynomially many unary electrons, and polynomial-length rational thresholds at separation b - a >= 1. Thus Z_max = 1 and c = 1; membership in QMA is not asserted; and neither b - a >= 1 nor the input gap 1/L constrains the molecular excitation gap. Each stage preserves the promise because every energy error is inverse-polynomially small relative to gamma, and gamma itself is inverse-polynomial in source length.
The proof starts from a QMA-hard problem about spins.
The spin source (Cubitt-Montanaro-Piddock, 2018/2019).
Consider a finite square-lattice graph (vertices are cells of Z^2, edges
join cells at distance 1), possibly with some signed real couplings J_e.
The Hamiltonian is
H_in = sum_e J_e h^s_{LR}, h^s_{LR} = sigma_L . sigma_R, (2)
where sigma = (sigma_x,sigma_y,sigma_z) are the Pauli matrices acting on
a spin-1/2 at each site, and sigma_L . sigma_R means
sigma_x_L sigma_x_R + sigma_y_L sigma_y_R + sigma_z_L sigma_z_R. There
are no one-site terms.
The cited result supplies three facts, not re-proved here: deciding the ground-energy threshold of (2) is QMA-hard with inverse-polynomial promise gap; the lattice has polynomially many vertices and edges and polynomially bounded couplings; and the couplings are real with polynomially many nonzero entries. The incoming instance is only reshaped to a clean rational format.
Identity of h^s. On two spins, h^s = sigma . sigma = 2 SWAP - I, where
swap exchanges the two spins. Its eigenvalues are
+1 on the three triplet states (total spin 1),
-3 on the singlet state (total spin 0).
So h^s favors singlets when positive, triplets when negative.
Normalizing the source to rationals.
Step 1: a rational lower bound on the gap. Because the source gap is
inverse polynomial in the input size, there is a rational delta > 0, known
in advance and inverse polynomial, with delta <= (source gap).
Step 2: rational couplings and thresholds. Let q be the number of
weighted edges. Set eps_0 = delta/100. Round each coupling J_e to a
rational within eps_0/max(1,q). The total change in Hamiltonian norm
(supremum over unit vectors of |change|) is at most q * (3 * eps_0/max(1,q))
= 3 eps_0, since |h^s_{LR}| has operator norm 3. Write the incoming
instance's scalar as c_0 I (c_0 = 0 if there is none) and its thresholds
as a_0 < b_0. Round c_0, a_0, b_0 within eps_0 to rationals c~, a~, b~.
The rounded Hamiltonian H~_in = sum J~_e h^s_{LR} obeys
|H~_in - (H_in - c_0 I)| <= 3 eps_0 . (3)
Define rational thresholds
a_s = a~ - c~ + 5 eps_0 , b_s = b~ - c~ - 5 eps_0 . (4)
Then YES/NO of the source stays intact: the 5-eps_0 shifts absorb the
3 eps_0 of Hamiltonian rounding plus eps_0 each for the scalar and the
threshold. And
b_s - a_s >= (b_0 - a_0) - 2 eps_0 - 10 eps_0 >= delta - 12 eps_0
= 0.88 delta > 0 . (5)
Step 3: geometric extent. A connected component of the square grid with V
vertices has coordinate diameter at most V-1 (walk along lattice edges).
Translating components into disjoint horizontal strips leaves all weighted
graphs unchanged and keeps the total coordinate range polynomial. So the
normalized instance lives in a polynomial-size bounding box.
Finally set
n = max(2, source input size), gamma = min(1, b_s - a_s). (6)
From now on, "the source instance" means the rational signed Heisenberg
Hamiltonian H_in with thresholds a_s < b_s and gap gamma. One may assume
there is at least one vertex (add a free spin if the source is empty) and
one may drop couplings of zero. All magnitudes and inverse gaps
encountered later are polynomially bounded in n.
The source couplings J_e can be positive or negative. The orbital
construction naturally produces positive exchange (the electron kinetic
energy prefers singlets). So first convert the signed
source into a positive-exchange spin model. The tool is a "singlet-pair
mediator": a fresh pair of spins forced into a singlet, with weak links to
the original spins. Its effective coupling to the two ends has the sign one chooses.
This is the mechanism of Piddock and Montanaro Lemma 5. The needed rule: eliminate excited states and keep second order. The quantitative statement, with fewer hypotheses than the usual Schrieffer-Wolff theorem and no requirement that the first-order block pbp vanish, follows.
Schrieffer, J. R. and Wolff, P. A., "Relation between the Anderson and
Kondo Hamiltonians", Phys. Rev. 149 (1966), 491-492.
https://doi.org/10.1103/PhysRev.149.491
Bravyi, S., DiVincenzo, D. P., and Loss, D., "Schrieffer-Wolff
transformation for many-body systems", Ann. Physics 326 (2011),
2793-2823. arXiv:1105.0675.
local: 1105.0675.pdf
Lemma 1 (second-order bottom energy).
Let A and B be self-adjoint operators on a finite-dimensional Hilbert
space. Suppose A has kernel projection P != 0, satisfies A >= g Q where
Q = I - P and g > 0, and B has norm ||B|| <= g/4. Then
E(A+B) = E( (P B P - P B A^{-1} B P) restricted to P ) + O(||B||^3/g^2),
where A^{-1} on the right means the inverse of A on the range of Q, with
value zero on the kernel of A.
Proof. First, |E(A+B)| <= ||B||: the bottom of A+B is at most <x,(A+B)x> at
any unit x, and at least -||B|| since A >= 0 everywhere on the space.
Let x be a nearly-ground unit vector for A+B with energy E = E(A+B). Write
x = p + q with p = Px and q = Qx. The eigenvalue equation (A+B)x = Ex,
projected onto the range of Q, reads
(A|_Q) q + QB(p+q) = E q .
So q = -( A|_Q + QBQ|_Q - E )^{-1} QBp; indeed p != 0, since p = 0 would
force energy >= g - ||B|| > -||B||, contradicting |E| <= ||B||. Because
|E| <= ||B|| and ||QBQ|| <= ||B|| <= g/4,
the inverse (A|_Q + QBQ|_Q - E)^{-1} differs from (A|_Q)^{-1} by O(||B||/g^2)
in operator norm, by the resolvent identity:
(A + C)^{-1} - A^{-1} = -A^{-1} C (A+C)^{-1} .
Substituting back gives, up to O(||B||^2/g) on q and O(||B||^3/g^2) on the
energy, q ~ -A^{-1} QBp. Project the eigenvalue equation onto P:
E p = P B p + P B q = P B p - P B A^{-1} B p + O(||B||^3/g^2) .
Taking <p, .> and dividing by ||p||^2 (which is 1 - ||q||^2 = 1 - O(||B||^2/g^2))
gives the claimed lower estimate.
For the upper estimate, choose a unit ground vector p of the compressed
operator W = (PBP - PBA^{-1}BP)|_P, and try x = p - A^{-1} QB p. Then
<x,(A+B)x> equals <p,Wp> plus O(||B||^3/g^2), and ||x||^2 = 1 + O(||B||^2/g^2);
dividing gives the same value, so E(A+B) <= E(W) + O(||B||^3/g^2). Together
with the lower estimate this proves the lemma. (The case Q = 0 is
immediate.) QED
For each surviving signed source edge e = LR of weight J != 0, add a
fresh pair of spins A_e, B_e, coupled as follows:
central bond: Delta * h^s_{A_e B_e} (Delta > 0 integer)
two "spokes": sqrt(Delta) * v_e * ( h^s_{L A_e} + h^s_{C_e R} )
where
v_e = sqrt(2 |J|),
C_e = B_e if J > 0, C_e = A_e if J < 0.
So a positive edge attaches L--A_e and B_e--R; a negative edge attaches
L--A_e and A_e--R (both spokes touch A_e, so the sign of the effective
LR coupling flips). The strength of the spokes is sqrt(Delta) v_e, which is
"small" compared to the central gap 4 Delta, so Lemma 1 applies.
The gadget as an effective LR coupling.
Step 1: make the fresh pair A_e B_e sit in its singlet. The operator
A_e, B_e's coupling Delta h^s prefers the singlet (h^s has eigenvalue -3
there, versus +1 on triplets). Add a fixed penalty
A = Delta * sum_e ( h^s_{A_e B_e} + 3 I )
so the singlet is the unique zero of A on the fresh pairs and
A >= 4 Delta on the Q space (each pair's triplet energy 1 minus (-3)
gives 4). A >= g Q with g = 4 Delta.
Step 2: the rest is B = H_in(rest) + spoke terms. Its norm is at most
poly(n)(1 + sqrt Delta) since the number of edges and the weights are
polynomial and sqrt(Delta) is the only new strong scale. If Delta is a
sufficiently large polynomial in n, the bound ||B|| <= g/4 = Delta holds.
Step 3: apply Lemma 1. The kernel P of A is spanned by product states in
which every fresh pair is a singlet; on that space the original signed
terms and the source-unaffected part act unchanged. The kernel restriction
of the first-order term P B P keeps the unreplaced part. The second-order
correction is -P B A^{-1} B P. Only two "hops" through B can act on the
same fresh pair, and they have energy denominator g = 4 Delta (the gap of
A on that pair): so A^{-1} on a double action on pair e is 1/(4 Delta)
times the pair identity.
Now use the Pauli algebra on a singlet. If |s> is the singlet on A_e B_e
and a Pauli operator sigma_mu sits on A_e only while sigma_nu sits on B_e
only, then
<s| sigma_mu^(A) sigma_nu^(B) |s> = - delta_{mu nu} (sum over mu),
<s| sigma_mu^(A) sigma_nu^(A) |s> = + delta_{mu nu} (sum over mu).
(I.e. the singlet correlation matrix is +delta for two Paulis on the same
member and -delta for one on each member.) Tracing the second-order term
over the singlet and over the intermediates gives, per edge e,
J * h^s_{LR} - (3/2) v_e^2 I = J h^s_{LR} - 3|J| I , (7)
because v_e^2 = 2|J|, hence (3/2) v_e^2 = 3|J|. The sign of the effective
LR term is the sign of J by construction (positive edge: L--A_e, B_e--R
distinct members of the singlet; negative edge: L--A_e, A_e--R, both
endpoints connect to the same member A_e, flipping the correlation via
the singlet relation above).
The correction to the ground energy from Lemma 1 is O(||B||^3/g^2) =
O( poly(n)(1 + sqrt Delta)^3 / Delta^2 ), which is O(poly(n)/sqrt(Delta));
choose Delta a large enough polynomial in n that this is <= gamma/100.
Step 4: all edges simultaneously. Because the operations commute across
distinct fresh pairs (each pair touches only its own two spokes), one can
replace every surviving signed edge at the same time, with the same large
Delta. The total shift is sum_e (3 Delta + (3/2) v_e^2) - (correction).
The singlet-pair gadget produces a positive spinner vertex (the fresh pair)
and two spokes of positive weight sqrt(Delta) v_e, one joining L to A_e and
one joining C_e to R. But the link strengths K_e must themselves be encoded in the distances between
orbitals. Having an intermediate A_e or C_e in the middle of a "spoke" is
fine, but a spoke of positive weight is still an h^s coupling, and the goal is
to realize all positive h^s couplings as orbital hoppings.
Subdivide each positive spoke into a path of length three using the same
gadget with positive sign (a fresh pair forced into a singlet, attached at
its two ends). Each subdivider replaces one h^s coupling K h^s_{XY} by a
chain X -- fresh pair -- Y, with two couplings sqrt(Delta') K h^s on each
side and a central bond Delta' h^s. Lemma 1 leaves the effective X-Y
coupling equal to (K/2) h^s up to a scalar shift; because one may rescale
the whole target (and absorb the 1/2 into K_e), the only care is that the
sign is preserved and the strength is positive.
Two rounds of subdivision (three per spoke, nine per original edge, because
each round triples the path) produce per original edge:
- one central bond A_e--B_e (positive), and
- two paths of length nine, joining L to A_e and C_e to R.
Each round uses its own polynomial Delta' (there are only three rounds
total, so there is no need to fit infinitely many Deltas), and each round
keeps errors O(poly(n)/sqrt(Delta')) = O(gamma/100). Summing the errors of
three rounds gives O(3 gamma/100) < gamma/25.
The final positive-exchange Hamiltonian is
H_+ = sum_{links} K_e h^s_{ij} , (8)
with every K_e positive, every K_e and 1/K_e polynomially bounded in n,
and the number of sites m (including the original vertices and every
fresh spacer vertex) polynomial. There is a polynomial-time computable
scalar shift C_s (sum of all the scalar shifts from the singlet-pair
gadgets and the subdivisions and their second-order corrections) with
| E(H_+) - E(H_in) + C_s | <= 4 gamma/100 . (9)
Every intermediate weight K_e and its shift are computable (as real
numbers) to inverse-polynomial absolute accuracy in polynomial time.
Square roots need not be exact rationals at this stage; they will be
absorbed into a final rational approximation.
The positive-exchange Heisenberg Hamiltonian H_+ = sum K_e h^s_{ij} has polynomial size m, every K_e and 1/K_e polynomially bounded, and a scalar shift C_s satisfying (9).
A Hamiltonian H_+ = sum K_e h^s_{ij} is just a weighted graph. In the
well construction each vertex is realized as a localized spatial orbital
and each link as
a distance between orbitals: the hopping strength t_ij drops off like
exp(-distance). So the link lengths must encode the K_e: each link has a
different length, chosen so its hopping equals the desired value. A
planar drawing of the modified graph is built in which every link
length can be adjusted independently, with nonlink distances at least 1.4.
Every original edge becomes a small "contact gadget", building on the
cmp (Cubitt-Montanaro-Piddock) lattice.
Place original lattice vertices at grid cells with spacing 17. For each
original edge e = LR, choose a local frame whose first coordinate runs
along the edge from L to R. In this frame:
positive edge (J > 0): center pair A_e = (8,0), B_e = (9,0).
negative edge (J < 0): center pair A_e = (8.5,0), B_e = (8.5,1).
Picture (topology only, not to scale):
J > 0: L ~~~ A_e === B_e ~~~ R
J < 0: L ~~~ A_e ~~~ R
|
B_e
~~~ = one bent nine-link path (length 8 or 8.5; drawn zig-zag so that
off-path sites stay >= 1.4 apart)
=== = the central unit-length bond A_e--B_e
L, R = the two original lattice vertices, 17 apart
A_e, B_e = the fresh singlet-pair sites from the positive-exchange conversion
The central bond A_e--B_e has length 1 in each case. The two spokes are
nine-link paths of length nine each:
positive edge: L to A_e and B_e to R;
negative edge: L to A_e and A_e to R (both touch A_e).
Each nine-link path joins two points 8 (positive) or 8.5 (negative) apart
on the axis. Lay it out along the axis with steps s_1,...,s_9 whose
projections on the axis sum to that span and whose transverse components
cancel: first two, middle (fifth), and last two steps are unit forward
steps (transverse 0). Steps 3, 4 slope forward with first component 3/4
(positive) or 7/8 (negative) and one sign of transverse offset; steps 6, 7
mirror those offsets with the opposite sign. This makes both paths lie
within a narrow strip of transverse offset < 1.324.
Nonlink separation. The distance between any two sites that are not
consecutive sites on a link is at least 1.4. Consecutive along-axis
separations are at least 3/4; nonconsecutive sites on the same path are at
least 1.5 apart; the extra negative-sign side site B_e is at least sqrt(2)
from everything except its central neighbor; distinct parallel lines never
interfere because the off-axis transverse offsets are < 1.324 and the
spacing is 17. Original vertices other than a gadget's two endpoints sit
far from all internal sites, by the 17-spacing. This verifies at once
that (with the size-one distances written above) nonlinks are at least
1.4 apart.
Now show that the length of each of the 18 links per gadget can be set
to any value in [1 - eps, 1 + eps] for some fixed small eps > 0, all at
once, without ever bringing nonlinks closer than 1.2. The original lattice
vertices stay fixed.
Fix one gadget. It has internal vertices (A_e, B_e and the path vertices)
and one central link plus the path links. To change the central A_e-B_e
distance, move B_e directly (a one-dimensional freedom). To change each of
the twelve path-link lengths, move the path's 8 internal vertices. This
gives for each path a map
d : (internal vertex displacements) -> (twelve length increments)
One must show d is surjective near zero. Look at the linearization. Let
e_1,...,e_9 be the unit step vectors of the path. A linear dependence with
coefficients a_1,...,a_9 among the rows of the fixed-end length
differential would give, at each internal vertex j (1 <= j <= 8),
a_j e_j - a_{j+1} e_{j+1} = 0 .
This forces all vectors a_j e_j to be equal. Since the path has some pair
of nonparallel steps (steps 3,4,6,7 have transverse components), all
coefficients must vanish. So the rows of the differential are linearly
independent, thus the whole gadget's link-length differential is
surjective on the internal vertex variables.
Change variables. Choose a fixed right inverse M for the (surjective)
differential restricted to a fixed small ball. In these new coordinates
the length increment map is
q |-> q + R(q) , with R(0) = 0, DR(0) = 0 .
On a sufficiently small ball DR(q) has ||DR(q)|| < 1/2. So for any small
target increment d (with ||d|| small), the fixed-point iteration
q_{k+1} = d - R(q_k)
is a contraction on that ball and converges to a solution of the prescribed
increments. Choose the ball small enough that nonlink separations stay
above 1.2, and then fix eps accordingly.
Because original lattice vertices stay fixed, each gadget can be solved
independently. Their number (polynomial in n) does not shrink the allowed
length interval [1 - eps, 1 + eps]: gadgets do not interfere, so one can
implement all eighteen independent adjustments at once.
Computability. Each such solve is a fixed-dimensional contraction with
algebraic initial data (rational or nested-square-root). Iterating with
guard precision computes every internal position to any inverse-polynomial
accuracy in polynomial time.
The unit-distance layout above is then rescaled by D = log(lambda), giving minimum spacing at least 0.95 D, nonlink spacing at least 1.2 D, and spatial extent poly(n) * O(D).
The orbit centers u_i are placed at the rescaled layout positions.
Every link is independently adjustable in [D(1 - eps), D(1 + eps)], and nonlinks remain farther apart than 1.2 D.
The sign conversion produced the target H_+ = sum K_e h^s_{ij} on m sites; the planar layout gives a placement with independent link lengths. Each site hosts one smooth localized spatial orbital, and the target hoppings are calibrated by choosing separations.
Because the eventual nuclei have charge 1, the physical scale is set by a
large number lambda. Work in dilated coordinates
y = lambda x = (r, z) in R^3, r in R^2 (plane), z in R (vertical).
After dividing energies by lambda^2:
- the kinetic operator becomes the usual -Delta/2;
- each unit nucleus contributes an attraction of coupling 1/lambda;
- electron repulsion has coefficient 1/lambda.
This "one lambda of dilution" is why the final electronic charges come out
as exactly one.
The planar mode. Take
f(r) = (1 - |r|^2)_+^16 , (10)
which is a smooth nonnegative bump supported in the unit disk, positive
inside, and vanishing to order 16 at the boundary. Define
phi = (-Delta_r + 1)^{-1} f , so -Delta_r phi + phi = f,
varphi = phi / ||phi||_2 ,
W(r) = - f(r) / (2 phi(r)) . (11)
Then phi > 0, radial, exponentially decaying, and W is nonpositive and
compactly supported (inside the unit disk). Check the eigen-relation:
-Delta phi + phi = f => -Delta phi / 2 = (f - phi)/2 ,
-Delta phi / 2 + W phi = (f - phi)/2 - f/2 = -phi/2 .
So (-Delta_r/2 + W) phi = -phi/2, i.e. varphi is the ground mode with energy
exactly -1/2. Derivatives of phi through order six are bounded. The
representation
phi(r) = int_0^\infty e^{-t} int_{R^2} (4 pi t)^{-1} e^{-|r-b|^2/(4t)}
f(b) db dt
(from the heat kernel; this is the standard Gaussian-heat formula for the
resolvent of -Delta + 1) gives, for |r| > 2, the tail bounds
C^{-1}(1+|r|)^{-C} e^{-|r|} <= |D^k phi(r)| <= C(1+|r|)^C e^{-|r|}
for the derivatives k one needs, using t + (|r|-1)^2/(4t) >= |r| - 1. The
lower bound comes from integrating the density over a small half-disk and a
band |r|/2 < t < |r|/2 + 1. In particular phi has a fixed positive lower
bound on the unit disk, so W and all its derivatives are bounded there, and
W is zero outside.
Isolated gap. The planar operator -Delta_r/2 + W has the unique ground state
varphi and a fixed positive gap above it. Proof: the shifted quadratic form
on a test function u is int phi^2 |nabla(u/phi)|^2 / 2 (nonnegativity by
closure and local positivity, and the kernel is exactly the multiples of
phi). If a unit sequence perpendicular to varphi had shifted form tending
to zero, it would be bounded in H^1 and have a weakly and locally L^2
convergent subsequence; by compact support of W the limit would have to lie
in the kernel (multiples of phi), so the limit would be perpendicular to and
parallel to varphi at once, hence zero; but then the 1/2 shift precludes the
form values from tending to zero. So a fixed gap c_pl > 0 separates the
ground mode varphi from the rest.
The vertical mode. For the z direction use the normalized Gaussian
g(z) = (omega/pi)^{1/4} e^{- omega z^2 / 2}, omega = sqrt(4 pi rho),
the ground mode of -d^2/dz^2/2 + 2 pi rho z^2 with energy omega/2. (rho > 0
is a large fixed integer chosen in the slab-and-wells density construction.) The oscillator has a fixed
gap omega above its ground mode.
The three-dimensional mode. Per site i, centered at planar position u_i
(from the contact-gadget layout, at scale D):
psi_i(y) = varphi(r - u_i) g(z), e_* = -1/2 + omega/2, (12)
the energy of the (planar ground) tensor (vertical ground). psi_i is
smooth, decaying exponentially at infinity with rate controlled by
phi and g.
Define the direct electron-electron coefficients
v_{ij} = iint psi_i(y)^2 psi_j(y')^2 / |y - y'| dy dy'
= v(|u_i - u_j|) , (13)
U = v(0) .
(U is the same-site coefficient; v_{ij} is the coefficient between sites i,j.)
They are uniformly bounded (the psi_i are uniformly bounded and decay), and
v is uniformly Lipschitz (translate the decaying density and differentiate).
Coercivity (a fixed lower bound on the whole matrix v). Now prove
(v_{ij}) >= c I, c > 0, uniformly over all separated geometries, once the
polynomial scale lambda is sufficiently large. For real coefficients a_i,
write F = sum_i a_i psi_i^2 (a smooth decaying density). Fix a fixed smooth
nonzero nonnegative test bump chi_0 supported near the origin, and let chi_i
be its disjoint translates to the centers (u_i, 0). Put
w = sum_i a_i chi_i, and Phi_F = the Coulomb potential of F, i.e. the function
with -Delta Phi_F = F.
Integration by parts + Cauchy-Schwarz:
int F w = int (-Delta Phi_F) w = int (nabla Phi_F) . (nabla w)
<= ||nabla Phi_F||_2 ||nabla w||_2 .
Now ||nabla Phi_F||_2^2 = int Phi_F (-Delta Phi_F) = int Phi_F F
= (1/(4 pi)) iint F(y) F(y') / |y - y'| dy dy' (this is exactly the Coulomb
double integral of F; it is finite for the smooth decaying F). Therefore
|int F w|^2 <= (1/(4 pi)) ( Coulomb integral of F ) * int |nabla w|^2 . (14)
Lower bound on the left. Define the matrix A_ij = int psi_i^2 chi_j.
Its diagonal entries are a fixed a_* > 0 (same localized overlap at each
site). Off-diagonal entries are exponentially small in D: each off-diagonal
requires psi_i^2 to reach chi_j at distance >= 0.95 D, costing exp(-c D).
Row and column sums of the off-diagonal entries are at most poly(m) e^{-c D};
choose lambda such that these are at most a_*/4 each. Then A = a_* I + R
with ||R|| <= a_*/2 (operator norm), and
int F w = sum_{i,j} a_i a_j A_ij = a^T A a >= (a_*/2) sum_i a_i^2 . (15)
Upper bound on ||nabla w||. The chi_i have disjoint supports and fixed
shapes, so
int |nabla w|^2 = sum_i a_i^2 int |nabla chi_0|^2 = O( sum_i a_i^2 ) . (16)
Combine (14), (15), and (16): the Coulomb double integral of F is at
least a fixed c times sum_i a_i^2. But F = sum_i a_i psi_i^2, so the Coulomb
integral is exactly sum_{i,j} a_i a_j v_{ij} = a^T V a with V = (v_{ij}).
Hence
a^T V a >= c sum_i a_i^2 for all real a , i.e. V >= c I . (17)
A direct consequence. Apply (17) to two sites with coefficients (1, -1):
(1,-1)^T V (1,-1) = v_{ii} + v_{jj} - 2 v_{ij} = 2U - 2 v(d)
which is >= 2c. Therefore
U - v(d) >= c > 0 for every distance d >= 0.95 D in use , (18)
and it is uniformly bounded above. So the superexchange denominator
U - v_{ij} below is uniformly positive and bounded.
Superexchange: two electrons on neighbouring sites momentarily pile onto one site, then split back; Pauli blocking restricts this to singlets, lowering their energy by roughly t^2/(U - v).
The fermion model. Consider 2m canonical fermion modes: at each site i, two
spin states sigma in {up,down}. Restrict to exactly m electrons (the
filling). Let n_i = sum_sigma c_{i sigma}^dagger c_{i sigma} be the total
occupancy at site i, so 0 <= n_i <= 2 and sum_i n_i = m. Multiply dilated
energies by lambda, so the Coulomb coefficients below carry no 1/lambda.
The direct Coulomb interaction, written in the mode basis, is
(U/2) sum_i n_i (n_i - 1) + sum_{i<j} v_{ij} n_i n_j . (19)
It is nonnegative and vanishes exactly when every site is singly occupied
(n_i = 1): then n_i(n_i-1) = 0 and the pairing term is zero at m = sum n_i.
Away from single occupancy it is positive.
One-site defect counterm. Later a diagonal one-body term
-sum_i s_i n_i will be needed, with s_i = sum_{j != i} v_{ij}. Using v_{ii} = U and
sum_i n_i = m, a short algebra shows
(U/2) sum_i n_i(n_i-1) + sum_{i<j} v_{ij} n_i n_j - sum_i s_i n_i
= (1/2) sum_{i,j} v_{ij} (n_i - 1)(n_j - 1) - sum_{i<j} v_{ij} . (C)
Proof of (C). Use symmetry v_{ij} = v_{ji} and s_i = sum_{j != i} v_{ij},
so the full row sum is sum_j v_{ij} = U + s_i. Write
M1 = sum_{i,j} v_{ij} n_i , M0 = sum_{i,j} v_{ij} .
By symmetry sum_{i,j} v_{ij} n_j = M1, and
M1 = sum_i n_i (U + s_i) = U m + sum_i s_i n_i (used sum_i n_i = m),
M0 = sum_i (U + s_i) = U m + 2 sum_{i<j} v_{ij} (used sum_i s_i = 2 sum_{i<j} v_{ij}).
Also, by symmetry again,
sum_{i,j} v_{ij} n_i n_j = U sum_i n_i^2 + 2 sum_{i<j} v_{ij} n_i n_j .
Now expand the right-hand side of (C):
(1/2) sum_{i,j} v_{ij}(n_i-1)(n_j-1) - sum_{i<j} v_{ij}
= (1/2)[ sum_{i,j} v_{ij} n_i n_j - 2 M1 + M0 ] - sum_{i<j} v_{ij}
= (1/2) U sum_i n_i^2 + sum_{i<j} v_{ij} n_i n_j
- U m - sum_i s_i n_i + (1/2) U m + sum_{i<j} v_{ij}
- sum_{i<j} v_{ij}
= (1/2) U sum_i n_i^2 + sum_{i<j} v_{ij} n_i n_j
- (1/2) U m - sum_i s_i n_i .
Since (U/2) sum_i n_i(n_i-1) = (U/2) sum_i n_i^2 - (U/2) m, this last line
is exactly the left-hand side of (C). QED
So define s_i = sum_{j != i} v_{ij}; the required diagonal one-body term
-sum_i s_i n_i will be produced by slightly deepening the wells.
The finite hamiltonian. The finite fermion model to be realized is
H_d(t) = (1/2) sum_{i,j} v_{ij} (n_i - 1)(n_j - 1) - sum_{i<j} v_{ij}
- sum_{ e={i,j} } t_e sum_{sigma} ( c_{i sigma}^dagger c_{j sigma}
+ c_{j sigma}^dagger c_{i sigma} ) . (20)
The quadratic penalty (1/2) sum v_{ij}(n_i-1)(n_j-1) is nonnegative and
its kernel is exactly the spin space with every site singly occupied (n_i =
1 all i). Its gap: a nonzero integer defect vector (n_i - 1)_i has sum zero
and hence squared Euclidean length at least 2; by coercivity (17) with
integer entries the penalty is at least c. So it has a fixed positive gap.
First order of hopping vanishes. On the singly-occupied spin space, a
single hop creates a doubly-occupied site and an empty site. That excitation
is penalized by the quadratic term, so it drops out of the first-order
compression onto the singly-occupied space.
Second order = superexchange. A single allowed hop on edge e = {i,j}
creates a doubly occupied site and an empty site. Its defect vector is
(1,-1) on that pair and zero elsewhere, so its excitation energy is exactly
U - v_{ij} (21)
(evaluated with the quadratic penalty plus scalar; derivation via (17)).
Only a hop on the same link can restore single occupancy in one more step.
For hopping strength t_e, the second-order (bottom-energy) operator is
( t_e^2 / (U - v_{ij}) ) ( h^s_{ij} - I ) . (22)
To verify the factor. On a singlet, the amplitude for one hop is set by
the fermion algebra; the doubly occupied intermediate has amplitude sqrt(2)
times t_e (two orthogonal doubly occupied states for the singlet), while
the triplet does not hop into a same-site pair. Squaring the amplitude and
dividing by the excitation energy U - v_{ij} gives the h^s-dependent term.
Physically this is Anderson's superexchange. Concretely, apply the
creation/annihilation operators to
(c_{i up}^dagger c_{j down}^dagger - c_{i down}^dagger c_{j up}^dagger)
|0> / sqrt(2) (the singlet on i,j)
and use that a second hop on the same link returns a spin state. Other
sites are placed after i,j in the ordering, which changes only the overall
sign on the spin space. The result matches the earlier statement.
Target hoppings. Introduce a common positive parameter tau (chosen below)
and set
t_e^tar = tau sqrt( K_e (U - v_{ij}) ) . (23)
Substituting into (22) and using (21) gives the effective spin term
tau^2 ( H_+ - sum_e K_e ) (24)
on the singly-occupied space. The hopping operator has norm at most
tau * poly(n) uniformly over admissible separated geometries. By Lemma 1
the energy remainder of the second-order expansion is at most
C tau^3 poly(n) (since the penalty gap is at least c, uniformly). Choose a
rational inverse-polynomial tau > 0 small enough that Lemma 1 applies
(tau poly(n) <= gap/4) and this remainder is at most gamma tau^2/100. This
choice of tau is independent of lambda and of the final centers.
Note: the occupancy penalty involves every pair {i,j} (through v_{ij}), but
the target hopping on a link e depends only on its own distance through
v_{ij} = v(d). So the calibration of one link is a one-dimensional problem.
The modes supply the hopping profile
T(d) = - int_{R^2} varphi(r) W(r) varphi(r - d e_1) dr
= ( 1 / (2 ||phi||_2^2) ) int f(r) phi(r - d e_1) dr . (25)
The second form follows from W = -f/(2 phi):
-varphi(r) W(r) = (phi(r)/||phi||) * f(r)/(2 phi(r)) = f(r)/(2 ||phi||).
Multiply by varphi(r - de_1) = phi(r - de_1)/||phi|| and integrate:
T(d) = (1/(2 ||phi||^2)) int f(r) phi(r - d e_1) dr .
This is manifestly positive (f, phi > 0), uniformly Lipschitz, and by the
tail bound exp(-|r|) of phi it satisfies, for some absolute C,
c (1+d)^{-C} e^{-d} <= T(d) <= C (1+d)^{C} e^{-d} . (26)
The hopping to be supplied by the continuum. The cubature replacement shows that the real
continuum orbital matrix, evaluated on the orbitals centered at u_i, has
hopping entries t_e = lambda T(|u_i - u_j|) in the units of H_d, up to an
error of order p(n,D) lambda^{-0.1}. So to realize the target t_e^tar of
(23), one must choose each link length d_e so that
lambda T(d_e) = tau sqrt( K_e (U - v(d_e)) ) . (27)
The occupancy penalty involves every pair (through v_{ij}), but on a link
the target depends only on d_e. So (27) is one scalar equation per link.
Define the residual
F_e(d) = lambda T(d) - tau sqrt( K_e ( U - v(d) ) ) . (28)
A root exists. Evaluate at d = (1 - eps/2) D and d = (1 + eps/2) D, both
inside the adjustable range of the contact-gadget layout:
T((1-eps/2)D) ~ exp(-(1-eps/2)D) => lambda T ~ lambda^{eps/2} up to poly(D),
T((1+eps/2)D) ~ exp(-(1+eps/2)D) => lambda T ~ lambda^{-eps/2} up to poly(D).
The target tau sqrt(K_e(U-v(d))) is bounded above by p(n) (tau and K_e are
polynomial, U-v is bounded) and below by a fixed inverse polynomial (tau
inverse poly, K_e inverse poly, U-v >= c). At the near endpoint the first
term dominates => F_e > 0. At the far endpoint the first term is dwarfed
=> F_e < 0. So by continuity F_e has a root in between.
Bisection with approximate signs. Evaluate F_e to inverse-polynomial
accuracy. If the sign is ambiguous at evaluation accuracy, the residual is
already small enough, so stop. Otherwise keep the sign-change bracket, and
continue until bracket width times a Lipschitz bound is within tolerance.
No monotonicity or derivative lower bound is needed: only need some root
with small residual, and any root works.
Resulting centers. The calibrated distances d_e, rescaled back, are the
planar separations |u_i - u_j| for link pairs. the contact-gadget layout's contraction then
produces internal-site coordinates realizing all these link lengths at once
with nonlinks separated; rational approximation of the resulting positions
gives the centers u_i.
From now on, H_d means H_d(t) at these final centers, with
t_e = lambda T(|u_i - u_j|). Choose the per-link calibration accuracy so
that replacing t_e^tar by t_e changes the total Hamiltonian norm by at most
gamma tau^2/100. Combining this and the second-order remainder from
the superexchange calculation (which is at most gamma tau^2/100) and the internal error of
Lemma 1, one gets
| E(H_d) + sum_{i<j} v_{ij} - tau^2 ( E(H_+) - sum_e K_e ) |
<= 2 gamma tau^2 / 100 . (29)
Line (29) says the bottom of H_d, shifted by the direct-Coulomb scalar
sum_{i<j} v_{ij}, equals tau^2 times the positive-exchange ground energy
H_+ minus a scalar, up to the stated error. This is the finite-fermion
target to be realized as a continuum molecule.
Every analytic error produced below has the form p(n,D) lambda^{-alpha}, with required smallness at worst inverse polynomial in n. The scales lambda, D, and the cubature cell h are then chosen so that all conditions hold at once.
Take k = ceil(C n^A), h = 1/k, lambda = 8 k^3/rho, with A, C fixed positive integers chosen large enough. The relation rho h^3/8 = 1/lambda gives each of the eight equal nodes of one cubature cell (side h) mass 1/lambda, so each node becomes a unit nucleus in physical coordinates.
Finally D = log lambda. Then D is polynomial in n, and "exp(-D) beats
poly(n)" and similar statements hold by taking lambda a large enough power
of n.
Now construct an attractive Coulomb potential whose low-energy matrix
realizes the prescribed modes and hoppings. The potential is the Coulomb
field of an honest nonnegative charge density: a broad constant-density
slab (supplying transverse confinement and a positive background) plus the
planar wells (supplying the orbitals). Adding the wells changes the density
by a bounded, compactly supported term, so the total density stays
nonnegative.
The construction has three parts: the slab and the well-depth
counterterm, formula (30); a fixed one-particle gap outside the modes
(Proposition 1); the compressed one-particle matrix on the modes.
Everything is in dilated coordinates y = (r,z); energies are the dilated
ones (divide final energies by lambda^2 later). The nuclei are not yet
points; they will be the density's field.
Choose integers
S with D <= S <= 10 D , H = ceil( lambda^5 ) , (32)
and the box Omega = [-H,H]^2 x [-S,S]. (The exact scale of S costs only
poly(D) factors; a binary bit-length scale works. H is huge so the
horizontal sides are far away.) Define the slab field and its value at 0:
F_H(y) = - rho int_Omega db / |y-b| , C_H = F_H(0) . (33)
The slab potential F_H is the attraction of a uniform density rho on the
box; because the box is tall and wide, near the planar wells F_H behaves
like a constant minus a vertical harmonic well centered on z=0, plus small
error (see (34)).
Why we need a well-depth adjustment. The occupancy identity (C) of the
superexchange calculation requires a one-body diagonal term -s_i/lambda on
site i, where
s_i = sum_{j!=i} v_{ij}. In the mode basis, the diagonal entry from
-s_i/lambda must be supplied by a small deepening of well i. Because
<psi_i, W_i psi_i> = int varphi^2 W < 0 (W nonpositive), increasing the
depth of well i lowers its diagonal energy, so the required sign is correct
and the size is tunable. Let
W_i(r) = W(r - u_i) (the planar well, translated to site i).
Fix a smooth cutoff zeta with 0 <= zeta <= 1, supported inside (-1,1),
equal to 1 for |z| <= 1/2, with piecewise polynomial rational data. Then
put
b_i = s_i / ( - int varphi^2 W ) >= 0 ,
V_ell(r,z) = zeta(z/S) sum_i (1 + b_i/lambda) W_i(r) . (35)
The b_i are nonnegative and polynomially bounded independently of D (they
depend only on the geometry of the wells through v_{ij} and the fixed
integral -int varphi^2 W). Hence 0 <= b_i/lambda < 1 at the scales. The
functions V_ell are supported in disjoint vertical cylinders around the u_i,
inside Omega.
Choosing rho. The derivatives of V_ell one needs are bounded independently
of S (the zeta factor's derivatives are controlled by S >= D and the well
supports are disjoint in the plane; the b_i/lambda terms are small). Use
1 + b_i/lambda <= 2 and the fixed bounds on W, zeta. Choose the fixed
integer rho so large that
| Delta V_ell / (4 pi) | <= rho / 2 everywhere. (30)
There is no circularity: the bound used to pick rho depends only on
2, W, zeta, not on the vertical Gaussian or on b_i. The Gaussian and the
exact b_i are then fixed afterwards as needed.
The density. Using the fundamental solution of the Laplacian
(-Delta_y (1/|y-b|) = 4 pi delta_b) and the compact support of V_ell inside
Omega, the whole potential F_H + V_ell is expressed as a Coulomb integral of a
nonnegative density:
V_ell(y) = - int_Omega ( Delta V_ell(b) / (4 pi) ) / |y - b| db ,
F_H(y) + V_ell(y)
= - int_Omega [ rho + Delta V_ell(b)/(4 pi) ] / |y-b| db . (36)
By (30) the bracketed density rho + Delta V_ell/(4 pi) lies between rho/2
and 3 rho/2 on Omega. Its total mass is exactly rho |Omega|, because
int Delta V_ell = 0 (V_ell is compactly supported and smooth). The centers
fit well inside the horizontal box at the scales. So (36) is a
nonnegative continuous charge density whose Coulomb field is the exact
one-particle potential use below (before the point-nucleus replacement).
Background estimates. The following hold with constants independent of
n and lambda:
(a) F_H(y) - C_H >= 2 pi rho min(z^2, S^2) - C S^3 / H ,
(b) | F_H(y) - C_H - 2 pi rho z^2 | <= C ( S^3 + S|r|^2 ) / H ,
for |r| <= H/2 and |z| <= S . (34)
Justification of (34). At every z, F_H is minimized over r at r=0, and is
also nondecreasing in |z| for fixed r: this follows one coordinate at a time
by integrating an even decreasing kernel over a centered interval. At r=0,
extending the horizontal integration to the plane in the difference between
z and 0 gives, for |z| <= S, exactly
2 pi rho int_{-S}^{S} ( |z-w| - |w| ) dw = 2 pi rho z^2
(integrate by parts / use polar integration around the axis). The omitted
difference from restricting the horizontal square costs at most
c s^3 / H, by inverse-cube decay of the kernel outside horizontal radius H.
For |r| <= H/2, translating the integration square gives a change at most
c s |r|^2 / H relative to the horizontal center. The first derivatives then
come from sides at distance >= H/2 (bounded by C S/H per unit distance),
and second derivatives are bounded by C S / H by integrating first
derivatives of the kernel on the sides. Evenness at r=0 gives (a) at z=0.
Finally, on the whole box, 0 >= F_H(y) >= C_H (F_H is negative and most
negative at y=0), and |C_H| <= C S H.
Together (36) and (34) say: near the wells the potential is the vertical
harmonic oscillator 2 pi rho z^2 (min energy omega/2, gap omega) plus the
planar multiwell sum_i W_i, plus perturbations. The precise statements follow.
The potential F_H + V_ell is now expressed as the Coulomb field of a
nonnegative density (36). Before replacing that density by point nuclei,
one needs two properties:
- a fixed energy gap above the modes, controlling arbitrary continuum
states;
- a small residual on the modes, controlling mixing with the mode space.
The shifted one-particle operator. Define
L_c = - Delta_y / 2 + F_H + V_ell - C_H - e_* . (37)
(The spin identity is suppressed.) It is the one-particle operator of the
continuum, shifted to put the mode bottoms near zero. Recall e_* =
-1/2 + omega/2.
The mode projection. Let P project onto the span of the psi_i, including
both spins when present. Let Gamma_ij = <psi_i, psi_j> be the spatial Gram
matrix. Exponential tails and the 0.95 D spacing give the bound
|| Gamma - I || <= p(n,D) lambda^{-0.8} . (38)
(Indeed each off-diagonal overlap is O(e^{-0.9 |u_i-u_j|}) <= e^{-0.9*0.95 D};
leave part of the decay as an integrable factor and sum poly many.) Use the
orthonormal modes
psihat_j = sum_i psi_i (Gamma^{-1/2})_{ij} ,
|| Gamma^{-1/2} - I || <= p(n,D) lambda^{-0.8} . (39)
Then unit vectors in Ran P have polynomially bounded one-particle kinetic
energy (they are finite combinations of the smooth modes with weights of
polynomial size).
Proposition 1 (continuous one-particle control).
With the density and separated centers constructed above, once the
polynomial scale is sufficiently large, there are fixed constants
c > 0 and C < infinity such that, as quadratic-form inequalities,
(I-P) L_c (I-P) >= c (I-P) on the complement,
L_c >= - Delta_y/2 - C (kinetic lower bound),
|| L_c P || <= p(n,D) lambda^{-0.9} (operator norm Ran P -> L^2).
(40)
All estimates hold with either spin included.
Proof.
Step 1: lower bound by a separable operator. By (34)(a), (36), and
zeta W_i >= W_i (zeta <= 1, W_i <= 0) with the b_i/lambda terms small and
bounded everywhere, L_c is at least the separable operator
A_sep = (-Delta_r/2 + sum_i W_i) + (-d^2/dz^2/2 + 2 pi rho min(z^2,S^2))
- e_* , (41)
minus an error O(S^3/H) + O(p(n,D)/lambda) that is uniform in all states
(as a form inequality). Call the first bracket (planar) A_pl and the second
(vertical) A_v, each shifted so its ground energy is near zero:
A_pl has ground energy -1/2 (mode varphi at each site) and A_v has ground
energy omega/2 (mode g).
Step 2: the planar gap. Use a quadratic partition of unity in the plane:
cutoffs theta_i about each site u_i (equal to 1 inside radius D/8, supported
at radius D/3) plus one exterior function theta_ext. Build them by smooth
angle interpolation between 0 and pi/2 on the transition annuli using sine
and cosine, so that sum_i |nabla theta_i|^2 + |nabla theta_ext|^2 <= C/D^2.
The IMS localization identity (product rule with kinetic coefficient 1/2)
then costs O(1/D^2) in the form. On each site term use the isolated ground
state and gap of -Delta_r/2 + W; on the exterior term use the
free kinetic lower bound and the fact that no well sits on its support. With
site cutoffs theta_i and varphi_i = varphi(.-u_i), the site term evaluated
on theta_i u is at least
(-1/2 + c_1) ||theta_i u||^2 - c_1 |<theta_i varphi_i, u>|^2 , (42)
where c_1 > 0 may be taken below both the isolated gap and 1/2. The exterior
term also absorbs its share of the (-1/2 + c_1) mass, because its support
carries no well. The sum of rank-one operators |theta_i varphi_i><...| over
i differs from P_r = sum_i |varphi_i><varphi_i| by at most
O(m e^{-c' D}) + O(p(n,D) lambda^{-0.8}) (mode tails and Gram correction).
Hence
A_pl >= -1/2 + c_1 (I - P_r) - o(1) I . (43)
Simon, B., "Semiclassical analysis of low lying eigenvalues. I,
non-degenerate minima: asymptotic expansions", Ann. Inst. H. Poincare A
38 (1983), 295-308. local: AIHPA_1983__38_3_295_0.pdf
Step 3: the vertical gap. Localize the z problem with an inner cutoff equal
to 1 on |z| <= S/2 (supported in |z| <= S) and an exterior function. On the
inner term use the harmonic oscillator (ground omega/2, gap omega); on the
exterior term use the large-z potential lower bound 2 pi rho min(z^2,S^2)
>= 2 pi rho (S/2)^2 for |z| >= S/2, which dwarfs the Gaussian tail loss. The
localization costs O(S^{-2}) in the form, and the rank-one mode is P_g =
|g><g| with a Gaussian-weight loss of order e^{-c S}. Hence
A_v >= omega/2 + c_2 (I - P_g) - o(1) I . (44)
Step 4: tensor the gaps. On the spatial space P = P_r (tensor) P_g. The
operator inequality
(I - P_r) (tensor) I + I (tensor) (I - P_g) >= I - P (45)
holds because the left side is the indicator of "at least one factor is
outside", which dominates "product is outside". Combining (43) and (44)
through (45), the separable operator A_sep has, on I - P, a fixed gap
A_sep - e_* >= (-1/2 + omega/2 - e_*) + c(I-P) - o(1) I
= c(I-P) - o(1) I , (46)
since e_* = -1/2 + omega/2 makes the constant term zero. Choosing the
polynomial scale so the total error O(S^3/H) + O(p(n,D)/lambda) + o(1) is
at most c/2 gives the first inequality in (40). The slowly small errors
o(d^{-2}), o(s^{-2}) are only needed to establish a fixed gap, not a per-
particle energy accuracy, so they do not blow any budget.
Step 5: kinetic lower bound and the mode-space norm. The kinetic lower bound
L_c >= -Delta_y/2 - C follows from (34) (F_H - C_H bounded below by
-c s^3/h >= -C on the useful region and 0 >= F_H >= C_H off it) together with
bounded V_ell and the constant -C_H - e_* absorbed into C. For the operator
norm ||L_c P||: (34) and the crude global background bound plus mode tails
give
|| ( F_H - C_H - 2 pi rho z^2 ) psi_i ||_2 <= p(n,D) lambda^{-2} .
(47)
Indeed only moments with |r| <= H/2 and |z| <= S matter there, and outside
that region the Gaussian in z or the planar exponential in |r| beats the
polynomial magnitudes, using S >= log lambda. The cutoff defect in V_ell
versus sum_i (1 + b_i/lambda) W_i acting on any mode is likewise negligible.
Now (-Delta_y/2 + 2 pi rho z^2 + W_j - e_*) psi_j = 0 exactly (tensor of the
planar ground eigenvalue -1/2 and the vertical ground eigenvalue omega/2,
shifted by e_*). So the residual L_c psi_j is (i) the other wells W_i,
i != j, acting on the small planar tails of psi_j, and (ii) the term
(b_j/lambda) W_j on psi_j, and (iii) the errors already bounded in (47).
All three are of size at most p(n,D) lambda^{-0.9} on unit mode vectors,
giving the last inequality in (40). QED
A fixed complementary gap is not enough: one needs the matrix on the modes to
a sharper scale, because the continuum energy will be compared to H_d after
multiplying by lambda. This subsection computes that matrix.
Compute the continuous one-body compression. On the nonorthogonal functions
psi_i (before Gram correction), the column on mode j comes, up to the
background and cutoff errors already bounded in the one-particle-gap step, from applying
( sum_{k != j} W_k + lambda^{-1} sum_k b_k W_k ) psi_j . (48)
Row j (diagonal). The depth-deepening b-term contributes diagonally as
follows: with I = <psi_j, W_j psi_j> = int varphi^2 W < 0, the j-th piece
of
lambda^{-1} sum_k b_k W_k on psi_j is lambda^{-1} b_j I = lambda^{-1}
( s_j / (-I) ) I = - s_j / lambda. Since s_j = sum_{i != j} v_{ij}, this is
the same number as -lambda^{-1} sum_{i != j} v_{ij} promised in the
construction. The other-well terms sum_{k != j} W_k acting on psi_j carry two
planar tail factors (spacing >= 0.95 D) and are exponentially small, so they
do not affect the diagonal. Hence A_jj = -s_j/lambda up to negligible error.
Row i != j (off-diagonal). In the first sum sum_{k != j} W_k, the term with
k = i gives -T(|u_i - u_j|) (this is exactly the hopping profile (25),
evaluated with the well at i acting on psi_j). The cases k != i, j are
products of two tail factors at a third site; by the spacing and compact
well support, each such integral is O(lambda^{-1.8}). The explicit
lambda^{-1} sum_k b_k W_k terms give only negligible off-diagonal errors (at
least one planar tail on each well support), bounded by p(n,D) lambda^{-1.9}
altogether per entry. Nonlink hoppings cost at most p(n,D) lambda^{-1.2} per
term, because two wells at a nonlink separation overlap by two tails.
Pass to orthonormal modes. Let A_ij = <psi_i, L_c psi_j> (the shifted
matrix: the constant -C_H - e_* is already subtracted inside L_c, which
removes the large background scalar; do not re-add it). The matrix in the
orthonormal modes is Gamma^{-1/2} A Gamma^{-1/2}. Using (38) and (40),
|| Gamma^{-1/2} A Gamma^{-1/2} - A ||
<= C || Gamma - I || * || A || <= p(n,D) lambda^{-1.7} . (49)
The correction therefore remains negligible after multiplication by lambda.
The spatial matrix. Define
M_ij = - s_i if i = j,
= - lambda T(|u_i - u_j|) if {i,j} is a link,
= 0 otherwise. (50)
Putting the row estimates and (49) together,
|| lambda Gamma^{-1/2} A Gamma^{-1/2} - M || <= p(n,D) lambda^{-0.1} .
(51)
with the identity on spin. This is the required one-particle compression:
the diagonal supplies the depth counterterm (the -s_i one-body term), and
the link entries supply the calibrated hopping -lambda T(d_e). Comparing
with H_d of (20), whose hopping is -t_e on links and whose penalty
produces -s_i as the linear counterterm, one sees the compressed one-particle
matrix is exactly the one-particle part of H_d, up to p(n,D) lambda^{-0.1}.
Where we are. There is a continuous charge density whose one-particle
operator L_c has (i) the modes as approximate low-energy states with a fixed
complementary gap, and (ii) a compressed matrix equal to the one-particle
part of H_d to accuracy p(n,D) lambda^{-0.1}. The cubature replacement replaces the
density by honest unit point nuclei; the comparison with H_d is then completed.
The input to the slab-and-wells density construction was a nonnegative charge density (36). One must now
replace its Coulomb field by finitely many point nuclei of equal mass
1/lambda, so that in physical coordinates each nucleus has charge exactly 1.
There are two distinct error requirements:
(i) a form bound on arbitrary continuum states (weak, but on all states);
(ii) a compression bound on the smooth localized modes only (stronger).
The comparison with H_d uses energies multiplied by lambda. So the
compressed error must be smaller than lambda^{-1} up to a polynomial margin.
The error coupling the mode space to its complement can be larger: the fixed
complementary gap makes its energy contribution quadratic in
the coupling, so a weaker bound suffices there.
Cubature side h will be chosen as
h = (8 / (rho lambda))^{1/3} , so rho h^3 / 8 = 1/lambda . (52)
A cubature rule exact through polynomial degree three will then give an
O(h^4) compressed error and an O(h^2) form error. Both resulting energy
errors are O(lambda^{-1/3}) in the units of H_d, up to polynomial factors.
The replacement uses Moser's deformation method (Moser 1965; Dacorogna-Moser 1990): tile the uniform slab with an evenly spaced cubature, then push nodes by a smooth compactly-supported flow so that the pushed-node density equals the target smear; the discrete Coulomb field mimics the continuum field on smooth test functions, and every node keeps common weight (charge one).
Define the interpolation path from the uniform slab density rho to the
target density,
q_t = rho + t Delta V_ell / (4 pi) , 0 <= t <= 1 , (53)
so q_0 = rho and q_1 = rho + Delta V_ell/(4 pi) is exactly the density of
(36). Since by (30) |Delta V_ell/(4 pi)| <= rho/2, there is q_t >= rho/2
for all t. On full space define the time-dependent vector field
V_t(b) = - nabla V_ell(b) / (4 pi q_t(b)) , 0 <= t <= 1 . (54)
Because V_ell is smooth, compactly supported inside the well cylinders
(strictly inside the box), and the denominator q_t >= rho/2, the field V_t
has uniformly bounded derivatives through order five and vanishes outside
those cylinders.
The flow. Let G_t be the flow of V_t starting at time zero, and put G = G_1.
Standard ODE theory (differentiate the equation, induct on the order) gives
bounded derivatives of G_t through order four and a bounded Lipschitz
inverse. The flow is the identity off the cylinders, so it preserves Omega.
Continuity equation. Directly from V_t = -nabla V_ell/(4 pi q_t) and
partial_t q_t = Delta V_ell/(4 pi),
partial_t q_t + div( q_t V_t ) = 0 . (55)
Integrating along the flow gives the transport identity
q_t(G_t(b)) det DG_t(b) = rho , and hence
int_Omega F(b) q_1(b) db = rho int_Omega F(G(b)) db (56)
for any integrable F supported in Omega. So the flow carries the uniform
measure rho db onto the target density q_1 db. The second line is the crucial
one for cubature: applying F to a transported point G(b) and keeping the
weight rho db reproduces int F q_1.
The cubature. Because H, S, k are integers (k = ceil(C n^A), h = 1/k),
cubes of side h tile Omega. On each cube Q use the eight tensor products of
two equal-weight Gauss points per coordinate, at the midpoint plus or minus
h/(2 sqrt(3)). This eight-point rule integrates polynomials of total degree
at most three exactly, by matching the moments on each interval. When mapped
these nodes by the flow G, their common weight stays rho h^3 / 8 = 1/lambda
by (52). (The rule's exactness is used in the original cube coordinates,
on integrands composed with G.)
The rational rounding. Each exact mapped point b_alpha (the image under G of
a cubature node) will be approximated within Euclidean error at most
delta_R <= (1+H)^{-3} (1+S)^{-2} lambda^{-3} (57)
by a rational point b~_alpha. The physical nuclear positions are then
x~_alpha = b~_alpha / lambda. Because the exact points have separation
controlled by node separation and the bi-Lipschitz bound of G, and delta_R
is far smaller than that separation, the rounded points stay distinct.
So the data fed to the physical problem are:
- polynomially many (namely 8 |Omega| h^{-3}) rational positions
x~_alpha = b~_alpha / lambda,
- every one of them a unit nucleus (the physical charge is then exactly 1),
- the electron number N = m,
- rational thresholds from the computability discussion.
The single-particle operator with the actual point nuclei is
L_p = - Delta_y/2 - sum_alpha (1/lambda) / | y - b~_alpha | - C_H - e_* ,
(58)
where b~_alpha are the rounded node positions and the common charge is
1/lambda in dilated units. The difference from L_c is exactly the difference
between the density's Coulomb field and the point-charge sum. Splitting it
into cubature error and rounding error:
L_p - L_c = [ rho int_Omega db / |y - G(b)|
- sum_alpha (1/lambda) / | y - b_alpha | ]
+ [ sum_alpha (1/lambda)/|y - b_alpha|
- sum_alpha (1/lambda)/|y - b~_alpha| ] . (59)
(The first bracket uses (56): rho int_Omega F(G(b)) db = int F q_1, with
F(x) = 1/|y-x|, so the density field equals rho int_Omega db/|y-G(b)|.)
Now bound each bracket as a form and as a compression.
Proposition 2 (unit-nucleus replacement).
For the transported degree-three cubature and rational rounding above,
with delta_R <= (1+H)^{-3}(1+S)^{-2}lambda^{-3}, the shifted continuous
and point-nucleus operators L_c, L_p satisfy
| <u, (L_p - L_c) u> | <= p(n,D) lambda^{-2/3} ( ||u||^2 + ||nabla u||^2 ),
|| P (L_p - L_c) P ||_P <= p(n,D) lambda^{-4/3} . (60)
The first estimate holds for every one-particle H^1 function u. The
second is the operator norm of the compressed form on the localized-mode
space Ran P. Both bounds are in dilated energy units.
Proof.
Step 1: the cubature error per cube. Start at the exact transported
nodes b_alpha = G(a_{Q,nu}). On an original grid cube Q with Gauss nodes
a_{Q,nu} (1 <= nu <= 8), Taylor's theorem plus exactness through degree three
give the common error estimate
| rho int_Q F(G(b)) db - (rho h^3/8) sum_{nu=1}^8 F(G(a_{Q,nu})) |
<= C h^7 sup_Q | D^4 (F o G) | , (61)
whenever F o G has bounded derivatives through order four on Q. (The degree-3 exactness
cancels the Taylor terms of order <= 3; the remainder is order 4, integrated
over a cube of volume h^3 with eight points.) Apply this (a) to the
singular Coulomb kernel away from its singularity, then (b) to the smooth
fields K_ij coming from the modes (Step 3).
Step 2: the general form estimate (first inequality in (60)). For fixed
evaluation point y, the charge integral is the kernel |y - G(b)|^{-1} against
uniform density in b. Split cubes Q into close ones (whose images touch a
ball of radius a fixed large multiple C' h about y) and far ones.
close cubes. Only O(1) cubes are close (the images are bi-Lipschitz and
h small). Their continuous contribution to the k=0 term is O(h^2) because
the transformed density (the Jacobian of G) is bounded. The discrete
contribution is sum over close nodes of C h^3 / |y - b_alpha| with
|y - b_alpha| <= C' h, which the Hardy estimate of Step 4 controls exactly
against ||u||^2 + ||nabla u||^2 on balls.
Hardy, G. H., Littlewood, J. E., and Polya, G., "Inequalities", Cambridge
University Press, 1934, Theorem 326. ISBN: 978-0521358804.
Lieb, E. H. and Loss, M., "Analysis", 2nd ed., Graduate Studies in
Mathematics 14, AMS, 2001. ISBN: 978-0821827833.
far cubes. The pure fourth derivative of the kernel is bounded by
C|y - G(b)|^{-5}. Within each far cube, distances to y are comparable, and
the bi-Lipschitz change of variables has bounded Jacobian. Summing the
pure fourth-derivative contributions therefore costs at most
C h^4 int_{|x-y| > ch} |x-y|^{-5} dx = O(h^2) . (62)
The remaining chain-rule terms in D^4(F o G) contain a derivative of G of
order at least two together with a kernel derivative bounded by
C|y - G(b)|^{-j}, 2 <= j <= 4. Such terms are nonzero only on cells meeting
the well cylinders (combined slab volume O(m(S+1))). Below distance one
they are bounded by the same inverse fifth power; beyond distance one their
total contribution is at most C h^4 m(S+1).
So the pointwise absolute error is at most C(h^2 + h^4 m(S+1)) plus the
close-node sum. In particular there is no factor the size of the large
horizontal area in this bound.
Step 3: the compressed estimate (second inequality). For compression use the
smooth Coulomb fields
K_ij(b) = int |y - b|^{-1} psi_i(y) psi_j(y) dy . (63)
Their derivatives through order four are uniformly bounded (differentiate
the translated densities). The fourth derivatives obey the pointwise decay
| D_b^4 K_ij(b) | <= C ( 1 + | b - (u_i,0) | )^{-5} . (64)
(Indeed derivatives of the product density psi_i psi_j decay exponentially
about the planar center (u_i,0), uniformly because the other factor has
bounded derivatives. At a large distance d = |b - (u_i,0)|, use a smooth
density cutoff supported within d/2 of that center and equal to 1 inside
d/3, holding the cutoff radius fixed while differentiating. There the pure
kernel derivative decays; on the remaining tail, differentiate the cutoff
density instead, giving exponentially small bounds.) Applying the same
degree-three cubature (61) to K_ij o G now costs, for each entry,
O( h^4 (1 + m(S+1)) ) , (65)
because pure fourth derivatives have summable volume bounds as in (62),
and the other chain-rule terms live only on the exceptional cells. Gram
correction (via (38)) and operator norms cost at most polynomial factors.
This proves the second inequality of (60) for the exact nodes.
Step 4: the close-node Hardy bound and rounding. The close-node balls have
bounded overlap also upon any fixed dilation. Around each node one has
int_{|y - b_alpha| <= C'h} |u(y)|^2 / |y - b_alpha| dy
<= C int_{|y - b_alpha| <= 2C'h} ( h^{-1} |u|^2 + h |nabla u|^2 ) .
(66)
This follows by a cutoff and the three-dimensional Hardy inequality
int |u|^2 / |y|^2 <= 4 int |nabla u|^2 (which follows from the divergence
formula for y/|y|^2 together with Cauchy-Schwarz). Summing (66) over the
bounded-overlap balls gives the general form bound, since h^2 = O(lambda^{-2/3}).
Finally the rounding (second bracket of (59)). A center shift of size
delta_R changes the inverse distance by at most
delta_R / ( |y - b_alpha| |y - b~_alpha| ) . (67)
Cauchy-Schwarz and translated Hardy inequalities bound its form by
4 delta_R ||nabla u||^2. Summing the masses gives the total bound
4 rho |Omega| delta_R ||nabla u||^2 , where
rho |Omega| delta_R <= C H^{-1} S^{-1} lambda^{-3} <= C lambda^{-8}
(68)
(using |Omega| = 8 H^2 S, delta_R <= (1+H)^{-3}(1+S)^{-2} lambda^{-3},
H = lambda^5, S >= log lambda). This is smaller than both required errors,
also after compression, because the modes have polynomial kinetic bounds
(39). This completes (60). QED
This yields the two replacement estimates (60) with constants
independent of the horizontal area of the slab; the full mass enters only
the rounding bound (68), where the chosen coordinate accuracy absorbs it.
These combine with the fixed complementary gap (40) and the compressed
matrix (51) to compare the full continuum with H_d.
The point nuclei are now fixed. One must compare the electronic
energy with the finite Hamiltonian H_d on all continuum states.
In dilated coordinates the full operator after dividing by lambda^2 and
shifting by -m(C_H + e_*) is
Hcal = sum_{l=1}^m (L_p)_l + lambda^{-1} R ,
R = sum_{1<=l<s<=m} |y_l - y_s|^{-1} , (69)
with L_p from (58). Work on the antisymmetric space with spins and its
H^1 form domain, exactly as in the physical model.
The form is closed and lower bounded. Fiberwise Hardy together with
Cauchy-Schwarz gives, for every eta > 0 and every fixed center b,
int |u|^2 / |y_l - b| <= 2 ||u||_2 ||nabla_l u||_2
<= eta ||nabla_l u||_2^2 + eta^{-1} ||u||_2^2 . (70)
The same estimate applies to a pair singularity by fixing the other
coordinates. For each finite nuclear configuration, summing (70) with a
sufficiently small eta gives a closed, lower-bounded form whose shifted form
norm is equivalent to the H^1 norm. Smooth compactly supported antisymmetric
functions form a core. So the spectral infimum is the infimum of the form
over normalized H^1 states, and every form estimate below extends by
continuity from the core. No ground eigenvector is assumed.
The projection Pi. Let Pi project onto the subspace in which all m particles
lie in P (the m-electron space built from the 2m one-particle spin modes).
Write ||u||_1^2 = ||u||^2 + sum_l ||nabla_l u||^2. Because the modes are
smooth and decaying, they lie in the form domain, so the compression
(Pi Hcal Pi)|_{Ran Pi} is well defined and equals, by (51), the finite
fermion operator below. The one-particle kinetic bound (39) gives
||w||_1 <= p(n,D) ||w|| for w in Ran Pi.
E_Pi = min spec( ( Pi Hcal Pi ) |_{Ran Pi} ) . (71)
This is a genuine finite-dimensional minimum, because Pi projects onto the
m-electron space of 2m one-particle spin modes.
The complement receives a uniform lower bound, and its coupling to the mode space enters the energy only quadratically by square completion.
Proposition 3 (eliminating the continuum complement).
For the constructed point nuclei with N = m and sufficiently large
polynomial scale lambda, the full antisymmetric m-electron operator in
dilated and shifted units satisfies
E_Pi - p(n,D) lambda^{-4/3} <= inf Spec Hcal <= E_Pi . (72)
The lower bound holds on every normalized antisymmetric H^1 state,
including every spin sector.
Proof.
Step 1: complement lower bound for the one-body part. Write
H_c^0 = sum_l (L_c)_l. By (40), as one-particle forms,
L_c >= c' (I - P) - p(n,D) lambda^{-0.9} I (73)
for some c' > 0 (the diagonal and cross terms involving P are controlled by
||L_c P||). Since sum_l (I - P)_l >= I - Pi (at least one particle outside
means the product is outside), summing (73) gives a fixed positive lower
bound for H_c^0 on Ran(I - Pi), once m p(n,D) lambda^{-0.9} is small.
Also valid is the kinetic lower bound H_c^0 >= -sum_l Delta_l/2 - C m of
(40). A
convex combination with a small weight of order 1/m on the kinetic bound keeps
a positive L^2 term and gives
<q, H_c^0 q> >= (c''/m) ||q||_1^2 , c'' > 0 , (74)
for q in Ran(I - Pi). By (60) on each coordinate and the positivity of R,
this estimate persists with half the coefficient for the full Hcal.
Step 2: the diagonal on the mode space and the mixed block. For
w in Ran Pi, (40) gives ||H_c^0 w|| <= p(n,D) lambda^{-0.9} ||w||. Hardy
applied one particle coordinate at a time gives ||R w|| <= p(n,D) ||w||.
Polarizing the form bound (60) and using the H^1 bound on w controls the
point-nucleus replacement in the mixed block. Together,
| <q, Hcal w> | <= p(n,D) lambda^{-2/3} ||w|| ||q||_1 , (75)
for w in Ran Pi and q in Ran(I - Pi). The same estimates show the compressed
form has norm less than c''/(4m) at sufficiently large scale.
Step 3: square completion. Set alpha = c''/(4m) and let
beta = p(n,D) lambda^{-2/3} bound (75). Decompose any normalized
form-domain vector as u = w + q, with w = Pi u and q = (I - Pi) u.
Subtracting E_Pi makes the w-diagonal form nonnegative and leaves at least
alpha ||q||_1^2 on the complement (74). Hence, all pairings interpreted as
forms,
<u, Hcal u> - E_Pi ||u||^2
>= alpha ||q||_1^2 - 2 beta ||w|| ||q||_1
>= -(beta^2/alpha) ||w||^2
>= -p(n,D) lambda^{-4/3} , (76)
using the elementary inequality a s^2 - 2 b s t >= -b^2/a t^2. This proves
the lower bound in (72); the upper bound follows by testing a minimizing
vector of the compressed operator (which lies in Ran Pi). The lower-bound
argument applies to every normalized antisymmetric state, in any spin sector,
without any binding or ground-eigenvector assumption. QED
It remains to evaluate lambda E_Pi and match it to E(H_d). the compression computation
computed the one-particle matrix M on the orthonormal modes to accuracy
p(n,D) lambda^{-0.1}. Multiplying the compressed replacement error (60) by
lambda costs only p(n,D) lambda^{-1/3}, so the same matrix accuracy
p(n,D) lambda^{-0.1} holds for the point nuclei in the units of H_d.
The pair matrix. For the repulsion R, the pair matrix on the original modes
(before Gram correction) has spatial entries from products of two transition
densities psi_i psi_k acting on psi_j psi_l. Keeping only psi_i^2, psi_j^2 as
the respective densities gives the direct coefficients v_{ij} of (13). Any
unequal indices within a density give an exponentially small error:
|| psi_i psi_k ||_{L^1} = O(lambda^{-0.8}) for i != k , (77)
and the Coulomb field of the other density product is uniformly bounded in
absolute value by orbital boundedness and decay. The matrix of direct
interactions in the product basis is bounded, and Gram correction changes
its entries negligibly (by (38)).
Dimension dependence. These entrywise estimates give a many-electron
operator-norm error of at most p(n,D) lambda^{-0.1}. To see the dimension
factor, write the one- and two-body matrices in the 2m orthonormal spin
modes. Each coefficient multiplies a monomial c_a^dagger c_b (one-body) or
c_a^dagger c_b^dagger c_d c_c (two-body), whose norm is at most one. There
are only O(m^2) and O(m^4) such monomials, so summing their coefficient
errors costs a polynomial factor over the entire fermion space.
Identifying the retained operator. In this representation, the retained
two-body term is
(U/2) sum_i n_i (n_i - 1) + sum_{i<j} v_{ij} n_i n_j ,
and the retained one-body diagonal is
- sum_i ( sum_{j != i} v_{ij} ) n_i = - sum_i s_i n_i .
At the filling sum_i n_i = m, the identity (C) of the superexchange calculation combines these
two lines into the first line of H_d (20), including its scalar
-sum_{i<j} v_{ij}. The hopping matrix M from (50) gives the remaining line
of the calibrated H_d. So the compressed operator is exactly H_d.
Proposition 4 (continuum-to-fermion energy comparison).
For the constructed unit nuclei with N = m, let Hcal be the dilated and
shifted continuum operator (69) and H_d the finite fermion Hamiltonian
(20). At sufficiently large polynomial scale,
| lambda inf Spec Hcal - E(H_d) | <= p(n,D) lambda^{-0.1} . (78)
Proof. The one- and two-body matrix estimates just obtained compare
lambda E_Pi with E(H_d) within p(n,D) lambda^{-0.1}: the one-body part
matches by (51), the pair part by (77) and (79). Multiplying (72)
by lambda contributes only p(n,D) lambda^{-1/3} (because lambda times
lambda^{-4/3} is lambda^{-1/3}), which fits within the bound. QED
Distinct estimates play distinct roles; in particular the arbitrary-state
form error cannot be used as a direct energy error after multiplication by
lambda.
+----------------------------------+----------------+---------------------+
| Estimate | Error scale | Role |
+----------------------------------+----------------+---------------------+
| General form (60) | lambda^{-2/3} | Mixed blocks, H^1 |
| Compressed form (60) | lambda^{-4/3} | Direct on modes |
| Complement (72) | lambda^{-4/3} | Square of mixed blk |
| Last two x lambda | lambda^{-1/3} | In units of H_d |
| Total comparison (78) | lambda^{-0.1} | + other mode errors |
+----------------------------------+----------------+---------------------+
Rows are in the units shown, up to polynomial factors p(n,D). The O(D^{-2}) and O(S^{-2}) localization losses only establish a fixed complementary gap and do not enter the energy-precision budget.
Choose the polynomial scale in lambda large enough that (78) costs at
most tau^2 gamma/100 and that
lambda tau^2 gamma >= 8 . (80)
All earlier smallness, separation, and calibration estimates hold
simultaneously by the growth of k. Lifting from the
dilated coordinates back to physical coordinates gives the physical
electronic ground energy exactly as
E_0 = lambda^2 ( m(C_H + e_*) + inf Spec Hcal ) . (81)
The affine map. Combine (9), (29), (78) and (81). Define
F(xi) = lambda^2 m (C_H + e_*)
+ lambda [ - sum_{i<j} v_{ij}
+ tau^2 ( xi - C_s - sum_e K_e ) ] . (82)
The total error is at most 7 lambda tau^2 gamma/100, and hence at most
lambda tau^2 gamma/8:
| E_0 - F( E(H_in) ) | <= lambda tau^2 gamma / 8 . (83)
The error decomposition: (9) contributes 4 gamma/100 in the
H_in variable (which maps to at most 4 lambda tau^2 gamma/100 through F);
(29) contributes 2 gamma tau^2/100; and (78) contributes at most
tau^2 gamma/100 by choice of lambda. Summing gives 7/100 < 1/8.
The output thresholds. Take rational numbers a and b approximating
a ~ F( a_s + gamma/4 ) , b ~ F( b_s - gamma/4 ) , (84)
each within lambda tau^2 gamma/32. The computability discussion shows these are computable
and have polynomial rational bit length.
Yes case. Suppose the source is YES: E(H_in) <= a_s. Then
E_0 <= F(a_s) + lambda tau^2 gamma/8
< F( a_s + gamma/4 ) - lambda tau^2 gamma/32
<= a . (85)
The middle step uses F(a_s + gamma/4) - F(a_s) = lambda tau^2 gamma/4, which
is bigger than lambda tau^2 gamma/8 + lambda tau^2 gamma/32.
No case. Suppose the source is NO: E(H_in) >= b_s. Then, symmetrically,
E_0 >= F(b_s) - lambda tau^2 gamma/8
> F( b_s - gamma/4 ) + lambda tau^2 gamma/32
>= b . (86)
Threshold separation. Since F is affine with slope lambda tau^2 and
b_s - a_s >= gamma,
b - a >= lambda tau^2 ( b_s - a_s - gamma/2 ) - lambda tau^2 gamma/16
>= lambda tau^2 (gamma/2) - lambda tau^2 gamma/16
= (7/16) lambda tau^2 gamma
>= 7/2 > 1 . (87)
using (80) in the last line. So both promises survive rational rounding,
and the reduction meets the required output separation b - a >= 1.
The configuration and thresholds are next certified to be computable in deterministic polynomial time with polynomial encoding length. The scale order is: (1) bounded well data (f, phi, W) and fixed integer rho from derivative bound (30); (2) source cutoff and spin-gadget penalties Delta, Delta', Delta'' from n and gamma by the normalization to rationals; (3) rational tau from the uniform Coulomb coercivity constant c and the polynomial size bounds; (4) integers k = ceil(C n^A), h = 1/k, lambda = 8 k^3/rho with A, C fixed absolute positive integers then large enough.
At stage 4, K_e, 1/K_e, the layout extent, and the b_i are already polynomial in n, so they do not affect k. Every remaining condition is p(n, log lambda) lambda^{-alpha} <= inverse polynomial in n, and the hopping endpoint brackets of (27) sit at the same fixed-power rate; the polynomials p(n,D) come from fixed-order derivatives and fixed-dimensional data, and the large box enters only through (34) and summable cubature derivatives. A large enough A meets every asymptotic requirement; C covers all n >= 2 plus the fixed conditions; hence H, |Omega|, and 1/delta_R are polynomial in n.
The large scalar shift. The dominating scalar is C_H times lambda^2 m; evaluating it within
tau^2 gamma / ( 128 lambda m ) (88)
uses at most a quarter of the allowed threshold error lambda tau^2 gamma/32. The other scalar terms fit the remaining budget the same way, their amplification factors and count being polynomially bounded.
The orbital integrals. The modes come from phi = (-Delta_r + 1)^{-1} f, so
in the Gaussian heat formula, substitute t = s^2 and rescale the Gaussian
integration variable. For a derivative D^alpha of any required order,
D^alpha phi(r) = (1/pi) int_0^\infty int_{R^2} 2s e^{-s^2} e^{-|a|^2}
D^alpha f(r + 2 s a) da ds . (89)
This integrand is regular at s = 0. The needed derivatives of f are bounded
piecewise polynomials with bounded first derivatives. Polynomial integration
cutoffs make the Gaussian tails smaller than any prescribed inverse-
polynomial error; on the remaining bounded domain the integrand and its
Lipschitz bound are polynomially bounded. Elementary grid sums in this fixed
dimension therefore suffice. Division by phi in W and its derivatives occurs
only on the unit disk, where phi has a fixed positive lower bound. The other
decaying orbital integrals (normalization, hopping integrals) can likewise be
truncated about their centers with polynomial cutoffs.
Coulomb integrations. For v_{ij} and C_H, replace |x|^{-1} inside radius
epsilon by epsilon^{-1}, with epsilon a sufficiently small inverse polynomial.
Integrating the local singularity against bounded densities gives an error of
order epsilon^2 times the relevant polynomial volume bounds. The capped kernel
has Lipschitz constant at most epsilon^{-2}, so grid sums again give the
required accuracy in polynomial time. All integrations have fixed dimension
and fixed nesting depth.
Elementary functions. Exponentials, logarithms, square roots, and other
elementary constants can be evaluated by rational approximations to the
required accuracy, since arguments and reciprocal tolerances are polynomially
bounded. Square roots at zero are handled by bisection. Guard precision for
the nested evaluations keeps accumulated errors within the allotted budgets.
Hopping bisection. The link equation (27) has strict endpoint brackets by
(26) and the polynomial bounds on its target. In approximate bisection,
evaluate the residual to a small fraction of its per-link tolerance. If the
sign is ambiguous, the residual is already small enough; otherwise keep the
sign-change bracket until its width times the Lipschitz bound is within
tolerance. No derivative lower bound is needed. The fixed layout contractions
of the contact-gadget layout and the rational center approximations can be computed more
accurately than these tolerances require.
Once the rational centers u_i are fixed, their exact real coefficients b_i
define a reference density and a reference flow G. The images are approximated
under that flow.
The flow lasts for time one, with uniformly bounded velocities and first
derivatives in space and time (the Moser step). Euler stepping with an
inverse-polynomial time mesh therefore suffices. For mesh nu (reciprocal of
the number of steps), initial error a_0, per-step arithmetic error a_1, and
velocity evaluation error eps_v, the final position error is
O( a_0 + nu + eps_v + a_1 / nu ) . (90)
Indeed a step multiplies the preceding error by at most 1 + C nu and adds
O(nu^2 + nu eps_v + a_1). Each parameter can be chosen inverse-polynomially
small so the final error is at most delta_R. Rounding intermediate positions
on sufficiently fine rational grids keeps bit lengths polynomial. Near the
support boundaries the specified smooth vanishing and derivative bounds give
the same error control. Numerical evaluation of b_i is part of the velocity
evaluation error, so all numerical and final rational-rounding errors can be
charged to the single displacement budget delta_R without changing the
reference density or any cell mass.
None of these operations uses a many-electron wavefunction: only single-particle quadrature, fixed-dimensional flow integration, and rational rounding. There are 8 |Omega| h^{-3} point nuclei, polynomially many in n, each of charge exactly one. Their rational physical positions x~_alpha = b~_alpha / lambda, the electron number N = m (in unary), and the rational thresholds a, b of (84) all have polynomial encoding length and are computable by the scale order, numerical evaluation, and flow evaluation. Reducing rational outputs to lowest terms gives the input format of the promise problem.
On any input, first verify source syntax and the fixed polynomial bounds used in the construction; if they fail, output a fixed instance (which obeys its own promise vacuously and does not affect either promised implication). Otherwise run the deterministic construction above.
The affine map F is (82). The estimates (9), (29), (78) and (81) chain to a single bound on the source-instance molecule,
| E_0 - F( E(H_in) ) | <= lambda tau^2 gamma / 8 . (91)
Each constant entering (91) is fixed by the source instance or by k, A, C, and is computable in deterministic polynomial time to the required inverse-polynomial accuracy. One source instance thus yields one molecule plus thresholds.
The output thresholds are rationals a < b approximating
F(a_s + gamma/4) and F(b_s - gamma/4), each within lambda tau^2 gamma/32.
If the source is YES (E(H_in) <= a_s), then (91) and the definition of a
give E_0 <= a (85). So the molecule is in the YES set.
If the source is NO (E(H_in) >= b_s), then (91) and the definition of b
give E_0 >= b (86). So the molecule is in the NO set.
By (87), b - a >= 7/2 > 1, so the produced instance meets its own input promise with room to spare; every energy error stays inside the built-in slack, so no YES instance becomes a NO instance and no NO becomes a YES.
Theorem (main).
The electronic promise problem for Hamiltonian (1) on the full spinful
fermionic Hilbert space H_N, with distinct rational nuclear positions, all
nuclear charges equal to one, unary electron number, and rational
thresholds separated by at least L^{-1}, is QMA-hard under deterministic
classical polynomial-time many-one reductions. The reduction produces
polynomially many nuclei and electrons, polynomial rational bit lengths,
and threshold separation at least one. The only output data are the
nuclear positions, electron number, and thresholds.
Proof. Given any signed Heisenberg instance, the construction above runs in six deterministic polynomial-time stages: (1) normalize to rational couplings and thresholds, obtaining H_in with gap gamma = min(1, b_s - a_s) > 0; (2) convert to positive exchange by singlet-pair mediators, giving H_+ = sum K_e h^s_{ij} with (9); (3) build the localized modes with direct Coulomb coefficients and calibrated hoppings, giving H_d with (29); (4) replace the smooth slab density by unit point nuclei via Moser transport and degree-three cubature, giving the form and compression bounds (60); (5) compare the full continuum infimum to E(H_d) using the complement gap and square completion, giving (72) and (78); (6) compute the physical thresholds through the affine map F, verifying (85)-(87) and the threshold gap > 1.
All six stages run in deterministic polynomial time; the outputs are polynomially many rational data plus a unary electron number, all at charge-one distinct rational positions. By (85) and (86) both promises are preserved, so the source's QMA-hardness transfers and the electronic problem is QMA-hard. QED