\documentclass[11pt]{article}
\usepackage[margin=1.1in]{geometry}
\usepackage{amsmath,amssymb,amsthm,mathtools,booktabs,longtable,array}
\usepackage{xcolor}
\usepackage[colorlinks=true,linkcolor=blue!50!black,citecolor=blue!50!black,urlcolor=blue!50!black]{hyperref}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{corollary}[theorem]{Corollary}
\theoremstyle{definition}\newtheorem{definition}[theorem]{Definition}
\theoremstyle{remark}\newtheorem{remark}[theorem]{Remark}
\newcommand{\R}{\mathbb R}\newcommand{\E}{\mathbb E}\newcommand{\abs}[1]{\lvert #1\rvert}\newcommand{\ip}[2]{\langle #1,#2\rangle}
\newcommand{\norm}[1]{\lVert #1\rVert}
\newcommand{\M}{\mathcal M}
\renewcommand{\Re}{\operatorname{Re}}
\newcommand{\spec}{\sigma}
\title{A Bound Helium Dianion Does Not Exist}
\author{Alec Shelley\thanks{Department of Applied Physics, Stanford University, Stanford, CA 94305, USA.}\ \ and\ \ Daniel M. Tartakovsky\thanks{Department of Energy Science and Engineering, Stanford University, Stanford, CA 94305, USA.}}
\date{October 2026}
\begin{document}\maketitle

\begin{abstract}
For a non-relativistic atom with $Z$ protons and $N$ bound electrons, Lieb proved in 1984 that $N<2Z+1$. Several improved bounds have been developed, but they all leave open the question of whether a bound helium dianion can exist. We prove that the answer is no. The ground-state energies of the four- and three-electron helium Hamiltonians coincide, $E_4(Z{=}2)=E_3(Z{=}2)$, so a helium nucleus binds at most three electrons. A weighted form of Lieb's identity reduces the problem to the nonnegativity of an explicit two-electron Schr\"odinger operator $m$, and Barta's inequality further reduces it to the pointwise nonnegativity of an explicit function of three variables. We prove the nonnegativity of this function by a combination of analytic estimates and a computer-assisted certification in interval arithmetic, and we publish the certificate with the paper.
\end{abstract}

\section{Introduction and main result}

Let $H_N(Z)$ be the nonrelativistic Hamiltonian of $N$ electrons and a fixed nucleus of charge $Z$, acting on antisymmetric functions of $N$ position--spin variables with two spin states, in atomic units:
\begin{equation}\label{eq:HN}
H_N(Z)=\sum_{i=1}^N\Big(-\tfrac12\Delta_i-\frac Z{r_i}\Big)+\sum_{1\le i<j\le N}\frac1{r_{ij}},\qquad r_i=\abs{x_i},\ r_{ij}=\abs{x_i-x_j}.
\end{equation}
 Let $E_N(Z)=\inf\spec H_N(Z)$. By the HVZ theorem \cite{RS4} $E_N\le E_{N-1}$, and $H_N$ has a bound state (an eigenvalue below its essential spectrum) if and only if $E_N<E_{N-1}$. The \emph{ionization problem} asks for the critical number $N_c(Z)$, the largest $N$ with $E_N<E_{N-1}$. Physically, no atom binds more than one extra electron \cite{BGS,Andersen}. Mathematically, a classic result of Lieb \cite{Lieb} shows that a bound state requires $N<2Z+1$, so $N_c(Z)\le2Z$; for hydrogen this gives $N_c(1)\le2$, which is sharp. Lieb's bound has since been improved to $N_c(Z)\le1.22Z+3Z^{1/3}$ by Nam \cite{Nam}, to $N_c(Z)\le1.5211Z+1+0.294Z^{1/3}$ for $Z\ge6$ by Benguria, Gonz\'alez-Brantes and Tubino \cite{BGT}, and, for $Z\ge4$, to $N_c(Z)\le1.1185Z+4Z^{1/3}$ by Hundertmark, Pattakos and Schulz \cite{HPS}. Recently an unrefereed, AI-generated manuscript released by OpenAI \cite{OAI} has claimed a proof of the ionization conjecture, that is that there exists a universal constant $C$ so that $N_c(Z)\le Z+C$, however their constant is not explicit. None of these results improves the $Z=2$ case, where it has been open for forty years whether or not bound dianionic helium exists. We prove that it does not:

\begin{theorem}\label{thm:main}
$E_4(2)=E_3(2)$. Consequently $N_c(2)\le3$: a helium nucleus binds at most three electrons.
\end{theorem}

Our proof is computer-assisted. Its mathematical part (Sections~\ref{sec:lieb}--\ref{sec:charts}) reduces the theorem to the nonnegativity of one explicit elementary function $q$ of three real variables on an explicit closed region, and establishes analytically the regularity needed to certify that nonnegativity by Taylor bounds. In Section~\ref{sec:cert} we prove that $q\ge0$ on that region in interval arithmetic with outward rounding. We follow the standards of computer-assisted proofs in analysis \cite{Hales,Tucker}. Rigorous computer-assisted lower bounds for atomic and ionic ground-state energies have recently been obtained by Hofacker \cite{Hofacker}; the question here is of a different kind, since no energy bound decides whether an eigenvalue exists. All floating-point operations that enter the proof are enclosed rigorously, the certificate is published, and the parts that cannot be certified by a computer are proved by hand.

\paragraph{Strategy.} Three ideas carry the proof.

\emph{(i) A reduction to a two-electron operator.} Lieb's argument multiplies the eigenvalue equation by $r_k$ and sums over $k$; the resulting identity has a pair term $\E\sum_{i<j}(r_i+r_j)/r_{ij}$ that is bounded below by $N-1$ pointwise, and $Z\ge N-1$ follows. With a one-body weight $w(r)$ growing linearly at infinity in place of $r$, the same identity acquires a kinetic term, and the pointwise bound on the pair term fails (a regular tetrahedron gives $1.837<2$, Section~\ref{sec:lieb}); but combined with the energy constraint through a Lagrange multiplier $\beta$, the identity becomes the statement that an explicit four-electron operator $\M$ is nonnegative. $\M$ contains one- and two-body terms only, so $\M=\sum_{i<j}m_{ij}$ for a two-electron operator $m$, and $m\ge0$ suffices. The choice of $w$ matters: with $w=r^2$ or with weights flat in the core, $m$ has a negative ``core-pair'' eigenvalue, because on $L^2(\R^6)$ nothing prevents both electrons of a pair from occupying the $1s$ shell. With $w(r)=r^4/(r^3+8)$ and $\beta=0.9$ the operator $m$ has a single bound state at $+0.0355$, below its thresholds $0.0468$ and $0.1022$.

\emph{(ii) A supersolution.} Barta's inequality (the ground-state transform) gives $m\ge\inf q$ for the local energy $q=m\varphi/\varphi$ of any positive $\varphi$. The exact ground state has $q\equiv0.0355$; an explicit $\varphi$ with exact nuclear and electron--electron cusps, $\sqrt r$ tails, three exponential terms and $400$ coefficients has $q\ge0.0102$ on $1.3\times10^6$ test points. Nonnegativity of $q$, a function of $(r_1,r_2,r_{12})$, is then the whole problem; no spectral computation is involved.

\emph{(iii) Certification.} The cusps cancel the Coulomb singularities of $q$ exactly, and we remove them analytically before evaluation. Barta's inequality is used only on the bounded region $\max(r_1,r_2)\le240$: outside it, $m$ is bounded below by elementary operator inequalities (one electron far from the nucleus sees the one-body ground-state energy plus the asymptotic pair term), and an IMS localization glues the two regions at the price of an explicit error $0.0145$, which the certification absorbs by proving $q\ge0.0145$ on $120\le\max(r_1,r_2)\le240$. On the bounded region, in three coordinate charts (interior, scaled coalescence, and a cone chart at the triple coalescence point) $q$ is smooth up to the boundary, and second-order Taylor bounds with interval-evaluated Hessians certify the required lower bound on boxes; branch-and-bound covers each chart.

Theorem~\ref{thm:main} excludes eigenvalues of $H_4(2)$ below the bottom $E_3(2)$ of its essential spectrum. It does not exclude eigenvalues embedded in the continuum and protected by symmetry (the four-electron analogues of the He$^-$ $1s2s2p\,{}^4P^o$ state), since the reduction of Theorem~\ref{thm:reduction} uses $E<E_3$; we return to this and to the other atoms in Section~\ref{sec:discussion}.

\paragraph{Organization.} Section~\ref{sec:lieb} proves the reduction theorem. Section~\ref{sec:m} studies the two-electron operator $m$ and fixes the parameters. Section~\ref{sec:barta} proves Barta's inequality in the form we need, with cusps. Section~\ref{sec:phi} defines $\varphi$ and derives the singularity-free form of its local energy. Section~\ref{sec:charts} proves the far-field theorem, which confines the pointwise problem to a bounded region, then defines the charts on that region and proves the regularity of $q$ in each. Section~\ref{sec:cert} describes the certified evaluation and states its result; Section~\ref{sec:proof} assembles the proof; Section~\ref{sec:discussion} discusses extensions. Appendix~\ref{app:coef} lists the basis and the coefficients of $\varphi$; Appendix~\ref{app:numerics} records the non-rigorous numerical computations that guided the choices (they are not part of the proof).

\section{Lieb's identity with a one-body weight, and the reduction}\label{sec:lieb}

Throughout, $\E f=\ip\Psi{f\Psi}$ for a normalized $\Psi$. Let $w:[0,\infty)\to[0,\infty)$ be $C^2$ with $w(r)\le C(1+r)$, $w'$ and $\Delta w:=w''+2w'/r$ bounded, and write $w_k=w(r_k)$.

\begin{lemma}[Weighted Lieb inequality]\label{lem:lieb}
Let $\Psi$ be a normalized eigenfunction of $H_N(Z)$ with eigenvalue $E<E_{N-1}(Z)$. Then
\begin{equation}\label{eq:lieb}
Z\,\E\sum_k\frac{w_k}{r_k}\;\ge\;\E\sum_{i<j}\frac{w_i+w_j}{r_{ij}}+\sum_k\Big[\tfrac12\int w_k\abs{\nabla_k\Psi}^2-\tfrac14\int(\Delta w)_k\abs\Psi^2\Big]+(E_{N-1}-E)\,\E\sum_kw_k.
\end{equation}
\end{lemma}
\begin{proof}
Since $E<E_{N-1}\le\inf\spec_{\rm ess}H_N$, $\Psi$ decays exponentially \cite{Agmon}, so $w_k\Psi$ and $w_k^{1/2}\Psi$ lie in the form domain of $H_N$ and $\ip{w_k\Psi}{(H_N-E)\Psi}=0$. Split $H_N=H^{(\hat k)}+h_k$ with $h_k=-\tfrac12\Delta_k-Z/r_k+\sum_{j\ne k}1/r_{kj}$ and $H^{(\hat k)}$ the Hamiltonian of the other $N-1$ electrons. As $w_k$ commutes with $H^{(\hat k)}$ and $w_k^{1/2}\Psi$ is antisymmetric in the other electrons, $\Re\ip{w_k\Psi}{H^{(\hat k)}\Psi}=\ip{w_k^{1/2}\Psi}{H^{(\hat k)}w_k^{1/2}\Psi}\ge E_{N-1}\int w_k\abs\Psi^2$. For the kinetic term, $\Re\ip{w\Psi}{-\Delta\Psi}=\int w\abs{\nabla\Psi}^2-\tfrac12\int\Delta w\abs\Psi^2$ by integration by parts, so $\Re\ip{w_k\Psi}{-\tfrac12\Delta_k\Psi}=\tfrac12\int w_k\abs{\nabla_k\Psi}^2-\tfrac14\int(\Delta w)_k\abs\Psi^2$. Collecting,
\[
0=\Re\ip{w_k\Psi}{(H_N-E)\Psi}\ \ge\ (E_{N-1}-E)\E w_k+\tfrac12\int w_k\abs{\nabla_k\Psi}^2-\tfrac14\int(\Delta w)_k\abs\Psi^2-Z\,\E\frac{w_k}{r_k}+\sum_{j\ne k}\E\frac{w_k}{r_{kj}} ,
\]
and summing over $k$ gives \eqref{eq:lieb}.
\end{proof}

With $w=r$ (or Lieb's $w=r/(1+\epsilon r)$, $\epsilon\to0$) the kinetic term is nonnegative \cite{Lieb}, and $(r_i+r_j)/r_{ij}\ge1$ pointwise, so \eqref{eq:lieb} gives $ZN=Z\,\E\sum_kr_k/r_k>\binom N2$, i.e.\ $N<2Z+1$ \cite{Lieb}. No pointwise bound of this type can reach $N-1\le Z$ for a dianion: for any weight, the ratio $\sum_{i<j}(w_i+w_j)/r_{ij}\big/\sum_kw_k/r_k$ equals $3\sqrt{3/8}=1.837<2$ on a regular tetrahedron centered at the nucleus, for $N=4$. The energy constraint is what closes the gap.

Fix $Z=2$, $N=4$, a number $\beta>0$ and a weight $w$ as above, and define, on antisymmetric four-electron functions,
\begin{equation}\label{eq:M}
\M=\sum_{k=1}^4\Big[-\tfrac12\nabla_k\!\cdot(\beta+w_k)\nabla_k-\frac{2(\beta+w_k)}{r_k}-\tfrac14(\Delta w)(r_k)\Big]+\sum_{i<j}\frac{\beta+w_i+w_j}{r_{ij}}+\beta\,\mathcal E ,
\end{equation}
where $\mathcal E>0$ is any number with $E_2(2)\le-\mathcal E$. The kinetic term is the quadratic form $\tfrac12\sum_k\int(\beta+w_k)\abs{\nabla_k\Psi}^2$, so $\M$ is the self-adjoint operator of a closed, bounded-below quadratic form on $\{\Psi:\sum_k\int(\beta+w_k)\abs{\nabla_k\Psi}^2<\infty\}$ (the Coulomb terms are form-bounded relative to the kinetic term with relative bound zero, by Hardy's inequality); only the quadratic form is used below.

\begin{theorem}[Reduction]\label{thm:reduction}
If $\M\ge0$ as a quadratic form on antisymmetric four-electron functions, then $E_4(2)=E_3(2)$.
\end{theorem}
\begin{proof}
Suppose $E_4<E_3$ and let $\Psi$ be a normalized ground state, so $\ip\Psi{H_4\Psi}=E_4$. Expanding \eqref{eq:M},
\begin{align*}
\ip\Psi{\M\Psi}=\beta\big(\ip\Psi{H_4\Psi}+\mathcal E\big)+\Big\{\sum_k\Big[\tfrac12\int w_k\abs{\nabla_k\Psi}^2&-\tfrac14\int(\Delta w)_k\abs\Psi^2\Big]\\
&+\E\sum_{i<j}\frac{w_i+w_j}{r_{ij}}-2\,\E\sum_k\frac{w_k}{r_k}\Big\}.
\end{align*}
By Lemma~\ref{lem:lieb} with $Z=2$, the curly bracket is $\le-(E_3-E_4)\E\sum_kw_k<0$. Since $E_4<E_3$ by hypothesis, $E_3\le E_2$ by HVZ and $E_2\le-\mathcal E$, $\beta(E_4+\mathcal E)<0$. Hence $\ip\Psi{\M\Psi}<0$, contradicting $\M\ge0$.
\end{proof}

We take $\mathcal E=2.9037$: the Rayleigh--Ritz value $-2.903724377$ of \cite{Nakashima} is an upper bound on $E_2(2)$, so $E_2\le-2.9037$. (Any cruder variational bound below $-2.9037$ would do, for instance Pekeris' \cite{Pekeris}.)

\begin{remark}
The multiplier can be split: since any eigenstate satisfies the virial relation $\langle2T+V\rangle=0$, one may replace $\beta(H_4+\mathcal E)$ by $\beta_TT+\beta_VV+(2\beta_V-\beta_T)\mathcal E$ with $\beta_T>0$, $2\beta_V\ge\beta_T$. Our final choice has $\beta_T=\beta_V$, and we do not use this freedom; it was useful in the search.
\end{remark}

\section{The two-electron operator}\label{sec:m}

Since $\M$ contains one- and two-body terms only, it decomposes exactly into six two-electron operators:
\begin{equation}\label{eq:pair}
\M=\sum_{i<j}m_{ij},\qquad m=\tfrac13(h_1+h_2)+\frac{\beta+w_1+w_2}{r_{12}}+\tfrac16\beta\mathcal E,\qquad h=-\tfrac12\nabla\!\cdot(\beta+w)\nabla-\frac{2(\beta+w)}{r}-\tfrac14\Delta w ,
\end{equation}
where $m_{ij}$ acts on the variables $x_i,x_j$. Each $m_{ij}$ is unitarily equivalent to $m$ on $L^2(\R^6)$, so $\M\ge6\inf\spec(m)$ and
\begin{equation}\label{eq:suff}
m\ge0\ \text{ on }L^2(\R^6)\quad\Longrightarrow\quad\M\ge0.
\end{equation}
No symmetry restriction is imposed on $m$: the decomposition is used only as an inequality.

From now on
\begin{equation}\label{eq:params}
\beta=\tfrac9{10},\qquad w(r)=\frac{r^4}{r^3+8},\qquad \mathcal E=2.9037,\qquad C:=\tfrac16\beta\mathcal E=0.435555.
\end{equation}
Then $w$ is smooth, $w(r)=r^4/8+O(r^7)$ at the origin, $w(r)=r-8/r^2+O(r^{-5})$ at infinity, $0\le w'\le\tfrac43$, and $\Delta w=\big(w''+2w'/r\big)$ is a bounded rational function, $\Delta w(0)=0$, $\Delta w(r)\sim2/r$. Writing $a(r)=\beta+w(r)$,
\begin{equation}\label{eq:m}
m=\sum_{k=1}^2\Big[-\tfrac16\nabla_k\!\cdot a(r_k)\nabla_k-\frac{2a(r_k)}{3r_k}-\tfrac1{12}(\Delta w)(r_k)\Big]+\frac{\beta+w(r_1)+w(r_2)}{r_{12}}+C.
\end{equation}

\paragraph{Why these parameters.} The thresholds of $m$ are explicit. With one electron in the ground state of $\tfrac13h$ (radial eigenvalue $\tfrac13e_0$) and the other at infinity, where $2a(r)/(3r)\to\tfrac23$ and $w_2/r_{12}\to1$, the energy is $t_a=\tfrac13e_0+\tfrac13+C$; with both at infinity it is at least $t_b=-\tfrac13+C$ (attained for antipodal electrons, where $(w_1+w_2)/r_{12}\to1$). Numerically $e_0=-2.16619$, $t_a=0.0468$, $t_b=0.1022$, and $m$ has exactly one eigenvalue below $t_a$, at $+0.0355$ (Appendix~\ref{app:numerics}). The parameters were chosen, among the family $w=r^{q+1}/(r^q+c^q)$, to maximize $\inf\spec(m)$ subject to both thresholds being comfortably positive; weights that stay flat in the core ($q\ge4$, $c\ge2.5$) produce a negative core-pair eigenvalue and are unusable. None of these numerical facts enters the proof; they explain why the supersolution of Section~\ref{sec:phi} can exist.

\begin{remark}[Perron--Frobenius]
$m$ is a Schr\"odinger operator with a scalar, strictly positive diffusion coefficient; its semigroup is positivity improving, so if $\inf\spec(m)<\inf\spec_{\rm ess}(m)$ the ground state is unique, positive, and therefore rotation-invariant and exchange-symmetric. This is why the numerics of Appendix~\ref{app:numerics} need only the $L=0$ symmetric sector, and why a supersolution of the symmetric form \eqref{eq:phi} is natural. The proof below does not use it.
\end{remark}

\section{Barta's inequality with cusps}\label{sec:barta}

Let $\Sigma=\{x\in\R^6:\ r_1r_2r_{12}=0\}$, a closed set of codimension three. Let $\varphi>0$ be continuous on $\R^6$, $C^2$ on $\R^6\setminus\Sigma$, locally Lipschitz on $\R^6$, and such that $q:=m\varphi/\varphi$ (computed on $\R^6\setminus\Sigma$) is locally bounded on $\R^6$.

\begin{lemma}[Ground-state transform]\label{lem:barta}
For every $\psi\in C_c^\infty(\R^6)$, writing $\psi=\varphi v$,
\begin{equation}\label{eq:gst}
\ip\psi{m\psi}=\tfrac16\sum_{k=1}^2\int a(r_k)\,\varphi^2\abs{\nabla_kv}^2+\int q\,\abs\psi^2.
\end{equation}
Consequently $m\ge\inf_{\R^6\setminus\Sigma}q$.
\end{lemma}
\begin{proof}
Let $\Sigma_\epsilon$ be the $\epsilon$-neighbourhood of $\Sigma$ and $\Omega_\epsilon=\R^6\setminus\Sigma_\epsilon$. On $\Omega_\epsilon$, $\varphi\in C^2$ and the pointwise identity $a\abs{\nabla_k(\varphi v)}^2=a\varphi^2\abs{\nabla_kv}^2+a\,\nabla_k(\varphi^2)\cdot v\nabla_kv+a v^2\abs{\nabla_k\varphi}^2=a\varphi^2\abs{\nabla_kv}^2+\tfrac12a\nabla_k(\varphi^2)\cdot\nabla_k(v^2)+av^2\abs{\nabla_k\varphi}^2$ holds; integrating by parts the middle term over $\Omega_\epsilon$,
\[
\tfrac16\sum_k\int_{\Omega_\epsilon}a\abs{\nabla_k\psi}^2=\tfrac16\sum_k\int_{\Omega_\epsilon}a\varphi^2\abs{\nabla_kv}^2+\int_{\Omega_\epsilon}\frac{-\tfrac16\sum_k\nabla_k\!\cdot(a\nabla_k\varphi)}{\varphi}\abs\psi^2+\tfrac1{12}\sum_k\int_{\partial\Omega_\epsilon}a\,v^2\,\nabla_k(\varphi^2)\cdot n.
\]
The boundary term is bounded by $C\epsilon^2\sup\abs{\nabla\varphi^2}\sup\abs v^2\to0$ since $\varphi$ is Lipschitz and $\partial\Omega_\epsilon\cap\operatorname{supp}\psi$ has area $O(\epsilon^2)$. Adding the potential terms, which are locally integrable (Hardy), and letting $\epsilon\to0$ (dominated convergence, using that $q$ is locally bounded and $\abs{\nabla\varphi}$ is bounded on compacts) gives \eqref{eq:gst} for the quadratic form of $m$ on $C_c^\infty$, which is a form core. The first term in \eqref{eq:gst} is nonnegative, so $\ip\psi{m\psi}\ge\inf q\,\norm\psi^2$.
\end{proof}

The hypotheses are exactly what cusps provide: a function of the form $e^{-S}$ with $S$ Lipschitz, smooth off $\Sigma$, has bounded $q$ near $\Sigma$ if and only if the Coulomb singularities of $m$ are cancelled by the $1/r$ singularities of $\Delta S$, which is the Kato cusp condition. For \eqref{eq:m} these conditions are
\begin{equation}\label{eq:cusps}
\partial_{r_k}S\to2\ \ (r_k\to0),\qquad \partial_{r_{12}}S\to-\kappa(r)\ \ (r_{12}\to0,\ r_1=r_2=r),\qquad \kappa(r)=\frac{3\,(\beta+2w(r))}{2\,(\beta+w(r))} ,
\end{equation}
the first because $-\tfrac16a\cdot\big(-2\partial_{r}S/r\big)$ must cancel $-\tfrac23a/r$, the second because $\tfrac16(a_1+a_2)\cdot\big(-2\partial_uS/u\big)$ must cancel $(\beta+2w)/u$. Note $\kappa(0)=\tfrac32$, $\kappa(\infty)=3$.

\section{The supersolution and its local energy}\label{sec:phi}

\subsection{Definition}
Set $u=r_{12}$, $t=\cos\theta_{12}=\hat x_1\cdot\hat x_2$, $\rho=\sqrt{(r_1^2+r_2^2)/2}$ and $b=2$. Let $f_0,\dots,f_{21}$ be the radial functions
\begin{gather*}
f_0=\frac{r^2}{1+r},\quad f_1=\frac{r^2}{(1+r)^2},\quad f_2=\frac{r^2}{(1+r)^3},\quad f_3=\frac{r^2}{4+r},\quad f_4=\frac{r^2}{(4+r)^2},\\
f_5=\sqrt{1+r^2}-1,\quad f_6=(1+r^2)^{1/4}-1,\quad f_7=(1+r^2/16)^{1/4}-1,\quad f_8=\frac{r^2}{(0.3+r)^2},\\
f_9=r^2e^{-r},\quad f_{10}=\frac{r^2}{(2+r)^2},\quad f_{11}=\frac{r^2}{1+r^2},\quad f_{12}=\frac{r^2}{(\tfrac12+r)^2},\quad f_{13}=\frac{r^2}{(2+r)^3},\\
f_{14}=\frac{r^2}{(2+r)^4},\quad f_{15}=\frac{r^2}{(6+r)^2},\quad f_{16}=\frac{r^2}{(6+r)^3},\quad f_{17}=r^2e^{-2r},\\
f_{18}=r^2e^{-r/2},\quad f_{19}=\frac{r^2}{(1+r^2)^2},\quad f_{20}=\frac{r^4}{(1+r^2)^3},\quad f_{21}=(1+r^2/4)^{-1/2}-1,
\end{gather*}
all with $f_j(0)=f_j'(0)=0$; $f_0,f_3,f_5$ grow linearly, $f_6,f_7$ like $\sqrt r$, the others are bounded. Let $\mathfrak B=\{1,2,4,8,10,11,12,13,15,17,18\}$. The basis $B_1,\dots,B_{200}$ consists of the \emph{even} functions
\begin{gather*}
P_j=f_j(r_1)+f_j(r_2)\ (0\le j\le21),\qquad PP_{jk}=f_j(r_1)f_k(r_2)+f_k(r_1)f_j(r_2)\ (j\le k\in\mathfrak B),\\
U_k=\frac{u^2}{(1+u)^k}\ (k=2,3,4),
\end{gather*}
\begin{gather*}
UP_j=\frac{u^2}{(1+u)^2}P_j\ (j\in\{1,2,4,11,12,13,15\}),\qquad VP_j=\frac{u^2}{(1+u)^3}P_j\ (j\in\{1,2,4\}),\\
UPP_j=\frac{u^2}{(1+u)^2}f_j(r_1)f_j(r_2)\ (j\in\{1,2,4\}),
\end{gather*}
\begin{gather*}
UPP_{jk}=\frac{u^2}{(1+u)^2}PP_{jk}\ ((j,k)\in\{(1,2),(1,4),(2,4)\}),\\
D_1=\frac{x_1\!\cdot\!x_2}{(1+r_1^2)(1+r_2^2)},\qquad D_2=\frac{D_1}{1+r_1+r_2},
\end{gather*}
and the \emph{odd} functions (antisymmetric under $x_1\leftrightarrow x_2$)
\begin{gather*}
Q_j=f_j(r_1)-f_j(r_2),\qquad QQ_{jk}=f_j(r_1)f_k(r_2)-f_k(r_1)f_j(r_2)\ (j<k\in\mathfrak B),\\
UQ_j=\frac{u^2}{(1+u)^2}Q_j,\qquad VQ_j=\frac{u^2}{(1+u)^3}Q_j,
\end{gather*}
\[
UQQ_{jk}=\frac{u^2}{(1+u)^2}QQ_{jk},\quad D_3=\frac{(x_1\!\cdot\!x_2)(r_1^2-r_2^2)}{(1+r_1^2)^2(1+r_2^2)^2},
\]
with the same index sets, and $x_1\cdot x_2=(r_1^2+r_2^2-u^2)/2$. The supersolution is
\begin{equation}\label{eq:phi}
\varphi=e^{-S^{(1)}}+e^{-S^{(2)}}+e^{-S^{(3)}},\qquad S^{(\tau)}=2(r_1+r_2)-\frac{\kappa(\rho)\,u}{1+bu}+\sum_{k=1}^{200}c^{(\tau)}_kB_k ,
\end{equation}
where $c^{(2)}_k=\pm c^{(1)}_k$ with the minus sign on odd functions (so $e^{-S^{(2)}}$ is the exchange of $e^{-S^{(1)}}$ and $\varphi$ is exchange-symmetric), $c^{(3)}_k=0$ for odd $k$, and the two independent coefficient vectors $c^{(1)},c^{(3)}\in\R^{200}$ are listed in Appendix~\ref{app:coef}. They satisfy the \emph{tail constraints}
\begin{equation}\label{eq:tails}
c^{(\tau)}_{P_0}+c^{(\tau)}_{P_3}+c^{(\tau)}_{P_5}=-2\ \ (\tau=1,3),\qquad c^{(1)}_{Q_0}+c^{(1)}_{Q_3}+c^{(1)}_{Q_5}=0 ,
\end{equation}
which remove the linear growth of every $S^{(\tau)}$ in both electrons, leaving $\sqrt r$ tails. (Why this is essential: if one electron of a term decays like $e^{-sr}$ while the other decays like $e^{-\alpha\sqrt r}$, that term still carries half the weight of $\varphi$ at the crossover radius $(\alpha/s)^2$, where its local energy is about $-\alpha^2/6$; the effect is invisible below $r\sim10^3$.)

\begin{lemma}\label{lem:phiprops}
$\varphi$ is continuous and positive on $\R^6$, real-analytic on $\R^6\setminus\Sigma$ except for the factors $r_k^2e^{-r_k}$, which are $C^\infty$, locally Lipschitz, and satisfies the cusp conditions \eqref{eq:cusps} for each term. Its local energy $q=m\varphi/\varphi$ is locally bounded on $\R^6$.
\end{lemma}
\begin{proof}
In the variables $(r_1,r_2,u)$, every $B_k$ and $\kappa(\rho)u/(1+bu)$ has $\partial_{r_1}B_k\to0$ as $r_1\to0$: the radial functions satisfy $f_j(0)=f_j'(0)=0$, the $u$-factors do not depend on $r_1$, $\rho$ depends on $r_1$ through $r_1^2$, and the dipole functions are $D=X\,g(r_1,r_2)$ with $X=x_1\!\cdot\!x_2=(r_1^2+r_2^2-u^2)/2$, so that $\partial_{r_1}D=r_1g+Xg_{r_1}$ with $X=r_1r_2\cos\theta_{12}\to0$ (this matters for $D_2$, whose factor $1/(1+r_1+r_2)$ has nonzero slope at $r_1=0$). Hence $\partial_{r_1}S^{(\tau)}\to2$; likewise at $r_2=0$. Near $u=0$ every $B_k$ is $O(u^2)$ or smooth in $x_1\cdot x_2$, so $\partial_uS^{(\tau)}\to-\kappa(\rho)=-\kappa(r)$ on $r_1=r_2=r$. Boundedness of $q$ then follows from the explicit regularized form \eqref{eq:qreg} below, whose every term is locally bounded; the weighted mean of bounded terms is bounded.
\end{proof}

\subsection{The local energy in Hylleraas coordinates}
For a single term $e^{-S}$ with $S=S(r_1,r_2,u)$, subscripts denoting partial derivatives, and
\begin{equation}\label{eq:cos}
\cos_1:=\hat x_1\cdot\frac{x_1-x_2}{u}=\frac{r_1-r_2t}{u}=\frac{r_1^2+u^2-r_2^2}{2r_1u},\qquad \cos_2:=\hat x_2\cdot\frac{x_2-x_1}{u}=\frac{r_2-r_1t}{u},
\end{equation}
one has $\abs{\nabla_1S}^2=S_1^2+S_u^2+2S_1S_u\cos_1$, $\Delta_1S=S_{11}+2S_1/r_1+S_{uu}+2S_u/u+2S_{1u}\cos_1$, $\hat x_1\cdot\nabla_1S=S_1+S_u\cos_1$, and from $-\nabla\cdot(a\nabla e^{-S})/e^{-S}=a(\abs{\nabla S}^2-\Delta S)-\nabla a\cdot\nabla S$,
\begin{equation}\label{eq:q}
q_S=C-\frac{\Delta w_1+\Delta w_2}{12}+\sum_{k=1}^2\Big[-\frac{a_k}6\big(\abs{\nabla_kS}^2-\Delta_kS\big)+\frac{w'_k}6\,\hat x_k\!\cdot\!\nabla_kS-\frac{2a_k}{3r_k}\Big]+\frac{\beta+w_1+w_2}{u}.
\end{equation}
For $\varphi=\sum_\tau e^{-S^{(\tau)}}$ the local energy is the weighted mean
\begin{equation}\label{eq:qmean}
q=\frac{\sum_\tau e^{-S^{(\tau)}}q_{S^{(\tau)}}}{\sum_\tau e^{-S^{(\tau)}}}.
\end{equation}

\subsection{Removing the singularities}
The terms of \eqref{eq:q} with $1/r_k$ and $1/u$ are individually singular; the cusps make their sums regular, and for a certified evaluation the cancellation must be performed symbolically. Write $J=-\kappa(\rho)u/(1+bu)$ and $\Sigma_c=\sum_kc_kB_k$ (for the term under consideration).

\begin{proposition}[Regularized local energy]\label{prop:qreg}
For each term, with
\[
E_k:=\frac{S_k-2}{r_k}=\frac{\partial_{r_k}(J+\Sigma_c)}{r_k},\qquad E_u^{b}:=\frac{\partial_u\Sigma_c}{u},\qquad g:=\frac{\beta\,\big(w_1+w_2-2w(\rho)\big)}{2\,(\beta+w(\rho))},
\]
\begin{align}
q_S\;=\;&C-\frac{\Delta w_1+\Delta w_2}{12}+\sum_{k=1}^2\Big[-\frac{a_k}6\big(g_k-l_k'\big)+\frac{w_k'}6d_k\Big]+\frac{a_1}3E_1+\frac{a_2}3E_2\notag\\
&+\frac{a_1+a_2}3\Big[\frac{\kappa(\rho)\,b\,(2+bu)}{(1+bu)^2}+E^b_u\Big]+\frac gu ,\label{eq:qreg}
\end{align}
where $g_k=S_k^2+S_u^2+2S_kS_u\cos_k$, $l_k'=S_{kk}+S_{uu}+2S_{ku}\cos_k$, $d_k=S_k+S_u\cos_k$. Each of $E_1,E_2,E^b_u$ is a finite sum of explicit functions that are real-analytic (resp.\ $C^\infty$ for the $e^{-r}$ terms) in the variables $(r_1,r_2,\cos\theta_{12})$ on $r_k\ge0$ including $r_k=0$, and in $u\ge0$ including $u=0$ (for the dipole terms $D=Xg$ one has $\partial_{r_1}D/r_1=g+r_2\cos\theta_{12}\,g_{r_1}$, which is regular in $(r_1,r_2,\cos\theta_{12})$ but not in $(r_1,r_2,u)$; for $J$, the quotients $u/\rho$ and $r_k/\rho$ that appear are bounded on the region but not continuous at the triple point $r_1=r_2=u=0$, and are regular functions of the chart coordinates of Section~\ref{sec:charts}); and $g$ vanishes to second order on $r_1=r_2$, with
\begin{equation}\label{eq:gbound}
\Big|\frac gu\Big|\le\frac{\beta\,M\,u}{16\,(\beta+w(\rho))},\qquad M=4\sup_{I}\abs{w''}+2\sup_I\Big|\frac{w'(x)}x\Big|,\quad I=[\min(r_1,r_2),\max(r_1,r_2)].
\end{equation}
\end{proposition}
\begin{proof}
The $1/r_k$ terms of \eqref{eq:q} are $\tfrac{a_k}3S_k/r_k-\tfrac23a_k/r_k=\tfrac{a_k}3(S_k-2)/r_k$, and $S_k-2=\partial_{r_k}(J+\Sigma_c)$; since every $B_k$ and $J$ depends on $r_k$ only through $r_k^2$ near $r_k=0$ (Lemma~\ref{lem:phiprops}), $\partial_{r_k}(\cdot)/r_k$ is regular, and the explicit expressions are obtained by dividing each term's derivative by $r_k$ and cancelling. The $1/u$ terms are $\tfrac{a_1+a_2}3S_u/u+(\beta+w_1+w_2)/u$; with $S_u=-\kappa(\rho)/(1+bu)^2+\partial_u\Sigma_c$,
\[
\frac{a_1+a_2}3\cdot\frac{S_u+\kappa(\rho)}u=\frac{a_1+a_2}3\Big[\kappa(\rho)\frac{1-(1+bu)^{-2}}u+E_u^b\Big]=\frac{a_1+a_2}3\Big[\frac{\kappa(\rho)b(2+bu)}{(1+bu)^2}+E_u^b\Big],
\]
and the remainder is $\big[\beta+w_1+w_2-\tfrac13(a_1+a_2)\kappa(\rho)\big]/u=g/u$, because with $a_k=\beta+w_k$ and $\kappa=3(\beta+2\bar w)/(2(\beta+\bar w))$, $\bar w=w(\rho)$,
\begin{align*}
\beta+w_1+w_2-\frac{(2\beta+w_1+w_2)(\beta+2\bar w)}{2(\beta+\bar w)}&=\frac{2(\beta+w_1+w_2)(\beta+\bar w)-(2\beta+w_1+w_2)(\beta+2\bar w)}{2(\beta+\bar w)}\\
&=\frac{\beta(w_1+w_2-2\bar w)}{2(\beta+\bar w)}.
\end{align*}
For \eqref{eq:gbound}: put $r=(r_1+r_2)/2$, $d=(r_1-r_2)/2$, $F(d)=w(r+d)+w(r-d)-2w(\sqrt{r^2+d^2})$. $F$ is even, $F(0)=0$, so $\abs{F(d)}\le\tfrac12d^2\sup_{[0,d]}\abs{F''}$, and $F''(d)=w''(r+d)+w''(r-d)-2\big[w''(\rho_d)d^2/\rho_d^2+w'(\rho_d)r^2/\rho_d^3\big]$ with $\rho_d=\sqrt{r^2+d^2}\ge r$, so $\abs{F''}\le M$ on $I$. Since $\abs{r_1-r_2}\le u$, $\abs g\le\beta u^2M/(16(\beta+\bar w))$.
\end{proof}
Away from coalescence $g/u$ is evaluated directly; the bound \eqref{eq:gbound} is used on boxes that approach $u=0$.

\section{The far field, and charts on a bounded region}\label{sec:charts}

Since $q$ depends on $(x_1,x_2)$ only through $(r_1,r_2,u)$ and is symmetric under $r_1\leftrightarrow r_2$, every pointwise statement about $q$ is made on the closed region $\{0\le r_1\le r_2,\ \abs{r_1-r_2}\le u\le r_1+r_2\}$. Barta's inequality is applied only on the bounded set $\{\rho_{\max}\le R'\}$, $\rho_{\max}:=\max(r_1,r_2)$; outside it $m$ is bounded below by elementary operator inequalities, and the two are glued by an IMS localization. Throughout,
\[
R=120,\qquad R'=2R=240.
\]


\subsection{The far field}
Write $m=\tfrac13h_1+\tfrac13h_2+V_{12}+C$ with $V_{12}=(\beta+w_1+w_2)/r_{12}$, and $\tfrac13h_k=\tfrac13T_k+p_k$, $T_k=-\tfrac12\nabla_k\!\cdot a\nabla_k\ge0$, $p_k=-\tfrac23a_k/r_k-\tfrac1{12}\Delta w_k$ (multiplication operators). We use the following elementary properties of $w(r)=r^4/(r^3+8)$: $w\le r$, so $a(r)/r\le1+\beta/r$; $r-w=8r/(r^3+8)\le1.0583$ (maximum at $r=4^{1/3}$) and $r-w\le8/r^2$; and $\Delta w=2/r-16/r^4+O(r^{-7})$ with $r\,\Delta w(r)\le2$ for $r\ge5$. Finally, the one-body operator $\tfrac13h$ on $L^2(\R^3)$ satisfies $\tfrac13h\ge\tfrac13e_0^{\rm rig}:=-0.7291$ (Remark~\ref{rem:checks}(iv): Barta's inequality with a certified radial supersolution; no symmetry assumption is needed, since the ground-state transform holds for every positive supersolution), hence $\tfrac13h_1\ge-0.7291$ on $L^2(\R^6)$.

\begin{lemma}[One electron far]\label{lem:far}
Let $t_a^{\rm rig}:=\tfrac13e_0^{\rm rig}+\tfrac13+C\ge0.0398$. For $\varrho\ge60$ and every $\psi$ in the form domain of $m$ supported in $\{r_2\ge\varrho\}$ (any $r_1$),
\[
\ip\psi{m\psi}\ \ge\ \big(t_a^{\rm rig}-0.93/\varrho\big)\norm\psi^2 ,
\]
and symmetrically for $\psi$ supported in $\{r_1\ge\varrho\}$. In particular the constant is $\ge0.0243$ for $\varrho=60$.
\end{lemma}
\begin{proof}
Dropping $T_2\ge0$ and using $\tfrac13h_1\ge\tfrac13e_0^{\rm rig}$, $m\ge\tfrac13e_0^{\rm rig}+C+p_2+V_{12}$ as quadratic forms. On $\{r_2\ge\varrho\}$: $-\tfrac23a_2/r_2\ge-\tfrac23-0.6/r_2$, $-\tfrac1{12}\Delta w_2\ge-0.167/r_2$, and, since $r_{12}\le r_1+r_2$,
\begin{align*}
V_{12}\ \ge\ \frac{\beta+w_1+w_2}{r_1+r_2}\ &=\ 1-\frac{(r_1-w_1)+(r_2-w_2)-\beta}{r_1+r_2}\\
&\ge\ 1-\frac{1.0583+8/r_2^2-0.9}{r_2}\ \ge\ 1-\frac{0.161}{r_2}\qquad(r_2\ge60),
\end{align*}
for every $r_1\ge0$. Hence $p_2+V_{12}\ge\tfrac13-0.93/\varrho$ on the support of $\psi$.
\end{proof}

\paragraph{The partition of unity.} Let $\lambda=\ln(r_1/r_2)$. Let $A(\rho)$ be the ramp equal to $0$ for $\rho\le R$, to $1$ for $\rho\ge R'$ and linear in $\ln\rho$ in between (slope $1/\ln2$), and $D(\lambda)$ the ramp equal to $0$ for $\lambda\le-\ln2$, to $1$ for $\lambda\ge\ln2$ and linear in between (slope $1/\ln4$). Define the Lipschitz functions
\begin{gather*}
\chi_N=\cos\big(\tfrac\pi2A(\rho_{\max})\big),\qquad \chi_2=\sin\big(\tfrac\pi2A(\rho_{\max})\big)\cos\big(\tfrac\pi2D(\lambda)\big),\\
\chi_1=\sin\big(\tfrac\pi2A(\rho_{\max})\big)\sin\big(\tfrac\pi2D(\lambda)\big),
\end{gather*}
so that $\chi_N^2+\chi_1^2+\chi_2^2=1$, $\operatorname{supp}\chi_N\subset\{\rho_{\max}\le R'\}$, $\operatorname{supp}\chi_2\subset\{\rho_{\max}\ge R,\ r_1\le2r_2\}\subset\{r_2\ge R/2\}$ and $\operatorname{supp}\chi_1\subset\{r_1\ge R/2\}$.

\begin{theorem}[Far field]\label{thm:far}
Suppose that the local energy $q$ of a function $\varphi$ satisfying the hypotheses of Lemma~\ref{lem:barta} obeys
\[
q\ge0\ \text{ on }\{\rho_{\max}\le R\}\setminus\Sigma,\qquad q\ge0.0145\ \text{ on }\{R\le\rho_{\max}\le R'\}\setminus\Sigma.
\]
Then $m\ge0$ on $L^2(\R^6)$.
\end{theorem}
\begin{proof}
For Lipschitz cutoffs $\chi_i$ with $\sum_i\chi_i^2=1$ and the divergence-form kinetic term, the IMS identity holds for every $\Psi$ in the form domain:
\[
\ip\Psi{m\Psi}=\sum_{i\in\{N,1,2\}}\ip{\chi_i\Psi}{m\chi_i\Psi}-\ip\Psi{W\Psi},\qquad W=\tfrac16\sum_{k=1,2}a(r_k)\sum_i\abs{\nabla_k\chi_i}^2
\]
(the cross terms cancel because $\sum_i\chi_i^2=1$; no smoothness beyond Lipschitz is needed). We bound $W$. A direct computation gives $\sum_i\abs{\nabla_k\chi_i}^2=(\tfrac\pi2)^2\big[\abs{\nabla_kA}^2+\sin^2(\tfrac\pi2A)\abs{\nabla_kD}^2\big]$ (the mixed terms cancel). Since $\abs{\nabla_k\rho_{\max}}\le1$ a.e.\ and only the electron with $r_k=\rho_{\max}$ contributes, the $A$-ramp contributes at most $\tfrac16a(\rho_{\max})(\tfrac\pi{2\ln2})^2\rho_{\max}^{-2}\le\tfrac16(1+\beta/R)\,5.136/R=0.0072$, for $\rho_{\max}\in[R,R']$ and zero elsewhere. The $D$-ramp has $\abs{\nabla_k\lambda}=1/r_k$, so it contributes at most $\tfrac16(\tfrac\pi{2\ln4})^2\sum_ka(r_k)/r_k^2\le0.214\sum_k(1/r_k+0.9/r_k^2)$, and on its support (where $\abs\lambda\le\ln2$ and $\rho_{\max}\ge R$) both $r_k\ge R/2=60$, giving at most $0.0073$. Hence $W$ is supported in $\{\rho_{\max}\ge R\}$ and $\sup W\le0.0145$.

For $\chi_N\Psi$, supported in $\{\rho_{\max}\le R'\}$, Lemma~\ref{lem:barta} gives $\ip{\chi_N\Psi}{m\chi_N\Psi}\ge\int q\,\chi_N^2\abs\Psi^2$. For $\chi_2\Psi$, supported in $\{r_2\ge60\}$, Lemma~\ref{lem:far} gives $\ip{\chi_2\Psi}{m\chi_2\Psi}\ge0.0243\norm{\chi_2\Psi}^2$, and likewise for $\chi_1\Psi$. Inserting into the IMS identity,
\[
\ip\Psi{m\Psi}\ \ge\ \int\Big[q\chi_N^2+0.0243\,(\chi_1^2+\chi_2^2)-W\Big]\abs\Psi^2\ \ge\ 0 ,
\]
since the bracket equals $q\ge0$ where $\rho_{\max}<R$ (there $W=0$ and $\chi_1=\chi_2=0$), is at least $\min(q,0.0243)-W\ge0.0145-0.0145=0$ where $R\le\rho_{\max}\le R'$, and equals $0.0243-W\ge0.0243-0.0073>0$ where $\rho_{\max}>R'$ ($\chi_N=0$ there, and only the $D$-ramp contributes to $W$).
\end{proof}

The choice $R=100$, $R'=200$ also works ($t_a^{\rm rig}-0.93/50=0.0212$ against $\sup W\le0.0179$) with a thinner margin. Note that the theorem uses nothing about the far field of $\varphi$: the tail constraints \eqref{eq:tails} are still what keeps $q$ bounded below at large $r$ on the bounded region, but no limit $r\to\infty$ of $q$ enters the proof, and the near-tie of the $\sqrt r$ coefficients of the three terms is irrelevant.

\subsection{Charts on the bounded region}
Fix $\nu_0=0.07$ and $\epsilon_0=2\times10^{-4}$. Write $x_1=\bar x+\tfrac u2n$, $x_2=\bar x-\tfrac u2n$ with $\abs n=1$, $R_c=\abs{\bar x}$, $\chi=\hat{\bar x}\cdot n\in[-1,1]$, and let $\nu=u/R_c$ be the \emph{scaled separation}. On the region $\{r_1\le r_2\le R'\}$ three charts are used; the coalescence set $u=0$ lies in the second only.

\paragraph{Interior chart $\mathcal C_{\rm int}$: $(r_1,r_2,t)$.} Domain $0\le r_1\le r_2\le R'$, $r_2\ge\epsilon_0$, $-1\le t\le1$, restricted to $\nu\ge\nu_0$. Here $u=\sqrt{(r_1-r_2)^2+2r_1r_2(1-t)}$, $\cos_1=(r_1-r_2t)/u$, $\cos_2=(r_2-r_1t)/u$, and
\begin{equation}\label{eq:chain}
\begin{gathered}
\frac{\partial u}{\partial r_1}=\cos_1,\qquad\frac{\partial u}{\partial r_2}=\cos_2,\qquad\frac{\partial u}{\partial t}=-\frac{r_1r_2}u=:u_t,\\
\frac{\partial\cos_1}{\partial r_1}=\frac{1-\cos_1^2}u,\qquad\frac{\partial\cos_1}{\partial r_2}=\frac{-t-\cos_1\cos_2}u,\qquad\frac{\partial\cos_1}{\partial t}=\frac{-r_2-\cos_1u_t}{u},
\end{gathered}
\end{equation}
and symmetrically for $\cos_2$; $\partial_tu_t=-r_1^2r_2^2/u^3$. On the domain $u>0$: if $u=0$ then $x_1=x_2$, so $R_c=r_2\ge\epsilon_0>0$ and $\nu=0<\nu_0$. In this chart $\nu=2u/\sqrt{(r_1-r_2)^2+2r_1r_2(1+t)}$, an expression without cancellation that is used to decide, box by box, whether a box lies in the domain.

\paragraph{Scaled coalescence chart $\mathcal C_{\rm coal}$: $(R_c,\nu,\chi)$.} Domain $0\le R_c\le R'$, $0\le\nu\le\nu_0$, $-1\le\chi\le1$. Then
\begin{equation}\label{eq:coal}
\begin{gathered}
r_1=R_c\sqrt{1+\tfrac{\nu^2}4+\nu\chi},\qquad r_2=R_c\sqrt{1+\tfrac{\nu^2}4-\nu\chi},\qquad u=R_c\nu,\\
\cos_1=\frac{2\chi+\nu}{2\sqrt{1+\nu^2/4+\nu\chi}},\qquad \cos_2=\frac{\nu-2\chi}{2\sqrt{1+\nu^2/4-\nu\chi}},
\end{gathered}
\end{equation}
and $\rho=R_c\sqrt{1+\nu^2/4}$, $x_1\!\cdot\!x_2=R_c^2(1-\nu^2/4)$, $\cos\theta_{12}=(1-\nu^2/4)/\sqrt{(1+\nu^2/4)^2-\nu^2\chi^2}$. The radicands are $\ge(1-\nu_0/2)^2>0$, so all of these are real-analytic on the closed domain, including $\nu=0$ (where $r_1=r_2=R_c$) and $R_c=0$ (the triple point). The chart is not injective at $R_c=0$ and not symmetric under $r_1\leftrightarrow r_2$ (which is $\chi\mapsto-\chi$); both are harmless, since the certification covers the whole domain. The domain contains the set $\{\nu\le\nu_0,\ \rho_{\max}\le R'\}$, and in particular every point with $t=\cos\theta_{12}>0.999$ and $r_1=r_2$ (there $u^2=2r^2(1-t)<0.002r^2$, $\nu<0.045$).

\paragraph{Cone chart $\mathcal C_{\rm cone}$: $(R,s,t)$, $R=r_2$, $s=r_1/r_2$.} Domain $0\le R\le\epsilon_0$, $0\le s\le1$, $-1\le t\le1$, restricted to $\nu\ge\nu_0$. Here $r_1=sR$, $r_2=R$, $u=R\,d(s,t)$ with $d=\sqrt{(1-s)^2+2s(1-t)}$, $\cos_1=(s-t)/d$, $\cos_2=(1-st)/d$, $\rho=R\sqrt{(1+s^2)/2}$, $\cos\theta_{12}=t$, and $\nu=2d/\sqrt{(1-s)^2+2s(1+t)}$ (independent of $R$). The direction cosines are singular only where $d=0$, i.e.\ on the coalescence line $s=t=1$, which $\nu\ge\nu_0$ excludes: $\nu\ge\nu_0$ forces $d\ge\tfrac{\nu_0}2\sqrt{(1-s)^2+2s(1+t)}$, which vanishes only at $s=1,t=-1$ where $d=2$. This chart is needed because the triple point $r_1=r_2=u=0$ is not a regular point of $(r_1,r_2,t)$ (there $\rho$, $u$ and the quotients $r_k/\rho$, $u/\rho$ are not smooth functions of $(r_1,r_2,t)$), whereas in $(R,s,t)$ every one of them is $R$ times an analytic function of $(s,t)$.

\begin{lemma}[Regularity in the charts]\label{lem:reg}
In each of the three charts, $q$ (defined by \eqref{eq:qmean}, \eqref{eq:qreg}) extends to a $C^\infty$ function on the closed chart domain, real-analytic except for the terms containing $e^{-r_k}$, which are $C^\infty$ with all derivatives bounded.
\end{lemma}
\begin{proof}
By Proposition~\ref{prop:qreg}, $q_S$ is a rational-exponential expression in the quantities $S_k$, $S_u$, $S_{kk}$, $S_{uu}$, $S_{ku}$ (regular wherever the basis functions are), $E_1$, $E_2$, $E_u^b$, the cusp term $\kappa(\rho)b(2+bu)/(1+bu)^2$, $g/u$, $\cos_k$, $w_k,w_k',\Delta w_k$, $a_k$, and the weights $e^{-(S^{(\tau)}-S^{(\tau')})}$ are smooth wherever the $S^{(\tau)}$ are. The radial functions $f_j$, $w$, $w'$, $\Delta w$, $\kappa$ are real-analytic on $[0,\infty)$ (rational, or algebraic with positive radicands; $e^{-r}$, $e^{-2r}$, $e^{-r/2}$ are entire), and $f_j(r)=r^2\tilde f_j(r)$ with $\tilde f_j$ analytic.

\emph{$\mathcal C_{\rm int}$.} On the domain $u>0$, so $u$, $\cos_k$ and the chain rule \eqref{eq:chain} are real-analytic in $(r_1,r_2,t)$, including at $r_1=0$ (where $u=r_2$). The quotients $E_k=\partial_{r_k}(J+\Sigma_c)/r_k$ are analytic in $(r_1,r_2,\cos\theta_{12})$ by Proposition~\ref{prop:qreg} ($\partial_{r_1}f_j(r_1)/r_1=2\tilde f_j+r_1\tilde f_j'$; for $D=Xg$, $\partial_{r_1}D/r_1=g+r_2\cos\theta_{12}\,g_{r_1}$; for $J$, $\rho\ge r_2/\sqrt2>0$ on the domain), $E_u^b$ is analytic in $u\ge0$ (every $u$-dependent basis function is $u^2$ times an analytic function of $u$, or depends on $u$ through $x_1\!\cdot\!x_2$), and $g/u$ is analytic where $u>0$.

\emph{$\mathcal C_{\rm coal}$.} By \eqref{eq:coal} all of $r_1,r_2,u,\rho,\cos_1,\cos_2,\cos\theta_{12},x_1\!\cdot\!x_2$ are real-analytic in $(R_c,\nu,\chi)$ on the closed domain, and $r_k=R_c\,\sigma_k(\nu,\chi)$, $\rho=R_c\,\sigma(\nu)$ with $\sigma_k,\sigma\ge1-\nu_0/2>0$. Hence every quotient in Proposition~\ref{prop:qreg} is regular: $\partial_{r_k}f_j(r_k)/r_k=2\tilde f_j(r_k)+r_k\tilde f_j'(r_k)$ is analytic in $r_k=R_c\sigma_k$; the $D$-terms give $g+r_{k'}\cos\theta_{12}\,g_{r_k}$; the derivatives of $J=-\kappa(\rho)u/(1+bu)$ are polynomials in $\kappa^{(n)}(\rho)$, $r_k/\rho=\sigma_k/\sigma$, $u/\rho=\nu/\sigma$ and $1/(1+bu)$, all analytic; $u/(1+bu)$, $u^2/(1+u)^k$ are analytic in $u=R_c\nu$; and $\kappa(\rho)b(2+bu)/(1+bu)^2$ is analytic. The term $g/u$: $g$ is a real-analytic function of $(r_1,r_2)$ vanishing to second order on the diagonal $r_1=r_2$, and in the chart $r_1-r_2=R_c\nu\,\delta(\nu,\chi)$ with $\delta$ analytic ($\delta=2\chi/(\sigma_1+\sigma_2)$), so $g=u^2\,G(R_c,\nu,\chi)$ with $G$ analytic and $g/u=u\,G$ extends analytically to $u=0$ and to $R_c=0$. (The certification bounds $g/u$ on boxes by \eqref{eq:gbound} rather than by Taylor expansion.) The weights $e^{-S^{(\tau)}}$ are analytic as compositions. Finally $\Delta w_k$, $w_k'$, $a_k$ are analytic in $r_k$.

\emph{$\mathcal C_{\rm cone}$.} On the domain $d>0$ (shown above), so $\cos_k$ and $\cos\theta_{12}=t$ are analytic in $(s,t)$, and $r_1=sR$, $r_2=R$, $u=Rd$, $\rho=R\sqrt{(1+s^2)/2}$ are $R$ times analytic functions of $(s,t)$ with $r_2/\rho,\,u/\rho,\,r_1/\rho=s\sqrt2/\sqrt{1+s^2}$ analytic. The same term-by-term inspection as for $\mathcal C_{\rm coal}$ applies, now including $s=0$ ($r_1=0$: $\partial_{r_1}f_j/r_1$ and $g+r_2\cos\theta_{12}g_{r_1}$ are analytic in $r_1=sR$) and $R=0$. The chain rule from $(r_1,r_2,u)$ to $(R,s,t)$ has the analytic Jacobian $\partial(r_1,r_2,u)/\partial(R,s,t)$ and Hessian, and $g/u=u\,G$ as before with $u=Rd\ge0$.
\end{proof}

\begin{remark}
In $\mathcal C_{\rm int}$ and $\mathcal C_{\rm cone}$ the condition $\nu\ge\nu_0$ is imposed box by box: a box whose enclosure of $\nu$ (by the cancellation-free expressions above) lies entirely in $[0,\nu_0]$ is \emph{deferred} to $\mathcal C_{\rm coal}$ and recorded, not evaluated; every other box is evaluated and must be accepted. The coverage statement of Theorem~\ref{thm:cert} is that the accepted boxes and the deferred boxes together tile the chart domains, and that every deferred box lies in the domain of $\mathcal C_{\rm coal}$, which is certified entirely.
\end{remark}

\section{The certified evaluation}\label{sec:cert}

\subsection{Rigorous enclosures}
All arithmetic in this section is interval arithmetic on IEEE double precision with outward rounding \cite{Moore,Tucker}: each elementary operation $+,-,\times,\div,\sqrt{\ }$ is correctly rounded by the IEEE standard, and $\exp$, $\log$ in the libraries used are accurate to less than one ulp; we round every result outward by two ulps, so that every computed interval encloses the exact real value for all inputs in the input intervals (the evaluator is a vectorized interval library on IEEE doubles written for this purpose; no multiprecision arithmetic is used). One step is not outward-rounded: the combination of the enclosures into the right side of \eqref{eq:taylor} below (the products of the gradient and Hessian bounds with the half-widths, and the final sum of at most twelve terms) is carried out in ordinary double precision. Each of these operations has relative error $2^{-53}$, so the computed value $\tilde b$ of the right side differs from its exact value $b$ by at most $24\cdot2^{-53}(\abs{q(m)}+L)\le3\times10^{-15}\,(\abs{q(m)}+L)$, where $L\ge0$ is the subtracted loss term; and on an accepted box $L=q(m)-b\le q(m)$, so $\abs{\tilde b-b}\le6\times10^{-15}\,q(m)$. A box is therefore correctly accepted whenever $\tilde b$ exceeds its threshold by more than $6\times10^{-15}q(m)$; the coverage audit recomputed $q(m)$ on the thousand boxes of smallest excess in each piece: the largest value among them is $5.0$ in the interior, $14.4$ in the shell, $1.1$ in the gap, $0.71$ in the sliver, $12.9$ in the ball and $245$ in the strip, so on these boxes the rounding error is at most $1.5\times10^{-12}$, against minimum excesses of $9.4\times10^{-10}$, $9.5\times10^{-8}$, $5.6\times10^{-8}$, $4.9\times10^{-4}$, $0.57$ and $3.9\times10^{-4}$ respectively (at least four orders of magnitude in every piece). On every other box the excess is at least as large, and since $q\le3.6\times10^{3}$ on the region (its maximum is on the coalescence line at $R_c=240$, Appendix~\ref{app:numerics}) the rounding error is below $2.2\times10^{-11}$. (A later version of the evaluator rounds this step outward as well.) The coefficients $c^{(\tau)}_k$ are the exact double-precision numbers of Appendix~\ref{app:coef} (printed to 17 significant digits); the function $\varphi$ \emph{is} the function with those coefficients, so no rounding of coefficients is involved.

\subsection{The bound on a box}
Let $y=(y_1,y_2,y_3)$ be the coordinates of a chart and $Q=\prod_i[m_i-h_i,m_i+h_i]$ a box in the chart domain. By Taylor's theorem with Lagrange remainder along the segment from $m$ to $y$, for every $y\in Q$,
\begin{equation}\label{eq:taylor}
q(y)\ \ge\ q(m)-\sum_i\abs{\partial_iq(m)}\,h_i-\tfrac12\sum_{i,j}\sup_Q\abs{\partial_i\partial_jq}\,h_ih_j.
\end{equation}
The value $q(m)$ and the gradient $\nabla q(m)$ are enclosed by evaluating exact expressions at the point $m$ (degenerate intervals), the Hessian by evaluating exact expressions on the interval box $Q$. The expressions are generated symbolically: with $S^{(\tau)}$ and the quantities $S_1$, $S_2$, $S_u$, $S_{11}$, $S_{22}$, $S_{uu}$, $S_{1u}$, $S_{2u}$, $E_1$, $E_2$, $E^b_u$ of each term, and $u$, treated as formal functions of $y$, $q$ is the rational-exponential expression \eqref{eq:qmean}--\eqref{eq:qreg} (with the weights written as $e^{-(S^{(\tau)}-S_{\rm ref})}$ for a per-box constant $S_{\rm ref}$), and $\partial_iq$, $\partial_i\partial_jq$ are computed by symbolic differentiation; the formal derivatives of the $S$-quantities are supplied from tables of all $(r_1,r_2,u)$-derivatives of $J$ and of the $200$ basis functions up to fourth order (and of the regularized $E$-forms up to second order), converted to chart coordinates by the chain rule \eqref{eq:chain} (resp.\ by the interval Jacobian and Hessian of the chart maps \eqref{eq:coal} and of $\mathcal C_{\rm cone}$, generated symbolically); the dipole terms $D_1,D_2,D_3$ and the cusp function $J$ are differentiated directly in the chart variables, where they are regular (Lemma~\ref{lem:reg}). Every table expression is evaluated in two algebraically equivalent forms (a partial-fraction form, tight at large $r$ and $u$, and the original product form, tight on intervals that are wide relative to their distance from $0$) and the two enclosures are intersected. Two devices are essential for the enclosures to be useful: the order-zero and order-one quantities are themselves enclosed in centered form, $D(Q)\subset D(m)+\sum_i\sup_Q\abs{\partial_iD}\,[-h_i,h_i]$, intersected with the naive enclosure (the coefficients $c_k$ cancel strongly, and a naive enclosure of $S$ over a box of relative size $10^{-3}$ is a hundred times too wide); and the weights are referenced to $S_{\rm ref}=\min_\tau S^{(\tau)}(m)$ (otherwise the symbolic Hessian of the softmax ratio contains astronomically cancelling terms). With these, the right side of \eqref{eq:taylor} is within a few $10^{-4}$ of $q(m)$ for boxes of relative size $10^{-3}$ in the interior, so boxes can be taken of the size $q/\abs{\nabla q}$, which is the optimum of any first-order method.

\subsection{Branch and bound}
Each chart domain is partitioned into cells (logarithmically in the radial coordinates), each cell into boxes; a box is accepted when the right side of \eqref{eq:taylor} exceeds the cell's threshold ($0$ for cells inside $\{\rho_{\max}\le R\}$, $0.0145$ for cells meeting $\{\rho_{\max}\ge R\}$, as required by Theorem~\ref{thm:far}), deferred when it lies in $\{\nu\le\nu_0\}$ (Remark after Lemma~\ref{lem:reg}; only in $\mathcal C_{\rm int}$ and $\mathcal C_{\rm cone}$), and otherwise bisected in the coordinate direction with the largest loss term in \eqref{eq:taylor}, down to a minimal relative size $\delta_{\min}=10^{-5}$. The procedure terminates with every box accepted or deferred; the set of boxes with their bounds is the certificate.

\begin{theorem}[Certified lower bounds]\label{thm:cert}
For the function $\varphi$ of \eqref{eq:phi} with the coefficients of Appendix~\ref{app:coef}, $q\ge0$ on $\{\rho_{\max}\le120\}\setminus\Sigma$ and $q\ge0.0145$ on $\{120\le\rho_{\max}\le240\}\setminus\Sigma$. More precisely, the branch-and-bound procedure terminated on each of the following pieces with every box accepted or deferred, the pieces cover $\{r_1\le r_2\le240\}$, and every deferred box lies in the domain of $\mathcal C_{\rm coal}$:
\begin{center}
\footnotesize
\setlength{\tabcolsep}{3pt}
\begin{tabular}{@{}llllrrl@{}}
\toprule
piece & chart & domain ($r_1\le r_2$ throughout) & thr. & boxes & evaluations & min.\ bound\\
\midrule
interior & $\mathcal C_{\rm int}$ & $10^{-4}\le r_1,r_2\le10$, $t\le0.999$ & $0$ & $49\,024\,580$ & $98\,040\,712$ & $9.4\times10^{-10}$\\
shell & $\mathcal C_{\rm int}$ & $10\le r_2\le240$, $r_1\ge10^{-4}$ & $0$ / $0.0145$ & $1\,034\,744$ & $1\,844\,706$ & $1.8\times10^{-7}$\\
gap & $\mathcal C_{\rm int}$ & $2\!\times\!10^{-4}\le r_2\le10$, $r_1\ge10^{-4}$, $t\ge0.999$ & $0$ & $231\,288$ & $454\,293$ & $5.6\times10^{-8}$\\
sliver & $\mathcal C_{\rm int}$ & $r_1\le10^{-4}$, $2\!\times\!10^{-4}\le r_2\le240$ & $0$ / $0.0145$ & $4\,730$ & $6\,772$ & $4.9\times10^{-4}$ / $0.062$\\
ball & $\mathcal C_{\rm cone}$ & $r_2\le2\!\times\!10^{-4}$, all $s,t$ & $0$ & $1\,694$ & $2\,470$ & $0.57$\\
strip & $\mathcal C_{\rm coal}$ & $R_c\le240$, $\nu\le0.07$, all $\chi$ & $0$ / $0.0145$ & $126\,932$ & $218\,024$ & $3.9\times10^{-4}$ / $0.029$\\
\bottomrule
\end{tabular}
\end{center}
(``thr.'' is the threshold; where two values are given, the second refers to the cells meeting $\{\rho_{\max}\ge R\}$, and ``all $t$'' is understood where no restriction on $t$ is shown.) In the shell, the gap and the ball, $21\,517$, $18\,889$ and $569$ boxes respectively were deferred to the strip; the sliver and the strip deferred none. No box failed in any piece; in the cells with threshold $0.0145$ the minimum certified bound exceeds the threshold by $9.5\times10^{-8}$ (shell), which is how branch-and-bound terminates. Three successive versions of the evaluator were used (for the interior, for the shell, and for the four remaining pieces), differing in the tightness of the table enclosures and in the charts supported, each validated by the tests of Remark~\ref{rem:checks} before use; every certificate file records the version that produced it. The computation took about $560$ core-hours on Stanford's Sherlock cluster (interior $\approx520$, shell $20$, the four remaining pieces $16$); the certificate ($570$\,MB: every box with its bound, the deferred boxes, the evaluators, the drivers, the logs and the coverage audit) and the code are archived at Zenodo, DOIs 10.5281/zenodo.23753889 (certificate) and 10.5281/zenodo.23754394 (code).
\end{theorem}

\begin{remark}[Checks on the evaluator]\label{rem:checks}
(i) On degenerate boxes the enclosure of $q$ agrees with an independent floating-point implementation of \eqref{eq:q}--\eqref{eq:qmean} (which does not use the regularized form) to a relative $6\times10^{-11}$ in $\mathcal C_{\rm int}$ ($3.4\times10^{-12}$ for the evaluator version used on the interior) and $10^{-10}$ through $\mathcal C_{\rm coal}$ and $\mathcal C_{\rm cone}$, at all test points with $r_1,u\ge10^{-4}$ (at points with $r_1$ or $u\sim10^{-6}$ the unregularized reference itself loses accuracy by cancellation, and the agreement is $10^{-6}$). (ii) On boxes of relative size $10^{-3}$ the right side of \eqref{eq:taylor} is within $0.2\%$ of $q(m)$ at all regular test points, and on boxes touching $r_1=0$, $R=0$ or $u=0$ the bounds are finite and positive (e.g.\ the whole neighbourhood $R\le2\times10^{-4}$, $s\le0.05$ of the triple point in one box of $\mathcal C_{\rm cone}$ gives $12.3$, against $q\approx13.4$). (iii) An earlier supersolution $\varphi_0$ in the sub-basis $f_0,\dots,f_{11}$ ($74$ functions, $148$ coefficients, floating-point minimum $q\approx0.0047$ at $(r_1,r_2,u)\approx(0.75,3.4,2.6)$), from which $\varphi$ was obtained by re-optimizing the short-range coefficients with the tail-active ones fixed, was certified independently on the whole interior chart $r_1\le r_2\le10$, $\cos\theta\le0.999$: $46\,434\,464$ boxes, $92\,860\,480$ evaluations, no failed box. Two different supersolutions certified on the same region is a check of the entire pipeline that does not depend on either function. (iv) The one-body analogue: for the radial operator $\tfrac13h$ the same construction with ten radial functions gives a supersolution whose local energy is certified to be $\ge-0.729046$ on $[0,10^6]$ (mean-value interval bounds on $117\,820$ adaptive segments) and $\ge-0.680636$ on $[10^6,\infty)$ (interval evaluation in the variable $s=r^{-1/2}\in[0,10^{-3}]$ on $2\,000$ segments, the exponentially decaying basis function bounded separately), hence $\ge-0.7291$ on $[0,\infty)$, against the numerical $\tfrac13e_0=-0.7221$; this is the constant $\tfrac13e_0^{\rm rig}$ of Lemma~\ref{lem:far}. (v) Re-evaluating random samples of certified boxes on other CPU models reproduces the recorded bounds to a relative $6\times10^{-13}$ but not bit for bit (the last bits of the vectorized $\exp$ and $\log$ differ between CPU models); this is covered by the two-ulp outward rounding, and it is why the certificate is a set of boxes with bounds rather than a hash.
\end{remark}

\section{Proof of Theorem~\ref{thm:main}}\label{sec:proof}
By Lemma~\ref{lem:phiprops}, $\varphi$ satisfies the hypotheses of Lemma~\ref{lem:barta}. By Theorem~\ref{thm:cert}, $q\ge0$ on $\{\rho_{\max}\le120\}\setminus\Sigma$ and $q\ge0.0145$ on $\{120\le\rho_{\max}\le240\}\setminus\Sigma$, so Theorem~\ref{thm:far} gives $m\ge0$ on $L^2(\R^6)$. By \eqref{eq:suff}, $\M\ge0$ on antisymmetric four-electron functions. By Theorem~\ref{thm:reduction}, $E_4(2)=E_3(2)$. Since Lieb's bound \cite{Lieb} gives $E_N=E_{N-1}$ for $N\ge5=2Z+1$, $N_c(2)\le3$. \qed

\section{Discussion}\label{sec:discussion}

\paragraph{Margins.} The numbers behind the proof are not tight: $\inf\spec(m)=0.0355$ while the supersolution reaches $0.0102$ pointwise, and $\M\ge6\inf\spec(m)$ loses the Pauli principle entirely (Yang's bound on the two-body density matrix \cite{Yang} would give $\M\ge2(\mu_0+\mu_1+\mu_2)\approx0.26$ instead of $0.21$, but requires certifying three eigenvalues). The reduction itself is far from sharp: on the trial states we examined, $\ip\Psi{\M\Psi}$ is of order one. A better $w$ or a configuration-dependent decomposition of $\M$ would widen every margin.

\paragraph{He$^-$.} For $N=3$ the mechanism gives no margin: the pair term tends to $2=Z$ as the outer electron recedes, consistent with He$^-$ being nearly bound. Proving $E_3(2)=E_2(2)$ would need an argument that sees the exchange correlation of the $1s^2$ core, which the one-body weight does not.

\paragraph{Heavier dianions.} Theorem~\ref{thm:reduction} and the pair decomposition hold for every $(Z,N)$; for Li$^{2-}$ ($Z=3,N=5$) the pair operator has $\tfrac14(h_1+h_2)$ and ten pairs, and the core-pair obstruction is more severe because two electrons genuinely occupy the core. Whether a supersolution exists there is a numerical question of the same kind; the certification cost is unchanged (three variables). A uniform treatment of all $Z$ would need the pair decomposition to be replaced by something that respects the shell structure.

\paragraph{Resonances.} The theorem leaves open symmetry-protected embedded eigenvalues. Excluding them in a given symmetry sector would require Lemma~\ref{lem:lieb} with $E_{N-1}$ replaced by the bottom of the spectrum of $H_{N-1}$ in the sector reached by removing one electron, and a separate operator $\M$ for each sector.

\paragraph{Acknowledgments.} This research was supported in part by the Air Force Office of Scientific Research under grants FA9550-21-1-0381 and FA9550-24-1-0237, and by the Office of Advanced Scientific Computing Research (ASCR) within the Department of Energy Office of Science under award number DE-SC0023163. The computations were performed on the Sherlock cluster; we thank Stanford University and the Stanford Research Computing Center for computational resources and support.

\paragraph{Data availability.} The certificate (every certified box with its bound, the deferred boxes, the coverage audit, the evaluators, drivers and logs as run) and the coefficient file defining $\varphi$ are archived at Zenodo, DOI \href{https://doi.org/10.5281/zenodo.23753889}{10.5281/zenodo.23753889} (CC BY 4.0); the code is archived at DOI \href{https://doi.org/10.5281/zenodo.23754394}{10.5281/zenodo.23754394} (MIT).

\paragraph{Conflict of interest.} The authors declare no conflict of interest.

\appendix
\section{Basis and coefficients}\label{app:coef}
The coefficients are the exact double-precision numbers of the file \texttt{barta5\_coefficients\_f2.pkl} (printed to 17 significant digits, which determines a double uniquely); the basis is indexed as in Section~\ref{sec:phi} in the order $P_0,\dots,P_{21}$, $Q_0,\dots,Q_{21}$, $PP/QQ$, $U$, $UP/UQ$, $VP/VQ/UPP$, $UPP/UQQ$, $D$.
\input{coefficients}

\section{Numerical computations (not part of the proof)}\label{app:numerics}
\paragraph{Spectrum of $m$.} The $L=0$ sector of $m$ was discretized on a cell-centered radial grid $r_{i+1/2}=R(i/n)^2$ with a finite-volume form of the variable-mass kinetic energy, partial waves $P_l(\cos\theta_{12})$, $l\le l_{\max}$, coupled through the multipole expansion of $(\beta+w_1+w_2)/r_{12}$, and the lowest eigenvalues found by shift-invert Lanczos with the shift below the spectrum. On helium ($\beta=1$, $w=0$) the scheme converges from below to the $l\le2$ configuration-interaction limit $-2.9005$ with errors $19,10,6,4$ mHa at $n=60,90,130,180$. For \eqref{eq:params}: $\mu_0=0.0336,0.0344,0.0350,0.0352,0.0354$ at $n=90,120,170,240,320$ ($R=40$, $l_{\max}=3$), extrapolating like $n^{-2}$ to $0.0355$; $l_{\max}=5$ changes the value by $-10^{-4}$. A general-$(L,\pi,\text{exchange})$ solver in the coupled basis, validated on the helium $1^1S,2^3S,2^3P,2^1P$ levels, finds no state below $t_a$ in any other sector (the $L=0$ antisymmetric sector starts at $0.113$, $L=1$ at $0.134$, $L=2$ at $0.17$).

\paragraph{Construction of $\varphi$.} The coefficients were found in two stages. First, in the sub-basis $f_0,\dots,f_{11}$ ($74$ basis functions), by maximizing a soft minimum of $q$ over a grid of $14\,225$ points in $(r_1,r_2,u)$ reaching $r=10^5$ and resolving $u$ down to $10^{-3}r$ at coalescence, by L-BFGS with analytic coefficient gradients, in stages of decreasing temperature; the two-term in--out ansatz saturates at $\min q\approx-0.012$, and the third, exchange-symmetric term removes the obstruction at the equal-radius crossing, giving $\varphi_0$ with $\min q=+0.0047$ on independent fine grids. Second, in the full basis, by re-optimizing with the \emph{tail-active} coefficients ($P_j,Q_j$ for $j\in\{5,6,7,15,16,18,21\}$, products involving $f_{15},f_{18}$, and all $UP,UQ,VP,VQ,UPP,UQQ,D$ terms) held at their $\varphi_0$ values, over $38\,000$ grid points including collinear and both-far configurations. Freezing the tails is what makes the far field of $\varphi$ provably that of $\varphi_0$ (same $\sqrt r$ coefficients $\alpha_\tau$ and the same large-$r$ behaviour of $q$, which keeps the margin on $120\le\rho_{\max}\le240$ required by Theorem~\ref{thm:far}); unconstrained re-optimization in the richer basis produced functions with larger interior margin and negative local energy at collinear configurations with $r\sim10$--$40$, undetected by the optimization grid. The result has $\min q=+0.0209$ on the optimization grid and, on independent fine grids ($651\,700$ interior points with $r\le60$, $u$ resolved to $10^{-4}r$ at coalescence and to $10^{-6}$ of its range at the collinear limits; $644\,000$ points to $r=10^4$ including $\cos\theta=\pm(1-10^{-6})$), $\min q=+0.0102$ with first percentile $+0.025$, the minimum at the compact configuration $(r_1,r_2,u)\approx(0.15,0.37,0.52)$; for $r_1,r_2\ge10$ the local energy rises from $0.02$ to $0.10\approx t_b$. The maximum of $q$ on the region $\rho_{\max}\le240$ is on the coalescence line $u=0$, where $q\approx5w(R_c)$ grows with the distance of the pair from the nucleus; at $R_c=240$ it is $3.5\times10^3$.

\begin{thebibliography}{99}
\bibitem{RS4} M. Reed, B. Simon, \emph{Methods of Modern Mathematical Physics IV: Analysis of Operators}, Academic Press (1978), Theorem XIII.17.
\bibitem{Lieb} E.~H. Lieb, Bound on the maximum negative ionization of atoms and molecules, Phys. Rev. A \textbf{29}, 3018--3028 (1984).
\bibitem{Nam} P.~T. Nam, New bounds on the maximum ionization of atoms, Commun. Math. Phys. \textbf{312}, 427--445 (2012).
\bibitem{BGT} R.~D. Benguria, J.~M. Gonz\'alez-Brantes, T. Tubino, New bounds on the excess charge for atomic systems, arXiv:2207.08328.
\bibitem{HPS} D. Hundertmark, N. Pattakos, M.~R. Schulz, On the excess charge problem of atoms, arXiv:2504.18487.
\bibitem{OAI} OpenAI, Uniform excess charge for Coulomb molecules and the outer radius of neutral atoms, OpenAI Math Release preprint, September 24, 2026, \url{https://github.com/openai/math} (AI-generated, unrefereed).
\bibitem{Hofacker} M. Hofacker, Rigorous lower bounds for atomic and ionic ground-state energies, arXiv:2610.34480.
\bibitem{BGS} A.~I. Boldyrev, M. Gutowski, J. Simons, Small multiply charged anions as building blocks in chemistry, Acc. Chem. Res. \textbf{29}, 497--502 (1996).
\bibitem{Andersen} T. Andersen, Atomic negative ions: structure, dynamics and collisions, Phys. Rep. \textbf{394}, 157--313 (2004).
\bibitem{Agmon} S. Agmon, \emph{Lectures on Exponential Decay of Solutions of Second-Order Elliptic Equations}, Princeton University Press (1982).
\bibitem{Nakashima} H. Nakashima, H. Nakatsuji, Solving the Schr\"odinger equation for helium atom and its isoelectronic ions with the free iterative complement interaction (ICI) method, J. Chem. Phys. \textbf{127}, 224104 (2007).
\bibitem{Pekeris} C.~L. Pekeris, Ground state of two-electron atoms, Phys. Rev. \textbf{112}, 1649--1658 (1958).
\bibitem{Yang} C.~N. Yang, Concept of off-diagonal long-range order and the quantum phases of liquid He and of superconductors, Rev. Mod. Phys. \textbf{34}, 694--704 (1962).
\bibitem{Hales} T. Hales et al., A formal proof of the Kepler conjecture, Forum Math. Pi \textbf{5}, e2 (2017).
\bibitem{Tucker} W. Tucker, \emph{Validated Numerics: A Short Introduction to Rigorous Computations}, Princeton University Press (2011).
\bibitem{Moore} R.~E. Moore, R.~B. Kearfott, M.~J. Cloud, \emph{Introduction to Interval Analysis}, SIAM (2009).
\end{thebibliography}
\end{document}