% A self-contained rewrite of "A sharp Fourier certificate for planar circle packing".
% Goal: clear, correct, concise, high-school approachable, full proof.
% Single-file version: everything (data, tables, lemmas) inlined.
\documentclass[11pt]{article}
\usepackage[OT1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage[margin=1in]{geometry}
\usepackage{amsmath,amssymb,amsthm,mathtools}
\usepackage{booktabs,longtable,array}
\usepackage{microtype}
\usepackage{xcolor}
\usepackage[numbers,sort&compress]{natbib}
\usepackage{hyperref}
\hypersetup{colorlinks=true,linkcolor=blue!45!black,citecolor=blue!45!black,urlcolor=blue!45!black}
\newcommand{\RR}{\mathbb R}
\newcommand{\CC}{\mathbb C}
\newcommand{\ZZ}{\mathbb Z}
\newcommand{\dd}{\mathop{}\!\mathrm d}
\newcommand{\norm}[1]{\lVert#1\rVert}
\newcommand{\abs}[1]{\lvert#1\rvert}
\DeclareMathOperator{\vol}{vol}
\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}
\numberwithin{equation}{section}
\setlength{\emergencystretch}{1.5em}
\title{A sharp Fourier certificate for planar circle packing}
\author{}
\date{}
\begin{document}
\maketitle
\begin{abstract}
We prove the two-dimensional sharpness conjecture of Cohn and Elkies:
their ``two-point'' linear programming bound with a suitable auxiliary
function attains the optimal circle-packing density $\pi/(2\sqrt3)$.
The construction gives a real radial Schwartz function $f:\RR^2\to\RR$
satisfying $\widehat f(0)=1$, $f(0)=2/\sqrt3$, $\widehat f\ge0$ everywhere,
and $f\le0$ outside the unit disk, so the Cohn--Elkies argument applies.
The same certificate recovers the classical fact that equality among
periodic packings forces the triangular lattice.
Every numerical assertion reduces to a finite list of exact rational
comparisons, checked by interval arithmetic and by an independent
rational algorithm whose scope and error bound we state.
\end{abstract}
\tableofcontents
\bigskip
%=====================================================================
\section{Introduction}
\label{sec:intro}
%=====================================================================
\subsection{The planar circle-packing problem}
\label{subsec:problem}
Place nonoverlapping disks of radius $1/2$ in the plane. If
$\mathcal C\subset\RR^2$ is the set of centers, the separation condition
is $\abs{x-y}\ge1$ for distinct $x,y\in\mathcal C$. For $R>0$ let $B_R$
be the disk of radius $R$ centered at origin. The covered-area upper
density is
\[
\overline\delta(\mathcal C)
=\limsup_{R\to\infty}
\frac{\vol\bigl(B_R\cap\bigcup_{c\in\mathcal C}B(c,1/2)\bigr)}
{\vol B_R}.
\]
The triangular lattice $\Lambda_0=\ZZ(1,0)+\ZZ(1/2,\sqrt3/2)$ has density
$\pi/(2\sqrt3)\approx0.9069$, and Thue's theorem says no packing can do
better. Hales \cite{Hales2000} presents an elementary geometry-only proof.
\subsection{The main theorem}
\label{subsec:main-thm}
We prove a completely different, purely Fourier-analytic upper bound that
is sharp. On $\RR^2$ use the convention
\[
\widehat f(\xi)=\int_{\RR^2}f(x)\,e^{-2\pi i x\cdot\xi}\,\dd x.
\]
For real radial integrable $f$ the transform is real and radial. A
\emph{Schwartz function} is smooth and its derivatives decay faster than
any inverse power of $\abs x$; this is stronger than the admissibility
needed below.
\begin{theorem}[Sharp Fourier certificate]
\label{thm:main}
There exists a real radial Schwartz function $f:\RR^2\to\RR$ with
\begin{equation}\label{eq:main-conditions}
\widehat f(0)=1,\quad f(0)=\frac2{\sqrt3},\quad
\widehat f(\xi)\ge0\ \text{for all }\xi,\quad
f(x)\le0\ \text{for }\abs x\ge1.
\end{equation}
\end{theorem}
Section~\ref{sec:cohnelkies} shows that such an $f$ gives the packing
bound $\overline\delta\le\frac{\pi}{4}\cdot\frac{f(0)}{\widehat f(0)}
=\frac{\pi}{2\sqrt3}$ for every packing. So Theorem~\ref{thm:main} proves
Thue's theorem again, without geometry. Section~\ref{subsec:strategy}
gives a roadmap; the body of the paper carries it out.
\subsection{Where the difficulty lies}
\label{subsec:consistency}
Cohn and Elkies \cite{CohnElkies2003} proved that \emph{any} $f$ with the
sign conditions \eqref{eq:main-conditions} gives the density bound
$f(0)/\widehat f(0)$, and asked (Conjecture~7.3 of \cite{CohnElkies2003})
whether a single $f$ attains equality in dimension two. In dimensions
$8$ and $24$ the answer is yes, by Viazovska \cite{Viazovska2017} and by
Cohn, Kumar, Miller, Radchenko, and Viazovska \cite{CKMRV2017}, who built
explicit candidates from modular forms. Those methods do not carry over
directly.
The difficulty is global. Theorem~\ref{thm:main} requires $f\le0$
\emph{everywhere} outside the unit disk, and $\widehat f\ge0$ everywhere.
A function that vanishes with the correct signs on an infinite periodic
grid is not obtained by local patching. We obtain such a function by
combining a finite number of exact inequalities (one for each bounded
interval between consecutive nodes) with an analytic error bound on the
unbounded region. The finite inequalities are checked exactly, by
interval arithmetic on rational inputs. No sign is decided by comparing
two approximate real numbers.
\subsection{Strategy of the proof}
\label{subsec:strategy}
The proof has four stages, matching the sections below.
\begin{enumerate}
\item \textbf{A linear family of candidates.} We build two-parameter
families of entire functions $p_1,p_2$ whose ``Gaussian modulations''
$F$, $\widehat F$ form a Fourier pair (Section~\ref{sec:build}).
These are trigonometric Hermite interpolants multiplied by a Gaussian;
the Gaussian gives decay, and the trigonometric polynomial gives exact
zeros at a chosen periodic set.
\item \textbf{A solvable interpolation system.} We impose double zeros at
each node of a periodic set that contains every positive value of the
triangular shell form, plus a simple zero at the origin. Solving these
conditions is a linear system on the coefficient lists; its kernel
operator factorizes into a finite block (solved exactly by interval
arithmetic on rational inputs) and a remainder whose operator norm is
small enough to be inverted by a Neumann series (Section~\ref{sec:solve}).
\item \textbf{Signs via Bernstein polynomials.} Near each prescribed zero,
we show that $p_i$ has the correct sign by checking finitely many
Bernstein coefficients of a truncated Taylor quotient. On the far tail,
a lower bound on the trigonometric factor $P$ and a uniform
second-derivative bound on the error term give the same conclusion
(Section~\ref{sec:signs}).
\item \textbf{Exact normalization.} The simple zero at $1$ carries the
prescribed slope $-Q_1$; Poisson summation on the triangular lattice
then fixes both origin values by an elementary calculation involving
only the six shortest lattice vectors (Section~\ref{sec:poisson}).
\end{enumerate}
All analytic steps of the proof appear in the main text; the only
external input is a list of exact rational inequalities, stated as
lemmas in Sections~\ref{sec:solve}--\ref{sec:signs} and reproduced by
the interval-arithmetic procedure of Appendix~\ref{app:moments}.
%=====================================================================
\section{The Cohn--Elkies linear programming bound}
\label{sec:cohnelkies}
%=====================================================================
We recall the Cohn--Elkies argument so the rest of the paper has a
self-contained logical base. Call a pair $(f,\widehat f)$ \emph{admissible}
if $f$ is real, radial, Schwartz, and both $f$ and $\widehat f$ are
integrable on $\RR^2$.
\begin{theorem}[Cohn--Elkies bound]
\label{thm:cohnelkies}
Let $f$ be admissible, radial, with $\widehat f\ge0$ everywhere and
$f(x)\le0$ for all $\abs x\ge1$. Then every packing of unit-diameter
disks in the plane has covered-area upper density
\[
\overline\delta(\mathcal C)
\le\frac{\pi}{4}\cdot\frac{f(0)}{\widehat f(0)}.
\]
\end{theorem}
We split the proof into the periodic case (elementary, transparent) and the
general case (quote Cohn--Elkies).
\subsection{Periodic case}
\label{subsec:periodic}
Let $\mathcal C=\bigcup_{j=1}^N(L+t_j)$ be a periodic packing of period
lattice $L$ of covolume $V$, with $t_j$ distinct modulo $L$ and any two
distinct centers at distance $\ge1$. The density is $N/V$.
Consider the absolutely convergent double sum
\[
S=\sum_{j,k=1}^N\sum_{v\in L}f(v+t_j-t_k).
\]
Since $f\le0$ outside the unit disk and $f$ is Schwartz, every summand
is $\le0$ except the $N$ diagonal contributions $f(0)$. Hence $S\le N f(0)$.
Poisson summation on each translated lattice gives
\[
S=\frac1V\sum_{w\in L^*}\widehat f(w)
\Bigl\lvert\sum_{j=1}^N e^{2\pi i w\cdot t_j}\Bigr\rvert^2.
\]
Because $\widehat f\ge0$, each term is nonnegative; the $w=0$ term alone
is $(N^2/V)\,\widehat f(0)$. So $S\ge N^2\widehat f(0)/V$.
Combining, $N^2\widehat f(0)/V\le N f(0)$, so $N/V\le f(0)/\widehat f(0)$.
Multiplying by the disk area $\pi/4$ gives the covered-area bound.
This dichotomy is standard: separation makes $f$ nonpositive on all
off-diagonal contributions; the transform sign makes every Poisson
coefficient nonnegative.
\subsection{General (aperiodic) packings}
\label{subsec:aperiodic}
Cohn and Elkies in \cite[Theorem~3.1]{CohnElkies2003} extended the argument
to arbitrary (not necessarily periodic) packings, using a centered
upper-density averaging argument rather than a fixed period lattice.
The three hypotheses they require ($f$ admissible, $\widehat f\ge0$,
$f\le0$ off the unit ball) hold because $f$ is Schwartz and because
Theorem~\ref{thm:main} arranges the signs. Their conclusion is the same
density upper bound; their center-density normalization, when multiplied
by the disk area $\pi$ (for unit-diameter disks), agrees with ours. So
Theorem~\ref{thm:cohnelkies} holds in complete generality, and the
remainder of this paper is devoted to constructing $f$ satisfying
\eqref{eq:main-conditions}.
%=====================================================================
\section{Construction: Gaussian--Hermite cardinal functions}
\label{sec:build}
%=====================================================================
We build $f$ in one shot from two coefficient lists. This section sets up
the machinery. A few constants appear repeatedly and are fixed now.
\subsection{The parameter choices}
\label{subsec:params}
Put
\begin{equation}\label{eq:params}
b=\frac{\sqrt3}{2},\quad h=\frac25,\quad B=b^{-2}=\frac43,\quad
s=b\abs x^2.
\end{equation}
The coordinate $s$ is used both on the physical side $x\in\RR^2$ and on the
Fourier side $\xi\in\RR^2$. The single-scale choice of $s$ is what makes
the physical-space and Fourier formulas symmetric below.
The relevant lattice is
\begin{equation}\label{eq:lattice}
\Lambda=\ZZ u+\ZZ v,\quad
u=b^{-1/2}(1,0),\quad v=b^{-1/2}(1/2,b).
\end{equation}
Its covolume is $1$, and its dual is $\Lambda^*=\mathcal R\Lambda$ where
$\mathcal R$ is $90^\circ$ rotation. For $a=ju+kv$ with $j,k\in\ZZ$ the
squared radial coordinate is
\begin{equation}\label{eq:shell}
s_a:=b\abs a^2=j^2+jk+k^2.
\end{equation}
So $\Lambda$ and $\Lambda^*$ have the same shell radii, a key symmetry
that turns the interpolation problem into a closed coupled system.
\subsection{The interpolation nodes}
\label{subsec:nodes}
Let
\begin{equation}\label{eq:nodes}
A=\{0,1,3,4,7,9\},\quad
\mathcal N=(12\ZZ+A)\cap(0,\infty).
\end{equation}
\emph{The set $\mathcal N$ contains every positive value of $j^2+jk+k^2$.}
Indeed, $j^2+jk+k^2\equiv(j-k)^2\pmod3$, so it is never $2\bmod 3$; and
mod $4$ it is $0$ when $j,k$ are both even, otherwise odd. Combining, only
residues in $A$ appear. The inclusion is strict: $15\in\mathcal N$ but
$15=j^2+jk+k^2$ would force $(2j+k)^2\equiv2k^2\pmod5$, so $5\mid k$ and
$5\mid j$, making $25\mid15$. So we add one extra periodic node, $15$,
and its translates.
We impose zeros at all of $\mathcal N$, not merely at the lattice shell
radii. This is deliberate: a periodic node pattern gives us a fixed
trigonometric vanishing polynomial.
\subsection{The vanishing polynomial}
\label{subsec:poly}
Define, for $s\in\CC$,
\begin{equation}\label{eq:poly}
P(s)=\prod_{a\in A}\Bigl(2\sin\frac{\pi(s-a)}{12}\Bigr)^2
=\sum_{j=-6}^{6}P_j\,e^{i\pi t_j s},\quad t_j=\frac j6.
\end{equation}
Because $P$ has period $12$ and its zeros are exactly the residues $A$ mod
$12$, each zero of $P$ on $\mathcal N$ is double. Coefficients in the
second expression are given explicitly by
\begin{equation}\label{eq:coeffs}
(P_0,P_1,P_2,P_3,P_4,P_5,P_6)
=\Bigl(5,\,-1+2ib,\,1+2ib,\,-2,\,-\frac12+ib,\,-1-2ib,\,1\Bigr),
\end{equation}
with $P_{-j}=\overline{P_j}$.
\begin{remark}
The derivation: putting $w=e^{i\pi s/6}$ and $\rho=e^{i\pi/6}$, the $a$-th
factor is $-\rho^{-a}w^{-1}(w-\rho^a)^2$. The six minus signs cancel because
$\rho^{-\sum A}=\rho^{-24}=1$. Then
$P(s)=w^{-6}\prod_{a\in A}(w-\rho^a)^2$; expanding the Laurent polynomial
gives \eqref{eq:coeffs}.
\end{remark}
\subsection{Node data $Q_n$, $D_n$}
\label{subsec:nodedata}
For each zero $n\in12\ZZ+A$ write
\begin{equation}\label{eq:node-expansion}
P(n+u)=Q_n u^2\bigl(1+D_nu+O(u^2)\bigr).
\end{equation}
Both $Q_n$ and $D_n$ are $12$-periodic in $n$. If $a\in A$ is the residue
of $n$ mod $12$, then
\begin{equation}\label{eq:node-data}
Q_n=\Bigl(\frac\pi6\Bigr)^2
\prod_{a'\in A,\,a'\ne a}
\Bigl(2\sin\frac{\pi(a-a')}{12}\Bigr)^2,
\qquad
D_n=\frac\pi6\sum_{a'\in A,\,a'\ne a}
\cot\frac{\pi(a-a')}{12}.
\end{equation}
The exponent $2$ in $P$ means a single double zero at $a$; logarithmic
differentiation of the other five factors gives $D_n$. The elementary
trigonometric values at multiples of $\pi/12$ reduce to
\begin{equation}\label{eq:node-table}
\begin{array}{c|c|c}
a & Q_n & D_n\\ \hline
0 & \pi^2/3 & -7\pi\sqrt3/18\\
1 & (2\pi^2/3)(2-\sqrt3) & \pi(3+\sqrt3)/18\\
3 & (2\pi^2/3)(2-\sqrt3) & -\pi(3+\sqrt3)/18\\
4 & \pi^2/3 & 7\pi\sqrt3/18\\
7 & (2\pi^2/3)(2+\sqrt3) & -\pi(3-\sqrt3)/18\\
9 & (2\pi^2/3)(2+\sqrt3) & \pi(3-\sqrt3)/18
\end{array}
\end{equation}
Uniform bounds follow immediately. Using $157/50<\pi<22/7$ and
$\sqrt3<26/15$, one checks
\begin{equation}\label{eq:node-bounds}
Q_n>\frac74,\qquad \abs{D_n}<\frac{53}{25}<\frac{17}{8}.
\end{equation}
The precise rational lower bound $Q_n>197192/112500$ exceeds $7/4$
by a positive margin; we will use this when $Q$ errors accumulate.
\subsection{The cardinal function and its measure representation}
\label{subsec:cardinal}
Given a constant $C\in\RR$ and a real bi-list $(c_n,d_n)_{n\in\{0\}\cup\mathcal N}$
satisfying
\begin{equation}\label{eq:l1}
\norm{(c,d)}_1:=\sum_{n\in\{0\}\cup\mathcal N}
(\abs{c_n}+\abs{d_n})<\infty,
\end{equation}
set, away from the nodes,
\begin{equation}\label{eq:cardinal}
p(s)=P(s)R(s),\quad
R(s)=C+\sum_{n\in\{0\}\cup\mathcal N}
\Bigl(\frac{c_n}{(s-n)^2}+\frac{d_n}{s-n}\Bigr).
\end{equation}
The double zero of $P$ cancels the double poles of each $1/(s-n)^2$ term
and also kills the simple poles of each $1/(s-n)$ term. So $p$ extends to
an entire function. The following lemma gives an explicit representation,
which will be the engine of everything downstream.
\begin{lemma}[Finite-measure representation]
\label{lem:measure}
Let $p$ be defined by \eqref{eq:cardinal} with $\norm{(c,d)}_1<\infty$.
Then there is a finite complex Borel measure $\mu$ on $[-1,1]$ with
$\mu(-E)=\overline{\mu(E)}$ for all Borel $E$, such that
\begin{equation}\label{eq:measure-rep}
p(s)=\int_{-1}^{1}e^{i\pi ts}\,\dd\mu(t)\qquad(s\in\CC),
\end{equation}
\begin{equation}\label{eq:TV}
\norm{\mu}_{TV}\le25\abs C+\frac{65\pi^2}{18}\sum_n\abs{c_n}
+\frac{32\pi}{3}\sum_n\abs{d_n},
\end{equation}
and $\mu$ has atoms of mass $CP_j$ at $t_j=j/6$ for $-6\le j\le6$. For
real $s$ the remaining part of $\mu$ can be folded onto $[0,1]$: on each
interval $t_{j-1}<t<t_j$, the folded density is
$\sum_n(c_n u_n(t)+d_n v_n(t))$, where
\begin{equation}\label{eq:folded}
u_n(t)=-2\pi^2e^{-i\pi tn}\sum_{l=j}^{6}(t_l-t)P_l e^{i\pi t_ln},
\quad
v_n(t)=2\pi i e^{-i\pi tn}\sum_{l=j}^{6}P_l e^{i\pi t_ln}.
\end{equation}
\end{lemma}
\begin{proof}
The division identities: for $\tau\in\RR$ and $u\in\CC$,
\begin{equation}\label{eq:division}
\frac{e^{i\pi\tau u}-1}{u}=i\pi\int_0^\tau e^{i\pi tu}\,\dd t,\quad
\frac{e^{i\pi\tau u}-1-i\pi\tau u}{u^2}
=-\pi^2\int_0^\tau(\tau-t)e^{i\pi tu}\,\dd t.
\end{equation}
They hold with the limiting values $i\pi\tau$, $-\pi^2\tau^2/2$ at $u=0$,
and extend to complex $u$. Now fix a node $n$ and write $s=n+u$. Inserting
$\tau=t_l$ into \eqref{eq:division}, multiplying by $P_l e^{i\pi t_l n}$,
and summing over $l$:
the subtracted constants sum to $P(n)P'\!(n)$ terms, which vanish because
$P$ has a double zero at $n$. So the direct series $P(s)/(s-n)^2$ and
$P(s)/(s-n)$ become oscillatory integrals. On the positive interval
$t_{j-1}<t<t_j$ only frequencies $l\ge j$ contribute; on the negative
interval only $l\le -j$ contribute. Combining gives the density and its
conjugate reflection; the factors of $2$ in \eqref{eq:folded} come from
adding the positive and negative pieces.
Total variation: because $\abs{e^{i\pi tn}}=1$ for real $t$ and
$\abs{e^{-i\pi tn}}=1$, and $\int_0^{t_l}(t_l-t)\dd t=t_l^2/2$ for $t>0$,
\[
\int_0^1\abs{u_n(t)}\dd t
\le2\pi^2\sum_{l=1}^6\abs{P_l}\int_0^{t_l}(t_l-t)\dd t
=\pi^2\sum_{l=1}^6\abs{P_l}t_l^2,
\]
\[
\int_0^1\abs{v_n(t)}\dd t
\le2\pi\sum_{l=1}^6\abs{P_l}t_l.
\]
Direct calculation with $\abs{P_l}$ from \eqref{eq:coeffs}: for
$l=1,\dots,6$, $\abs{P_l}=(2,2,2,1,2,1)$; and
$\sum_{l=1}^6\abs{P_l}t_l^2=65/18$,
$\sum_{l=1}^6\abs{P_l}t_l=32/6=16/3$. So the bounds become
$65\pi^2/18$ and $32\pi/3$. The constant part contributes atoms at
$t_j=-6,\dots,6$ with masses $CP_j$, and $\sum_{j=-6}^6\abs{P_j}=25$.
Entireness: for $s$ in a compact subset of $\CC$, the integrand in
\eqref{eq:measure-rep} and each $s$-derivative is bounded uniformly in
$t\in[-1,1]$. So the integral is entire and locally uniform limits of
finite cardinal sums exist. Away from nodes the series
\eqref{eq:cardinal} converges locally uniformly to the same limit.
Conjugate symmetry of the measure makes $p$ real on the real axis.
\end{proof}
Consequently for every integer $k\ge0$,
\begin{equation}\label{eq:deriv-bound}
\abs{p^{(k)}(s)}\le\pi^k\norm{\mu}_{TV}
\qquad(s\in\RR).
\end{equation}
\subsection{The node values}
\label{subsec:nodevalues}
From \eqref{eq:cardinal} and \eqref{eq:node-expansion}, each non-own
summand vanishes to second order at $n$. Hence
\begin{equation}\label{eq:node-values}
p(n)=Q_nc_n,\qquad
p'(n)=Q_n(d_n+D_nc_n)\qquad(n\in\{0\}\cup\mathcal N).
\end{equation}
For finite input lists this follows by direct algebra; the general
$\ell^1$ case follows from Lemma~\ref{lem:measure} together with
locally uniform convergence.
\subsection{The Gaussian modulation and Fourier formula}
\label{subsec:fouriertransform}
Define
\begin{equation}\label{eq:q}
q_p(x)=e^{-\pi hs}p(s),\qquad s=b\abs x^2.
\end{equation}
Because $e^{-\pi h s}$ decays like a Gaussian in $\abs x^2$ and $p$ is
polynomially bounded, $q_p$ is Schwartz (indeed real, radial, Schwartz, as
the reader can check by the derivative bound \eqref{eq:deriv-bound}).
Introduce
\begin{equation}\label{eq:lam-z}
\lambda(t)=\frac{i}{b(t+ih)},\qquad
z(t)=-\frac{B}{t+ih}-ih\qquad(t\in[-1,1]),
\end{equation}
and, using the measure $\mu$ of Lemma~\ref{lem:measure},
\begin{equation}\label{eq:kernel}
K_p(s)=\int_{-1}^{1}\lambda(t)\,e^{i\pi z(t)s}\,\dd\mu(t).
\end{equation}
\begin{proposition}[Fourier pairing]
\label{prop:fourier}
The $\widehat{q_p}$ is real, radial, Schwartz, and
\[
\widehat{q_p}(\xi)=e^{-\pi hs}K_p(s),
\qquad s=b\abs\xi^2.
\]
The function $K_p$ is entire and real on the real axis. Its Laplace abscissa
in the positive $s$ direction is strictly negative:
\begin{equation}\label{eq:damping}
\operatorname{Im}z(t)=h\Bigl(\frac B{t^2+h^2}-1\Bigr)
\ge\frac{26}{435}>0\qquad(-1\le t\le1).
\end{equation}
Hence every derivative of $K_p(s)$ decays exponentially as $s\to+\infty$.
\end{proposition}
\begin{proof}
Using \eqref{eq:measure-rep}:
$q_p(x)=\int_{-1}^1e^{-\pi b(h-it)\abs x^2}\dd\mu(t)$. Each spatial
derivative is a polynomial times $e^{-\pi b h\abs x^2}$; the measure has
finite variation. Differentiate under the integral: the result is
$C^\infty$ with Gaussian decay.
Two-dimensional Gaussian transform: for $\operatorname{Re}a>0$,
\begin{equation}\label{eq:gauss-T}
\widehat{e^{-\pi a\abs x^2}}(\xi)
=a^{-1}e^{-\pi\abs\xi^2/a}.
\end{equation}
Apply with $a=b(h-it)$, $t\in[-1,1]$. As $\operatorname{Re}a=bh=(\sqrt3/2)(2/5)>0$,
Fubini applies with the common Gaussian majorant $e^{-\pi b h\abs x^2}$.
Compute
$a^{-1}=\lambda(t)$, and $-\pi\abs\xi^2/a+\pi hs=i\pi z(t)s$ by direct
algebra, giving \eqref{eq:kernel}.
For \eqref{eq:damping}, compute
$\operatorname{Im}z(t)=h(B/(t^2+h^2)-1)$; since $t^2\le1$ and $h^2=4/25$
and $B=4/3$, the minimum is at $\abs t=1$ and equals
$h(B/(1+h^2)-1)=(2/5)((4/3)/(29/25)-1)=(2/5)(100/87-1)=26/435$.
Decay: for $s\ge0$,
$\abs{K_p^{(k)}(s)}\le\norm\mu_{TV}\max_{\abs t\le1}
(\abs{\lambda(t)}(\pi\abs{z(t)})^k)e^{-26\pi s/435}$.
\end{proof}
\begin{remark}
A standard polar-coordinate calculation gives the Hankel form
$K_p(s)=e^{\pi hs}\frac\pi b\int_0^\infty e^{-\pi hu}p(u)
J_0(2\pi\sqrt{us}/b)\dd u$ for real $s\ge0$. This is not needed below;
it identifies $K_p$ as the radial Fourier transform of $q_p$.
\end{remark}
\subsection{The coupled interpolation conditions}
\label{subsec:interp-cond}
Let $p_1,p_2$ have coefficient lists $(c_{i,n},d_{i,n})_{n}$, constants
$C_i$, and let $K_i=K_{p_i}$. Define combined functions
\begin{equation}\label{eq:H}
H_i=p_i+K_{3-i},\qquad i\in\{1,2\}.
\end{equation}
We impose that the two functions
\begin{equation}\label{eq:pair}
F(x)=e^{-\pi hs}H_1(s),\qquad
\widehat F(\xi)=e^{-\pi hs}H_2(s),\qquad s=b\abs{x}^2=b\abs{\xi}^2
\end{equation}
form a Fourier pair. Indeed $F=q_{p_1}+\widehat{q_{p_2}}$, and
Fourier inversion gives $\widehat F=\widehat{q_{p_1}}+q_{p_2}$, matching
because everything involved is even.
The \emph{interpolation conditions} are
\begin{equation}\label{eq:interp}
H_1(n)=H_2(n)=0,\quad
H_2'(n)=0,\quad
H_1'(n)=-Q_n\mathbf1_{\{n=1\}}
\qquad(n\in\mathcal N).
\end{equation}
So both $H_1$ and $H_2$ vanish doubly at every positive node, except
$H_1$ has a \emph{simple} zero at $1$. This single simple zero is what
we will use to set the density normalization.
\begin{lemma}[Coefficient form]
\label{lem:coeffform}
For any two real $\ell^1$ lists, \eqref{eq:interp} is equivalent to
\begin{equation}\label{eq:interp-coeff}
c_{i,n}=-\frac{K_{3-i}(n)}{Q_n},
\qquad
d_{i,n}=\frac{D_nK_{3-i}(n)-K_{3-i}'(n)}{Q_n}
-\mathbf1_{\{i=1,n=1\}}
\end{equation}
for $i\in\{1,2\}$ and all $n\in\mathcal N$.
\end{lemma}
\begin{proof}
By \eqref{eq:node-values},
$H_i(n)=Q_nc_{i,n}+K_{3-i}(n)$ and
$H_i'(n)=Q_n(d_{i,n}+D_nc_{i,n})+K_{3-i}'(n)$.
Setting $H_i(n)=0$ gives the first equation; substituting into the
derivative identity and setting $H_i'(n)=-\mathbf1_{\{i=1,n=1\}}Q_n$ gives
the second. The division by $Q_n$ is legal because $Q_n>7/4>0$.
\end{proof}
\subsection{The fixed data at zero}
\label{subsec:fixed}
We need to fix the coefficient data at the node $0$, and the constants
$C_i$, to make interpolation conditions at $1$ (\emph{simple-zero
contribution}) and at $0$ consistent with the required physical and
Fourier normalization. Precisely, take
\begin{equation}\label{eq:fixed}
(c_{1,0},d_{1,0},C_1)=(1,\,0.44,\,-0.013),\quad
(c_{2,0},d_{2,0},C_2)=(0,\,-0.368,\,0.017).
\end{equation}
The finite decimals are exact rationals. These four numbers plus the two
constants are the only free parameters of the construction; everything
else is determined by \eqref{eq:interp-coeff}.
Remark: the choice of $(0.44,-0.013,-0.368,0.017)$ is not forced by the
mathematics; it is a search outcome. Any coefficient list satisfying
\eqref{eq:interp} with the derivative sign displayed extends to the same
argument; only the numeric certificate sizes depend on the choice. We use
these values because they make the finite certificates check out.
Section~\ref{sec:solve} solves \eqref{eq:interp-coeff} for the remaining
positive-node coefficients. Section~\ref{sec:signs} then proves the sign
conditions of Theorem~\ref{thm:main}. Section~\ref{sec:poisson} fixes
the origin value and dilates to the physical exclusion radius. Section
\ref{sec:packing} completes the proof.
%=====================================================================
\section{Solving the coupled interpolation equations}
\label{sec:solve}
%=====================================================================
We must find two real $\ell^1$ coefficient lists
$x_i=((c_{i,n},d_{i,n}))_{n\in\mathcal N}$ satisfying
\eqref{eq:interp-coeff}, with the fixed data \eqref{eq:fixed} at node $0$.
The two lists are coupled because $c_{i,n}$ uses $K_{3-i}(n)$. Our plan:
\begin{enumerate}
\item Rewrite the coupled system as a pair of scalar functional equations
involving a single real linear operator $S$ on $\ell^1$.
\item Split $S$ into a finite $6\times6$ block (the small nodes
$J=\{1,3,4\}$), a finite intermediate block ($7\le n\le100$),
and an infinite tail ($n>100$).
\item Certify numerically, exactly, a bound $\norm{S_{>x,*}}<\delta_x$
on the tail, and a bound $\norm{(I-\eta D)^{-1}}<32$ on the finite
block, for $\eta=\pm1$.
\item Combine these via the Schur complement to get a bounded inverse
of $I-\eta S$ on all of $\ell^1$, and use it to prove both existence
and quantitative enclosure of the solution.
\end{enumerate}
Steps 1--3 are structural; step 4 is where the inequality certificates
enter. We explain each.
\subsection{The operator $S$}
\label{subsec:S-def}
Let
\begin{equation}\label{eq:X}
X=\ell^1(\mathcal N;\RR^2),\qquad
\norm{x}_1=\sum_{n\in\mathcal N}(\abs{c_n}+\abs{d_n})
\quad\text{for }x=((c_n,d_n)).
\end{equation}
For $x\in X$, let $p_x$ be \eqref{eq:cardinal} with positive-node
coefficients given by $x$ and with $C=c_0=d_0=0$. Define $S:X\to
\RR^{\mathcal N\times\{c,d\}}$ coordinatewise by
\begin{equation}\label{eq:S}
(Sx)_{c,n}=-\frac{K_{p_x}(n)}{Q_n},\qquad
(Sx)_{d,n}=\frac{D_nK_{p_x}(n)-K_{p_x}'(n)}{Q_n}.
\end{equation}
By Proposition~\ref{prop:fourier} plus the uniform node bounds
\eqref{eq:node-bounds} and the damping \eqref{eq:damping}, $S$ is bounded
on $X$ (Section~\ref{subsec:tailenv} makes this quantitative).
Lemma~\ref{lem:coeffform} rewrites the coupled system as
\begin{equation}\label{eq:scalar-pair}
x_1 - Sx_2 = g_1,\qquad
x_2 - Sx_1 = g_2,
\end{equation}
where
$g_i \in X$ collects the contribution of the fixed node-0 data and the
constant $C_i$, plus the single nonzero entry coming from the
simple-zero marker at $n=1$ in the second equation. Concretely,
$g_i=((g_{i,c,n},g_{i,d,n}))_n$ with
\begin{equation}\label{eq:g}
g_{i,c,n}=-\frac{K_{p^0_{3-i}}(n)}{Q_n},\quad
g_{i,d,n}=
\frac{D_nK_{p^0_{3-i}}(n)-K_{p^0_{3-i}}'(n)}{Q_n}
-\mathbf1_{\{i=1\}}\mathbf1_{\{n=1\}},
\end{equation}
where $p^0_i$ is the cardinal function \eqref{eq:cardinal} corresponding
to the fixed node-0 and constant data of $\widetilde p_i$, with all
positive-node entries zero.
Adding and subtracting \eqref{eq:scalar-pair},
\begin{equation}\label{eq:plusminus}
(I-S)(x_1-x_2)=g_1-g_2,\qquad
(I+S)(x_1+x_2)=g_1+g_2.
\end{equation}
So we need to invert $I-S$ and $I+S$ on $X$. Once that holds, $x_1$ and
$x_2$ are determined uniquely. The rest of this section is exactly that
inversion.
\subsection{The block decomposition}
\label{subsec:blocks}
Split $\mathcal N$ into
\begin{equation}\label{eq:JEB}
J=\{1,3,4\},\qquad E=\mathcal N\cap[7,100],\qquad
I_f=J\cup E=\mathcal N\cap[1,100],
\end{equation}
and the infinite tail $T=\mathcal N\cap(100,\infty)$. We write
$X=X_f\oplus X_t$ accordingly. On $X_f$ we further split $X_f=X_J\oplus X_E$.
Write $D=S_{J,J}\in\RR^{6\times 6}$, a square block on $J$, and
$L=S_{E,J}$, $U=S_{J,E}$, $V=S_{E,E}$. Then
$S_{I_f,I_f}$ is the $2\times2$ block matrix
$\begin{pmatrix}D&U\\L&V\end{pmatrix}$.
\subsection{Uniform tail envelopes}
\label{subsec:tailenv}
For the operator $S$ restricted to output nodes $>x$ (any input column),
we use the following. It comes from the damping of the kernel plus the
periodicity of $D_n,Q_n$ and the folded-mass bounds.
\begin{lemma}[Tail envelope]
\label{lem:tail}
Let $\mathcal N$ be \eqref{eq:nodes}. Let $b_j^*, a_j^*, g_j, M_j,
\beta_j$ be the envelope constants in Table~\ref{tab:env}. For any
$x\ge0$, the operator norm on $X$ of the output-nodes-$>x$ restriction is
\begin{equation}\label{eq:tail-bound}
\norm{S_{>x,*}}
\le\sum_{n\in\mathcal N\cap(x,x+12]}
\sum_{j=1}^6\frac{1+\abs{D_n}+\pi a_{j-1}^*}{Q_n}
\frac{M_j L_j(n)}{1-e^{-12\pi g_j}},
\end{equation}
where
$L_j(s)=b_{j-1}^*e^{-\pi g_js}\min\{1/6,1/(\pi\beta_js)\}$.
For $x=0,4,16,100$, the right side is respectively at most
$27.618723$, $.269094$, $.014122$, $4.617010\cdot10^{-10}$. We use the
slightly weaker bounds
\begin{equation}\label{eq:tails}
\norm{S_{>0,*}}<28,\quad \norm{S_{>4,*}}<.28,\quad
\norm{S_{>16,*}}<.015,\quad \norm{S_{>100,*}}<5\cdot10^{-10}.
\end{equation}
For input $C=1$ (the atom-only column) and $x=100$, the corresponding
bound is $<3.131186\cdot10^{-8}$, which we round to
\begin{equation}\label{eq:atomtail}
\norm{(S_{\text{atom}})_{>100,*}}<3.2\cdot10^{-8}.
\end{equation}
\end{lemma}
\begin{table}[htbp]
\centering
\begin{tabular}{c|rrr|rr}
$j$ & $b_j^*$ & $a_j^*$ & $g_j$ & $M_j$ & $\beta_j$\\ \hline
0 & 2.887 & 2.934 & 2.933 & --- & ---\\
1 & 2.665 & 2.713 & 2.440 & 56 & 0\\
2 & 2.218 & 2.269 & 1.567 & 48 & 4.8\\
3 & 1.804 & 1.859 & .900 & 35 & 3.17\\
4 & 1.486 & 1.548 & .482 & 25 & 1.946\\
5 & 1.250 & 1.320 & .224 & 19 & 1.217\\
6 & 1.073 & 1.151 & .0597 & 6.284 & .792
\end{tabular}
\caption{Envelope constants for $|\lambda(t)|$, $|z(t)|$,
$\operatorname{Im}z(t)$, and folded densities on each interval
$t_{j-1}\le t\le t_j$.}
\label{tab:env}
\end{table}
\begin{proof}[Proof sketch]
On each piece $[t_{j-1},t_j]$:
$|\lambda(t)|\le b_{j-1}^*$, $|z(t)|\le a_{j-1}^*$,
$\operatorname{Im}z(t)\ge g_j+\beta_j(t_j-t)$,
and the folded integrands satisfy
$|u_n(t)|,|v_n(t)|\le M_j$. These are direct checks from the closed
formulas \eqref{eq:lam-z} and \eqref{eq:folded}: $|\lambda|$ and
$|z|$ decrease in $t$, $\operatorname{Im}z$ is minimized at piece
endpoints, and the density bound uses periodicity of $P_l e^{i\pi t_l n}$
in $n$ mod $12$. Integrating $\lambda(t)e^{i\pi z(t)s}\dd t$ against
$e^{-\pi g_j s}$, times the $1/6$ piece length, produces the
$L_j(s)$ factor. Derivatives add $(\pi a_{j-1}^*)^\ell$. Summing over
$n> x$ and using periodicity of $D_n,Q_n$ and geometric decay of the
$L_j(n+12k)/L_j(n)$ factor gives \eqref{eq:tail-bound}.
\end{proof}
Two corollaries. First, $S$ is compact on $X$; it is the norm limit of
its finite-rank truncations $P_MS$ (the right side of
\eqref{eq:tail-bound} is $O(e^{-cM})$). Second, after building the
finite-block inverse, we can use \eqref{eq:tails} directly on the
Schur complement.
\subsection{The finite certificates}
\label{subsec:finitecert}
We need quantitative bounds on the finite $6\times6$ block $D=S_{J,J}$
of the operator $S$. Let $W_\pm$ be the explicit rational matrices
\begin{equation}\label{eq:Wplus}
W_+=10^{-4}\begin{pmatrix}
4159&4452&-2003&-10467&-2421&-1683\\
2714&6573&1689&8258&104& -335\\
-4675&4785&7003&-13100&-1355&-1163\\
36853&-37142&27514&105973&1100&2012\\
57838&-56941&36296&158623&29097&17293\\
-4705&3589&-3311&-11932&-2562&6919
\end{pmatrix},
\end{equation}
\begin{equation}\label{eq:Wminus}
W_-=10^{-4}\begin{pmatrix}
13265&-1358&-385&-381&815&-164\\
616&10426&686&-484&182&775\\
-146&4&9842&215&267&179\\
-2761&2483&-958&5174&-328&206\\
-615&-861&730&-1388&5619&-4370\\
874&820&1148&-431&905&11927
\end{pmatrix}.
\end{equation}
The rows/columns are ordered so they correspond to the $J$-nodes and
$c$-$d$ split: $(c_1,c_3,c_4,d_1,d_3,d_4)$.
\begin{lemma}[Finite certificate]
\label{lem:finitecert}
For each $\eta\in\{1,-1\}$ the finite matrix $I-\eta D$ is invertible.
Further,
\begin{equation}\label{eq:finite-bounds}
\norm{(I-\eta D)^{-1}}<32,\quad
\norm{(I-\eta D)^{-1}S_{J,E}}<3.7,\quad
\norm{S_{E\cap[7,16],J}}<.09.
\end{equation}
\end{lemma}
\begin{proof}[Proof sketch]
The proof uses only two things:
(a) $\norm{W_\eta(I-\eta D)}$ is small,
(b) $\norm{W_\eta}$ is bounded.
Specifically, exact interval arithmetic gives
\begin{equation}\label{eq:defect}
\begin{aligned}
\norm{W_+(I-D)}&<.001,\quad \norm{W_-(I+D)}<.001,\\
\norm{W_+}&<30.84,\quad \norm{W_-}&<1.83.
\end{aligned}
\end{equation}
Because $\norm{I-W_\eta(I-\eta D)}<.001<1$, the matrix $W_\eta(I-\eta D)$
is invertible by the geometric series, and
$(I-\eta D)^{-1}=(W_\eta(I-\eta D))^{-1}W_\eta$, so
\[
\norm{(I-\eta D)^{-1}}
\le\norm{W_\eta}\big/(1-\norm{I-W_\eta(I-\eta D)}).
\]
For $\eta=+1$: $30.84/(1-.001)<31$, well under $32$. For $\eta=-1$:
$1.83/(1-.001)<2$. The $\norm{(I-\eta D)^{-1}S_{J,E}}$ bounds follow from
direct interval columns of $S_{J,E}$ combined with $\norm{(I-\eta D)^{-1}}$.
The exterior bound $\norm{S_{E\cap[7,16],J}}<.09$ is one more interval
column max. All are strict inequalities on exact rational numbers; see
Appendix~\ref{app:rational} for the exact algorithm.
\end{proof}
\begin{lemma}[Residual certificate]
\label{lem:residualcert}
Let $\widetilde c_{i,n},\widetilde d_{i,n}$ be the table stored in
Table~\ref{tab:coeffs}, including the fixed node-0 data
\eqref{eq:fixed}. Then for $i\in\{1,2\}$ and $n\in I_f$,
\begin{equation}\label{eq:residuals}
\begin{aligned}
\abs{\widetilde c_{i,n}+\frac{\widetilde K_{3-i}(n)}{Q_n}}
&<3\cdot10^{-10},\\
\abs{\widetilde d_{i,n}+\mathbf1_{\{i=1,n=1\}}
+\frac{\widetilde K_{3-i}'(n)-D_n\widetilde K_{3-i}(n)}{Q_n}}
&<3\cdot10^{-10}.
\end{aligned}
\end{equation}
\end{lemma}
\begin{proof}[Proof sketch]
Direct interval evaluation of the moments
$U_{\widetilde p_i,l}(m,1)$ and $X_{\widetilde p_{3-i},l}(m,1)$, via
Fej\'er quadrature with $N=256$, combining exact rational coefficient
inputs and analytic moment remainder bound $<10^{-22}$; see
Appendix~\ref{app:rational}. The certificate is $3\cdot10^{-10}$, much
looser than the moment remainder.
\end{proof}
\begin{lemma}[List-norm certificate]
\label{lem:listnorm}
For each $i$, $\norm{\widetilde x_i}_1<2.29$ (including node zero).
The tail components satisfy
$\norm{(\widetilde x_i)_{n>40}}_1<.00002$.
\end{lemma}
\begin{proof}[Proof sketch]
Direct sum of the ledgered coefficients from Table~\ref{tab:coeffs}, plus
the node-0 and constant entries.
\end{proof}
\subsection{Inversion on all of $X$}
\label{subsec:inversion}
\begin{lemma}[Full inverse]
\label{lem:inverse}
For $\eta\in\{1,-1\}$ the operator $I-\eta S$ has a bounded inverse on $X$.
If $r=r_f+r_t$ with $r_f$ supported on $I_f$ and $r_t$ on
$\mathcal N\cap(100,\infty)$, then
\begin{equation}\label{eq:inverse}
\norm{(I-\eta S)^{-1}r}_1
\le91\norm{r_f}_1+2540\norm{r_t}_1.
\end{equation}
\end{lemma}
\begin{proof}
Step 1: invert on the finite rank block $I_f$.
Let $V=(I-\eta D)^{-1}$, $L=S_{E,J}$, $U=S_{J,E}$. Then
\begin{equation}\label{eq:Lbound}
\norm L\le\norm{S_{E\cap[7,16],J}}+\norm{S_{E\cap(16,100],J}}
<.09+.015<.11.
\end{equation}
Eliminating the $J$-coordinates from $(I-\eta S_{I_f,I_f})x_f=r_f$
reduces to the Schur complement $I-\eta S_{E,E}-LVU$ on $X_E$, where
$V=(I-\eta D)^{-1}$, $L=S_{E,J}$, $U=S_{J,E}$. Its difference from the
identity has norm
\[
\norm{\eta S_{E,E}+LVU}
\le\norm{S_{E,E}}+\norm L\norm{VU}
<.28+.11\cdot 3.7=.687<1.
\]
(The strict bound $\norm{VU}=\norm{(I-\eta D)^{-1}S_{J,E}}<3.7$ is part
of Lemma~\ref{lem:finitecert}.) So the Schur complement is invertible by
a Neumann series. Solving the
$J$-equations then gives a bounded inverse $V_f$ on $I_f$ with
\[
\norm{V_f}
\le\max\Bigl\{
32+\frac{(1+3.7)\cdot.11\cdot 32}{1-.687},\quad
\frac{1+3.7}{1-.687}\Bigr\}<90.
\]
Step 2: invert on the infinite tail through another Schur complement.
Write $S$ as the block matrix
$\begin{pmatrix}S_{ff}&S_{ft}\\S_{tf}&S_{tt}\end{pmatrix}$ w.r.t.
$X=X_f\oplus X_t$. The equation $(I-\eta S)x=r$ becomes, after
elimination of $x_f$,
\begin{equation}\label{eq:tailschur}
\bigl(I-\eta S_{t,t}-S_{t,f}V_fS_{f,t}\bigr)x_t
=r_t+\eta S_{t,f}V_fr_f.
\end{equation}
By Lemma~\ref{lem:tail} with $\delta=5\cdot10^{-10}$:
$\norm{S_{t,t}}<\delta$, $\norm{S_{t,f}}<\delta$, $\norm{S_{f,t}}<28$.
Hence
$\norm{\eta S_{t,t}+S_{t,f}V_fS_{f,t}}<\delta(1+90\cdot28)<1$.
Neumann series. Solving then gives
\[
\norm{x_t}_1\le\frac{\norm{r_t}_1+90\delta\norm{r_f}_1}{1-2521\delta},
\quad
\norm{x_f}_1\le90\norm{r_f}_1+2520\norm{x_t}_1.
\]
Adding, we get $\norm{x}_1\le91\norm{r_f}_1+2540\norm{r_t}_1$.
\end{proof}
\subsection{Existence and enclosure}
\label{subsec:existence}
\begin{proposition}[Exact solution and enclosure]
\label{prop:existence}
There are unique real lists $x_1,x_2\in X$ satisfying
\eqref{eq:interp-coeff} with the fixed data \eqref{eq:fixed}. For each
$i\in\{1,2\}$,
\begin{equation}\label{eq:enclosure}
\sum_{n\in\mathcal N}
\bigl(\abs{c_{i,n}-\widetilde c_{i,n}}
+\abs{d_{i,n}-\widetilde d_{i,n}}\bigr)<.00003.
\end{equation}
The total $L^1$ norm of each exact list, including its node-0 entries,
is $<2.3$.
\end{proposition}
\begin{proof}
Use Lemma~\ref{lem:inverse} on the scalar pair \eqref{eq:plusminus}.
Define $g_i$ by \eqref{eq:g} and set
$\rho_i=\widetilde x_i-S\widetilde x_{3-i}-g_i$. This is the ``residual''
of the tabulated lists against the equation.
Split $\rho_i=r_{i,f}+r_{i,t}$:
\begin{itemize}
\item Finite part: Lemma~\ref{lem:residualcert} gives
$\norm{r_{i,f}}_1<1.1\cdot10^{-7}$ (102 entries each
$<3\cdot10^{-10}$; double for the two equations).
\item Tail part: since both the tabulated positive-node coefficients
for $n>100$ are zero, we use the atom-tail transfer. Let us estimate.
$\widetilde x_i$ on the tail is zero; so is $g_i$ on tail except
$p^0_i$'s contribution. The all-node norm $<2.29$ from Lemma~\ref{lem:listnorm}
bounds $\norm{\widetilde p_i}_1$. Then Lemma~\ref{lem:tail} plus
\eqref{eq:atomtail} gives
$\norm{r_{i,t}}_1<2.29\cdot5\cdot10^{-10}+.017\cdot3.2\cdot10^{-8}$.
\end{itemize}
Then
$\norm{x_i-\widetilde x_i}_1$ is bounded by taking half the sum/difference
of the solutions with residuals $-(\rho_1\pm\rho_2)$, and applying
Lemma~\ref{lem:inverse} twice:
\[
\norm{x_i-\widetilde x_i}_1
\le2\bigl(91\cdot1.1\cdot10^{-7}
+2540(2.29\cdot5\cdot10^{-10}+.017\cdot3.2\cdot10^{-8})\bigr)
<.00003.
\]
Uses the fact that the fixed node-0 data and constants agree between
$\widetilde p_i$ and $p_i$. This also shows
$\norm{x_i}_1<2.29+.00003<2.3$. Finally, uniqueness in $X$ follows from
invertibility of $I-\eta S$; and $x_i$ satisfies \eqref{eq:interp-coeff}
by design, hence \eqref{eq:interp} by Lemma~\ref{lem:coeffform}.
\end{proof}
%=====================================================================
\section{Proving the signs}
\label{sec:signs}
%=====================================================================
We now prove that
\[
H_1(s)<0\ (s\ge1,\ s\notin\mathcal N),\qquad
H_2(s)>0\ (s\ge0,\ s\notin\mathcal N),
\]
with exact zero sets at $\mathcal N$ and correct multiplicities. Write
$\sigma_1=-1$, $\sigma_2=1$, so the sign target is
$\sigma_iH_i(s)>0$ off the nodes. Also set
\[
\delta=3\cdot10^{-5}
\]
as the $\ell^1$ distance between exact and tabulated coefficient lists
(marginally weaker than Proposition~\ref{prop:existence}). Write
$\widetilde H_i=\widetilde p_i+K_{\widetilde p_{3-i}}$ for the tabulated
pair.
The proof has three parts: (i) a finite Taylor + Bernstein argument on
each half-gap in $(0,41.5]$, (ii) a quadratic lower bound for the
vanishing polynomial $P$, and (iii) a second-derivative bound for the
error between the exact and tabulated $H_i$ that lets the P-lower-bound
dominate on $(41.5,\infty)$.
\subsection{Reduction to Taylor quotients}
\label{subsec:taylor-reductions}
Let $\mathcal N_0=\{0\}\cup\mathcal N$. For each $m\in\mathcal N_0\cap
[0,40]$ let $m^+$ be its successor in $\mathcal N_0$ and include
$y=(m^+-m)/2$; if $m>0$ let $m^-$ be its predecessor and include
$y=(m^--m)/2$. This $y$ is signed; the corresponding half-gap is
$I_{m,y}=\{m+yv:0\le v\le1\}$.
Orders:
\begin{itemize}
\item For index $2$, $\nu=0$ at $m=0$ and $\nu=2$ at every other $m$.
\item For index $1$, exclude any interval not inside $[1,\infty)$;
$\nu=1$ at $m=1$ (only the rightward half-gap) and $\nu=2$
elsewhere.
\end{itemize}
There are $22$ nodes in $\mathcal N_0\cap[0,40]$, $43$ index-2 half-gaps,
$41$ index-1 half-gaps; total $84$ intervals.
Fix a triple $(i,m,y)$ from the list. For $s=m+u$ with $u$ in the
interval, we want $\sigma_iH_i(m+u)/u^\nu>0$. (At $u=0$ the quotient is
interpreted by continuity; the lower-order derivatives of $H_i$ vanish
because of the interpolation conditions.)
The tabulated function $\widetilde H_i$ approximates $H_i$ but its
Taylor jets are not exact. Write $\Delta_i=H_i-\widetilde H_i$. The
key error-cleaning identity is
\begin{align}\label{eq:deleted-jets}
&\frac{\sigma_i H_i(m+u)}{u^\nu}
-\sigma_i\sum_{r=0}^{\infty}
\frac{\widetilde H_i^{(r+\nu)}(m)}{(r+\nu)!}u^r
\nonumber\\
&=
\frac{\sigma_i}{u^\nu}
\Bigl(\Delta_i(m+u)
-\sum_{j=0}^{\nu-1}\frac{\Delta_i^{(j)}(m)}{j!}u^j\Bigr).
\end{align}
For $\nu=0$ the inner sum is empty. For $\nu>0$, Taylor's integral
remainder gives the bound
\[
\Bigl|\frac{\Delta_i(m+u)-\sum_{j<\nu}\Delta_i^{(j)}(m)u^j/j!}{u^\nu}\Bigr|
\le\frac1{\nu!}\sup_{s\in I_{m,y}}\abs{\Delta_i^{(\nu)}(s)}.
\]
The point: \emph{we subtract the noise's lower Taylor terms}
before dividing by $u^\nu$, otherwise that division would be singular
at $u=0$.
The finite cardinal series on the LHS is a polynomial in $v=u/y$ of
degree $\le 28$ (truncate \eqref{eq:deleted-jets} at $r=28$;
$\widetilde H_i$ and its derivatives are the sum of the moment
integrals, which are analytic so their Taylor series has full domain).
The truncation error is an infinite tail that we control separately.
\subsection{Bernstein range enclosure}
\label{subsec:bernstein}
For a polynomial
$Q(v)=\sum_{r=0}^{28}T_r v^r$ of degree $28$, its \emph{Bernstein form}
on $[0,1]$ is
\[
Q(v)=\sum_{k=0}^{28}B_k\binom{28}{k}v^k(1-v)^{28-k},
\qquad
B_k=\sum_{r=0}^{k}\frac{\binom{k}{r}}{\binom{28}{r}}T_r.
\]
The functions $v\mapsto\binom{28}{k}v^k(1-v)^{28-k}$ are nonnegative and
sum to $1$ on $[0,1]$. Hence, if every $B_k\ge \beta$, then $Q\ge\beta$
on all of $[0,1]$.
This is what we use: verify the $B_k$ lower bound by exact rational
arithmetic for $29$ values of $k$ on $84$ intervals, then the sign
quotient is $>\beta-(0.001)$ after unit scaling, and greater than $0$
if $\beta>0.001$.
\begin{lemma}[Finite Bernstein certificate]
\label{lem:bernstein}
For every specified triple $(i,m,y)$, put
\begin{equation}\label{eq:Tr}
T_r=\sigma_i y^r
\frac{\widetilde H_i^{(r+\nu)}(m)}{(r+\nu)!},
\qquad 0\le r\le 28.
\end{equation}
Then for every $0\le k\le28$,
\begin{equation}\label{eq:Bk}
B_k=\sum_{r=0}^k
\frac{\binom{k}{r}}{\binom{28}{r}}T_r>\begin{cases}
0.74,&m=0,\\
0.18,&m=1,\\
0.02,&m\in\{3,4\},\\
0.009,&4<m\le40.
\end{cases}
\end{equation}
\end{lemma}
\begin{proof}[Proof sketch]
This is a finite interval comparison: $2436$ inequalities
($29$ values of $k$ times $84$ half-gaps). Table~\ref{tab:bernstein}
records a strict lower bound on every $B_k$ across the entire half-gap
set, appearing as an entry indexed by $(m,i)$; the resulting lower bounds
are much stronger than the thresholds in \eqref{eq:Bk}. The verification
uses the moment quadrature from Appendix~\ref{app:rational} with
$N=256$, propagating $10^{-22}$ moment remainder into a
$<5.2\cdot10^{-17}$ Bernstein remainder. See Appendix~\ref{app:rational}
for details.
\end{proof}
\begin{table}[htbp]
\centering
\begin{tabular}{c|cc}
$m$&$i=1$&$i=2$\\ \hline
$0$&$\text{---}$&$.762$\\
$1$&$.314$&$.191$\\
$3,4$&$.028$&$.024$\\
$7,9$&$.114$&$.112$\\
$12,13,15,16$&$.0105$&$.0094$\\
$19,21$&$.123$&$.113$\\
$24,25,27,28$&$.0111$&$.0113$\\
$31,33$&$.121$&$.126$\\
$36,37,39,40$&$.0107$&$.0115$
\end{tabular}
\caption{Uniform lower bounds on $B_k$ over all $k$ and both half-gap
directions. For example, the $(2,16)$-entry is contained in
$[.00951875955680,.00951875955684]$ at $k=28$.}
\label{tab:bernstein}
\end{table}
\subsection{Bounding the two errors}
\label{subsec:twoerrors}
There are two error sources to control: (a) the coefficient error between
$H_i$ and $\widetilde H_i$, which we bound on the whole half-gap via
\eqref{eq:deleted-jets}; (b) the Taylor truncation error at degree $28$.
\paragraph{Coefficient error.}
The density bounds of Table~\ref{tab:env}, together with the moment
bound proved in Appendix~\ref{app:moments} and the $\ell^1$ distance $\delta$, give for
all $s\ge s_*\ge0$ and $\nu\in\{1,2\}$:
\begin{equation}\label{eq:coef-envelope}
\abs{\Delta_i^{(\nu)}(s)}
\le\delta\sum_{j=1}^6 M_j
\Bigl(\frac{(\pi t_j)^\nu}{6}+L_j(s_*)(\pi a_{j-1}^*)^\nu\Bigr).
\end{equation}
There is no atom contribution because the constants $C_i$ did not
change; the node-0 data did not change either. Direct substitution gives
\begin{equation}\label{eq:coef-caps}
\frac1{\nu!}\sup_{I_{m,y}}\abs{\Delta_i^{(\nu)}}<\begin{cases}
0.004,&m=0,\\
0.011,&m=1,\\
0.002,&m\ge3.
\end{cases}
\end{equation}
Here one uses $s_*=0$ for $m=0$, $s_*=1/2$ for $m=1$, and $s_*=2$ for
$m\ge3$.
\paragraph{Truncation error.}
Let $Y=\abs y$ and $q=29+\nu$. The tabulated cardinal measure has
total variation $<80$ on $[-1,1]$ (from Table~\ref{tab:env} and the
tabulated norm $<2.3$; not from the atom-only column). Bounding
$e^{i\pi t s}$ by its Taylor series in $s$ and absorbing $u^{-\nu}$:
\begin{equation}\label{eq:taylor-bound}
\begin{aligned}
\mathcal E(m,Y,\nu)
&:=\frac{Y^{-\nu}}{q!}\frac{1}{1-\pi a_0^*Y/(q+1)}
\Bigl(80(\pi Y)^q\\
&\qquad+2.3\sum_{j=1}^6M_jL_j(m)(\pi a_{j-1}^*Y)^q\\
&\qquad+0.017\sum_{j=0}^6L_j^P(m)(\pi a_j^*Y)^q\Bigr).
\end{aligned}
\end{equation}
The bound \eqref{eq:taylor-bound} works uniformly for positive or negative
$y$.
Five representative substitutions suffice:
\[
(m,Y,\nu)\in\{(0,1/2,0),(1,1,1),(1,1,2),(3,1,2),(4,3/2,2)\}.
\]
Numeric evaluation gives
\begin{equation}\label{eq:taylor-cases}
\mathcal E\le\begin{array}{cc}
1.578\cdot10^{-10}&(0,1/2,0)\\
2.007\cdot10^{-5}&(1,1,1)\\
5.735\cdot10^{-6}&(1,1,2)\\
9.581\cdot10^{-11}&(3,1,2)\\
6.133\cdot10^{-7}&(4,3/2,2)
\end{array}
\end{equation}
all less than $0.001$. Increasing $m$ or decreasing $Y$ can only decrease
the bound (each $L_j$ decreases in its argument), so these five cases
cover every required half-gap.
\paragraph{Combining.}
Combining Bernstein lower bounds, coefficient error caps, Taylor truncation
bounds, we obtain
\begin{equation}\label{eq:quotient}
\frac{\sigma_i H_i(s)}{(s-m)^\nu}>\begin{cases}
.735,&m=0,\\
.168,&m=1,\\
.017,&m\in\{3,4\},\\
.006,&4<m\le40.
\end{cases}
\qquad(s\in I_{m,y}).
\end{equation}
This is strictly positive, so $\sigma_iH_i(s)>0$ on each half-gap
(except at $s=m$ where both sides vanish). The union of all half-gaps is
$[0,41.5]$ for index $2$ and $[1,41.5]$ for index $1$.
\subsection{The unbounded region}
\label{subsec:unbounded}
For $s\ge s_0:=41.5$, split
\begin{equation}\label{eq:tailsplit}
p_i=P R_{i,0}+p_{i,1},\quad
R_{i,0}(s)=C_i+\sum_{n\le40}
\Bigl(\frac{c_{i,n}}{(s-n)^2}+\frac{d_{i,n}}{s-n}\Bigr),
\end{equation}
where $p_{i,1}$ is the cardinal contribution from exact coefficients at
nodes $>40$. Also define $\widetilde R_{i,0}$ from the tabulated
coefficients. All denominators in both $\widetilde R_{i,0}$ and $R_{i,0}$
are $\ge s_0-40=1.5$ on $[s_0,\infty)$.
Write $E_i:=p_{i,1}+K_{p_{3-i}}$ and $H_i=P R_{i,0}+E_i$. At each
$m\in\mathcal N$ with $m\ge43$ both $H_i$ and $P R_{i,0}$ vanish to
order $\ge2$ by the interpolation conditions and by the double zero of
$P$; hence $E_i$ does too. We will bound $|E_i''(s)|$ and $|\sigma_iR_{i,0}(s)|$
separately.
\paragraph{Rational-part lower bound.}
For the tabulated coefficients put $u_n^*=\sigma_i\widetilde c_{i,n}$,
$v_n^*=\sigma_i\widetilde d_{i,n}$, $[w]_-=\min\{w,0\}$. Then
\begin{equation}\label{eq:rational-lower}
\sigma_i\widetilde R_{i,0}(s)\ge
\sigma_iC_i+\frac{[\sum_{n\le40}v_n^*]_-}{s_0}
+\sum_{n\le40}\Bigl(
\frac{[u_n^*]_-}{(s_0-n)^2}
+\frac{[v_n^*]_- n}{s_0(s_0-n)}\Bigr)
\end{equation}
uniformly for $s\ge s_0$. Both displayed sums include $n=0$. The
inequality holds because each term with $[\,]_-$ is at its worst value at
$s=s_0$; discarded positive contributions are truly positive.
Exact rational evaluation of \eqref{eq:rational-lower} gives
$>.01206834$ for index $1$ and $>.01656468$ for index $2$. To pass to the
exact list we use
\begin{equation}\label{eq:rat-error}
\abs{R_{i,0}(s)-\widetilde R_{i,0}(s)}
\le\sum_{n\le40}\Bigl(
\frac{\abs{c_{i,n}-\widetilde c_{i,n}}}{(s-n)^2}
+\frac{\abs{d_{i,n}-\widetilde d_{i,n}}}{s-n}\Bigr)
<\frac{\delta}{s_0-40}=2\cdot10^{-5}
\end{equation}
for $s\ge s_0$. Therefore
$\sigma_i R_{i,0}(s)>0.011$ on $[s_0,\infty)$.
\paragraph{Second-derivative bound.}
The tabulated list has norm $<2.29$ on all nodes, and $<.00002$ on
$n>40$. Since the exact and tabulated coefficients are within $\delta$
in $\ell^1$, and $x_i$ has positive norm $<2.3$ (from
Proposition~\ref{prop:existence}):
\begin{align}\label{eq:second}
\sup_{s\ge s_0}\abs{E_i''(s)}
&\le .00005\cdot 77
+2.3\sum_{j=1}^6M_jL_j(s_0)(\pi a_{j-1}^*)^2
\nonumber\\
&\qquad+0.017\sum_{j=0}^6L_j^P(s_0)(\pi a_j^*)^2
<0.006.
\end{align}
The first term is the tail-from-$40$ direct part, using
$\sum_{j=1}^6M_j(\pi t_j)^2/6<77$ from Table~\ref{tab:env}; the other two
are the transformed kernel tails.
\paragraph{The P-barrier.}
\begin{lemma}
\label{lem:Plower}
For $s\in\RR$ let $m$ be a nearest zero of $P$. Then
\[
P(s)\ge0.68(s-m)^2.
\]
\end{lemma}
\begin{proof}
Consider $P(m+u)/u^2$ on each half-gap; extend at $u=0$ by
$Q_m>1.75$. Its logarithm is concave ($a=\pi/12$):
\begin{itemize}
\item Second derivative of $\log$ of $|2\sin(a(m+u-a'))|^2$ for $a'\ne m$
is $-2a^2\csc^2(a(m+u-a'))<0$.
\item Second derivative of $\log|\sin(au)/u|$ is
$1/u^2-a^2/\sin^2(au)\le0$ (since $|\sin(au)|\le|au|$).
\end{itemize}
A concave quantity is minimized at an endpoint. So on each half-gap the
minimum of $P(m+u)/u^2$ is $\ge\min\{Q_m,\,\text{midpoint value}\}$.
Put $\theta=\pi(s-2)/12$ and $c=\cos(2\theta)=\cos(\pi(s-2)/6)$. Algebraic
identities reduce each pair of nonvanishing factors to
$(1-2c)$, $(\sqrt3-2c)$, $(\sqrt3+2c)$ according to the pair of residues.
Reducing the whole product gives
\[
P(s)=(4c^2-3)^2(1-2c)^2,
\]
from which at each gap midpoint the quotient $P(s)/(s-m)^2$ is:
$12-8\sqrt2$ on $[0,1]$ and $[3,4]$;
$1$ on $[1,3]$;
$9$ on $[7,9]$;
$\frac89(3+\sqrt2+\sqrt3+\sqrt6)$ on $[4,7]$ and $[9,12]$.
The smallest is $12-8\sqrt2>0.68$: use $\sqrt2<283/200$ (since
$(283/200)^2>2$). Others $\ge1$.
\end{proof}
\paragraph{Closing.}
For $s\ge41.5$ choose a nearest node $m\ge43$ (at $s=41.5$, choose $m=43$
over the tied $m=40$). At $m\ge43$, both $H_i$ and $PR_{i,0}$ vanish
to order $\ge2$. So $E_i=H_i-PR_{i,0}$ vanishes to order $\ge2$ at $m$,
and $|s-m|\le s$ so the segment between $s$ and $m$ stays in
$[41.5,\infty)$. Hence by \eqref{eq:second},
$\abs{E_i(s)}\le\frac{0.006}{2}(s-m)^2=0.003(s-m)^2$.
Combining with Lemma~\ref{lem:Plower} and the rational lower bound:
\begin{equation}\label{eq:tail-margin}
\sigma_iH_i(s)\ge
\bigl(0.011\cdot 0.68-0.003\bigr)(s-m)^2
=0.00448(s-m)^2>0.
\end{equation}
This covers $[41.5,\infty)$ for both $i$.
\subsection{Global signs, zeros, and multiplicities}
\label{subsec:global-signs}
\begin{proposition}[Global signs]
\label{prop:global-signs}
The exact functions $H_1,H_2$ satisfy
\[
\{s\ge1:H_1(s)=0\}=\mathcal N,\quad
\{s\ge0:H_2(s)=0\}=\mathcal N,
\]
and strict inequalities off those sets. Multiplicities: the zero of
$H_1$ at $1$ is simple; its zeros at $n\in\mathcal N$, $n>1$, and the
zeros of $H_2$ at every node, are order exactly two (as zeros in $s$).
\end{proposition}
\begin{proof}
Combine the half-gap bounds \eqref{eq:quotient} on $[0,41.5]$ (or
$[1,41.5]$ for index $1$) with \eqref{eq:tail-margin} on $[41.5,\infty)$.
At each node $n\in\mathcal N$ the interpolation conditions give
$H_i(n)=0$. Between nodes the quotient is positive, so strictly positive
by continuity. Zero sets align.
Multiplicities. At $1$, $H_1'(1)=-Q_1\ne0$; simple. For a node
$m\ge1$, $m\ne1$: from \eqref{eq:quotient} with $\nu=2$, the endpoint value
$\sigma_iH_i''(m)/2>0$ (finite region, $m\le40$) or from \eqref{eq:tail-margin}
with $m\ge43$ and $u\to0$:
$\sigma_iH_i''(m)/2=\lim_{u\to0}\sigma_iH_i(m+u)/u^2\ge0.00448>0$.
So the zero has order exactly two. This handles both $i$.
\end{proof}
%=====================================================================
\section{Poisson normalization}
\label{sec:poisson}
%=====================================================================
We still need the origin values $F(0)=\widehat F(0)$. The prescribed
simple zero at $n=1$ of $H_1$ gives a nonzero slope there; Poisson
summation on the physical lattice then forces both origin values equal to
an explicit positive number. This is purely a normalization step: no
new inequality appears.
Throughout let
\[
G_i(s)=e^{-\pi h s}H_i(s),\quad F(x)=G_1(b\abs x^2),
\quad\widehat F(\xi)=G_2(b\abs\xi^2).
\]
The last two identifications are \eqref{eq:pair}. In terms of the
interpolation conditions:
\begin{equation}\label{eq:jets}
G_1(n)=G_2(n)=0,\quad G_2'(n)=0,\quad
G_1'(n)=-Q_1e^{-\pi h}\mathbf1_{\{n=1\}}\quad(n\in\mathcal N).
\end{equation}
Indeed $H_i(n)=0$ gives $G_i(n)=0$; differentiating
$G_i=e^{-\pi hs}H_i$ and using $H_i(n)=0$ gives
$G_i'(n)=e^{-\pi hn}H_i'(n)$, which is $-Q_1e^{-\pi h}$ at $n=1$
(from $H_1'(1)=-Q_1$), and zero at all other $n$ for $G_1$ and for
all $n$ for $G_2$.
\subsection{The dilation identity}
\label{subsec:dilate}
Recall $\Lambda$ from \eqref{eq:lattice}: covolume $1$, dual
$\Lambda^*=\mathcal R\Lambda$, and shell coordinate
$s_a=b\abs a^2=j^2+jk+k^2$. Every nonzero $s_a$ is a positive element of
$\mathcal N$ by the modular arithmetic of \eqref{eq:nodes}. The shell
$s_a=1$ has exactly $6$ lattice vectors:
\[
(j,k)=(1,0),(-1,0),(0,1),(0,-1),(1,-1),(-1,1).
\]
Observation: $4(j^2+jk+k^2)=(2j+k)^2+3k^2$, so $4s_a=4$ iff
$(2j+k)^2+3k^2=4$, whose solutions are exactly the six above.
Apply Poisson summation on $\Lambda$ to the rescaled $F_t(x)=F(\sqrt t x)$
for $t>0$:
\begin{equation}\label{eq:poisson-raw}
\sum_{a\in\Lambda}G_1(t s_a)
=\frac1t\sum_{a\in\Lambda}G_2(s_a/t).
\end{equation}
The derivation: $\widehat{F_t}(\xi)=t^{-1}\widehat F(\xi/\sqrt t)$
(the 2D Fourier scaling rule). $\Lambda^*=\mathcal R\Lambda$ has the same
$s(\cdot)$-values as $\Lambda$ by radiality, so summing $\widehat F$ over
$\Lambda^*$ is the same as summing over $\Lambda$. Poisson summation is
$\sum_{\Lambda}F_t=\frac1{\mathrm{covol}}\sum_{\Lambda^*}\widehat{F_t}$,
and the covolume is $1$.
Both sides are locally uniformly convergent in $t>0$, together with their
first two $t$-derivatives, because $F$ is Schwartz. The needed bound is
straightforward: on a compact $t$-interval,
$\sum_{a\in\Lambda}(1+\abs a)^{-M}$ converges if $M>2$ (number of lattice
points in a disk of radius $R$ grows as $O(R^2)$), and the first two
$t$-derivatives of $F(\sqrt t x)$ are bounded by $C_M(1+\abs x)^{-M}$.
Hence two termwise differentiations are legitimate.
\subsection{Proof of the origin value}
\label{subsec:origin}
\begin{lemma}[Origin value]
\label{lem:origin}
$F(0)=\widehat F(0)=G_1(0)=G_2(0)=6Q_1e^{-\pi h}>0$.
\end{lemma}
\begin{proof}
At $t=1$, \eqref{eq:poisson-raw} reads
$\sum_{a\in\Lambda}G_1(s_a)=\sum_{a\in\Lambda}G_2(s_a)$.
Every $a\ne0$ has $s_a\in\mathcal N$, where $G_i(s_a)=0$. So both sides
reduce to the $a=0$ term: $G_1(0)=G_2(0)$.
Differentiate \eqref{eq:poisson-raw} in $t$ and set $t=1$.
The left side is
$\sum_a s_a G_1'(s_a)$. Only $a$ with $s_a=1$ contribute, because at any
other nonzero node $s_a\in\mathcal N$ we have $G_1'(s_a)=0$ by
\eqref{eq:jets}; at $a=0$ the factor $s_a=0$. Six shell vectors with
$s_a=1$ give $-Q_1e^{-\pi h}$ each, total $6(-Q_1e^{-\pi h})$.
The right side, differentiated, is
\[
-\frac1{t^2}\sum_a G_2(s_a/t)
-\frac1{t^3}\sum_a s_aG_2'(s_a/t).
\]
At $t=1$: the first summand contributes $-G_2(0)$ (all other $a$ give
$G_2(s_a)=0$). The second summand at $t=1$ contributes
$-\sum_a s_aG_2'(s_a)$; each term with $s_a=0$ fetches an explicit $0$,
and each term with $s_a>0$ vanishes because $G_2'(s_a)=0$. So the second
summand is $0$.
Hence $-6Q_1e^{-\pi h}=-G_2(0)$, i.e. $G_2(0)=6Q_1e^{-\pi h}$. Positivity
follows from $Q_1>7/4$ in \eqref{eq:node-bounds}. Combining with
$G_1(0)=G_2(0)$ gives the result.
\end{proof}
\subsection{The final dilation and the main theorem}
\label{subsec:final-dilation}
Define
\begin{equation}\label{eq:scaling}
f(x)=\frac{F(x/\sqrt b)}{bF(0)},\quad
\widehat f(\xi)=\frac{\widehat F(\sqrt b\,\xi)}{F(0)}.
\end{equation}
The second line is the two-dimensional Fourier scaling rule applied to
the first, plus Lemma~\ref{lem:origin}. Since $F(0)>0$ and $F$,
$\widehat F$ are real radial Schwartz (Proposition~\ref{prop:fourier}),
so is $f$. And $f$'s transform is real radial Schwartz as well.
\begin{proof}[Proof of Theorem~\ref{thm:main}]
Take $f$ defined by \eqref{eq:scaling}.
\emph{Origin values.}
By Lemma~\ref{lem:origin},
$F(0)=\widehat F(0)=6Q_1e^{-\pi h}>0$. Then
\[
\widehat f(0)=\frac{\widehat F(0)}{F(0)}=1,
\qquad
f(0)=\frac{F(0)}{bF(0)}=\frac1b=\frac2{\sqrt3}.
\]
\emph{Fourier sign.}
By Proposition~\ref{prop:global-signs}, $H_2(s)>0$ for all $s\ge0$ off
the nodes, and $H_2(n)=0$ at the nodes. Since $e^{-\pi hs}>0$,
$G_2(s)\ge0$ for all $s\ge0$. That is, $\widehat F(\xi)\ge0$ for all
$\xi$. Dividing by $F(0)>0$ preserves sign:
$\widehat f\ge0$ everywhere.
\emph{Physical sign.}
For $\abs x\ge1$,
$b\abs{x/\sqrt b}^2=\abs x^2\ge1$. By Proposition~\ref{prop:global-signs},
$H_1(s)<0$ for $s\ge1$ off the nodes, i.e. $G_1(s)\le0$ (equality at
the nodes). So $F(x)=G_1(\abs x^2)\le0$ for $\abs x\ge1$, and dividing by
$bF(0)>0$ gives $f(x)\le0$.
\emph{Regularity.} $F$ is Schwartz by Proposition~\ref{prop:fourier};
scaling preserves Schwartz; so $f$ is Schwartz.
Together \eqref{eq:main-conditions} holds.
\end{proof}
\begin{remark}[Zero sets in final coordinates]
Writing $f_{\mathrm{rad}}(r)=f(r,0)$ and
$\widehat f_{\mathrm{rad}}(\rho)=\widehat f(\rho,0)$:
\[
f_{\mathrm{rad}}(r)=\frac{G_1(r^2)}{bF(0)},\qquad
\widehat f_{\mathrm{rad}}(\rho)=\frac{G_2(b^2\rho^2)}{F(0)}.
\]
Proposition~\ref{prop:global-signs} then identifies the exact zero sets:
\[
\{r^2:r\ge1,\ f_{\mathrm{rad}}(r)=0\}=\mathcal N,\qquad
\{\rho^2:\rho\ge0,\ \widehat f_{\mathrm{rad}}(\rho)=0\}
=b^{-2}\mathcal N=\tfrac43\mathcal N.
\]
Multiplicities transfer directly: the factor $e^{-\pi hs}$ is
nonvanishing, and $s=r^2$ and $s=b^2\rho^2$ have nonzero derivative at
positive radii. So the physical zero at $r=1$ is simple, every other
listed physical zero is exactly double, every listed Fourier zero is
exactly double. The physical radial slope at $r=1$ is
\[
f_{\mathrm{rad}}'(1)
=\frac{2G_1'(1)}{bF(0)}
=\frac{2\cdot(-Q_1e^{-\pi h})}{b\cdot6Q_1e^{-\pi h}}
=-\frac{1}{3b}=-\frac2{3\sqrt3}.
\]
\end{remark}
%=====================================================================
\section{From the certificate to the packing bound}
\label{sec:packing}
%=====================================================================
Theorem~\ref{thm:main} produces the auxiliary function required by
Theorem~\ref{thm:cohnelkies}. Therefore every packing of unit-diameter
disks in the plane has covered-area upper density
$\overline\delta\le\frac\pi4\cdot\frac{f(0)}{\widehat f(0)}
=\frac\pi4\cdot\frac{2/\sqrt3}{1}=\frac{\pi}{2\sqrt3}$.
This is the triangular-lattice value, established without geometry.
\subsection{Recovering the classical equality case}
\label{subsec:equality}
The certificate also recovers a classical uniqueness statement among
periodic packings, using only the physical zero set. Everything follows
from the fact that all physical zeros of $f$ outside the unit disk have
\emph{integer} squared radius.
\begin{corollary}[Periodic equality case]
\label{cor:uniqueness}
Suppose a periodic packing of unit-diameter disks attains covered-area
density $\pi/(2\sqrt3)$. Then its centers are a Euclidean isometry of
the triangular lattice $\ZZ(1,0)+\ZZ(1/2,b)$.
\end{corollary}
\begin{proof}
Write the packing as $\mathcal C=\bigcup_{j=1}^N(L+t_j)$ with period
lattice $L$ of covolume $V$. Equality of densities gives $N/V=2/\sqrt3$,
i.e. $f(0)/\widehat f(0)$.
Redo the Poisson comparison of \S\ref{subsec:periodic}. Both inequalities
used there are equalities:
\[
S=\sum_{j,k}\sum_v f(v+t_j-t_k)\le Nf(0), \quad
S\ge\frac{N^2}{V}\widehat f(0)=N\cdot\frac{2}{\sqrt3}\cdot1=Nf(0).
\]
For absolute convergence decompose the physical defect:
\[
0=Nf(0)-S=\sum_{(j,v)\ \text{off-diagonal}}(-f(v+t_j-t_k)).
\]
Every summand is nonnegative (all off-diagonal $f\le0$), so every one is
zero. Any difference $x-y$ of two distinct centers occurs as such a
summand, so
\[
\abs{x-y}^2\in\mathcal N\subset\ZZ,
\qquad x,y\in\mathcal C,\ x\ne y.
\]
Now dilate to convert integer squared distances into even integers.
Translate one center to $0$ and dilate by $\sqrt2$: put
$L_0=\sqrt2 L$, $\mathcal C_0=\sqrt2(\mathcal C-t_1)$, and set
$\Gamma=\mathrm{span}_{\ZZ}\mathcal C_0$. Since $0\in\mathcal C_0$ and all squared
distances in $\mathcal C_0$ are even integers, polarization gives
\[
\langle x,y\rangle
=\frac{\abs x^2+\abs y^2-\abs{x-y}^2}{2}\in\ZZ\quad(x,y\in\mathcal C_0).
\]
Integer linear combinations preserve integrality (cross terms get a
factor $2$), so $\Gamma$ inherits: integral pairings and even squared
norms. In particular pairing with an $L_0$-basis gives
$L_0\subseteq\Gamma\subseteq L_0^*$, so $\Gamma$ is a finite-index
sublattice of $L_0$.
Write $d=[\Gamma:L_0]$ and pick a basis $(w_1,w_2)$ of $\Gamma$. Its
Gram matrix is
$\begin{pmatrix}a&m\\m&c\end{pmatrix}$ with $a,c\in2\ZZ_{>0}$,
$m\in\ZZ$, and determinant $ac-m^2=\mathrm{covol}(\Gamma)^2>0$.
Since $a,c$ are even and $m^2\equiv0$ or $1$ mod $4$, $ac-m^2
\in\{0,3,7,\dots\}\pmod 4$. The positive value $\mathrm{covol}^2>0$
therefore satisfies $\mathrm{covol}(\Gamma)^2\ge3$, i.e.
$\mathrm{covol}(\Gamma)\ge\sqrt3$.
On the other hand $\mathcal C_0$ occupies $N$ distinct $L_0$-cosets of the
$d$ cosets of $\Gamma$, so $d\ge N$. Since $\mathrm{covol}(L_0)=2V$ and
$\mathrm{covol}(\Gamma)=2V/d$,
\[
\sqrt3\le\mathrm{covol}(\Gamma)
=\frac{2V}{d}\le\frac{2V}{N}=\sqrt3.
\]
The last equality uses $N/V=2/\sqrt3$, i.e.\ $2V/N=\sqrt3$. So all
inequalities are equalities, and $d=N$.
So $\mathcal C_0=\Gamma$.
It remains to identify $\Gamma$. Pick a shortest nonzero $w_1\in\Gamma$;
a shortest vector is primitive, so extend to a $\Gamma$-basis
$(w_1,w_2)$. Replace $w_2$ with an integer translate along $w_1$ so that
$m=\langle w_1,w_2\rangle$ satisfies $\abs m\le a/2$ where $a=\abs{w_1}^2$.
Then $c=\abs{w_2}^2\ge a$ (shortest vector) and $a\ge2$ (even).
Calculate
\[
3=ac-m^2\ge ac-a^2/4\ge\frac34a^2.
\]
So $a^2\le4$, i.e. $a=2$. Then $2c-m^2=3$ with $c$ even, $\abs m\le1$.
The only solution is $m=\pm1$, $c=2$. That is the Gram matrix of a
triangular lattice with shortest vector length $\sqrt2$; back-scaling by
$1/\sqrt2$ identifies $\mathcal C$ as an isometric image of
$\ZZ(1,0)+\ZZ(1/2,b)$.
\end{proof}
%=====================================================================
\appendix
\section{The finite coefficient table}
\label{app:coeffs}
%=====================================================================
Table~\ref{tab:coeffs} records the tabulated coefficient lists used
throughout, in the notation of Proposition~\ref{prop:existence}. Each
column is multiplied by $10^{10}$: an integer $a$ in the table means the
coefficient $a\cdot10^{-10}$. At every positive node $n>100$ both
$c_{i,n}=d_{i,n}=0$.
The node-zero data and constants are stated separately:
\begin{equation}\label{eq:coeffs-at-0}
(c_{1,0},d_{1,0},C_1)=(1,\,0.44,\,-0.013),\qquad
(c_{2,0},d_{2,0},C_2)=(0,\,-0.368,\,0.017),
\end{equation}
where the finite decimals are exact rational numbers. These already
appear in \eqref{eq:fixed}; we restate them for locality.
The uniform $\ell^1$ norms of the tabulated lists, including node zero,
satisfy (from integer addition of absolute values):
\begin{equation}\label{eq:table-norms}
\begin{aligned}
\sum_{\{0\}\cup\mathcal N}(\abs{c_{1,n}}+\abs{d_{1,n}})
&=2.2516666247<2.29,\\
\sum_{\{0\}\cup\mathcal N}(\abs{c_{2,n}}+\abs{d_{2,n}})
&=0.8927477314<2.29,
\end{aligned}
\end{equation}
and their tails beyond $40$ satisfy
\begin{equation}\label{eq:table-tails}
\begin{aligned}
\sum_{n>40}(\abs{c_{1,n}}+\abs{d_{1,n}})
&=0.0000197471<0.00002,\\
\sum_{n>40}(\abs{c_{2,n}}+\abs{d_{2,n}})
&=0.0000161623<0.00002.
\end{aligned}
\end{equation}
\begingroup
\small
\begin{longtable}{r r r r r}
\caption{Coefficient lists, scaled by $10^{10}$.}
\label{tab:coeffs}\\
\toprule
$n$ & $10^{10}c_{1,n}$ & $10^{10}d_{1,n}$
& $10^{10}c_{2,n}$ & $10^{10}d_{2,n}$\\
\midrule
\endfirsthead
\toprule
$n$ & $10^{10}c_{1,n}$ & $10^{10}d_{1,n}$
& $10^{10}c_{2,n}$ & $10^{10}d_{2,n}$\\
\midrule
\endhead
\bottomrule
\endlastfoot
1 & -670689128 & -6398628799 & -250018302 & 4107732082 \\
3 & 93524468 & 564897270 & -186215974 & -65204793 \\
4 & -64146977 & 121133118 & 82882809 & -327822296 \\
7 & -1358193 & -16329465 & 1979670 & 21637038 \\
9 & -2932910 & -940875 & 3636920 & 2145806 \\
12 & 3646952 & -34846939 & -6005382 & 34369042 \\
13 & 3183411 & 63030995 & -574123 & -75237042 \\
15 & 11390977 & 27049982 & -11312377 & -34985096 \\
16 & -5593246 & 13405071 & 5915710 & -10980954 \\
19 & 54905 & -1515007 & -120136 & 1498351 \\
21 & -210208 & -634307 & 177992 & 747946 \\
24 & 976077 & -804949 & -1037705 & 105643 \\
25 & -824472 & 7003719 & 1011732 & -6510015 \\
27 & 611891 & 4536095 & -428528 & -4440927 \\
28 & -472882 & -107613 & 399655 & 403751 \\
31 & 30744 & -112653 & -31773 & 90515 \\
33 & -4974 & -105568 & 1240 & 99081 \\
36 & 128313 & 246018 & -114876 & -267388 \\
37 & -176685 & 505187 & 166591 & -396209 \\
39 & -21638 & 472037 & 32128 & -402283 \\
40 & -18461 & -165597 & 10836 & 164974 \\
43 & 4565 & -2296 & -4050 & 547 \\
45 & 1477 & -10734 & -1524 & 8978 \\
48 & 10748 & 50050 & -8611 & -45048 \\
49 & -20506 & 7495 & 17354 & -266 \\
51 & -10361 & 29035 & 9641 & -21534 \\
52 & 2121 & -26523 & -2269 & 23159 \\
55 & 425 & 858 & -344 & -839 \\
57 & 302 & -632 & -266 & 463 \\
60 & 382 & 5658 & -231 & -4688 \\
61 & -1481 & -4320 & 1140 & 4098 \\
63 & -1494 & -352 & 1265 & 647 \\
64 & 558 & -2609 & -491 & 2114 \\
67 & 21 & 169 & -14 & -145 \\
69 & 34 & 13 & -28 & -18 \\
72 & -55 & 394 & 53 & -301 \\
73 & -21 & -805 & 5 & 685 \\
75 & -136 & -379 & 108 & 338 \\
76 & 70 & -144 & -58 & 104 \\
79 & -1 & 19 & 1 & -15 \\
81 & 2 & 9 & -2 & -8 \\
84 & -13 & 3 & 11 & 0 \\
85 & 13 & -87 & -11 & 70 \\
87 & -6 & -60 & 4 & 49 \\
88 & 6 & 5 & -4 & -5 \\
91 & 0 & 1 & 0 & -1 \\
93 & 0 & 1 & 0 & -1 \\
96 & -2 & -4 & 1 & 3 \\
97 & 2 & -6 & -2 & 4 \\
99 & 0 & -6 & 0 & 5 \\
100 & 0 & 2 & 0 & -2 \\
\end{longtable}
\endgroup
%=====================================================================
\section{Moment estimates and the finite certificate}
\label{app:moments}
\label{app:rational}
%=====================================================================
This appendix lays out the two tools used in Sections~\ref{sec:solve}
and~\ref{sec:signs}: the moment bound \eqref{eq:moment-bound-full}, and
the exact finite verification of the certificates via Fej\'er
quadrature. Both are standard; we state them explicitly so no step in
the main argument relies on an opaque numerical procedure.
\subsection{The moment bound}
\label{sec:moment-bound}
For a cardinal function
$p=P(C+\sum_n(c_n/(s-n)^2+d_n/(s-n)))$ with $w=\norm{(c,d)}_1<\infty$,
the transformed kernel \eqref{eq:kernel} has derivative
\[
K_p^{(\ell)}(s)
=\int_{-1}^1\lambda(t)(i\pi z(t))^\ell e^{i\pi z(t)s}\dd\mu(t).
\]
Splitting the folded measure onto the six pieces $t_{j-1}\le t\le t_j$
and bounding each factor by the entries of Table~\ref{tab:env} gives
the moment bound.
\begin{lemma}[Moment bound]
\label{lem:moment-bound}
Let $p$ be a cardinal function with coefficient norm $w<\infty$ and
constant $C\in\RR$. For every $s\ge0$ and every integer $\ell\ge0$:
\begin{equation}\label{eq:moment-bound-full}
\abs{K_p^{(\ell)}(s)}
\le w\sum_{j=1}^6M_jL_j(s)(\pi a_{j-1}^*)^\ell
+\abs C\sum_{j=0}^6L_j^P(s)(\pi a_j^*)^\ell,
\end{equation}
where
$L_j(s)=b_{j-1}^*e^{-\pi g_js}\min\{1/6,1/(\pi\beta_js)\}$,
$L_j^P(s)=m_j^*b_j^*e^{-\pi g_js}$, and the six tuple
$(b_j^*,a_j^*,g_j,M_j,\beta_j)$ is as in Table~\ref{tab:env}.
\end{lemma}
\begin{proof}
Starting from the folded measure of Lemma~\ref{lem:measure}: on piece
$j$, the integrand has factor $\lambda(t)(i\pi z(t))^\ell e^{i\pi z(t)s}$
times the folded density $\sum_n(c_nu_n(t)+d_nv_n(t))$. Bound each factor
by the corresponding entry of Table~\ref{tab:env}:
\[
\abs{\lambda(t)(i\pi z(t))^\ell e^{i\pi z(t)s}}
\le b_{j-1}^*(\pi a_{j-1}^*)^\ell e^{-\pi g_js}
e^{-\pi\beta_j(t_j-t)s}.
\]
The folded density has $\abs{u_n(t)},\abs{v_n(t)}\le M_j$, and the
triangle inequality over $n$ gives $w\cdot\max_j M_j$ overall. Integrate
$e^{-\pi\beta_j(t_j-t)s}$ over the piece of length $1/6$. The result is
$\min(1/6,1/(\pi\beta_js))$, and multiplying by
$b_{j-1}^*(\pi a_{j-1}^*)^\ell e^{-\pi g_js}$ gives $L_j(s)$.
The atom-only column (from $C$) uses the atom masses
$m^*=(5,4,4,4,2,4,2)$ and the factor $b_j^*(\pi a_j^*)^\ell e^{-\pi g_js}$,
i.e.\ $L_j^P(s)$. Sum over $j$.
\end{proof}
\subsection{Fej\'er quadrature}
\label{subsec:fejer}
The finite verification of Lemma~\ref{lem:residualcert} and
Lemma~\ref{lem:bernstein} rests on the following quadrature rule.
Set $N=256$. Define the cosine-sum weights
\begin{equation}\label{eq:fejer}
\alpha_v=\frac{v+1/2}{N},\quad
o_v=\frac1N\Bigl(1-2\sum_{a=1}^{N/2-1}
\frac{\cos(2a\pi\alpha_v)}{4a^2-1}\Bigr),
\quad
\tau_{jv}=\frac{2j-1+\cos(\pi\alpha_v)}{12}.
\end{equation}
Integrals over $(t_{j-1},t_j)$ are approximated by
$(1/6)\sum_{v=0}^{N-1}o_v f(\tau_{jv})$.
Properties (all exact and elementary):
\begin{enumerate}
\item Exact through degree $N-1=255$. Chebyshev orthogonality at the
discrete nodes $\theta_v=\pi\alpha_v$: substituting $x=\cos\theta$,
$\int_{-1}^1T_k(x)\dd x$ equals the discrete sum
$\sum_v o_vT_k(\cos\theta_v)$ for $0\le k\le255$. This follows because
the discrete cosine sums satisfy
$\sum_vo_v\cos(2k\pi\alpha_v)=1$ if $k=0$ and $0$ if $1\le k\le N-1$.
\item Positive weights: $o_v\ge\frac1N(1-\sum_{a=1}^{N/2-1}2/(4a^2-1))
=\frac1{N(N-1)}>0$.
\item Total mass: $\sum_vo_v=1$.
\end{enumerate}
\begin{lemma}[Quadrature moment error]
\label{lem:quad}
Let $p$ have a real coefficient list supported on $\{0\}\cup
\mathcal N\cap[1,100]$, with $\norm{(c,d)}_1\le3$, $|C|\le1$. Let
$m,y\in\RR$ and $\ell\in\{0,\ldots,30\}$ satisfy $0\le m\le100$ and
$|y|\le3/2$. Then the Fej\'er rule \eqref{eq:fejer}, applied to the
moments
\[
U_{p,\ell}(m,y)=\frac{y^\ell p^{(\ell)}(m)}{\ell!},\quad
X_{p,\ell}(m,y)=\frac{y^\ell K_p^{(\ell)}(m)}{\ell!},
\]
with exact evaluations at its nodes and at the atoms, commits error less
than $10^{-22}$ in each moment.
\end{lemma}
\begin{proof}[Sketch]
Each piece-$j$ density integrand, and each factor
$\lambda(t)e^{i\pi z(t)m}(i\pi z(t)y)^\ell/\ell!$, is analytic on the
disk of radius $\rho=1/6$ centered at $(2j-1)/12$, because the pole
$-ih$ has $\operatorname{Im}(-ih)=-2/5$, farther than $1/6$ from
$[0,1]$. A rational arithmetic computation bounds the moduli of the
integrand and its disk Cauchy-seminorms by $M_0<10^{53}$. The degree-255
Taylor tail is then $M_0 2^{-255}<3.455\cdot10^{-24}$. Adding the six
piece errors gives $<10^{-22}$. The atoms are exact.
\end{proof}
\subsection{Interval arithmetic verification}
\label{subsec:verifier}
The primary verifier \path{verification/verify_all.py} chains
\path{interval_certificate.py} and \path{interval_bounds.py} and writes
\path{finite.json} and \path{scalar.json} into \path{results/}. Each
script uses Arb ball arithmetic (via \path{python-flint}) at $192$ bits
of precision.
The comparisons made are:
\begin{itemize}
\item \emph{Finite certificate} (Lemma~\ref{lem:finitecert}): intervals
for $\norm{W_\eta(I-\eta D)}$ and $\norm{W_\eta}$, for
$\eta=\pm1$. Every enclosing ball must lie strictly below the
rational threshold ($10^{-3}$ and $32$).
\item \emph{Residual certificate} (Lemma~\ref{lem:residualcert}): 102
residual enclosures at every $n\in I_f$, threshold $3\cdot10^{-10}$.
\item \emph{Bernstein certificate} (Lemma~\ref{lem:bernstein}): 84
half-gaps x 29 $k$-values = 2436 individual Bernstein values, each
strictly above its threshold from \eqref{eq:Bk}. Sum total $2.32\cdot
10^{-12}$ is the propagated moment error.
\item \emph{Scalar substitutions}: 22 named rational comparisons: the
row-envelope cuts, the atom tail, the sign-table rows, the Taylor
representative cases, the disk integrand and quadrature bounds
(Tables~\ref{tab:env} and \ref{tab:bernstein}). All are
strict-rational-outward comparisons.
\end{itemize}
No step depends on \emph{observing} agreement between two floats. Each
acceptance is a machine-checked rational inequality; the ball-arithmetic
enclosure is the sole witness, and the algorithm stops when every
enclosure is strictly on the asserted side.
\subsection{The rational $N=64$ alternative}
\label{subsec:rational-small}
An independent verifier
\path{verification/check_sjj_rational.py} computes the $6\times6$ block
on $J=\{1,3,4\}$ with the Fej\'er rule at $N=64$, using exact rational
arithmetic throughout. Constants $\pi$ and $b=\sqrt3/2$ are externally
approximated by rationals $p_*$ and $b_*$ accurate to $10^{-36}$ (checked
via rational arctangent series and bisection of $\sqrt{3/4}$). The
rounded exponential $e_*(u)$ is constructed via
$v_0=[\sum_{q\le32}(ip_*u/1024)^q/q!]_{40}$, $v_{a+1}=[v_a^2]_{40}$
($a<10$), and error $|e_*(u)-e^{ip_*u}|<3.542\cdot10^{-34}$ whenever
$|p_*u|\le1000$ and $\Re(ip_*u)\le0$.
The rule is exact through degree $63$; the disk Cauchy bound gives
Taylor tail $M_0 2^{-62}<5.2\cdot10^{-15}$ in the value moment and
$M_1 2^{-62}<1.02\cdot10^{-13}$ in the derivative moment. The
propagated matrix-entry error is $<7\cdot10^{-14}<10^{-12}$ after
accounting for the rounding of the quadrature nodes and the moment
division by $Q_m>1.75$. This is checked as a separate independent path,
not as a substitute for the main Arb verifier.
%=====================================================================
\begin{thebibliography}{9}
%=====================================================================
\bibitem[Cohn--Elkies(2003)]{CohnElkies2003}
H.~Cohn and N.~Elkies.
\newblock New upper bounds on sphere packings {I}.
\newblock \emph{Annals of Mathematics}, 157(2):689--714, 2003.
\bibitem[Hales(2000)]{Hales2000}
T.~C. Hales.
\newblock Cannonballs and honeycombs.
\newblock \emph{Notices of the American Mathematical Society},
47(4):440--449, 2000.
\bibitem[Viazovska(2017)]{Viazovska2017}
M.~S. Viazovska.
\newblock The sphere packing problem in dimension 8.
\newblock \emph{Annals of Mathematics}, 185(3):991--1015, 2017.
\bibitem[Cohn--Kumar--Miller--Radchenko--Viazovska(2017)]{CKMRV2017}
H.~Cohn, A.~Kumar, S.~D. Miller, D.~Radchenko, and M.~Viazovska.
\newblock The sphere packing problem in dimension 24.
\newblock \emph{Annals of Mathematics}, 185(3):1017--1033, 2017.
\end{thebibliography}
\end{document}