\documentclass[11pt,letterpaper]{article}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage[textwidth=6.22in,top=.99in,bottom=.98in,headheight=15pt,headsep=17pt]{geometry}
\usepackage{amsmath,amsthm,newpxtext,newpxmath,microtype,booktabs,fancyhdr,needspace}
\usepackage[hidelinks]{hyperref}
\hypersetup{pdftitle={Quartic Fekete points on every real simplex},pdfauthor={Yunpeng Li},
  pdfsubject={Quartic Fekete configurations and cardinal polynomial inequalities on simplices},
  pdfkeywords={Fekete points, Fejér points, multivariate polynomial interpolation, optimal design}}
\pagestyle{fancy}
\fancyhf{}
\fancyhead[L]{\small Yunpeng Li}
\fancyhead[R]{\small Quartic Fekete points on every real simplex}
\fancyfoot[C]{\thepage}
\renewcommand{\headrulewidth}{.25pt}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Corollary}
\numberwithin{equation}{section}
\newcommand{\F}{\mathcal F}
\newcommand{\Fset}{\F_4}
\newcommand{\R}{\mathbb R}
\newcommand{\Q}{\mathbb Q}
\newcommand{\T}{\mathsf T}
\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator{\partn}{part}
\DeclareMathOperator{\supp}{supp}
\setlength{\emergencystretch}{1em}
\setlength{\abovecaptionskip}{4pt}
\setlength{\intextsep}{9pt plus 2pt minus 2pt}
\raggedbottom

\begin{document}
\thispagestyle{plain}
\begin{center}
  {\LARGE\bfseries Quartic Fekete points on every real simplex\par}
  \vspace{.8em}
  {\large\renewcommand{\thefootnote}{\fnsymbol{footnote}}%
   Yunpeng Li\footnote[1]{The Chinese University of Hong Kong, Shenzhen.
     \href{mailto:liyunpeng@cuhk.edu.cn}{\texttt{liyunpeng@cuhk.edu.cn}}}\par}
\end{center}

\begin{abstract}
We prove that on every real simplex Bos's quartic interpolation nodes
uniquely maximize the absolute Vandermonde determinant, up to ordering,
among all configurations of the prescribed size. In dimension four the squared Euclidean norm
of the cardinal vector exceeds one at the centroid. We bound its excess
by products of the cardinal polynomials for tetrahedral face centers.
These products also measure correlations between columns of the evaluation
matrix, which give the required loss in the determinant. An elementary
consequence of the same corrected squared-norm estimate is Bos's
fourth-power bound in every dimension, with equality exactly at the nodes.
The corrected estimate follows from cubic perturbations, reduction to
at most two distinct positive coordinate values, and coefficient monotonicity in their
multiplicities. All polynomial bounds are given explicitly. The common
base case is the tetrahedral Fej\'er bound previously stated by Panda and
Pal in equivalent optimal-design form. We give an independent proof of
that bound and determine all its equality points.
\end{abstract}

\noindent{\small\emph{2020 Mathematics Subject Classification.}
Primary 41A05. Secondary 41A63, 62K05.\par}
\noindent{\small\emph{Keywords.} Fekete points, Fej\'er points, polynomial interpolation,
simplex, optimal design.\par}

\section{The extremal configuration}\label{sec:intro}
The choice of interpolation nodes determines how an error in one datum
affects the interpolant elsewhere. Let $\ell_i$ be the cardinal polynomial
that is one at the $i$th node and zero at the others. An error $\varepsilon$
at that node changes the interpolated value at $x$ by $\varepsilon\ell_i(x)$.
There is a determinant interpretation of the same factor: replacing the
$i$th node by $x$ multiplies the Vandermonde determinant by $\ell_i(x)$.
Hence on a \emph{Fekete set}, one that maximizes the absolute determinant,
each $|\ell_i(x)|$ is at most one throughout the domain
\cite[Remark~0.8]{bos2023}.

On $[-1,1]$ the absolute determinant in the monomial basis is simply
$\prod_{i<j}|y_i-y_j|$. For degree $n\ge1$ its maximizing nodes are the
endpoints and the interior zeros of $P_n'$, where $P_n$ is the Legendre
polynomial of degree $n$ \cite[\S2]{taylorwingatevincent}.
Fej\'er proved that at these nodes the stronger
inequality $\sum_i\ell_i(x)^2\le1$ holds on the whole interval
\cite[pp.~264--265]{fejer1932}. This squared sum has a statistical
meaning: with independent errors of common variance $\sigma^2$ at the
nodes, the variance of the interpolated value at $x$ is
$\sigma^2\sum_i\ell_i(x)^2$.

On a simplex, Bos constructed families of nodes that agree with the
lower-dimensional construction on each face \cite{bos1983,bos2023}.
For quadratics the vertices and edge midpoints satisfy the Fej\'er
inequality in every dimension \cite[Proposition~2.1]{bos2023}.
For higher degrees, Taylor, Wingate and Vincent computed candidate
Fekete sets on the triangle through degree~19 \cite{taylorwingatevincent}.
Here we fix degree four and allow the dimension to vary. We prove that
Bos's nodes uniquely maximize the determinant even among configurations
with no prescribed symmetry or face structure.
Section~\ref{sec:higher-degree} examines which steps extend in the
polynomial degree and identifies the interior determinant to maximize
on a triangle with prescribed boundary nodes.

An invertible affine map preserves the extremal problem, so we use
barycentric coordinates on a real $d$-simplex and write
\[
 m=d+1,\qquad s=\sum_{i=1}^m x_i,\qquad
 \Delta_d=\{x\in\R^m:x_i\ge0,\ s=1\},\qquad N=\binom{m+3}{4}.
\]
Let $\mathcal P_4(\Delta_d)$ be the $N$-dimensional space of real
polynomials of total degree at most four on $\Delta_d$. For a basis
$b=(b_1,\ldots,b_N)$ of this space, let
$V_b(Y)=(b_j(y_i))_{i,j=1}^N$ for $Y=(y_1,\ldots,y_N)\in\Delta_d^N$.
A change of basis multiplies $|\det V_b(Y)|$ by a constant independent of $Y$.

Put $r=\sqrt5$, $\eta=\sqrt{3/7}$, $\alpha=(4+r)/11$ and
$\beta=(7-r)/22$. Bos's set $\F_m$ is obtained from the following
representatives:
\begin{equation}\label{eq:nodes}
 (1),\quad (1/2,1/2),\quad
 \bigl((1+\eta)/2,(1-\eta)/2\bigr),\quad
 (\alpha,\beta,\beta),\quad (1/4,1/4,1/4,1/4).
\end{equation}
Retain the vectors of length at most $m$, append zeros to length $m$,
and take all coordinate permutations.
The vertices, edge nodes, triangular face nodes, and tetrahedral face centers
give $m+3\binom m2+3\binom m3+\binom m4=N$ points in all.

\begin{theorem}[Quartic Fekete theorem]\label{thm:fekete}
For every $d\ge1$, every basis $b$ of $\mathcal P_4(\Delta_d)$, and every
$Y\in\Delta_d^N$,
\[
 |\det V_b(Y)|\le |\det V_b(\F_{d+1})|,
\]
where the nodes on the right may be ordered arbitrarily. Equality holds
if and only if $Y$ is an ordering of $\F_{d+1}$.
\end{theorem}

The maximality assertion is Bos's Conjecture~4.3 \cite[p.~285]{bos2023}.
For $z\in\F_m$ write
$L_{m,z}$ for the quartic cardinal polynomial, equal to one at $z$ and
zero at every other node. Lemma~\ref{lem:restriction} proves unisolvence.
In this basis the Vandermonde matrix at the nodes is the identity.
A sufficient condition for determinant maximality is the Fej\'er bound
$S_m(x)=\sum_{z\in\F_m} L_{m,z}(x)^2\le1$, studied in several variables
by Bos \cite{bos1990}. Indeed, Hadamard's determinant inequality bounds
the absolute determinant of a square matrix by the product of the
Euclidean lengths of its rows \cite[pp.~505--506]{hornjohnson}.
For the matrix $M=(L_{m,z}(y_i))_{i,z}$ this gives
\[
 |\det M|\le\prod_{i=1}^{N}\sqrt{S_m(y_i)}.
\]
For nonzero rows, equality in Hadamard's inequality holds precisely when
the rows are pairwise orthogonal. If $S_m\le1$ throughout
the simplex, the nodal determinant is maximal, and $|\det M|<1$ whenever
some $S_m(y_i)<1$. Bos proved the quartic bound on the triangle
\cite{bos1983} and posed the tetrahedral case as Conjecture~4.2
\cite[p.~285]{bos2023}. Panda and Pal subsequently stated the equivalent
four-ingredient design bound \cite[\S3, Remark~2.1]{pandapal}.

\begin{theorem}[Tetrahedral Fej\'er bound and its equality cases]\label{thm:tetra}
On $\Delta_3$ one has
\begin{equation}\label{eq:main}
 S_4(x)=\sum_{z\in\F_4}L_{4,z}(x)^2\le1 .
\end{equation}
Equality holds precisely at the 35 points of $\F_4$.
\end{theorem}

In the mixture model a design is a probability measure $\mu$, and
$D$-optimality means maximizing the determinant of its information
matrix $G_\mu=\int b(x)^{\T}b(x)\,d\mu(x)$, with $b$ viewed as a row.
For equal masses on $N$ points $Y$, its determinant is
$N^{-N}|\det V_b(Y)|^2$. Thus Fekete maximality compares these
equally weighted configurations, while $D$-optimality compares all
probability measures \cite[Propositions~0.5--0.6]{bos2023}.
Panda and Pal prove the three-ingredient case in Theorem~3.1 of
\cite{pandapal}. In the four-ingredient case they specify equal masses
$1/35$ at precisely the nodes $\F_4$. Their edge parameter is
$(7-\sqrt{21})/14=(1-\eta)/2$ and their face parameters are
$\beta,\beta,1-2\beta=\alpha$. They state $D$-optimality for this design
and omit its proof \cite[\S3, Remark~2.1]{pandapal}. In the cardinal
basis its information matrix is $I_{35}/35$. The Kiefer--Wolfowitz
equivalence theorem says that a design with nonsingular $G_\mu$ is
$D$-optimal precisely when its sensitivity $b(x)G_\mu^{-1}b(x)^{\T}$
is at most $N$ at every $x$ \cite{kieferwolfowitz}. Here the sensitivity
is $35S_4(x)$ and $N=35$, giving exactly \eqref{eq:main}. Section~\ref{sec:tetra}
gives an independent proof and determines the equality set, from which
uniqueness of the optimal measure follows.

The first obstruction to the same determinant argument occurs in
dimension four. At the centroid $c=(1/5,\ldots,1/5)$, Bos's calculation
\cite[\S4]{bos2023} gives
\begin{equation}\label{eq:centroid}
 S_5(c)=\frac{9298842+80243\sqrt5}{9453125}>1.
\end{equation}
Thus the uniform design on these 70 nodes fails the equivalence
criterion in the five-ingredient case proposed in
\cite[Remark~2.2 and Table~2]{pandapal}. Theorem~\ref{thm:fekete} gives a unique
maximizer among all 70-point configurations. The proof bounds the excess
$S_m-1$ pointwise by products of distinct face-center cardinals. The
same products force a loss in the determinant of their column Gram
matrix (Section~\ref{sec:det}). Section~\ref{sec:certificate} proves
the pointwise bound in all dimensions by cubic perturbations and uniform
algebraic inequalities in the coordinate multiplicities.

Bos also conjectured a uniform pointwise bound for the fourth powers.

\begin{theorem}[Fourth-power bound]\label{thm:four}
For every $d\ge1$, $m=d+1$, and $x\in\Delta_d$,
\[
 \sum_{z\in\F_m}L_{m,z}(x)^4\le1.
\]
Equality holds if and only if $x\in\F_m$.
\end{theorem}

This is Bos's Conjecture~4.6 \cite[p.~286]{bos2023}, with its equality cases.
Section~\ref{sec:four} derives it from the same corrected squared-norm
bound used in the determinant proof. Retaining the center cross terms
in $S_m^2$ absorbs the correction and gives all equality cases.
The fourth-power bound alone does not give the determinant
estimate. For instance, the matrix of order three with diagonal $-1/2$
and off-diagonal $3/4$ has row sums one and fourth-power row sums
$89/128$, yet its determinant is $25/16$.

\section{From cardinal polynomials to determinants}\label{sec:det}
\subsection{Restriction to faces and the row bound}
We represent each polynomial of degree at most four on $s=1$ by its
unique homogeneous quartic on $\R^m$. The monomials $x^\nu$ with
$|\nu|=4$ form a basis. We use the same notation $L_{m,z}$ for the
homogeneous representative.

\begin{lemma}\label{lem:restriction}
The set $\F_m$ is unisolvent. On a coordinate face the cardinals for
nodes on that face restrict to its cardinal basis, and the other
cardinals restrict to zero.
\end{lemma}
\begin{proof}
Suppose a homogeneous quartic vanishes at every node. Its restriction to
an edge vanishes, since the edge contains five distinct nodes. On a
triangular face the three coordinates divide the restricted polynomial,
so their product divides it. The remaining linear factor vanishes at
the three interior nodes. These are noncollinear because $\alpha\ne\beta$,
and therefore the restriction vanishes on the whole triangle. On a
tetrahedral face all four coordinates divide the restriction, which is
therefore a constant times their product. The centroid makes that
constant zero. These divisibility statements use homogeneity to extend
vanishing on a face to its coordinate subspace.

Every quartic monomial uses at most four coordinates. Its coefficient
can therefore be read on a coordinate subspace where the polynomial
vanishes by the preceding argument. All coefficients are zero. The node
count in Section~\ref{sec:intro} now gives unisolvence. On a face the
nodal values determine each restriction uniquely, and homogeneity extends
this rule to its entire coordinate subspace.
\end{proof}

For a four-element set $J$, Bos gives the cardinal at the center of the
corresponding face as $q_J(x)=256\prod_{i\in J}x_i$
\cite[Proposition~4.4]{bos2023}. Its nodal values also follow by direct
substitution. Put
\[
 S_m=\sum_{z\in\F_m}L_{m,z}^2,\qquad
 G_m=s^8-S_m,\qquad
 E_m=\sum_{J\ne K}q_J^3 q_K^3 ,
\]
where the last sum is over \emph{ordered} pairs of four-element sets.
The center products in $E_m$ vanish whenever at most four coordinates
are positive. At other points they bound the excess of $S_m$ above one
and enter the determinant estimate with the opposite sign.

\begin{lemma}\label{lem:centers}
For $x\in\Delta_{m-1}$, $\sum_J q_J(x)^2\le1$.
Equality holds precisely at tetrahedral face centers.
\end{lemma}
For $m\ge4$ this is a special case of the uniform-profile theorem of
Orlitsky, Santhanam, Viswanathan and Zhang \cite[Theorem~14]{orlitsky}.
In eight independent draws with outcome probabilities $x_i$, the
probability of the pattern $11223344$ (four distinct outcomes, each
drawn twice in the indicated order) is $4!\sum_{|J|=4}\prod_{i\in J}x_i^2$.
Their theorem gives the uniform distribution on four outcomes as maximizer,
since its support-size criterion holds already at four, $(5/4)^8/5>1$.
We give a direct proof in barycentric coordinates.
\begin{proof}
For $m<4$ the sum is empty. For $m\ge4$ it equals
$4^8 e_4(x_1^2,\ldots,x_m^2)$, where $e_k$ denotes the $k$th
elementary symmetric polynomial. If its maximum is positive, at least
four coordinates are positive there. Fix two positive coordinates $a,b$,
their sum $v$, and the remaining coordinates $w$. With $u=ab$ the
unscaled expression is
$e_4(w^2)+(v^2-2u)e_3(w^2)+u^2 e_2(w^2)$.
It is strictly convex on $0\le u\le v^2/4$ since $e_2(w^2)>0$.
Unequal positive $a,b$ would give an interior maximum. Thus the $k\ge4$
positive coordinates all equal $1/k$. The resulting value is
$c_k=\binom{k}{4}(4/k)^8$, with $c_4=1$ and
$c_{k+1}/c_k=k^8/((k-3)(k+1)^7)<1$. The ratio bound is immediate at
$k=4,5$ and follows for $k\ge6$ from
$(1-3/k)(1+1/k)^7\ge(1-3/k)(1+7/k)>1$. Thus equality requires four
equal positive coordinates, and each such point gives equality.
\end{proof}

\begin{lemma}\label{lem:correction}
On $\Delta_{m-1}$ one has $S_m\le1+E_m/16$. The inequality is strict away
from $\F_m$.
\end{lemma}
Section~\ref{sec:certificate} proves the lemma by cubic perturbations
and explicit algebraic inequalities, uniformly in the dimension. On faces of
dimension at most three there is at most one nonzero $q_J$, and the lemma
reduces by restriction to the tetrahedral Fej\'er bound.

\subsection{The 70-point case}
There are five four-element sets when $m=5$. Set
$T(x)=\sum_{J\ne K}q_J(x)q_K(x)$. Cauchy--Schwarz gives
$T=(\sum_J q_J)^2-\sum_J q_J^2\le4\sum_J q_J^2\le4$ by
Lemma~\ref{lem:centers}. Also $E_m\le T^2/4$.
Indeed, each product $q_Jq_K$ is at most $1/2$, and the sum of the
squares of the ordered products is at most $T^2/2$ because every
unordered product occurs twice. Lemma~\ref{lem:correction} gives
$S_m\le1+T^2/64$.

For an arbitrary 70-point configuration let $M$ denote its evaluation
matrix in the cardinal basis, with the five center columns last. Take
$R=\max_i T(y_i)$ and $t=R/32$. Use the inner product on row vectors
given by $A_t=\operatorname{diag}(I_{65},K_t)$, where $K_t$ has diagonal one and
off-diagonal $-t$. The $i$th diagonal entry of $MA_tM^{\mathsf T}$ is
\[
 S_m(y_i)-tT(y_i)\le1+\frac{T(y_i)(T(y_i)-2R)}{64}\le1.
\]
For a row attaining $R$, the middle expression equals $1-16t^2$. The eigenvalues
of $K_t$ are $1+t$ with multiplicity four and $1-4t$. As $0\le t\le1/8$,
\[
 \det K_t=(1+t)^4(1-4t)
 =1-t^2(10+20t+15t^2+4t^3)\ge1-13t^2 .
\]
Thus $A_t$ is positive definite. In its Gram-matrix form, Hadamard's
inequality gives $\det G\le\prod_i G_{ii}$ for every positive
semidefinite $G$ \cite[p.~505]{hornjohnson}. Applying it to
$G=MA_tM^{\T}$ gives
$|\det M|^2\det K_t=\det(MA_tM^{\mathsf T})\le1-16t^2$, proving
$|\det M|<1$ if $R>0$. If $R=0$, every point has a zero coordinate,
since all five $q_J$ are positive at every interior point. By face
restriction and Theorem~\ref{thm:tetra}, all rows then have norm at
most one. A point outside $\F_5$ gives a row of norm strictly below one,
so the product of the row norms is less than one. Equality therefore
requires 70 distinct nodes.

\subsection{A determinant estimate in arbitrary dimension}
For an $n\times n$ real matrix $M=[U\ Q]$ let $q_{ij}$ be the entries of
$Q$, with $p$ columns, and write $M_{i*}$ for row $i$. The sums over
$j\ne k$ below count ordered pairs of columns of $Q$.
\begin{lemma}\label{lem:matrix}
Suppose $q_{ij}\ge0$ and $\sum_jq_{ij}^2\le1$ for every row. If $M$ is
nonsingular, then
\begin{equation}\label{eq:det}
 \log |\det M|^2\le
 \sum_i\left(\|M_{i*}\|_2^2-1
       -\frac1{16}\sum_{j\ne k}q_{ij}^3q_{ik}^3\right).
\end{equation}
\end{lemma}
\begin{proof}
For a positive definite $k\times k$ matrix $X$, summing
$\log v\le v-1$ over its eigenvalues gives $\log\det X\le\tr X-k$.
This proves the case $p=0$ by taking $X=M^{\T}M$.
For $p>0$, $H=Q^{\T}Q$ is positive definite since $M$ is nonsingular.
Put $D(H)=\tr H-p-\log\det H$. This is Stein's covariance loss with the
identity as reference matrix \cite[eq.~(72)]{jamesstein}. Let
$\Pi=QH^{-1}Q^{\mathsf T}$ be the orthogonal projection onto the columns
of $Q$. Eliminating the $H$ block from $M^{\T}M$ leaves the Schur
complement $U^{\T}U-U^{\T}QH^{-1}Q^{\T}U=U^{\T}(I-\Pi)U$.
The block determinant formula \cite[eq.~(0.8.5.1)]{hornjohnson} gives
$|\det M|^2=\det H\,\det(U^{\T}(I-\Pi)U)$, with the empty
determinant equal to one if $p=n$. Projecting the columns of $U$ with
$I-\Pi$ can only shorten them, so the trace bound on the Schur
complement gives
\[
 \log|\det M|^2\le \tr(U^{\mathsf T}U)-(n-p)+\log\det H
 =\sum_i(\|M_{i*}\|_2^2-1)-D(H).
\]
Write $h_j=H_{jj}>0$ and $C_{jk}=H_{jk}/\sqrt{h_jh_k}$. The matrix $C$
is positive definite with diagonal one. Normalizing the columns of $Q$
gives the identity
\begin{equation}\label{eq:loss}
 D(H)=\sum_j(h_j-1-\log h_j)-\log\det C .
\end{equation}
The scalar inequality gives $h_j-1-\log h_j\ge0$. Hadamard's inequality
gives $\det C\le1$ because $C_{jj}=1$. We use these two nonnegative
terms to bound sums over disjoint sets of ordered pairs.

For $j\ne k$ let $w_{jk}=\sum_i q_{ij}^2q_{ik}^2$ and
$B_{jk}=w_{jk}/\sqrt{h_jh_k}$. Set $w_{jj}=B_{jj}=0$. The row hypothesis
gives $\sum_k w_{jk}=\sum_i q_{ij}^2\sum_{k\ne j}q_{ik}^2\le h_j$.
For every $v\in\R^p$,
\[
 |v^{\mathsf T}Bv|
 \le\sum_{j<k}w_{jk}(v_j^2/h_j+v_k^2/h_k)\le\|v\|_2^2 .
\]
Thus the spectrum of $B$ lies in $[-1,1]$. Also $0\le B_{jk}\le C_{jk}$
entrywise because $0\le q_{ij}q_{ik}\le1$. Hence
$\tr(B^2)\le\tr(BC)$ by summing the entrywise products. For every
positive definite $T$, applying $\log\det X\le\tr X-p$ to
$X=C^{1/2}TC^{1/2}$ gives
\[
 -\log\det C\ge p-\tr(CT)+\log\det T .
\]
Equality holds at $T=C^{-1}$, giving the variational formula for the
log determinant \cite[Example~3.23]{boydvandenberghe}. Choose
$T=I-B/2$, which is positive definite by the spectral bound on $B$.
Since $\tr C=p$,
\(\tr B=0\), and $\log(1-u)\ge-u-u^2$ for $|u|\le1/2$, we obtain
\begin{equation}\label{eq:tracebound}
 -\log\det C\ge\tfrac12\tr(BC)+\log\det(I-B/2)
 \ge\tfrac12\tr(BC)-\tfrac14\tr(B^2)\ge\tfrac14\tr(BC).
\end{equation}

Write $\sum_i\sum_{j\ne k}q_{ij}^3q_{ik}^3=E_{\rm lo}+E_{\rm hi}$. In
$E_{\rm lo}$ both $h_j,h_k$ are at most two. For each such pair,
\[
 B_{jk}C_{jk}
 =\frac{(\sum_i q_{ij}^2q_{ik}^2)(\sum_i q_{ij}q_{ik})}{h_jh_k}
 \ge\tfrac14\sum_i q_{ij}^3q_{ik}^3 .
\]
All cross terms are nonnegative. Summing and using \eqref{eq:tracebound}
gives $-\log\det C\ge E_{\rm lo}/16$.

For $E_{\rm hi}$, use $v-1-\log v\ge v/8$ for $v\ge2$, by
differentiation and $\log2<3/4$. For each row and $j$,
\[
 q_{ij}\sum_{k\ne j}q_{ik}^3
 \le q_{ij}(1-q_{ij}^2)^{3/2}\le\tfrac12.
\]
The indicator that $h_j>2$ or $h_k>2$ is at most
$\mathbf1_{h_j>2}+\mathbf1_{h_k>2}$. By symmetry in $j,k$,
\[
 E_{\rm hi}\le
  2\sum_{j:h_j>2}\sum_i q_{ij}^3\sum_{k\ne j}q_{ik}^3
  \le\sum_{j:h_j>2}h_j\le8\sum_j(h_j-1-\log h_j).
\]
Substitution into \eqref{eq:loss} gives
$D(H)\ge E_{\rm lo}/16+E_{\rm hi}/8\ge(E_{\rm lo}+E_{\rm hi})/16$,
as required.
\end{proof}

\begin{proof}[Proof of Theorem~\ref{thm:fekete}]
Fix any $Y=(y_1,\ldots,y_N)$ and set $M_{iz}=L_{m,z}(y_i)$. Interpolation
of each $b_j$ gives $V_b(Y)=M V_b(\F_m)$. For $m\le4$,
Theorem~\ref{thm:tetra}, face restriction and Hadamard's inequality
give $|\det M|\le1$, strictly if any point lies outside $\F_m$.
For $m\ge5$ take $Q$ to consist of the tetrahedral face-center columns. Singular
matrices have zero determinant. For any nonsingular matrix the row
hypothesis of Lemma~\ref{lem:matrix} follows from Lemma~\ref{lem:centers}.
Lemma~\ref{lem:correction} makes every summand in \eqref{eq:det}
nonpositive, strictly so for any point outside $\F_m$. Thus
$|\det M|\le1$. The nodal configuration gives the identity matrix.
Equality requires distinct nodal rows, hence precisely the prescribed set.
\end{proof}

\section{The corrected squared-norm bound}\label{sec:certificate}
Homogenize the assertion of Lemma~\ref{lem:correction} by setting
\[
 \Phi_m=s^{16}G_m+E_m/16,\qquad
 D_m=\sum_{i=1}^m\partial_{x_i},\qquad
 R_m=sD_mG_m+16mG_m .
\]
The degree of $\Phi_m$ is 24 and that of $R_m$ is eight. For $x\ge0$ with
$s>0$, homogeneity gives
\[
 \Phi_m(x)=s^{24}\bigl[1-S_m(x/s)+E_m(x/s)/16\bigr].
\]
Thus Lemma~\ref{lem:correction} is equivalent to nonnegativity of $\Phi_m$
on the orthant together with its stated equality set on $s=1$. If
$x_m=0$, the cardinal restriction rule gives $G_m(x',0)=G_{m-1}(x')$,
and every term of $E_m$ using the last coordinate vanishes. In particular
$\Phi_m(x',0)=\Phi_{m-1}(x')$. We prove strict positivity in the open
orthant for $m\ge5$ and recover the remaining points by face restriction.
Our induction follows segments from a facet to the centroid after
normalization to $s=1$. Bos used the same segments and raised the
degree of a polynomial to obtain nonnegative coefficients. This bounded
the partial sum omitting the tetrahedral face centers in dimensions
$3\le d\le7$ \cite[proof of Proposition~4.4]{bos2023}. Here we control
the derivative of the corrected polynomial $\Phi_m$ along the segment.

\paragraph{The input polynomial.}
For a nonnegative integer vector $v$ let $\partn(v)$ be its positive
entries in decreasing order. Symmetry and face restriction show that
the coefficient of $x^v$ in $363G_m$, for $|v|=8$, depends only on
$\lambda=\partn(v)$. Table~\ref{tab:input} records this ordinary
monomial coefficient $g_\lambda=a+b\sqrt5$. Multiplicities are
accounted for by summing over distinct monomials, with no additional
factor on an individual coefficient.

\begin{table}[!ht]
\centering\small
\caption{The 22 coefficients of $363G_m$. Digits concatenate the parts of
$\lambda$. Each pair represents $a+b\sqrt5$.}\label{tab:input}\smallskip
\begin{tabular}{@{}lrr@{\hskip 2.4em}lrr@{}}
\toprule
$\lambda$ & $a$ & $b$ & $\lambda$ & $a$ & $b$ \\
\midrule
8     &0&0                 &41111    &25065&$-174285$\\
71    &7260&0              &332      &4016166&$-246554$\\
62    &$-87120$&0          &3311     &547658&607874\\
611   &$-14421$&$-1881$    &3221     &2858902&$-208510$\\
53    &348480&0            &32111    &2810280&$-23430$\\
521   &211772&70004        &311111   &3671370&421110\\
5111  &285549&35541        &2222     &$-33498204$&732492\\
44    &$-537240$&0         &22211    &$-4607622$&$-66090$\\
431   &192445&$-115995$    &221111   &956922&$-383826$\\
422   &$-4545497$&313841   &2111111  &6116790&$-396270$\\
4211  &$-701079$&$-141739$ &11111111 &14636160&0\\
\bottomrule
\end{tabular}
\end{table}

The explicit cardinals in Appendix~\ref{sec:cardinals} give the table by
squaring and collecting monomials. The two non-midpoint edge cardinals
are conjugate under $\sqrt{21}\mapsto-\sqrt{21}$, so their squared sum
lies in $\Q(\sqrt5)$. This construction uses the five cardinal formulas
and face restriction in every dimension.

We also need the octic in elementary symmetric functions. Write
$e_j=e_j(x)$, $e_0=1$, and $e_1=s$. Let
\[
 \mathcal H(s,e_2,e_3,e_4)
   =s^8F(e_2/s^2,e_3/s^3,e_4/s^4),
\]
where $F$ is the explicit tetrahedral polynomial in \eqref{eq:F}; the
right side is polynomial after cancellation. The stable polynomial
$\mathcal G$ representing $G_m$ for every $m$ is
\begin{align}
363(\mathcal G-\mathcal H)={}&(-3885987+728963r)s^3e_5\notag\\
&+(17210374-2817998r)se_2e_5
 +(19321502+3933518r)e_3e_5\notag\\
&+(-8761898+3618862r)s^2e_6
 +(-4646668-8052364r)e_2e_6\notag\\
&+(7995164+1003412r)se_7
 +(-10666112+2914432r)e_8.\label{eq:stable}
\end{align}
Equation~\eqref{eq:stable} follows by expanding the elementary functions
and comparing with Table~\ref{tab:input}. Set $e_j=0$ for $j>m$ when
evaluating in $m$ coordinates. Since $D_me_j=(m-j+1)e_{j-1}$, define
\begin{equation}\label{eq:Rformal}
 \mathcal R_m=s\sum_{j=1}^8(m-j+1)e_{j-1}
                         \frac{\partial\mathcal G}{\partial e_j}
                    +16m\mathcal G,
 \qquad \mathcal B=\partial_m\mathcal R_m.
\end{equation}
Substitution of the actual elementary functions gives $R_m$.
When a coordinate vector is an argument of $\mathcal R_m$ or
$\mathcal B$, its elementary functions are understood to be substituted.
The formula is affine in the formal parameter $m$. On a face, this
parameter remains the ambient number of coordinates; it is not replaced
by the size of the positive support.

\paragraph{Varying three coordinates.}\label{par:cubic-variation}
Take three distinct positive coordinates with sum $u$, pairwise sum $v$
and product $k$. They are the roots of $z^3-uz^2+vz-k$. Choose disjoint
intervals about the roots, each with positive endpoints at which the
cubic has opposite signs. These signs persist for all sufficiently
small changes of $k$, so the perturbed cubic has one root in each
interval. The intervals can be arbitrarily small. Thus we can vary $k$
in either direction, keeping $u,v$ and all other coordinates fixed,
within the same face. This is the real-root perturbation used by Riener
\cite[Corollary~3.2(4) and proof of Proposition~3.3]{riener}.

\paragraph{Five barycentric coordinates.}
The correction in this case has the useful factorization
\begin{equation}\label{eq:E-five}
 E_5/16=2^{45}(x_1x_2x_3x_4x_5)^3e_3(x_1^3,\ldots,x_5^3).
\end{equation}
Indeed, two different four-element sets omit two different coordinates;
each unordered pair occurs twice in $E_5$.

\begin{lemma}\label{lem:five-analytic}
The polynomial $\Phi_5$ is positive at every point of the open orthant.
\end{lemma}
\begin{proof}
Work on $s=1$. We first exclude a local minimum having three distinct
positive coordinates. Hold the other two coordinates fixed, and write
their sum and product as $a,b$. For the selected three coordinates write
$u,v,k$ for their sum, pairwise sum and product. Thus $u=1-a$.
We vary $k$ in the two-sided interval just constructed.

Under the variation $e_2$ is fixed, while the derivatives of
$e_3,e_4,e_5$ with respect to $k$ are $1,a,b$. The terms quadratic in
these variables in $G_5$, obtained from Table~\ref{tab:input}, are
\[
 A(e_2)e_3^2+B e_3e_4+C e_4^2+D e_3e_5,
\]
where $A,B,C$ are given in \eqref{eq:F} and
$D=(19321502+3933518r)/363$. All remaining terms are constant or
linear along the variation. Consequently
$G_5''/2=A(e_2)+aB+a^2C+bD$. Since $A'>0$, $A'+D>0$,
$e_2=v+au+b\le u^2/3+au+b$, and $b\le a^2/4$, we obtain
\begin{align*}
 \tfrac12G_5''&\le
 \frac{-100806811+14769605r}{2178}a^2
 +\frac{17876302-8609402r}{1089}a
 +\frac{-376852+90950r}{99}\\
 &<-31000a^2-1200a-1700.
\end{align*}
The last bounds use $223/100<r<9/4$.

It remains to control the second derivative of the correction.
Set $d_3=a^3-3ab\ge0$. The sum of cubes and the sum of the pairwise
products of cubes of the selected triple are respectively
$u^3-3uv+3k$ and $v^3-3uvk+3k^2$. Equation~\eqref{eq:E-five} gives
\begin{align*}
 (E_5/16)''=2^{45}b^3\bigl[30k^4
 &+d_3(60k^3-36uvk^2+6v^3k)\\
 &+b^3(36k^2-18uvk+6u^3k)\bigr].
\end{align*}
Use $uv\ge9k$, $v\le u^2/3$ and $k\le u^3/27$. The term
$d_3(60k^3-36uvk^2)$ is nonpositive, and
$36k^2-18uvk+6u^3k\le-126k^2+6u^3k\le u^6/14$.
Moreover $b^3\le a^6/64$, $b^6\le a^{12}/4096$, and
$b^3d_3\le a^9/256$. For the last bound, the function
$b^3(a^3-3ab)$ increases for $0\le b\le a^2/4$.
It follows that
\[
 (E_5/16)''\le2^{45}\left[
 \frac{10a^6u^{12}}{64\cdot177147}
 +\frac{2a^9u^9}{256\cdot243}
 +\frac{a^{12}u^6}{4096\cdot14}\right]<39000a^2
\]
when $a>0$; the derivative is zero when $a=0$. To check the numerical
bound, divide by $a^2$ and maximize each product $a^j(1-a)^{16-j}$
at $a=j/16$. The three resulting upper bounds are
\[
 3840,\qquad \frac{1312993546389}{67108864}<19566,
 \qquad \frac{7119140625}{458752}<15519,
\]
whose sum is less than $39000$. Thus in all cases
\begin{equation}\label{eq:five-curvature}
 \Phi_5''<-23000a^2-2400a-3400<0.
\end{equation}

Every boundary point satisfies $\Phi_5\ge0$ by face restriction and
Theorem~\ref{thm:tetra}. A negative global minimum on the compact simplex
would therefore be interior. By \eqref{eq:five-curvature}, its coordinates
could have at most two distinct values. Such a point, including the
all-equal case, is a permutation of
$(\underbrace{A,\ldots,A}_{p},\underbrace{A+T,\ldots,A+T}_{5-p})$
with $A>0$, $T\ge0$, and $1\le p\le4$.
Appendix~\ref{sec:five-binary} gives explicit positive coefficient bounds
for all four binary forms, contradicting a negative minimum.
This proves nonnegativity everywhere. An interior zero would now be a
local minimum and would be excluded by the same argument. Homogeneity
extends strict positivity to the open orthant.
\end{proof}

\paragraph{A two-value reduction for the induction polynomial.}
\begin{lemma}\label{lem:two-global}
For $m\ge6$, every local minimum of $R_m$ relative to the interior of
any face of $\Delta_{m-1}$ has at most two distinct positive coordinate
values.
\end{lemma}
\begin{proof}
Choose three distinct positive coordinates, if they exist. Keep their
sum $u$ and pairwise sum $v$ fixed and vary their product $k$, as in
Lemma~\ref{lem:five-analytic}. This gives a two-sided curve in the same
face. Let $a,b$ be the elementary functions of degrees one and two in
the other coordinates. On $s=1$, $u=1-a$ and $e_2=v+au+b$ is fixed,
while the derivatives of $e_3,e_4,e_5,\ldots$ are $1,a,b,\ldots$.
Only the terms containing $e_3^2,e_3e_4,e_4^2,e_3e_5$ in
\eqref{eq:Rformal} contribute a second derivative. Their coefficients
have the form $A_*(e_2),B_*,C_*,D_*$. Thus
\[
 \tfrac12R_m''=A_*(e_2)+aB_*+a^2C_*+bD_*.
\]
Here $A_*'>0$ and $A_*'+D_*>0$. Using
$e_2\le(1-a)^2/3+a(1-a)+b$ and $b\le a^2/2$ gives
\begin{align}
\tfrac12R_m''\le{}&
 \frac{-462502064+160573264r}{1089}ma^2\notag\\
&+\left(\frac{246084376-126288176r}{1089}m
       +\frac{29150056-23606360r}{363}\right)a\notag\\
&+\frac{-56740394+9398698r}{1089}m
 +\frac{-17876302+8609402r}{363}\label{eq:global-curvature}\\
<{}&-90000ma^2-(30000m+60000)a-32000m+5000<0.\notag
\end{align}
The coefficient bounds use $223/100<r<9/4$; the last inequality holds
for $m\ge6$. A strictly negative second derivative contradicts a local
minimum on the admissible curve.
\end{proof}

Write $P\succeq Q$ when $P-Q$ has nonnegative coefficients in the
displayed polynomial variables. Such an inequality implies $P\ge Q$
when those variables are nonnegative. In a vector argument, $A^p$
will mean $p$ repeated entries equal to $A$.

\paragraph{One larger coordinate.}
\begin{lemma}\label{lem:one-large}
Let $p\ge1$ and $m\ge\max(6,p+1)$. For $A,T\ge0$ with $A+T>0$,
\[
 R_m(A^p,A+T,0^{m-p-1})>0.
\]
\end{lemma}
\begin{proof}
The coefficient bounds in Appendix~\ref{sec:uniform-algebra} give
\begin{align}
 \frac{363}{p}\mathcal R_{p+1}(A^p,A+T)
 &\ge10^5pAT^3\Psi(pA,T)+8712A^8+7260T^8 &&(p\ge5),\notag\\
 363\mathcal B(A^p,A+T)
 &\ge10^5pAT^3\Psi(pA,T)+8712A^8+7260T^8 &&(p\ge3),
 \label{eq:one-lower}
\end{align}
where
\[
 \Psi(Z,T)=Z^2(5Z^2-7ZT+3T^2)+T^2(7Z^2-ZT+T^2).
\]
Both quadratic forms are positive definite, with discriminants $-11$
and $-27$. The right sides of \eqref{eq:one-lower} are therefore positive.
For the remaining multiplicities, put $\sigma=(p+1)A+T$. We multiply
by a power of this positive linear form to test the coefficients, as
in P\'olya's theorem \cite{powersreznick}. Appendix~\ref{sec:uniform-algebra}
checks the exponent three and gives
\begin{align}
 \sigma^3\mathcal R_6(A^p,A+T)&\succeq4\sigma^{11}
                                      &&(p=1,2,3,4),\notag\\
 \sigma^3\mathcal B(A^p,A+T)&\succeq6\sigma^{11}
                                      &&(p=1,2).\label{eq:one-small}
\end{align}
For $p\ge5$, use
$\mathcal R_m=\mathcal R_{p+1}+(m-p-1)\mathcal B$.
For $p\le4$, use $\mathcal R_m=\mathcal R_6+(m-6)\mathcal B$.
The preceding positive bounds prove both cases.
\end{proof}

\paragraph{At least two larger coordinates.}
We interpolate the multiplicities by polynomials. For real parameters
$n,q$ define
\begin{equation}\label{eq:counts}
 E_j(n,q;A,T)=\sum_{h=0}^j\binom qh\binom{n-h}{j-h}A^{j-h}T^h,
 \qquad E_0=1.
\end{equation}
The binomial coefficients here are falling-factorial polynomials.
For integers $n\ge q\ge0$, $E_j$ is the elementary function on
$(A^{n-q},(A+T)^q)$. This follows by choosing the $h$ factors contributing
$T$ from the $q$ larger coordinates, and the remaining $j-h$ factors
contributing $A$ from the other $n-h$ coordinates. Set
\begin{equation}\label{eq:CB}
 C(n,q;A,T)=363\mathcal R_n(E_1,\ldots,E_8),\qquad
 B(n,q;A,T)=363\mathcal B(E_1,\ldots,E_8).
\end{equation}
The notation on the right means substitution into the elementary
variables of \eqref{eq:Rformal}. In particular, if $p+q=n\le m$,
\begin{equation}\label{eq:ambient}
 363R_m(A^p,(A+T)^q,0^{m-n})=C(n,q;A,T)+(m-n)B(n,q;A,T).
\end{equation}
There are two different derivatives in this argument. The derivative
defining $\mathcal B$ in \eqref{eq:Rformal} holds the elementary
variables fixed. In contrast, $\partial_n C$ and $\partial_n B$ below
differentiate the complete polynomials after the substitution
\eqref{eq:counts}, keeping $q,A,T$ fixed. Thus all dependence of the
elementary functions on $n$ is included.

\begin{lemma}\label{lem:two-large}
For integers $p\ge1$, $q\ge2$, $m\ge\max(6,p+q)$, and $A,T\ge0$
with $A+T>0$, the left side of \eqref{eq:ambient} is positive.
\end{lemma}
\begin{proof}
Appendix~\ref{sec:uniform-algebra} proves the following coefficient
bounds. For $u,v\ge0$, put $Z=(1+u)A+(1+v)T$. Then
\begin{align}
 \partial_n C(6+u,2+v;A,T)&\succeq4000Z^8,\notag\\
 \partial_n B(6+u,2+v;A,T)&\succeq69696AZ^7.
 \label{eq:monotone-counts}
\end{align}
The coefficient order is in $u,v,A,T$. For fixed $q\ge2$ and
nonnegative $A,T$, both $C(n,q;A,T)$ and $B(n,q;A,T)$ are consequently
nondecreasing in the real variable $n$ on $n\ge6$. More explicitly,
for $n\ge n_0\ge6$ and $P=C$ or $B$,
\begin{equation}\label{eq:count-integration}
 P(n,q;A,T)-P(n_0,q;A,T)
   =\int_{n_0}^n\partial_nP(t,q;A,T)\,dt\ge0.
\end{equation}
Only the polynomial extension in \eqref{eq:counts} is used at
noninteger $t$; no simplex of noninteger dimension is invoked.

The boundary family $n=q+1$ satisfies, for $w\ge0$,
\begin{align}
 C(6+w,5+w;A,T)&\succeq8712(1+w)^9(A+T)^8,\notag\\
 B(6+w,5+w;A,T)&\succeq8712(1+w)^8(A+T)^8.
 \label{eq:count-boundary}
\end{align}
Both right sides are positive. Finally, for
\[
 (n,q)\in\{(3,2),(4,2),(4,3),(5,2),(5,3),(5,4),
                         (6,2),(6,3),(6,4)\},
\]
the explicit bounds in the appendix give
\begin{align}
 C(n,q;A,T)+(6-n)B(n,q;A,T)&\succeq2000(nA+qT)^8,\notag\\
 B(n,q;A,T)&\succeq2000(nA+qT)^8.\label{eq:count-small}
\end{align}

Let $n=p+q$. If $n\le5$, write the right side of \eqref{eq:ambient} as
$C+(6-n)B+(m-6)B$ and apply \eqref{eq:count-small}.
Suppose $n\ge6$. If $q=2,3,4$, take $n_0=6$ in
\eqref{eq:count-integration} and use \eqref{eq:count-small}.
If $q\ge5$, then $p\ge1$ implies $n\ge q+1\ge6$; take
$n_0=q+1$ and use \eqref{eq:count-boundary} with $w=q-5$.
In this integration $q$ remains fixed: the boundary substitution
$n=6+w$, $q=5+w$ only supplies the initial value.
Both initial values are strictly positive whenever $A+T>0$,
including $A=0$ or $T=0$. By \eqref{eq:count-integration} they stay
positive even if a derivative vanishes. Hence $C,B>0$, and $m\ge n$ makes
\eqref{eq:ambient} positive.
This treats both unbounded multiplicities and the ambient dimension.
\end{proof}

We can now conclude that $R_m(x)>0$ for every $x\ge0$, $x\ne0$, and
$m\ge6$. Normalize to $s=1$ and take a minimizer. It lies in the
relative interior of the face determined by its positive coordinates, so
Lemma~\ref{lem:two-global} leaves at most two distinct positive values.
A support of size one gives $R_m(1,0,\ldots,0)=20(m-1)>0$.
All other possibilities are covered by Lemmas~\ref{lem:one-large}
and \ref{lem:two-large}; equal positive coordinates correspond to $T=0$.
Homogeneity removes the normalization $s=1$.

\begin{proof}[Completion of the proof of Lemma~\ref{lem:correction}]
For $m\le4$ the face restriction and Theorem~\ref{thm:tetra} give
$\Phi_m\ge0$, with equality on $s=1$ exactly at the nodes. The five-coordinate
argument in Lemma~\ref{lem:five-analytic} gives the next case. For $m\ge6$ we have
\[
 D_m\Phi_m=s^{15}R_m+\tfrac1{16}D_mE_m>0
 \qquad(0\ne x\ge0),
\]
since $E_m$ has nonnegative monomial coefficients in the original
variables. A zero coordinate reduces the assertion to $m-1$ coordinates.
For a point with all coordinates positive, put $t=\min_i x_i$,
$\mathbf1=(1,\ldots,1)$, and $y=x-t\mathbf1$. The face identity gives
$\Phi_m(y)\ge0$ by induction,
with $\Phi_m(0)=0$ if all coordinates of $x$ coincide. For $0<v\le t$ the
point $y+v\mathbf1$ is nonzero, so
\[
 \Phi_m(x)-\Phi_m(y)=\int_0^t D_m\Phi_m(y+v\mathbf1)\,dv>0 .
\]
Deleting zero coordinates reduces a point on $s=1$ to the positive
orthant on its support. Supports of size at least five give strict
positivity by the preceding argument. On smaller supports the
tetrahedral equality characterization gives exactly the nodes.
\end{proof}

\section{The tetrahedral Fej\'er bound}\label{sec:tetra}
We now prove Theorem~\ref{thm:tetra}, the common base case for the two
higher-dimensional results. Timofte's half-degree principle tests a
symmetric octic on the nonnegative orthant at points with at most four
distinct nonzero coordinates \cite[Corollary~2.1(1)]{timofte}. In
four variables this leaves the whole tetrahedron. For the specific octic
$G_4$, varying the product of three distinct positive coordinates while
fixing their sum and their pairwise sum gives strictly negative
curvature, as the Hessian below shows.

\subsection{The octic in elementary symmetric functions}
Enumerate $\Fset=\{\xi_1,\ldots,\xi_{35}\}$ and put $L_j=L_{4,\xi_j}$.
Let $v(x)=(x^\nu)_{|\nu|=4}$ be a row vector, $V$ the matrix with rows
$v(\xi_i)$, and $W=V^{-1}$. The entries of $v(x)W$ are the homogeneous
cardinal polynomials. Define the homogeneous deficit
\begin{equation}\label{eq:deficit}
  H(x):=G_4(x)=s^8-v(x)WW^{\T}v(x)^{\T}.
\end{equation}
On $\Delta_3$ this is $1-\sum_jL_j^2$. Permuting the coordinates permutes the
nodes and hence the cardinal polynomials, so $H$ is symmetric. We will
also evaluate \eqref{eq:deficit} at unnormalized coordinates. Its sign is
unchanged by positive rescaling.


Let $e_k$ denote the $k$th elementary symmetric polynomial in $x$,
and write $e=e_2$, $p=e_3$, $q=e_4$. On $s=1$ we have $H(x)=F(e,p,q)$ with
\begin{equation}\label{eq:F}
 F(e,p,q)=U(e)+P(e)p+R(e)q+A(e)p^2+Bpq+Cq^2.
\end{equation}
The coefficients in \eqref{eq:F} are
\begin{align*}
 U(e)&=20e-360e^2+2100e^3-3920e^4,\\
 P(e)&=\frac{4623-57r}{11}+\frac{-125333+7219r}{33}e
               +\left(\frac{29764}{3}-924r\right)e^2,\\
 R(e)&=\frac{16458-5230r}{11}+\frac{4491922-33070r}{363}e
               +\frac{-1742572+221332r}{33}e^2,\\
 A(e)&=\frac{-136228+14134r}{11}+\frac{77200-3296r}{3}e,\\
 B&=\frac{2845034-2736862r}{363},\qquad
 C=\frac{-17739344+1312048r}{363}.
\end{align*}
Equation~\eqref{eq:F} follows by substituting the cardinals in
Appendix~\ref{sec:cardinals} into \eqref{eq:deficit}. Equivalently,
homogenize each monomial $e^ip^jq^k$ in \eqref{eq:F} by multiplying it
by $s^{8-2i-3j-4k}$, expand the elementary functions, and compare with
the entries of Table~\ref{tab:input} supported on at most four
coordinates. The paired edge formulas explain the cancellation of
$\sqrt{21}$.

\begin{lemma}\label{lem:hessian}
For every $0\le e\le 3/8$, the Hessian of $F$ in $(p,q)$ at fixed $e$,
\[
 M(e)=\begin{pmatrix}2A(e)&B\\ B&2C\end{pmatrix},
\]
is negative definite. This interval contains $e_2(\Delta_3)$.
\end{lemma}
\begin{proof}
Since $2e=1-\sum_i x_i^2$ on $\Delta_3$ and $\sum_i x_i^2\geq1/4$, we have
$0\le e\le 3/8$. The slope of $A$ is positive, $C<0$, and
$A(3/8)=(-30078+9602r)/11<0$. The determinant $4A(e)C-B^2$ therefore
decreases with $e$. At the right endpoint it equals
\[
 \det M(3/8)=
 \frac{8(4149878507601-1515033492901r)}{131769}>0 .
\]
The signs follow from $2<r<9/4$. The two leading principal minors now
give negative definiteness.
\end{proof}

\subsection{A minimum has at most two distinct positive coordinates}

\begin{lemma}\label{lem:min}
Every local minimum of $H$ in the interior of $\Delta_3$ has at most two
distinct coordinate values. If the restriction of $H$ to a triangular
face has a local minimum in its relative interior, at least two of the
three positive coordinates coincide.
\end{lemma}
\begin{proof}
Suppose $x\in \Delta_3$ has three pairwise distinct positive coordinates
$a,b,c$. Write $d$ for the fourth coordinate and, after a permutation,
take $x=(a,b,c,d)$. Hold $d$ fixed. For both signs of $t$ sufficiently
close to zero, the polynomial
\[
 K(z)-t,\qquad K(z)=(z-a)(z-b)(z-c),
\]
has three distinct positive roots $a(t),b(t),c(t)$ tending to $a,b,c$,
respectively, by the \hyperref[par:cubic-variation]{sign-change argument
for three simple roots}. The sum of these roots stays $a+b+c$, the negative
of the cubic's quadratic coefficient. Thus $x(t)=(a(t),b(t),c(t),d)$
belongs to $\Delta_3$ and $x(t)\to x$.

Let $e,p,q$ be the invariants at $x$. The polynomial with roots the
coordinates of $x(t)$ is
\begin{equation}\label{eq:cubicpath}
 (z-d)(K(z)-t)
 =z^4-z^3+ez^2-(p+t)z+(q+dt).
\end{equation}
In particular, $e_2$ stays fixed and $H(x(t))=F(e,p+t,q+dt)$.
This is a quadratic in $t$. Using its Hessian from
Lemma~\ref{lem:hessian}, for every sufficiently small $t>0$ we obtain
\begin{equation}\label{eq:second}
 H(x(t))+H(x(-t))-2H(x)
 = t^2(1,d)M(e)(1,d)^{\T}<0.
\end{equation}
At least one of the two nearby points has smaller deficit than $x$.

Every interior point with three or four distinct values contains such a
triple. In particular, at $(u,u,v,w)$ choose $a=u$, $b=v$, $c=w$ and
$d=u$. Although the full quartic has a repeated root, the three roots
of $K$ are simple. For a face with three distinct positive coordinates
take $d=0$. The two nearby points then lie on that same face, so
\eqref{eq:second} also excludes a relative minimum there.
\end{proof}

\subsection{The remaining one-parameter families}

By Lemma~\ref{lem:min}, it suffices to study binary homogeneous
restrictions of $H$. We give their factorizations in nonnegative
variables $a,b$. Each identity can be obtained by substituting into
the homogeneous form of \eqref{eq:F}.

\paragraph{Edges and triangular faces.}
On an edge we have
\begin{equation}\label{eq:edge}
 H(a,b,0,0)=20ab(a-b)^2(a^2-5ab+b^2)^2.
\end{equation}
When $a+b=1$, the quadratic factor is $1-7ab$, whose roots give
$a=(1\pm\eta)/2$. Together with the endpoints and midpoint these are
the five zeros on the edge.

For a face point with a repeated coordinate the factorization is
\begin{equation}\label{eq:face}
 H(a,a,b,0)=ab\left(b-\frac{3+r}{2}a\right)^2 T(a,b).
\end{equation}
The remaining quartic can be written as
\begin{equation}\label{eq:sos}
 T(a,b)=40(b^2-\lambda ab+a^2)^2+
          a^2(A_0b^2+B_0ab+C_0a^2).
\end{equation}
The coefficients in this decomposition are
\[
 \begin{aligned}
 \lambda&=\frac{4397-383r}{880},&
 A_0&=\frac{23382299+5935673r}{29040},\\
 B_0&=\frac{-229907+88485r}{33},&
 C_0&=\frac{316632-138752r}{33}.
 \end{aligned}
\]
Here $A_0,C_0>0$, and
\[
 4A_0C_0-B_0^2
 =\frac{-4238873363448+1896536114072r}{59895}>0 .
\]
For these signs it suffices to use $559/250<r<9/4$, which follows
by squaring. Thus the quadratic form in \eqref{eq:sos} is positive
definite. When $a=0$, the first square is $40b^4$, so $T(a,b)>0$
whenever $(a,b)\ne(0,0)$.

Each triangular face is compact. Its minimum occurs on an edge or,
by Lemma~\ref{lem:min}, at a repeated-coordinate point. Equations
\eqref{eq:edge} and \eqref{eq:face} prove nonnegativity on the whole
boundary. Any zero in a face interior is therefore a relative minimum.
Up to permutation it has the form $(a,a,b,0)$ with $a,b>0$. Solving
$b=(3+r)a/2$ and $2a+b=1$ gives $a=\beta$, $b=\alpha$.

\paragraph{Two pairs and a triple.}
An interior point with at most two coordinate values has multiplicities
$2+2$ or $3+1$, allowing $a=b$ in either expression. For two pairs,
direct substitution gives
\begin{equation}\label{eq:pairs}
 H(a,a,b,b)=ab(a-b)^2
 \bigl[c(a-b)^4+d\,ab(a-b)^2+f\,a^2b^2\bigr],
\end{equation}
where
\[
 c=\frac{144384-17408r}{33},\quad
 d=\frac{-433948+285060r}{363},\quad
 f=\frac{2147592+1441384r}{363}.
\]
All three coefficients are positive by the same bounds on $r$. Thus
\eqref{eq:pairs} is strictly positive for distinct positive $a,b$.

\Needspace{8\baselineskip}
For a triple the analogous factorization is
\begin{equation}\label{eq:triple}
 H(a,a,a,b)=a(a-b)^2Q(a,b).
\end{equation}
Here $Q$ is the binary quintic given by
\begin{equation}\label{eq:quintic}
\begin{aligned}
 121Q(a,b)={}&(264682-26576r)a^5+(3606167+2531r)a^4b\\
 &+(-3940095+238179r)a^3b^2+(685905+148093r)a^2b^3\\
 &-(87021+1881r)ab^4+7260b^5 .
\end{aligned}
\end{equation}
The mixed signs in \eqref{eq:quintic} can be resolved by bounding
$g(t)=Q(t,1-t)$ on subintervals of $0\le t=a/(a+b)\le1$. We use the
Bernstein coefficient bounds described by Garloff \cite[\S\S1.1--1.2]{garloff}.
For each interval $[\ell,h]$ in Table~\ref{tab:quintic}, expand
$g(\ell+(h-\ell)z)=\sum_{k=0}^5 d_kz^k$. Its Bernstein coefficients
on that interval are
\begin{equation}\label{eq:bern}
 b_j=\sum_{k=0}^j d_k\frac{\binom{j}{k}}{\binom{5}{k}},
 \qquad
 g(\ell+(h-\ell)z)=\sum_{j=0}^5 b_j\binom{5}{j}z^j(1-z)^{5-j}.
\end{equation}
For a coefficient $b_j=u+vr$ with $u,v\in\Q$, a rational lower bound
is $u+559v/250$ if $v\geq0$ and $u+9v/4$ if $v<0$.
The table gives the integer parts of these bounds. For example, on
the first interval $b_1=(1989-171r)/220\ge6417/880>7$.

\begin{table}[htbp]
 \centering
 \caption{Integer lower bounds for the Bernstein coefficients of
 $Q(t,1-t)$.}\label{tab:quintic}
 \smallskip
 \begin{tabular}{@{}lrrrrrr@{}}
 \toprule
 Interval & $b_0$ & $b_1$ & $b_2$ & $b_3$ & $b_4$ & $b_5$\\
 \midrule
 $[0,1/4]$   & 60 & 7 & 29 & 35 & 23 & 17\\
 $[1/4,1/2]$ & 17 & 11 & 11 & 40 & 148 & 346\\
 $[1/2,1]$   & 346 & 743 & 1499 & 2704 & 3831 & 1693\\
 \bottomrule
 \end{tabular}
\end{table}

The Bernstein basis functions are nonnegative and sum to one on
$0\leq z\leq1$. It follows that $g>0$ throughout $[0,1]$.
By homogeneity $Q(a,b)=(a+b)^5g(a/(a+b))>0$ for all nonnegative
$a,b$ with $a+b>0$. Thus \eqref{eq:triple} is strictly positive at
every interior point of this family with $a\ne b$.

\begin{proof}[Proof of Theorem~\ref{thm:tetra}]
The continuous function $H$ attains its minimum on $\Delta_3$. The boundary
calculation excludes a negative boundary minimum. An interior minimum
has at most two distinct coordinate values by Lemma~\ref{lem:min}, and
\eqref{eq:pairs} and \eqref{eq:triple} exclude a negative value there.
This proves \eqref{eq:main}.

Every zero of $H$ is now a global minimum. In the interior the two
factorizations force all four coordinates to coincide, giving the
centroid. On a face interior the zeros are the permutations of
$(\alpha,\beta,\beta,0)$, and \eqref{eq:edge} gives precisely the five
nodes on each edge. These exhaust the orbits in \eqref{eq:nodes}. Conversely, the
cardinal conditions give equality at every listed point.
\end{proof}

\begin{corollary}\label{cor:design}
The uniform probability measure on $\F_4$ is the unique quartic
$D$-optimal design on $\Delta_3$.
\end{corollary}
\begin{proof}
Work in the cardinal basis $L=(L_1,\ldots,L_{35})^{\T}$ and put
$G_\mu=\int LL^{\T}\,d\mu$. By Theorem~\ref{thm:tetra},
$\tr G_\mu=\int\sum_jL_j^2\,d\mu\leq1$. The arithmetic--geometric
mean inequality for the eigenvalues of $G_\mu$ gives
\[
  \det G_\mu\leq(\tr G_\mu/35)^{35}\leq 35^{-35}.
\]
The uniform measure on $\Fset$ attains the bound since its matrix is
$I/35$. Equality forces $\tr G_\mu=1$ and every eigenvalue to equal
$1/35$. The trace condition and the equality set in
Theorem~\ref{thm:tetra} force $\mu$ to be supported on $\Fset$.
In that case $G_\mu$ is diagonal with the node masses on its diagonal,
so all masses are $1/35$.
\end{proof}

\section{The fourth-power bound and local squared norms}\label{sec:four}
The corrected squared-norm estimate also proves the fourth-power theorem.
The cross terms between distinct tetrahedral face-center cardinals
absorb the correction in Lemma~\ref{lem:correction}.

\begin{proof}[Proof of Theorem~\ref{thm:four}]
For $x\in\Delta_{m-1}$ set
\[
 A_m(x)=\sum_{J\ne K}q_J(x)^2q_K(x)^2,
\]
with ordered pairs as in the definition of $E_m$. By
Lemma~\ref{lem:centers},
\[
 0\le A_m\le\left(\sum_J q_J^2\right)^2\le1,
 \qquad
 q_Jq_K\le\frac{q_J^2+q_K^2}{2}\le\frac12\quad(J\ne K).
\]
Consequently $E_m\le A_m/2$. Expanding the square of $S_m$ and retaining
only the cross terms between center cardinals gives
\[
 \sum_{z\in\F_m}L_{m,z}^4\le S_m^2-A_m.
\]
Lemma~\ref{lem:correction} now yields
\begin{align*}
 \sum_{z\in\F_m}L_{m,z}^4
 &\le \left(1+\frac{E_m}{16}\right)^2-A_m\\
 &\le \left(1+\frac{A_m}{32}\right)^2-A_m
 =1-A_m\left(\frac{15}{16}-\frac{A_m}{1024}\right)\le1.
\end{align*}
If $x\notin\F_m$, Lemma~\ref{lem:correction} is strict. Since $S_m\ge0$
and $1+E_m/16>0$, the first inequality in this display is then strict.
The cardinal conditions give equality at each node. The argument also
covers $m<4$, when both center sums are empty.
\end{proof}

\subsection{The squared norm near a node}
Although $S_5$ exceeds one at the centroid of $\Delta_4$, its behavior
near each interpolation node persists in every dimension.

\begin{corollary}\label{cor:local}
For every $m\ge2$ and $z\in\F_m$ there is a relative neighborhood of $z$
in $\Delta_{m-1}$ in which $S_m(x)<1$ for $x\ne z$.
\end{corollary}
\begin{proof}
For $m\le4$ the assertion follows from Theorem~\ref{thm:tetra} and face
restriction. Suppose $m\ge5$, let $I=\supp z$, and take $j\notin I$.
The cardinal conditions and $s(z)=1$ give
\begin{equation}\label{eq:normal}
 \partial_jG_m(z)=8-2\partial_jL_{m,z}(z).
\end{equation}
For nodes with $|I|\le3$, restrict to the coordinate subspace indexed by
$I\cup\{j\}$. The derivative is therefore determined by a tetrahedral
cardinal polynomial. For $|I|=4$ the cardinal is $q_I$, which is
independent of $x_j$. The values in Table~\ref{tab:normal} follow.

\begin{table}[!ht]
\centering
\caption{Derivatives in a coordinate outside the support of a node.
Here $r=\sqrt5$ and $G_m=s^8-S_m$.}\label{tab:normal}\smallskip
\begin{tabular}{@{}lcc@{}}
\toprule
Node type & $\partial_jL_{m,z}(z)$ & $\partial_jG_m(z)$\\
\midrule
Vertex & $-6$ & $20$\\
Edge midpoint & $(-150+34r)/33$ & $(564-68r)/33$\\
Other edge node & $(-57-17r)/33$ & $(378+34r)/33$\\
Triangular face node & $-1$ & $10$\\
Tetrahedral face center & $0$ & $8$\\
\bottomrule
\end{tabular}
\end{table}

For completeness, the values can also be obtained directly from
Table~\ref{tab:input}. If the positive coordinates of $z$ are
$z_1,\ldots,z_k$, then
\[
 \partial_jG_m(z)=\frac1{363}
 \sum_{\substack{\nu\in\mathbb Z_{\ge0}^k\\|\nu|=7}}
       g_{\partn(\nu,1)}z_1^{\nu_1}\cdots z_k^{\nu_k}.
\]
The five substitutions are the representatives in \eqref{eq:nodes}.
All the derivatives in the last column of Table~\ref{tab:normal} are
positive, using $\sqrt5<9/4$ for the midpoint. By continuity there is a
ball about $z$ in which $\partial_jG_m>0$ for every $j\notin I$.

For $x$ in that ball with $x\ge0$, let $x_I$ be obtained by setting the
coordinates outside $I$ to zero. The segment $[x_I,x]$ stays in the
ball. Integration gives $G_m(x)>G_m(x_I)\ge0$ whenever $x\ne x_I$.
The last inequality follows by homogeneity from Theorem~\ref{thm:tetra}
on the face with support $I$. On that face and $s=1$ there are finitely
many nodal equality points. Shrinking the ball to exclude the other
nodes proves the assertion.
\end{proof}

\section{Extension to higher polynomial degree}\label{sec:higher-degree}
Fix a simplex and increase the degree from four to $n$. The determinant
estimate of Lemma~\ref{lem:matrix} still applies. For a unisolvent set
$X$ with cardinals $\ell_z$, put $S_X=\sum_{z\in X}\ell_z^2$.
A sufficient condition for $X$ to be Fekete is that a fixed subset
$q_1,\ldots,q_p$ of the cardinals satisfy, throughout the simplex,
\begin{equation}\label{eq:general-correction}
 q_j\ge0,\qquad \sum_jq_j^2\le1,\qquad
 S_X\le1+\tfrac1{16}\sum_{j\ne k}q_j^3q_k^3 .
\end{equation}
Indeed, applying the lemma to the matrix $(\ell_z(y_i))_{i,z}$ bounds
the determinant ratio for every competing configuration. The issue in
higher degree is how large a squared row norm this estimate can absorb.
If $a_j=q_j^2$ and $h=\sum_j a_j\le1$, then
\[
 \sum_{j\ne k}(a_ja_k)^{3/2}
 \le \sum_j a_j^2(h-a_j)\le h^3/4 .
\]
The first inequality uses $\sqrt{a_ja_k}\le(a_j+a_k)/2$ and the second
uses $a_j(h-a_j)\le h^2/4$. Consequently \eqref{eq:general-correction}
requires $S_X\le65/64$.

\paragraph{An obstruction on the quintic triangle.}
At the centroid $c$ of the triangle Bos reports $S_X(c)\approx1.2246$
for his 21 nodes \cite[Numerical Example~5.2, p.~289]{bos2023}.
One quintic already gives an exact separation. On $\sum x_i=1$ take
$f=27x_1x_2x_3(15e_2-4)$. It vanishes at the 15 boundary nodes and
$f(c)=1$. The six interior nodes form two three-point orbits of
$(t,t,1-2t)$, with $0.1480<u<0.1481$ and $0.4208<v<0.4209$ as
their parameters (Appendix~\ref{sec:quintic-triangle}). On such an orbit
$f=27t^2(1-2t)[1-5(1-3t)^2]$. Bounding the two factors separately
on the displayed intervals gives $|f|<1/4$ on the first orbit and
$|f|<1/2$ on the second. Interpolation at $c$ and Cauchy--Schwarz give
\[
 1\le S_X(c)\sum_{z\in X}f(z)^2
   < S_X(c)\,3\bigl((1/4)^2+(1/2)^2\bigr).
\]
Thus $S_X(c)>16/15>65/64$. Every choice of columns in
\eqref{eq:general-correction} encounters the same obstruction.

\paragraph{Two steps that depend on the degree.}
For a theorem uniform over simplices, the quartic face-center bound
survives in degree five. For $m\ge n$
and $n\ge2$, consider the natural products $q_J=n^n\prod_{i\in J}x_i$,
$|J|=n$. For $n=5$ they are the cardinals at the centers of Bos's
four-dimensional faces \cite[\S5, p.~291]{bos2023}. In general,
\[
 \max_{x\in\Delta_{m-1}}\sum_Jq_J(x)^2
   =\max_{n\le k\le m}\binom{k}{n}(n/k)^{2n}.
\]
The related maximization over all support sizes is the uniform-profile
theorem \cite[Theorem~14]{orlitsky}. Here the sum is
$n^{2n}e_n(x_1^2,\ldots,x_m^2)$. At a positive maximum, choose two
positive coordinates and fix their sum and all other coordinates $w$.
As in Lemma~\ref{lem:centers}, the expression is strictly convex in
their product: the coefficient of its square is $n^{2n}e_{n-2}(w^2)>0$.
Hence all $k$ positive coordinates are equal. For $n=5$ the ratio of
successive terms on the right is $1/g(1/k)$ with
$g(t)=(1-4t)(1+t)^9$. Its derivative has the sign of $1-8t$, and
$g(1/5)>g(0)=1$. The maximum is therefore one, attained at $k=5$.
For $n\ge6$, already $k=n+1$ gives $(n+1)(n/(n+1))^{2n}>1$.
At $n=6$ this is $6^{12}/7^{11}>1$, and for $n\ge7$ the value
exceeds $(n+1)e^{-2}>1$. Thus these products first violate the
uniform bound in degree six, on a six-dimensional simplex.

The reduction to two coordinate values also uses the degree.
Along the variation in Section~\ref{sec:certificate}, the sum and
pairwise sum of three coordinates are fixed. A symmetric polynomial
of degree $2n$ has degree at most $\lfloor 2n/3\rfloor$ in their
product. For the quartic squared norm the dependence is at most
quadratic. At degree five the term $e_1e_3^3$ is already possible,
so the curvature argument would require control of a cubic term.

\paragraph{A fixed triangle in degrees five and six.}
There is another way to approach higher degrees on the same domain.
In degree $n\ge3$, prescribe the three vertices and $n-1$ further
distinct points in the relative interior of each edge. The $3n$ boundary nodes leave
$r_n=\binom{n-1}{2}$ points $Z$ in the interior. A polynomial of degree
at most $n$ vanishing at the boundary nodes has $n+1$ zeros on each edge
and hence vanishes on all three edges,
so it is $w p$ with $w=x_1x_2x_3$ and $\deg p\le n-3$.
Bos studied this separation of boundary and interior nodes in
\cite{bos1991} and applied it to the quintic triangle in
\cite[\S5, pp.~287--288]{bos2023}. With the boundary fixed, the
quantity to maximize over $Z$ is
\begin{equation}\label{eq:weighted-triangle}
 \left(\prod_{z\in Z} w(z)\right)
 \left|\det\bigl(p_j(z)\bigr)_{z\in Z,\,1\le j\le r_n}\right|,
\end{equation}
where the $p_j$ form any basis of the polynomials of degree at most
$n-3$. To see this, complete the functions $wp_j$ to a degree-$n$
basis. Put the boundary rows first and the $wp_j$ columns last.
The upper right block vanishes, leaving the product of the two
diagonal determinants.

For quintics, \eqref{eq:weighted-triangle} is a weighted quadratic
determinant at six points. Bos computed it for two symmetric triples.
Appendix~\ref{sec:quintic-triangle} proves that his choice uniquely
maximizes it among permutation-invariant sets of six interior points.
For sextics the corresponding problem is a weighted cubic determinant
at ten points. An unrestricted quintic proof along this route would
have to maximize \eqref{eq:weighted-triangle} over arbitrary positions
of the six interior points and then allow the number and locations of
boundary nodes to vary. The unrestricted quintic case remains open here.

\appendix
\section{Explicit homogeneous cardinal polynomials}\label{sec:cardinals}
Bos gives barycentric formulas for the quartic cardinals
\cite{bos1983} and lists them by node type in \cite[\S4]{bos2023}.
We use homogeneous expressions in the coordinates on and off the
support of a node. The nodal evaluations below verify their coefficients.
For a fixed support $I$, take all sums below over its complement and put
\[
 \sigma_j=\sum_{i\notin I}x_i^j,\qquad
 \varepsilon_j=e_j((x_i)_{i\notin I}),\qquad
 M_{21}=\sigma_1\sigma_2-\sigma_3.
\]
Empty sums and positive-degree elementary functions of an empty set are
zero. The notation $a,b,c$ denotes coordinates in $I$.

At a vertex, with $I=\{i\}$ and $a=x_i$, the cardinal is
\begin{align}
132L_i={}&132a^4-792a^3\sigma_1+792a^2\sigma_2
 +(1566+342r)a^2\varepsilon_2\notag\\
&-132a\sigma_3-(480+120r)aM_{21}
 -(1554+306r)a\varepsilon_3.\label{eq:card-vertex}
\end{align}
For the midpoint on an edge with support $I=\{i,j\}$, write
$a=x_i,b=x_j$. Then
\begin{align}
132L_{ij}^{\rm mid}=ab\bigl[&-704(a^2+b^2)+3520ab
 +(-2400+544r)(a+b)\sigma_1\notag\\
&+(2528-864r)\sigma_2+(2432-448r)\varepsilon_2\bigr].
\label{eq:card-mid}
\end{align}
For the other two edge nodes, put $h=\sqrt{21}$ and
\begin{align*}
 U={}&1078(a-b)^2-(441+245r)(a+b)\sigma_1
       +(98+294r)\sigma_2+(1568+392r)\varepsilon_2,\\
 V={}&154(a^2-b^2)-(385+77r)(a-b)\sigma_1.
\end{align*}
Their cardinals are
\begin{equation}\label{eq:card-edge}
 132L_{ij}^{\pm}=ab(U\pm hV).
\end{equation}
The plus sign belongs to the node $((7+h)/14,(7-h)/14)$ in the
ordered coordinates $(i,j)$, and the minus sign belongs to its reversal.
In particular their squared sum is $2a^2b^2(U^2+21V^2)/132^2$.

Two details matter when comparing the edge formulas with those printed
in \cite[\S4, p.~284]{bos2023}. On $s=1$, the printed midpoint expression has
value $(14r-43)/33$ at the midpoint. Adding $76-14r$ inside its braces
gives \eqref{eq:card-mid}. In the oriented-edge expression, replacing
the repeated $\lambda_i$ in the second linear term by $\lambda_j$
gives $L_{ij}^{-}$ in \eqref{eq:card-edge} for $a=\lambda_i$, $b=\lambda_j$.

For a triangular node with value $\alpha$ at $i$ and values $\beta$ at
$j,k$, let $I=\{i,j,k\}$ and $(a,b,c)=(x_i,x_j,x_k)$. Then
\begin{equation}\label{eq:card-triangle}
 132L_{i\mid jk}=abc\bigl[(330+4818r)a
                    +(2244-2376r)(b+c)-(4818+66r)\sigma_1\bigr].
\end{equation}
Finally, the cardinal at a tetrahedral face center is
\begin{equation}\label{eq:card-center}
 L_I=256\prod_{i\in I}x_i\qquad(|I|=4).
\end{equation}

Each displayed polynomial contains the product of the coordinates in
its support. It therefore vanishes at nodes whose supports do not
contain that support. On the containing supports, substitution of the
five patterns in \eqref{eq:nodes} gives one at its own node and zero at
all the others. For example, on an edge the vertex expression is
$a(a-b)((a+b)^2-7ab)$, and \eqref{eq:card-edge} has value one at its
assigned non-midpoint node because $ab=1/7$ and
$a-b=\pm\sqrt{21}/7$. At every tetrahedral center containing the
support, the bracket in each of \eqref{eq:card-mid}--\eqref{eq:card-triangle}
vanishes. The remaining triangular substitutions use
$\alpha+2\beta=1$, $\alpha=(4+r)/11$, $\beta=(7-r)/22$, and $r^2=5$.
Unisolvence identifies these polynomials with the cardinal basis.

Squaring \eqref{eq:card-vertex}--\eqref{eq:card-center} and summing over
supports gives Table~\ref{tab:input} and \eqref{eq:stable}.
Only support sizes at most eight can occur. The cancellation in the
squared pair \eqref{eq:card-edge} explains why the resulting octic has
coefficients in $\Q(\sqrt5)$.

\section{The four binary forms in five coordinates}\label{sec:five-binary}

This appendix supplies the finite algebra in
Lemma~\ref{lem:five-analytic}. All coefficients and bounds needed for the
argument are given here. For $1\le p\le4$, put $q=5-p$ and define
\[
 H_p(z)=363G_5(\underbrace{1,\ldots,1}_{p},
                       \underbrace{z,\ldots,z}_{q})
       =\sum_{j=0}^8h_{p,j}z^j .
\]
Table~\ref{tab:five-octics} gives $H_1,H_2$ directly from
Table~\ref{tab:input}. Symmetry and homogeneity give
$H_3(z)=z^8H_2(1/z)$ and $H_4(z)=z^8H_1(1/z)$; these are polynomial
identities, including at $z=0$ after cancellation.

\begin{table}[!ht]
\centering\small
\caption{Coefficients $h_{p,j}=a_{p,j}+b_{p,j}\sqrt5$ of the two octics.}
\label{tab:five-octics}\smallskip
\begin{tabular}{@{}rrrrr@{}}
\toprule
$j$&$a_{1,j}$&$b_{1,j}$&$a_{2,j}$&$b_{2,j}$\\
\midrule
0&0&0&0&0\\
1&29040&0&2382336&$-287232$\\
2&$-435006$&$-11286$&$-15413304$&2738556\\
3&5077380&982212&62019474&$-5862582$\\
4&$-35500485$&$-1384047$&$-101452455$&6537219\\
5&104017092&348060&89414028&$-6274488$\\
6&$-118452912$&$-1138584$&$-53890698$&1690512\\
7&47579136&0&18460818&334098\\
8&0&0&794046&$-79728$\\
\bottomrule
\end{tabular}
\end{table}

For $A,T\ge0$ write the degree-24 binary form as
\begin{equation}\label{eq:five-binary}
 \Phi_5(A^p,(A+T)^q)
 =\sum_{j=0}^{24}\binom{24}{j}\beta_{p,j}(5A)^{24-j}(qT)^j,
\end{equation}
where powers inside the argument denote repeated coordinates.
Here is an explicit formula for every coefficient. Set
\[
 g_{p,\ell}=\sum_{h=\ell}^8\binom h\ell h_{p,h}.
\]
The coefficient $c_{p,j}=\binom{24}{j}5^{24-j}q^j\beta_{p,j}$ is
\begin{equation}\label{eq:five-binary-coeff}
 \begin{split}
 c_{p,j}={}&\frac1{363}\sum_{\substack{0\le\ell\le8\\0\le j-\ell\le16}}
 \binom{16}{j-\ell}5^{16-j+\ell}q^{j-\ell}g_{p,\ell}\\
 &+2^{45}\sum_{\ell=0}^3
 \binom p\ell\binom q{3-\ell}\binom{3q+9-3\ell}{j}.
 \end{split}
\end{equation}
Binomial coefficients outside their usual nonnegative range are zero.
The first sum expands $(5+qT)^{16}H_p(1+T)/363$.
The second follows from \eqref{eq:E-five}, since
\[
 E_5(1^p,(1+T)^q)/16
 =2^{45}\sum_{\ell=0}^3\binom p\ell\binom q{3-\ell}
                       (1+T)^{3q+9-3\ell}.
\]

\begin{table}[!ht]
\centering
\caption{Integer lower bounds for $1000\beta_{p,j}$.
Each entry applies to every index in the indicated range.}
\label{tab:five-bounds}\smallskip
\begin{tabular}{@{}lrrrr@{}}
\toprule
$j$&$p=1$&$p=2$&$p=3$&$p=4$\\
\midrule
$0,1$&3&3&3&3\\
$2$--$5$&4&6&9&21\\
$6$--$9$&22&50&91&103\\
$10$--$13$&65&153&213&193\\
$14$--$17$&125&319&309&453\\
$18$--$21$&144&513&416&575\\
$22,23$&66&425&312&666\\
$24$&0&258&0&0\\
\bottomrule
\end{tabular}
\end{table}

Substitution in \eqref{eq:five-binary-coeff} gives the rational lower
bounds in Table~\ref{tab:five-bounds}. They use only
$2236/1000<\sqrt5<2237/1000$: in each coefficient $a+b\sqrt5$, replace
the radical by $2236/1000$ for $b\ge0$ and by $2237/1000$ for $b<0$.
Thus every coefficient in \eqref{eq:five-binary} is nonnegative, and
$\beta_{p,0}\ge3/1000$ for every $p$. In particular
\[
 \Phi_5(A^p,(A+T)^q)\ge\frac3{1000}(5A)^{24}>0
 \qquad(A>0,\ T\ge0),
\]
as used in Lemma~\ref{lem:five-analytic}.

\section{Uniform algebra for the two-value families}\label{sec:uniform-algebra}
This appendix proves the coefficient bounds used in
Lemmas~\ref{lem:one-large} and \ref{lem:two-large}. All input polynomials
are specified by \eqref{eq:F}, \eqref{eq:stable}, and \eqref{eq:Rformal}.
The comparisons are finite rational inequalities after using
\begin{equation}\label{eq:radical-lower}
 a+br\ge\ell(a+br):=
 \begin{cases}
 a+2236b/1000,&b\ge0,\\
 a+2237b/1000,&b<0.
 \end{cases}
\end{equation}
Here $2236/1000<r<2237/1000$, as follows by squaring.
For coefficient comparisons we apply $\ell$ separately to each
coefficient after collecting it in $\Q(r)$. Equivalently, $\ell(a+br)$
is the minimum of the two rational endpoint values. We do not apply
$\ell$ separately to factors or treat it as a ring homomorphism.

\subsection{A direct formula for every coefficient}
We give a second way to obtain the binary coefficients, directly from
Table~\ref{tab:input}. It also specifies the finite arithmetic in this
appendix without a multivariate elementary-function expansion.
For a partition $\lambda=(\lambda_1,\ldots,\lambda_l)$, put
$d_\lambda=\prod_{j\ge1}m_j(\lambda)!$, where $m_j(\lambda)$ is the
multiplicity of the part $j$, and write
$(z)_{\underline h}=\prod_{j=0}^{h-1}(z-j)$. Define the nonnegative integers
$\kappa_{\lambda;i,h}$ and the polynomials $M_{\lambda,i}$ by
\begin{equation}\label{eq:kappa}
 \begin{aligned}
 \kappa_{\lambda;i,h}
   &=[t^iz^h]\prod_{a=1}^l\bigl(1+z((1+t)^{\lambda_a}-1)\bigr),\\
 M_{\lambda,i}(n,q)&=\frac1{d_\lambda}
   \sum_{h=0}^l\kappa_{\lambda;i,h}
       (q)_{\underline h}(n-h)_{\underline{l-h}}.
 \end{aligned}
\end{equation}
Indices outside their nonnegative ranges give zero. If $m_\lambda$
denotes the sum of distinct monomials with exponent partition
$\lambda$, then
\begin{equation}\label{eq:monomial-counts}
 m_\lambda(A^{n-q},(A+T)^q)
   =\sum_{i=0}^{|\lambda|}M_{\lambda,i}(n,q)A^{|\lambda|-i}T^i.
\end{equation}
To prove this identity, first label the $l$ parts. Choose the $h$
parts that contribute a positive power of $T$; the binomial
coefficients of those powers are counted by $\kappa$.
Their coordinate indices can be assigned in $(q)_{\underline h}$
ways. The other labeled parts can use any of the remaining $n-h$
coordinates, in $(n-h)_{\underline{l-h}}$ ways. Dividing by
$d_\lambda$ removes the labelings of equal parts. This proves
\eqref{eq:monomial-counts} for integer $n\ge q\ge0$ and specifies
its polynomial extension.

With $g_\lambda$ from Table~\ref{tab:input}, set
\begin{equation}\label{eq:HJ}
 H_i=\sum_{\lambda\vdash8}g_\lambda M_{\lambda,i},
 \qquad
 J_i=\sum_{\mu\vdash7}g_{\mu\cup(1)}M_{\mu,i}.
\end{equation}
The appended partition is sorted, and $H_i=0$ outside $0\le i\le8$,
$J_i=0$ outside $0\le i\le7$.
Write the normalized binary coefficients as
\begin{equation}\label{eq:normalized-CB}
 C=\sum_{i=0}^8\binom8i\widehat C_i(n,q)A^{8-i}T^i,
 \qquad
 B=\sum_{i=0}^8\binom8i\widehat B_i(n,q)A^{8-i}T^i.
\end{equation}
Then the explicit formulas are
\begin{equation}\label{eq:CB-direct}
 \boxed{\begin{aligned}
 \widehat C_i&=\frac{(24-i)nH_i+(9-i)qH_{i-1}}{\binom8i},\\
 \widehat B_i&=\frac{16H_i+nJ_i+qJ_{i-1}}{\binom8i}.
 \end{aligned}}
\end{equation}
Indeed, $H(A,T)=\sum_iH_iA^{8-i}T^i$ is $363G_n$ on the two-value
family. If one adds a coordinate $y$, the coefficient of $y$ in
$363G_{n+1}$ is $J(A,T)=\sum_iJ_iA^{7-i}T^i$: only partitions with
a part equal to one contribute. A common increase in the $n$
displayed coordinates differentiates $H$ with respect to $A$ while
holding $T$ fixed. Each of the $m-n$ zero coordinates contributes $J$.
Consequently, with $s=nA+qT$,
\[
 363R_m=s\partial_AH+16nH+(m-n)(sJ+16H).
\]
Comparing with \eqref{eq:ambient} and collecting $T^i$ proves
\eqref{eq:CB-direct} for integer $n\ge q\ge0$. It is also an identity
in the formal parameters: for each nonnegative integer $q$ it holds
for infinitely many $n$, and then each coefficient in $n$ agrees
for infinitely many $q$.

For clarity, the integers in \eqref{eq:kappa} can be generated by the
short recurrence
\[
 \kappa_{\varnothing;0,0}=1,\qquad
 \kappa_{\lambda\cup(a);i,h}
  =\kappa_{\lambda;i,h}
    +\sum_{b=1}^a\binom ab\kappa_{\lambda;i-b,h-1}.
\]
All other initial entries are zero. Thus every comparison below uses
only this integer recurrence, falling factorials of order at most
eight, and the two rational endpoints in \eqref{eq:radical-lower}.

\subsection{One larger coordinate}
Define
\[
 \sum_{i=0}^8c_i(p)T^i=\frac{363}{p}\mathcal R_{p+1}(1^p,1+T),
 \qquad
 \sum_{i=0}^8d_i(p)T^i=363\mathcal B(1^p,1+T).
\]
These polynomials are obtained by the substitution
\begin{equation}\label{eq:one-sub}
 e_j=\binom pj+(1+T)\binom p{j-1},\qquad 1\le j\le8.
\end{equation}
Equivalently, \eqref{eq:CB-direct} gives
\[
 c_i(p)=\frac{\binom8i}{p}\widehat C_i(p+1,1),\qquad
 d_i(p)=\binom8i\widehat B_i(p+1,1).
\]
The numerator defining $c_i$ vanishes at $p=0$, since
$\mathcal R_1(1+T)=0$, and is therefore divisible by $p$.
Both $c_i,d_i$ have degree at most $8-i$ in $p$.
Table~\ref{tab:one-minorants} specifies targets $\lambda_i(p)$ and
positive integers $\rho_i$.

\begin{table}[!ht]
\centering
\caption{Targets and residual bounds for the one-larger-coordinate family.}
\label{tab:one-minorants}\smallskip
\begin{tabular}{@{}rcr@{}}
\toprule
$i$&$\lambda_i(p)$&$\rho_i$\\
\midrule
0&0&8712\\
1&0&69696\\
2&0&180143\\
3&$500000p^5$&350859\\
4&$-700000p^4$&588215\\
5&$1000000p^3$&95119\\
6&$-100000p^2$&114867\\
7&$100000p$&4791\\
8&0&7260\\
\bottomrule
\end{tabular}
\end{table}

Substitution of \eqref{eq:one-sub} into \eqref{eq:Rformal}, followed by
collection in $T$ and expansion in $w$, gives
\begin{align}
 c_i(5+w)-\lambda_i(5+w)&\succeq\rho_i(1+w)^{8-i},\notag\\
 d_i(3+w)-\lambda_i(3+w)&\succeq\rho_i(1+w)^{8-i}.
 \label{eq:one-minorants}
\end{align}
More explicitly, for each $0\le j\le8-i$, applying $\ell$ to the
coefficient of $w^j$ on either left side gives at least
$\rho_i\binom{8-i}{j}$. This gives all coefficient comparisons in
\eqref{eq:one-minorants}. For example, $c_8=d_8=7260$, while
\[
 c_7(p)=(108999-1881r)p+65241+1881r.
\]
The other rows follow by the same falling-factorial expansion in
\eqref{eq:one-sub}; no coefficient of degree above $8-i$ occurs.

The five nonzero targets give
\[
 \sum_i\lambda_i(p)A^{8-i}T^i
 =10^5pAT^3\bigl[5(pA)^4-7(pA)^3T+10(pA)^2T^2-pAT^3+T^4\bigr].
\]
The bracket is $\Psi(pA,T)$ from \eqref{eq:one-lower}. The residual
constant and leading coefficients in \eqref{eq:one-minorants} are at
least $8712$ and $7260$, respectively. Dropping the other nonnegative
residuals proves \eqref{eq:one-lower}.

For \eqref{eq:one-small}, use the same substitution with the indicated
fixed $p$ and multiply by $((p+1)A+T)^3$. Normalize the coefficient
of $A^{11-j}T^j$ by $\binom{11}{j}(p+1)^{11-j}$.
For $\mathcal R_6$ at $p=1,2,3,4$, the lower bounds from
\eqref{eq:radical-lower}, for every $0\le j\le11$, are respectively
\[
 \frac{16791}{1000},\quad\frac{6359}{1000},\quad
 \frac{4339}{1000},\quad\frac{5690}{1000}.
\]
For $\mathcal B$ at $p=1,2$ they are $7485/1000$ and $6187/1000$.
These bounds exceed four and six, respectively, proving
\eqref{eq:one-small}.

\subsection{Two or more larger coordinates}
Use the normalized coefficients in \eqref{eq:normalized-CB}.
The substitution \eqref{eq:counts} gives degree at most $9-i$ in $n$
and at most $i$ in $q$ for $\widehat C_i$; for $\widehat B_i$ the
corresponding bounds are $8-i$ and $i$. Indeed, a term contributing
$T^h$ from $E_j$ has degree at most $j-h$ in $n$ and $h$ in $q$.
Weighted homogeneity of \eqref{eq:Rformal} then gives these bounds,
including the one additional factor of $n$ in $C$.

The following coefficient inequalities give the full arithmetic needed
for \eqref{eq:monotone-counts} and \eqref{eq:count-boundary}:
\begin{align}
 \ell\!\left([u^av^b]\,\partial_n\widehat C_i(6+u,2+v)\right)
 &\ge\mu_i\binom{8-i}{a}\binom ib,
 &&\substack{0\le i\le8,\ 0\le a\le8-i,\\0\le b\le i},\notag\\
 \ell\!\left([u^av^b]\,\partial_n\widehat B_i(6+u,2+v)\right)
 &\ge\nu_i\binom{7-i}{a}\binom ib,
 &&\substack{0\le i\le7,\ 0\le a\le7-i,\\0\le b\le i},
 \label{eq:derivative-coeff-bounds}\\
 \ell\!\left([w^a]\widehat C_i(6+w,5+w)\right)
 &\ge8712\binom9a,
 &&0\le i\le8,\ 0\le a\le9,\notag\\
 \ell\!\left([w^a]\widehat B_i(6+w,5+w)\right)
 &\ge8712\binom8a,
 &&0\le i\le8,\ 0\le a\le8.\notag
\end{align}
Here $[u^av^b]$ and $[w^a]$ denote ordinary coefficient extraction,
and the constants are given in Table~\ref{tab:count-minorants}.
Also $\partial_n\widehat B_8=0$. Each bound in
\eqref{eq:derivative-coeff-bounds} follows by inserting
\eqref{eq:kappa}--\eqref{eq:HJ} into \eqref{eq:CB-direct},
differentiating where indicated, expanding the shifted falling
factorials, and applying \eqref{eq:radical-lower} to each collected
coefficient. This gives respectively $165$, $120$, $90$, and $81$
rational comparisons in the four lines. The indicated degree bounds
include every coefficient. For the boundary substitutions, the total
degrees are at most nine and eight, respectively.

\begin{table}[!ht]
\centering
\caption{Rational lower bounds in \eqref{eq:derivative-coeff-bounds}.
For $i\le7$, $\mu_i=8712(9-i)$ and $\nu_i=8712(8-i)$.}
\label{tab:count-minorants}\smallskip
\begin{tabular}{@{}rrr@{\hskip2em}rrr@{}}
\toprule
$i$&$\mu_i$&$\nu_i$&$i$&$\mu_i$&$\nu_i$\\
\midrule
0&78408&69696&5&34848&26136\\
1&69696&60984&6&26136&17424\\
2&60984&52272&7&17424&8712\\
3&52272&43560&8&4101&0\\
4&43560&34848&&&\\
\bottomrule
\end{tabular}
\end{table}

For example, $\partial_n\widehat C_8(n,q)=\widehat B_8(q)$ is independent
of $n$. The lower bounds for the coefficient of $v^a$ in
$\widehat B_8(2+v)$, divided by $\binom8a$, for $a=0,\ldots,8$, are
\[
 (1159865,1542714,622638,119284,7399,4101,7668,9449,8712).
\]
This gives the exceptional value $\mu_8=4101$. The remaining rows use
the stated linear formulas for $\mu_i,\nu_i$. More precisely, for
$i\le7$ the normalized minimum in the first line of
\eqref{eq:derivative-coeff-bounds} occurs at $(a,b)=(8-i,i)$ and
equals $\mu_i$. In the second line it occurs at $(a,b)=(7-i,i)$
and equals $\nu_i$. The minima in the boundary lines occur at
$a=9$ and $a=8$, respectively, and equal $8712$ for every $i$.
In the exceptional row the smallest normalized lower bound occurs
at $a=0$, $b=5$. Its exact calculation from \eqref{eq:CB-direct} is
\[
 \frac{[v^5]\widehat B_8(2+v)}{\binom85}
   =\frac{6980759}{672}-\frac{629525}{224}r,
 \qquad
 \ell\left(\frac{6980759}{672}-\frac{629525}{224}r\right)
   =\frac{15748667}{3840}>4101.
\]

Summing the first line of \eqref{eq:derivative-coeff-bounds} over $a,b$
gives
\[
 \partial_n\widehat C_i(6+u,2+v)
 \succeq\mu_i(1+u)^{8-i}(1+v)^i.
\]
Since every $\mu_i\ge4000$, substitution into
\eqref{eq:normalized-CB} proves the first inequality in
\eqref{eq:monotone-counts}. For the second, use
$\nu_i\binom8i=69696\binom7i$ for $i\le7$, and
$\partial_n\widehat B_8=0$. The last two lines of
\eqref{eq:derivative-coeff-bounds}, followed by the binomial theorem in
$A,T$, give \eqref{eq:count-boundary}.

For the small supports put $C_6=C+(6-n)B$. Table~\ref{tab:small-counts}
bounds each coefficient of $C_6$ and $B$ from below after division by the
corresponding coefficient of $(nA+qT)^8$. The entries follow from
\eqref{eq:CB-direct} and \eqref{eq:radical-lower}; all exceed $2000$,
proving \eqref{eq:count-small}.

\begin{table}[!ht]
\centering
\caption{Lower bounds for all normalized coefficients and listed $q$ at each support size.}
\label{tab:small-counts}\smallskip
\begin{tabular}{@{}rcrr@{}}
\toprule
$n$&$q$&$C_6$&$B$\\
\midrule
3&2&4695&2584\\
4&2,3&5808&2904\\
5&2,3,4&2065&2181\\
6&2,3,4&3360&2270\\
\bottomrule
\end{tabular}
\end{table}

\clearpage
\section{The symmetric quintic triangle}\label{sec:quintic-triangle}
Write $O_t$ for the orbit of $(t,t,1-2t)$. In Bos's construction
\cite[\S5, p.~288]{bos2023}, $u,v$ are the roots in
$(0.1480,0.1481)$ and $(0.4208,0.4209)$ of the polynomial
\begin{align*}
 P(t)={}&2737800t^{10}-8660340t^9+11981963t^8-9523289t^7+4800577t^6\\
        &-1598001t^5+354286t^4-51415t^3+4649t^2-235t+5 .
\end{align*}
All ten roots are isolated by the alternating signs of $P(10^{-4}j)$
at $j=0$, $1000$, $1480$, $1481$, $2550$, $3000$, $4208$, $4209$,
$4580$, $5000$, $6000$.

\begin{proposition}\label{prop:quintic-symmetric}
Among 21-point configurations consisting of the three vertices, four
further points on each edge, and six interior points invariant as a set
under all permutations of the vertices, Bos's quintic nodes uniquely
maximize the absolute Vandermonde determinant, up to ordering.
\end{proposition}
\begin{proof}
For the boundary factor preceding \eqref{eq:weighted-triangle}, use the
three vertex coordinates and the four functions $x_ix_jx_i^k(x_i+x_j)^{3-k}$,
$0\le k\le3$, on each edge $ij$. Write $\tau_1,\ldots,\tau_4$ for the
$x_i$-coordinates of its four interior edge nodes. The edge block has
absolute determinant $\prod_a\tau_a(1-\tau_a)\prod_{a<b}|\tau_a-\tau_b|$.
This univariate Vandermonde product with endpoints $0,1$ is uniquely
maximized by the degree-five Gauss--Lobatto set
\cite[\S2]{taylorwingatevincent}.

An interior permutation orbit has size one, three or six, with only
the centroid giving size one. On a single six-point orbit $e_2$ is
constant, so the weighted quadratic determinant vanishes. Every
nonsingular symmetric choice therefore has the form $O_s\cup O_t$,
$0<s<t<1/2$. Write $w_t=t^2(1-2t)$. Bos's factorization
\cite[\S5, pp.~287--288]{bos2023} gives \eqref{eq:weighted-triangle}
as a fixed positive constant times
\[
 K(s,t)=\left|[w_sw_t(t-s)]^3(1-3s)^2(1-3t)^2(2-3s-3t)\right|.
\]
Set $\sigma=s+t$ and $\tau=st$. The boundary of the closed
parameter triangle has $K=0$. At a positive maximum the stationary
equations are $G(s,t)=G(t,s)=0$, where
\[
 G(s,t)=-\tfrac16s(1-2s)(s-t)(3s-1)(3\sigma-2)\,\partial_s\log K(s,t).
\]
Taking their difference divided by $s-t$, then substituting in their
sum, gives
\begin{align*}
 \tau&=\frac{((45\sigma-59)\sigma+27)\sigma-5}{84\sigma-55},\\
 J(\sigma)&=52650\sigma^5-166545\sigma^4+207899\sigma^3
             -127879\sigma^2+38700\sigma-4600=0 .
\end{align*}
The numerator in the first line is $3205/197568$ at $\sigma=55/84$.
Substitution in the sum after multiplication by $(84\sigma-55)^2$
gives $-2(3\sigma-2)J(\sigma)$, and $3\sigma=2$ gives $K=0$.
The value to compare is
\[
 K^2=(\sigma^2-4\tau)^3\tau^{12}(1-2\sigma+4\tau)^6
              (1-3\sigma+9\tau)^4(3\sigma-2)^2 .
\]
The intervals $10^{-6}(j,j+1)$ for
$j=398396,568845,650697,705430,839877$ each have opposite endpoint
signs for $J$ and account for all its roots. Substitution in the
formulas for $\tau$ and $K^2$ gives $10^{28}K^2<1/16$ in the
first, fourth and fifth intervals. In the third,
$(1-2s)(1-2t)=1-2\sigma+4\tau<0$, so this interval is inadmissible.
Since $10^{28}K(3/20,21/50)^2>1$, the maximum occurs in the second
interval. There $x^2-\sigma x+\tau$ has signs $+,-,-,+$ at
$0.1480,0.1481,0.4208,0.4209$, respectively. Eliminating $s$ from the
two stationary polynomials gives
\[
 \operatorname{Res}_{s}\bigl(G(s,t),G(t,s)\bigr)
    =3t(2t-1)(3t-1)^4P(t).
\]
At a positive maximum none of the factors preceding $P(t)$ vanishes.
Thus $P(t)=0$ and, by symmetry, $P(s)=0$. The isolating intervals
for Bos's decic now identify the unique maximizing pair as $(u,v)$.
\end{proof}

\clearpage
\begingroup
\small
\begin{thebibliography}{99}
\setlength{\itemsep}{0pt}
\raggedright
\bibitem{bos1983}
L.~Bos, Bounding the Lebesgue function for Lagrange interpolation in a
simplex, \emph{J. Approx. Theory} \textbf{38} (1983), 43--59.
\href{https://doi.org/10.1016/0021-9045(83)90140-5}{doi:10.1016/0021-9045(83)90140-5}.

\bibitem{bos1990}
L.~Bos, Some remarks on the Fej\'er problem for Lagrange interpolation in
several variables, \emph{J. Approx. Theory} \textbf{60} (1990), 133--140.
\href{https://doi.org/10.1016/0021-9045(90)90078-5}{doi:10.1016/0021-9045(90)90078-5}.

\bibitem{bos1991}
L.~Bos, On certain configurations of points in $\R^n$ which are
unisolvent for polynomial interpolation, \emph{J. Approx. Theory}
\textbf{64} (1991), 271--280.
\href{https://doi.org/10.1016/0021-9045(91)90063-G}{doi:10.1016/0021-9045(91)90063-G}.

\bibitem{bos2023}
L.~Bos, On Fekete points for a real simplex, \emph{Indag. Math. (N.S.)}
\textbf{34} (2023), 274--293.
\href{https://doi.org/10.1016/j.indag.2022.11.003}{doi:10.1016/j.indag.2022.11.003}.
Preprint version: \href{https://arxiv.org/html/2205.06498v1}{arXiv:2205.06498v1}.

\bibitem{boydvandenberghe}
S.~Boyd and L.~Vandenberghe, \emph{Convex Optimization},
Cambridge University Press, Cambridge, 2004.
\href{https://stanford.edu/~boyd/cvxbook/}{Authors' online edition}.

\bibitem{fejer1932}
L.~Fej\'er, Bestimmung derjenigen Abszissen eines Intervalles, f\"ur
welche die Quadratsumme der Grundfunktionen der Lagrangeschen
Interpolation im Intervalle ein m\"oglichst kleines Maximum besitzt,
\emph{Ann. Scuola Norm. Super. Pisa Cl. Sci.} (2) \textbf{1} (1932),
263--276. \href{https://www.numdam.org/article/ASNSP_1932_2_1_3_263_0.pdf}{Numdam}.

\bibitem{garloff}
J.~Garloff, The Bernstein algorithm, \emph{Interval Comput.}
(1993), no.~2, 154--168.
\href{https://interval.louisiana.edu/reliable-computing-journal/1993/interval-computations-1993-2-pp-154-168.pdf}{Journal PDF}.

\bibitem{hornjohnson}
R.~A.~Horn and C.~R.~Johnson, \emph{Matrix Analysis}, 2nd ed.,
Cambridge University Press, Cambridge, 2013.
\href{https://doi.org/10.1017/CBO9781139020411}{doi:10.1017/CBO9781139020411}.

\bibitem{jamesstein}
W.~James and C.~Stein, Estimation with quadratic loss, in
\emph{Proc. Fourth Berkeley Symp. Math. Statist. Probab.}, vol.~1,
University of California Press, Berkeley, 1961, 361--379.
\href{https://www.math.tau.ac.il/~yekutiel/eBayes/james_stein_1961.pdf}{Paper}.

\bibitem{kieferwolfowitz}
J.~Kiefer and J.~Wolfowitz, The equivalence of two extremum problems,
\emph{Canad. J. Math.} \textbf{12} (1960), 363--366.
\href{https://doi.org/10.4153/CJM-1960-030-4}{doi:10.4153/CJM-1960-030-4}.

\bibitem{orlitsky}
A.~Orlitsky, N.~P.~Santhanam, K.~Viswanathan and J.~Zhang,
On modeling profiles instead of values, in \emph{Proceedings of the
Twentieth Conference on Uncertainty in Artificial Intelligence},
AUAI Press, 2004, 426--435.
\href{https://arxiv.org/abs/1207.4175}{arXiv:1207.4175}.

\bibitem{pandapal}
M.~K.~Panda and M.~Pal, D-optimal designs for quartic polynomial model
for mixture experiment with three and four ingredients,
\emph{Calcutta Stat. Assoc. Bull.} \textbf{78} (2026), no.~1, 79--88.
Published online 15 September 2025.
\href{https://doi.org/10.1177/00080683251366118}{doi:10.1177/00080683251366118}.

\bibitem{powersreznick}
V.~Powers and B.~Reznick, A new bound for P\'olya's theorem with
applications to polynomials positive on polyhedra,
\emph{J. Pure Appl. Algebra} \textbf{164} (2001), 221--229.
\href{https://doi.org/10.1016/S0022-4049(00)00155-9}{doi:10.1016/S0022-4049(00)00155-9}.

\bibitem{riener}
C.~Riener, On the degree and half-degree principle for symmetric
polynomials, \emph{J. Pure Appl. Algebra} \textbf{216} (2012), 850--856.
\href{https://doi.org/10.1016/j.jpaa.2011.08.012}{doi:10.1016/j.jpaa.2011.08.012}.
\href{https://arxiv.org/html/1001.4464v2}{Preprint with the cited numbering}.

\bibitem{taylorwingatevincent}
M.~A.~Taylor, B.~A.~Wingate and R.~E.~Vincent, An algorithm for computing
Fekete points in the triangle, \emph{SIAM J. Numer. Anal.} \textbf{38}
(2000), 1707--1720.
\href{https://doi.org/10.1137/S0036142998337247}{doi:10.1137/S0036142998337247}.

\bibitem{timofte}
V.~Timofte, On the positivity of symmetric polynomial functions.
Part I: General results, \emph{J. Math. Anal. Appl.} \textbf{284}
(2003), 174--190.
\href{https://doi.org/10.1016/S0022-247X(03)00301-9}{doi:10.1016/S0022-247X(03)00301-9}.
\end{thebibliography}
\endgroup
\end{document}