The energy and the main results

For \(R>1\), let \(K_R=[-R,R]^2\). Fix smooth cutoffs \(\chi_R\in C_c^\infty(\mathbb R^2)\) satisfying \[ 0\leq\chi_R\leq1,\qquad \mathop{\mathrm{supp}}\chi_R\subset K_R,\qquad \chi_R=1\text{ on }K_{R-1},\qquad \sup_{R>1}\|\nabla\chi_R\|_\infty<\infty. \tag{1}\] Let \(\Lambda\subset\mathbb R^2\) be locally finite and simple, so every point has multiplicity one, and put \(\nu_\Lambda=\sum_{p\in\Lambda}\delta_p\). The admissible class \(\mathcal A_1\) consists of locally integrable vector fields \(E\) such that \[ \operatorname{div}E=2\pi(\nu_\Lambda-1),\qquad \mathop{\mathrm{curl}}E=0,\qquad \sup_{R>1}\frac{\nu_\Lambda(B_R)}{|B_R|}<\infty, \tag{2}\] where the equations hold in distributions and \(B_R\) is the disk of radius \(R\) centered at the origin. These conditions imply \(E\in L^q_{\mathrm{loc}}(\mathbb R^2;\mathbb R^2)\) for every \(1\leq q<2\): locally subtract the finitely many Coulomb singularities and a smooth field with the prescribed constant divergence; the remaining field has harmonic components. Section [sec:reduction] supplies the details. Consequently the same class results if local \(L^q\) integrability is required for any one, or for all, of these exponents.

For \(E\in\mathcal A_1\) define \[\begin{align*} W(E,\chi_R) &=\lim_{\eta\downarrow0}\left[ \frac12\int_{\mathbb R^2\setminus\bigcup_{p\in\Lambda}B(p,\eta)} \chi_R|E|^2\,\mathrm dx +\pi\log\eta\sum_{p\in\Lambda}\chi_R(p)\right],\tag{3}\\ W(E)&=\limsup_{R\to\infty}\frac{W(E,\chi_R)}{|K_R|}. \tag{4}\end{align*}\] The local puncture limit exists and is finite. The outer limsup is initially allowed to take either infinite value. The puncture limit is taken first, and the cutoff layer has the fixed width in (1) throughout.

Identify \(\mathbb R^2\) with \(\mathbb C\) when convenient. The triangular lattice of covolume one is \[ \Lambda_\triangle =\sqrt{\frac{2}{\sqrt3}} \left(\mathbb Z(1,0)+\mathbb Z\left(\frac12,\frac{\sqrt3}{2}\right)\right). \tag{5}\] Let \(h_\triangle\) be the mean-zero periodic distributional solution of \(-\Delta h_\triangle=2\pi(\delta_0-1)\) on \(\mathbb R^2/\Lambda_\triangle\), and set \(E_\triangle=-\nabla h_\triangle\).

In particular \(W(E)=-\infty\) is impossible in this class. The inequality includes fields of energy \(+\infty\). No separation, periodicity, or Voronoi-cell bound is imposed on an arbitrary competitor.

The finite-volume statement used to prove Theorem 1 has its own exact normalization. For an integer \(n\geq2\), write \(T_n=\mathbb R^2/(\sqrt n\,\mathbb Z^2)\). Given distinct points \(a_1,\ldots,a_n\in T_n\), let \(h\) have mean zero and satisfy \(-\Delta h=2\pi(\sum_{i=1}^n\delta_{a_i}-1)\) on \(T_n\). Define \[ \mathcal W_n(h)=\lim_{\eta\downarrow0} \left[\frac12\int_{T_n\setminus\bigcup_{i=1}^n B(a_i,\eta)} |\nabla h|^2\,\mathrm dx+\pi n\log\eta\right]. \tag{6}\]

This is a lower bound on square tori. It does not assert equality on any particular square torus. Constant harmonic fields increase the energy by \(n|c|^2/2\); the required orthogonality follows by integrating the periodic gradient, as proved in (20).

The planar minimum also determines the linear term in optimal logarithmic energy for finitely many points on the sphere.

This proves the \(d=2\) case of the logarithmic-energy conjecture of Brauchart, Hardin, and Saff (Brauchart et al. 2012, Conjecture 4). Bétermin and Sandier had already proved the existence of the linear coefficient and shown that it equals (8) precisely when the density-one triangular lattice minimizes their planar renormalized energy (Bétermin and Sandier 2018, Theorem 1.5). Section 8 checks the normalization and applies that equivalence to Theorem 1. The conclusion is an asymptotic value, not a construction of near-optimal point sets or an algorithm for Smale’s seventh problem.

Context and prior work

The finite-vortex renormalized energy of Bethuel, Brezis, and Hélein (Bethuel et al. 1994) is an antecedent of the infinite-system energy considered here. Sandier and Serfaty introduced the planar functional in their analysis of microscopic vortex patterns in the Ginzburg–Landau model. They proved existence of a minimum, approximation of its value by square-periodic fields, and triangular optimality among Bravais lattices (Sandier and Serfaty 2012, Theorems 1 and 2). Their Conjecture 1 asks whether the triangular value is also the minimum over the full admissible class. Theorem 1 resolves this conjecture positively. The current normalization in that paper has one point per area \(2\pi/m\); our density-one convention corresponds to its background \(m=2\pi\).

The lattice comparison belongs to a longer history. The minimization of convergent Epstein zeta sums at fixed covolume was developed by Rankin (Rankin 1953) and Cassels (Cassels 1959), with an error noted in (Cassels 1963) and a repair supplied by Ennola (Ennola 1964) and simplified by Diananda (Diananda 1964). Montgomery proved triangular minimality of the lattice theta function for every positive Gaussian parameter (Montgomery 1988, Theorem 1). These are comparisons within the lattice class; the issue here is the lower bound for arbitrary admissible configurations.

For the confining-potential class in their Theorem 2, Sandier and Serfaty show that microscopic limits of finite Coulomb minimizers concentrate on minimizers of the planar renormalized energy (Sandier and Serfaty 2015, Theorem 2). Its minimum also determines a next-order coefficient for optimal logarithmic energy on the sphere (Bétermin and Sandier 2018, Theorem 1.5). Petrache and Serfaty prove a conditional implication from universal optimality to minimality for Coulomb and Riesz renormalized energies in their configuration-energy normalization (Petrache and Serfaty 2020, Theorem 2 and Corollary 4). Their periodic comparison uses tori compatible with the reference lattice, not arbitrary square tori, and their planar conclusion is conditional on the Cohn–Kumar conjecture. The companion universal-optimality theorem supplies the Gaussian comparison used in that implication (OpenAI 2026, Theorem 1.1 and Corollary 1.3). The present proof treats the logarithmic functional directly and uses no universal-optimality result as an input.

Two geometric precedents help explain the finite-torus argument. Screening a point charge by a uniform disk yields a separation argument for planar confined logarithmic minimizers; see Lieb, Rougerie, and Yngvason (Lieb et al. 2018, Lemma 3 and Appendix A). Voronoi moment inequalities combine sharp cell comparisons with Euler’s bound on the average number of sides, as in Gruber’s proof of a theorem of Fejes Tóth (Gruber 1999). Affine bounds on cell energies, sharp at the regular hexagon, also appear in the Wasserstein model of Bourne, Peletier, and Theil for cells with at least three sides and mass \(v\ge1.5\cdot10^{-4}\) (Bourne et al. 2014, Lemma 12). Here the cell comparison must arise from a trial potential that is continuous off the sites and has the prescribed logarithmic singularities. Its nonradial correction and the two scalar estimates of Section [sec:calibration] address that requirement, including sectors whose perpendicular foot lies outside the corresponding edge.

The whole-plane reduction below uses Sandier and Serfaty’s local energy estimates and periodic approximation. The finite-torus argument provides the comparison with the triangular value. Throughout, the conclusion concerns energy per unit area: compact changes to a configuration can preserve this value, so the minimum statement alone does not classify all minimizing configurations.

How the comparison works

The proof separates a finite-torus comparison from the passage to arbitrary fields. The finite comparison identifies the triangular value. Sandier and Serfaty’s approximation theorem supplies periodic fields approaching the unknown full-plane minimum; it does not identify that minimum.

Section [sec:reduction] checks the conversion from that approximation theorem to our normalization and fixed-width square cutoffs. The main point is to control the boundary error for an arbitrary field whose energy limsup is bounded above, without assuming separation or density one in advance. Local energy estimates give the needed charge density and boundary-strip bound. Proposition 7 then transfers any common lower bound on square tori to every field in \(\mathcal A_1\), including the possible infinite-energy cases.

We next prove the finite comparison. On each square torus the point energy attains a minimum. A screened logarithmic potential and the minimum principle separate the minimizing points, as shown in Section [sec:geometry]. Separation is used only for that minimum; once it is bounded below, every other configuration is bounded below by minimality. This is why the main theorems impose no separation condition on their competitors.

The lifted minimizer has convex Voronoi cells, the regions closer to one site than to any other. Each edge and its incident cell determine a sector joining that edge to the site. Counting incidences on the torus gives at most six sectors per point on average. The comparison is not an estimate of the electrostatic energy of an isolated cell. Instead, we construct a single trial potential \(H\) on the torus. Its sector formulas have matching traces and satisfy \(H=-\log r+O(r^2)\) near each site. Thus \(H\) and the Green potential \(h\) have matching logarithmic singularities, while the constant term is normalized to zero in \(H\). Proposition 12 gives the exact identity \[\mathcal W_n(h) =\sum_S J(S)+\frac12\int_{T_n}|\nabla(h-H)|^2.\] Here \(J(S)\) is the trial potential’s sector contribution: its linear background-interaction term minus one half of its renormalized Dirichlet integral, as defined in (31). The identity cancels the matching logarithmic singularities and retains the nonnegative square exactly.

The construction of \(H\) uses radial functions common to every sector. At a shared Voronoi edge its trace depends only on the equal distances to the two neighboring sites; along a ray to a vertex it depends only on the vertex radius. These observations glue the local formulas continuously off the sites and give an \(H^1\) function on every punctured torus. The identity is proved for any compatible flat torus whose area equals its number of sites, not just for square tori. A separate one-site triangular torus supplies six regular sectors; there the chosen potential has exactly the triangular Green gradient, so the square term vanishes.

Section [sec:calibration] constructs the radial functions from the triangular Green expansion and subtracts fixed multiples of sector area and angle. The desired bound has the form \[J(S)\ge K+\lambda A(S)+\mu\theta(S),\qquad K<0.\] Areas and angles sum exactly to \(n\) and \(2\pi n\) on a square torus. The incidence count then gives \(n(6K+\lambda+2\pi\mu)\), and exact calibration on the triangular torus identifies this with \(nW(E_\triangle)\).

The perpendicular from a site to a Voronoi edge need not lie on the edge, so one endpoint angle of a sector can be negative. A lower bound for a scalar primitive alone does not control the resulting difference of two primitive values. A second estimate bounds the total negative part of its derivative. Together these estimates handle both ordinary and signed sectors without imposing a regular-polygon hypothesis.

Section [sec:tails] proves summable analytic bounds for all omitted lattice modes. Section 7 first verifies the retained finite core and then uses the analytic error bounds to enclose the full functions and their derivatives. Midpoint and interpolation errors cover a compact rectangle between mesh points. At the exact triangular point, analytic stationarity and a positive Hessian bound establish the non-strict minimum; an exact formula handles the unbounded exterior. Section 8 combines these results with the geometry and the whole-plane reduction, then derives Corollary 3. Appendix 9 records two complementary methods: an independent coefficient recurrence and sharper exact tail bounds. They preserve additional quantitative information without adding premises to the main comparison.

Local regularity and the puncture limit

The local integrability convention causes no restriction on the competitors: under the distributional equations, membership in \(L^1_{\mathrm{loc}}\) already implies membership in every \(L^q_{\mathrm{loc}}\), \(1\le q<2\).

In particular the only possible infinite values of \(W(E)\) arise in the large-square limsup. No puncture limit is interchanged with that limsup below.

Approximation and local energy estimates

We distinguish the current normalization from the electric-field normalization of \(\mathcal A_1\). For \(b>0\), write \(\mathcal C_b\) for currents \(j\in L^1_{\mathrm{loc}}(\mathbb R^2;\mathbb R^2)\) with \[\operatorname{div}j=0,\qquad \operatorname{curl}j=\rho-b,\qquad \rho=2\pi\sum_{v\in\Lambda}\delta_v, \qquad \rho(B_s)=O(s^2),\] where \(\Lambda\) is simple and locally finite. Denote their local puncture functional by \(\mathsf W(j,\xi)\); its formula is the same as that of \(W(E,\xi)\). Put \(Q_s=(-s,s)^2\).

For such a current, the inverse component rotation \(F=(j_2,-j_1)\) satisfies \(\operatorname{div}F=\rho-b\) and \(\operatorname{curl}F=0\). In the proof of Lemma 4, use the local model \(\nabla\log|x-p|-(b/2)(x-p)\) near a charge and \(-(b/2)x\) off the charges. This gives the radial singular expansion for \(F\) and its rotated expansion for \(j\). Orthogonality preserves the squared norm, so the same finite local puncture conclusion holds for \(j\) and every smooth compactly supported cutoff. This local argument uses simplicity and local finiteness, but not the charge-growth bound.

For balls and for \(Q_s\), choose cutoffs \(0\le\zeta_{U_s}\le1\), compactly supported in \(U_s\), equal to one at distance at least one from its boundary, and with uniformly bounded gradients. Write \[\mathsf W_U(j)=\limsup_{s\to\infty} \frac{\mathsf W(j,\zeta_{U_s})}{|U_s|}.\] Theorem 1(1),(3),(4) of (Sandier and Serfaty 2012) makes the cutoff choice immaterial within each fixed shape family, gives a common finite minimum \[ M=\min_{\mathcal C_1}\mathsf W_{\mathrm{ball}} =\min_{\mathcal C_1}\mathsf W_Q\in\mathbb R, \tag{10}\] and provides square-periodic \(j_k\in\mathcal C_1\) with \(\mathsf W_{\mathrm{ball}}(j_k)\to M\). Their square periods have side lengths \(\sqrt{2\pi n_k}\) for positive integers \(n_k\). The common minimum does not assert equality of the ball and square values of an arbitrary fixed field. Their conditions (1.4)–(1.5) require increasing bounded open sets with \(\bigcap_{s>0}\overline U_s=\{0\}\), left-continuous volumes, \(|U_s-U_s|\le C|U_s|\), and, for every fixed \(a\in\mathbb R^2\), \[\frac{|(a+U_s)\mathbin\triangle U_s|}{|U_s|}\longrightarrow0.\] They further require constants \(C\) and \(\theta<2\) such that \[U_s+B_1\subset U_{s+C},\qquad U_{s+1}\subset U_s+B_C,\qquad |U_{s+1}\setminus U_s|=O(s^\theta).\] Balls and squares satisfy these conditions with \(\theta=1\) and \(|U_s-U_s|=4|U_s|\).

We also use the following nonsharp consequences of Lemma 4.7 and Proposition 4.9 of (Sandier and Serfaty 2012). The latter is an application of the mass-displacement method developed for Ginzburg–Landau energies in (Sandier and Serfaty 2011). In the proof of their Lemma 4.7, Sandier and Serfaty credit the dyadic-shell argument for the local \(L^r\) estimate to Struwe (Struwe 1994); the precise weighted estimate used here is theirs. Let \(U\) be bounded and open, \(\widehat U=U+B_1\), and suppose \(\operatorname{div}j=0\) and \(\operatorname{curl}j=\rho-1\) in \(\widehat U\). The number \(n\) of simple charges there is finite, and \(j\) is locally square-integrable off those charges. For \(1<r<2\), \(0\le\xi\le1\) compactly supported in \(U\), and \(\|\nabla\xi\|_\infty\le C_0\), Lemma 4.7 implies \[ \int_U\xi^{r/2}|j|^r \le C_{r,C_0}(|U|+1)^{1-r/2} \bigl(\mathsf W(j,\xi)_++n\log(2+n)+1\bigr)^{r/2}, \tag{11}\] where \(a_+=\max\{a,0\}\). For \(n\ge1\), this follows by replacing the energy by its positive part, enlarging the logarithmic count, and absorbing the fixed cutoff bounds into the constant. The statement of Lemma 4.7 and the final display of its proof differ by a factor two on the energy term; this weakening covers either form. If \(n=0\), there is no puncture term and Hölder’s inequality gives directly \[\int_U\xi^{r/2}|j|^r \le |U|^{1-r/2}\left(\int_U\xi|j|^2\right)^{r/2} =|U|^{1-r/2}\bigl(2\mathsf W(j,\xi)\bigr)^{r/2},\] which also implies (11). Proposition 4.9 provides one measure \(g_U\), supported in \(\widehat U\) and satisfying \(g_U\ge-C\,\,\mathrm dx\), such that for every smooth compactly supported \(\xi\) in \(U\), \[ \left|\mathsf W(j,\xi)-\int\xi\,\,\mathrm dg_U\right| \le C n_\xi\bigl(\log(2+n_\xi)+1\bigr) \|\nabla\xi\|_\infty, \tag{12}\] where \[n_\xi=\#\{v\in\Lambda: B(v,C)\cap\operatorname{supp}\nabla\xi\ne\varnothing\}.\] The constants in (12) are universal. In particular the same \(g_U\) works for two different cutoffs. The local smoothness needed in these results follows as in Lemma 4, after rotating the components.

Scaling and comparison of fixed-width cutoffs

Periodic averaging and the harmonic component

For a periodic field \(F\) with one unit logarithmic singularity at each of \(N\) charges in a period cell \(P\), define its cell energy by \[\mathcal E_P(F)=\lim_{\eta\downarrow0} \left[\frac12\int_{P\setminus\bigcup_v B(v,\eta)}|F|^2\,\,\mathrm dx +\pi N\log\eta\right].\] The integral is interpreted on the flat torus, so disks meeting a chosen cell boundary are identified periodically.

On \(T_n=\mathbb R^2/(\sqrt n\,\mathbb Z^2)\), let \(h\) be the mean-zero solution of \(-\Delta h=2\pi(\sum_{i=1}^n\delta_{v_i}-1)\) for distinct points \(v_1,\ldots,v_n\), with energy defined by (6). The compatibility condition holds because \(|T_n|=n\); the mean-zero solution exists and is unique by the Fourier series for the Laplacian on the torus. Any periodic locally integrable curl-free field \(E\) with divergence \(-\Delta h\) differs from \(-\nabla h\) by a smooth field of zero divergence and curl. Each component is harmonic on the torus, hence constant. Therefore \(E=-\nabla h+c\), with \(c\in\mathbb R^2\). As a distributional derivative of a periodic function, \(\nabla h\) has integral zero; it is in \(L^1(T_n)\) by Lemma 4. The cross integral on the punctured torus thus tends to zero, giving the exact identity \[ \mathcal E_{T_n}(E)=\mathcal W_n(h)+\frac n2|c|^2. \tag{20}\] Applying the same fixed orthogonal rotation to the fields preserves this energy identity.

Green representation and existence

By (20), every admissible periodic field with these charges differs from \(-\nabla h\) by a constant vector \(c\), which adds \(n|c|^2/2\) to the energy. It therefore suffices to minimize over the mean-zero potentials considered here.

Separation of minimizers

The argument uses a screened logarithmic potential and the minimum principle. A closely related exclusion argument for minimizing configurations of a planar confined system, with the same single-charge screening radius, appears in Lieb, Rougerie, and Yngvason (Lieb et al. 2018, Lemma 3 and Appendix A). The following proof treats the torus Green function and its background directly.

Counting sectors on the torus

For \(d>0\) and real \(\alpha,\beta\) with \(|\alpha|,|\beta|<\pi/2\) and \(\alpha+\beta>0\), define the sector \[ S(d,\alpha,\beta)= \{(r,\varphi):-\alpha<\varphi<\beta,\ 0<r<d\sec\varphi\}. \tag{22}\] Its polar axis is the normal from its site toward the supporting edge.

For reference, the area and angle of (22) are \[ A(S)=\frac12\int_{-\alpha}^{\beta}d^2\sec^2\varphi\,\,\mathrm d\varphi, \qquad \theta(S)=\alpha+\beta. \tag{23}\] In particular the sector has positive angle even when one endpoint parameter is negative.

Compatible flat-torus data and trace matching

Let \(\Gamma\subset\mathbb R^2\) be a rank-two lattice and put \(\mathbb T=\mathbb R^2/\Gamma\). Let \(N\geq1\) be an integer, let \(v_1,\ldots,v_N\) be distinct points of \(\mathbb T\), and let \(P\) be their full \(\Gamma\)-periodic lift. Suppose that the Voronoi cells of \(P\) are bounded convex polygons with finitely many face–edge incidences modulo \(\Gamma\). Assume that each incidence has the sector form (22) and that, for one named constant \(d_0>0\), every sector has \(d\geq d_0\). We call these compatible periodic Voronoi data with lower edge distance \(d_0\). This is conditional geometric data, not a separation assertion about an arbitrary competitor.

The sectors partition the \(N\) cells modulo \(\Gamma\), up to their boundary rays. Their total area is \(|\mathbb T|\), and their total angle is \(2\pi N\): each convex cell has its site in its interior because \(d_0>0\), so its sectors have total angle \(2\pi\). These facts do not require an incidence bound.

Call a radius \(\eta>0\) admissible if it is smaller than the injectivity radius of \(\mathbb T\) and, when \(N\geq2\), smaller than half the minimum torus distance between distinct sites. Thus the torus disks \(B(v_i,\eta)\) are embedded and pairwise disjoint; for \(N=1\) there is no pairwise-distance condition. Write \[\Omega_\eta=\mathbb T\setminus\bigcup_{i=1}^N B(v_i,\eta).\] In puncture limits we may further take \(\eta<d_0\), so every deleted disk lies in its Voronoi cell. This restriction does not change a limit as \(\eta\downarrow0\).

Fix compatible periodic Voronoi data with lower edge distance \(d_0\). Let \(I\) be a finite or countable index set, and choose integers \(k_i\geq6\), \(i\in I\); different indices may have the same integer. Suppose that \[B\in C^1((0,\infty)),\qquad D_i\in C^1((0,\infty))\] satisfy the following assumptions: \[ B(r)=-\log r+O(r^2),\qquad B'(r)=-1/r+O(r) \quad(r\downarrow0), \tag{24}\] and, for every \(T_0>d_0\), \[ \sum_{i\in I}\ \sup_{d_0\leq t\leq T_0} \bigl((1+k_i)|D_i(t)|+t|D_i'(t)|\bigr)<\infty. \tag{25}\] These are the exact inputs needed in the gluing and duality argument. The later radial construction will verify them directly.

For a sector (22), define \(t(\varphi)=d\sec\varphi\) and \[ H(r,\varphi)=B(r)+ \sum_{i\in I}\left(\frac r{t(\varphi)}\right)^{k_i} D_i(t(\varphi)). \tag{26}\] The same radial functions are used in every sector.

The common notation \(D_i=b_i-u_i\), \(B=\sum_i u_i\) is one way to supply these data. In that notation the common edge trace is simply \(\sum_i b_i(r)\). The gluing argument itself only uses \(B\) and \(D_i\).

The exact dual identity

Assume now that \(|\mathbb T|=N\). Let \(h\) be the mean-zero distributional solution on \(\mathbb T\) of \[ -\Delta h=2\pi\left(\sum_{i=1}^N\delta_{v_i}-1\right). \tag{28}\] The area condition is exactly the compatibility condition for this equation. With the convention \[\Gamma^*=\{k\in\mathbb R^2:k\cdot\gamma\in\mathbb Z \text{ for every }\gamma\in\Gamma\},\] the Fourier mode \(e^{2\pi i k\cdot x}\) has eigenvalue \(4\pi^2|k|^2\) under \(-\Delta\). The constant Fourier coefficient of the right-hand side vanishes. Dividing its coefficient at each nonzero \(k\in\Gamma^*\) by \(4\pi^2|k|^2\), and setting the coefficient at zero to zero, gives a distributional solution. It is unique because the kernel of the Laplacian on the torus consists of constants. Define \[ \mathcal W_{\mathbb T}(h)=\lim_{\eta\downarrow0} \left[ \frac12\int_{\Omega_\eta}|\nabla h|^2\,\,\mathrm dx+\pi N\log\eta \right], \tag{29}\] where the limit is taken through admissible radii. The local calculation in Lemma 4, applied in a torus chart, gives the finite puncture limit. When \(\mathbb T=T_n\) and \(N=n\), this is exactly \(\mathcal W_n(h)\) in (6).

For a sector \(S=S(d,\alpha,\beta)\), define \[ \int_S^{\mathrm{ren}}|\nabla H|^2 =\lim_{\eta\downarrow0}\left[ \int_{S\cap\{r>\eta\}}|\nabla H|^2\,\,\mathrm dx +\theta(S)\log\eta\right]. \tag{30}\] Lemma 11 shows that this limit is finite. Its counterterm uses the positive angle \(\theta(S)\), also for a sector with a negative endpoint parameter. Define \[ J(S)=-2\pi\int_S H\,\,\mathrm dx -\frac12\int_S^{\mathrm{ren}}|\nabla H|^2. \tag{31}\] The sectorwise estimates in the proof of Lemma 11 also show that these expressions are well defined on any individual sector with \(d\geq d_0\) and the endpoint conditions in (22), without first choosing a periodic tiling.

If the gradient hypotheses were unchanged but the value expansion were instead \(H=-\log r+c_i+o(1)\) at \(v_i\), the right-hand side of the last identity would contain the additional term \(2\pi\sum_i c_i\). This explains why the zero constant term in the radial construction is part of the normalization.

For a square-torus minimizer, Lemma 10 supplies the compatible data with \(\Gamma=\sqrt n\,\mathbb Z^2\), \(N=n\), and \(d_0=d_*=r_*/2=1/(2\sqrt\pi)\). Its torus has area \(n\), so Proposition 12 applies with \(\mathcal W_{\mathbb T}=\mathcal W_n\). The triangular application in the next section supplies its own period and sectors and does not use a square-torus minimizer.

From scalar estimates to configuration bounds

The two scalar bounds below control ordinary and signed sectors by different arguments. Once they yield an affine lower bound for each sector, the dual identity and the incidence count give the energy bound. We formulate this implication for compatible flat tori first; separation of a square-torus minimizer then removes the geometric hypotheses from the finite-volume comparison.

The next section chooses particular radial data and exact constants \(K<0,\lambda,\mu\). For these data, the primitive identity (52) and the two scalar bounds yield \[J(S)\geq K+\lambda A(S)+\mu\theta(S)\] for every relevant sector. Proposition 18 applies the compatible-torus conclusion to this construction and identifies \(6K+\lambda+2\pi\mu\) with the triangular energy. Section 8 uses the square-torus conclusion, whose minimizer argument has already removed the separation hypothesis.

The local triangular Green function

Set \[ a=\sqrt{\frac{2}{\sqrt3}},\qquad p=\frac a2,\qquad R_0=\frac{2p}{\sqrt3},\qquad x_0=\frac{R_0}{2}. \tag{38}\] Thus \(p\) and \(R_0\) are the inradius and circumradius of the covolume-one triangular Voronoi hexagon, and \(x_0=p\tan(\pi/6)\). Identify the plane with \(\mathbb C\), and write \(\partial_z=(\partial_x-i\partial_y)/2\), with a nearest lattice point on the positive real axis. For \(m\ge1\) put \[ k_m=6m,\qquad c_m=\frac1{k_m}\sum_{z\in\Lambda_\triangle\setminus\{0\}}z^{-k_m}. \tag{39}\] These sums converge absolutely and are real by reflection symmetry. We also introduce an auxiliary index \(0\), with \(k_0=6\). It is distinct from index \(1\): both exponents equal \(6\), but they will carry different radial functions.

Matching the regular hexagon

Write \(\gamma(r)=-\log r+\pi r^2/2\) for \(r>0\). We first choose the radial data on the regular hexagon, where the test potential must equal \(h_\triangle-C_\triangle\). In a sector whose outward normal is the positive real axis, the edge has distance \(p\) and the endpoint angles are \(\pm\pi/6\). At an edge point of radius \(t=p\sec\varphi\), \[p+i\sqrt{t^2-p^2}=te^{i|\varphi|}.\] The boundary value of the \(m\)th Green-function mode is therefore \(c_mt^{k_m}\cos(k_m\varphi)\). Expressing it as a function of the radius gives the polynomial \[ V_m(r)=c_m\sum_{j=0}^{k_m/2}(-1)^j\binom{k_m}{2j} p^{k_m-2j}(r^2-p^2)^j, \qquad p\le r\le R_0. \tag{41}\] On either vertex ray, the same mode equals \((-1)^mc_mr^{k_m}\), since \(k_m\pi/6=m\pi\). This is the radial baseline we use for that mode inside the hexagon. The correction in (26) then recovers the angular dependence exactly: \[(-1)^mc_mr^{k_m} +(r/t)^{k_m}\bigl(V_m(t)-(-1)^mc_mt^{k_m}\bigr) =c_mr^{k_m}\cos(k_m\varphi).\] Thus the edge data \(b_m=V_m\) and the vertex-ray baseline determine the desired potential on each regular sector. Its remaining part \(\gamma\) is already radial and needs no correction there.

Continuation to arbitrary sectors

Other Voronoi sectors can have edge radii below \(p\) or above \(R_0\). We now extend the radial data to those radii, keeping two derivatives and making the correction vanish beyond a fixed radius. The assembled potential is required to agree with \(h_\triangle-C_\triangle\) only inside the regular hexagon; the continuation is chosen for the scalar estimates, not as an extension of that Green function.

Define the clipped polynomial \[S(y)=\begin{cases} 0,&y\le0,\\ 10y^3-15y^4+6y^5,&0<y<1,\\ 1,&y\ge1. \end{cases}\] Let \(h_s=3/40\) and \(c=3/50\). For \(0\le z\le h_s\) put \[P(z)=z-\frac{6z^3}{h_s^2}+\frac{8z^4}{h_s^3}-\frac{3z^5}{h_s^4},\qquad Q(z)=\frac{z^2}{2}-\frac{3z^3}{2h_s} +\frac{3z^4}{2h_s^2}-\frac{z^5}{2h_s^3}.\] For a function \(V\) with two left derivatives at \(R_0\), its continuation is defined for \(r\ge R_0\) by \[ (\mathcal HV)(r)= \begin{cases} V(R_0)+P(r-R_0)V'(R_0)+Q(r-R_0)V''(R_0),&r\le R_0+h_s,\\ V(R_0),&r>R_0+h_s. \end{cases} \tag{42}\] Indeed \((P,P',P'')(0)=(0,1,0)\) and \((Q,Q',Q'')(0)=(0,0,1)\); all six values at \(h_s\) vanish. Thus this operation preserves two derivatives at \(R_0\) and joins a constant with two derivatives at \(R_0+h_s\).

Define \(b_0=\gamma\) up to \(R_0\), continued by (42). For \(m\ge1\), set \(b_m=V_m\) on \([p,R_0]\) and use \(\mathcal HV_m\) above \(R_0\). Below \(p\) put \[ b_m(r)=\begin{cases} 0,&0<r\le p-c,\\ S((r-p+c)/c)\bigl(T_{0m}+T_{1m}(r-p)+T_{2m}(r-p)^2\bigr),&p-c<r<p, \end{cases} \tag{43}\] where \[\begin{split} T_{0m}&=c_mp^{k_m},\\ T_{1m}&=-\frac{k_m(k_m-1)}pT_{0m},\\ T_{2m}&=\frac{k_m(k_m-1)}{p^2} \left(\frac{(k_m-2)(k_m-3)}6-\frac12\right)T_{0m}. \end{split}\] Direct differentiation of (41) gives \(T_{0m}=V_m(p)\), \(T_{1m}=V_m'(p)\), and \(2T_{2m}=V_m''(p)\). Consequently every \(b_i\) is \(C^2\) away from the origin.

Now define \[U_m(r)=(-1)^mc_mr^{k_m}\bigl(1-S((r-R_0)/(1/10))\bigr),\qquad \tau(r)=S\left(\frac{r-R_0-1/1000}{3/40}\right),\] and set \[ \begin{split} u_0&=(1-\tau)\gamma+\tau b_0,\qquad u_m=(1-\tau)U_m+\tau b_m\quad(m\ge1),\\ B&=\sum_{i\ge0}u_i,\qquad D_i=b_i-u_i,\qquad l=R_0+\frac{19}{250}. \end{split} \tag{44}\] Here the auxiliary mode \(0\) lets the radial part \(\gamma\) join its constant continuation by the same correction mechanism as the other modes. Its exponent \(k_0=6\) makes this correction vanish at the pole without changing the singular normalization. The coefficient estimates proved below justify the sums and their first three derivatives on every open spline interval. In particular, \(B,D_i\) are \(C^2\) and piecewise \(C^3\) away from zero. Near zero, \(B=\gamma+\sum_{m\ge1}(-1)^mc_mr^{k_m}\); normal convergence of (40) therefore gives \[ \begin{gathered} B(r)=-\log r+O(r^2),\qquad B'(r)=-1/r+O(r)\quad(r\downarrow0),\\ D_i(r)=0,\qquad B(r)=B(l)\quad(r\ge l). \end{gathered} \tag{45}\] The fixed radial data also satisfy (25) for every \(d_0>0\). Indeed, there are only finitely many retained modes, while for the omitted modes (59) and the vanishing of \(D_i\) beyond \(l\) bound the summand on \([d_0,T_0]\) by \((1+k_i+T_0)w_{k_i}\). The explicit geometric terms in (58) make this series summable. This conclusion uses the analytic mode bounds, not the later finite-core estimates.

For a sector with edge distance \(d\) and endpoint angles \(-\alpha<\beta\), write \(t(\varphi)=d\sec\varphi\) and use the test potential \[ H(r,\varphi)=B(r)+\sum_{i\ge0} \left(\frac r{t(\varphi)}\right)^{k_i}D_i(t(\varphi)), \qquad 0<r<t(\varphi),\quad -\alpha<\varphi<\beta. \tag{46}\] At the edge its trace is \(\sum_i b_i(t)\); at a ray ending at a vertex of radius \(r_v\) its trace is \(B(r)+\sum_i(r/r_v)^{k_i}D_i(r_v)\). These expressions explain both trace matching conditions in the glued Voronoi construction.

The sector integral

The following formulas separate radial integration from sector geometry. For \(t>0\) define \[ \begin{split} B_1(t)&=B'(t)+1/t,\qquad S_1(t)=\int_0^t B(r)r\,dr,\\ S_2(t)&=\int_0^t\bigl(-2B_1(r)+rB_1(r)^2\bigr)\,dr,\\ H_k(t)&=t^{-k}\int_0^tB_1(r)r^k\,dr,\qquad Q_i(t)=tD_i'(t)-k_iD_i(t). \end{split} \tag{47}\] Here \(H_k\) is indexed by the exponent, so the two modes with exponent six use the same \(H_6\). They remain separate in every sum involving \(D_i\). Put \[ \begin{split} f_0(t)={}&-2\pi\left(S_1(t)+t^2\sum_i\frac{D_i(t)}{k_i+2}\right)\\ &-\frac12\Biggl(\log t+S_2(t) +\sum_i\bigl(-2D_i(t)+2k_iD_i(t)H_{k_i}(t)\bigr)\\ &\qquad+\sum_{i,j}\frac{k_ik_j}{k_i+k_j}D_i(t)D_j(t)\Biggr),\\ e(t)={}&\frac12\sum_{i,j}\frac{Q_i(t)Q_j(t)}{k_i+k_j}. \end{split} \tag{48}\] All sums here and below include index \(0\). The angular coefficient is nonnegative, since \[ e(t)=\frac12\int_0^1 \left(\sum_iQ_i(t)v^{k_i}\right)^2\frac{dv}{v}\ge0. \tag{49}\]

Let \(J\) denote the sector contribution to the dual energy, \[J=-2\pi\int_{\rm sector}H -\frac12\int_{\rm sector}^{\rm ren}|\nabla H|^2,\] where the renormalized integral adds \((\alpha+\beta)\log\eta\) before letting the radial exclusion \(\eta\) tend to zero. Then \[ J=\int_{-\alpha}^{\beta} \bigl(f_0(d\sec\varphi)-\tan^2\varphi\,e(d\sec\varphi)\bigr) \,d\varphi. \tag{50}\] To verify this identity, differentiate (46): \[H_r=-\frac1r+B_1(r)+\sum_i\frac{k_i}{r}(r/t)^{k_i}D_i(t),\qquad H_\varphi=\tan\varphi\sum_i(r/t)^{k_i}Q_i(t).\] The source integral gives \(S_1+t^2\sum_iD_i/(k_i+2)\). The radial square contributes \(\log t+S_2\), its cross terms contribute \(-2D_i+2k_iD_iH_{k_i}\), and its quadratic terms have coefficient \(k_ik_j/(k_i+k_j)\). Integrating the angular square gives \(\tan^2\varphi\sum Q_iQ_j/(k_i+k_j)\). These are precisely (48). The estimates below permit termwise integration.

Subtracting area and angle

Given real constants \(\lambda,\mu\), define \[ \begin{split} f(t)&=f_0(t)-\lambda t^2/2-\mu,\qquad t=\sqrt{d^2+s^2},\\ G(d,s)&=\frac d{t^2}f(t)-\frac{s^2}{dt^2}e(t),\qquad g_d(x)=\int_0^xG(d,s)\,ds\quad(d>0,\ x\in\mathbb R). \end{split} \tag{51}\] The evenness of \(G\) in \(s\) makes \(g_d\) odd. If \(A\) is the sector area and \(\theta=\alpha+\beta\), substitution \(s=d\tan\varphi\) in (50) gives \[ J-\lambda A-\mu\theta =g_d(d\tan\alpha)+g_d(d\tan\beta). \tag{52}\] This identity holds for signed endpoint angles; only \(|\alpha|,|\beta|<\pi/2\) and \(\alpha+\beta>0\) are required.

We choose \(\lambda,\mu\) so that the adjusted primitive is stationary at the regular sector \((d,x)=(p,x_0)\), and set \(K\) to twice its value there. The exact choices are \[ \begin{split} \lambda=\frac1{px_0}\biggl[&\frac p{x_0}\bigl(f_0(R_0)-f_0(p)\bigr) -\frac{x_0}{p}e(R_0)\\ &+p\int_0^{x_0}\frac{f_0(\sqrt{p^2+s^2})-f_0(p)}{s^2}\,ds +\frac1p\int_0^{x_0}e(\sqrt{p^2+s^2})\,ds\biggr],\\ \mu={}&f_0(R_0)-\frac13e(R_0)-\lambda R_0^2/2,\qquad K=2g_p(x_0). \end{split} \tag{53}\] The first integrand extends continuously to \(s=0\) with value \(f_0'(p)/(2p)\).

The two scalar estimates needed by the geometry

The remaining analytic task is to establish \[ \begin{gathered} K<0,\qquad 2g_d(x)\ge K\quad(d\ge b,\ x\ge0),\qquad b=\frac{28209}{100000},\\ \int_0^\infty\max\{0,-G(d,s)\}\,ds\le-K\quad(d\ge b). \end{gathered} \tag{55}\] The rational endpoint \(b=28209/100000\) is below the geometric threshold \(d_*=1/(2\sqrt\pi)\). Thus the domain \(d\ge b\) covers all edge distances supplied by the separated-minimizer geometry; the elementary comparison is justified in Section 8.

The second estimate is essential when the perpendicular from a point to a Voronoi edge falls outside that edge. Its role is made precise in the following implication, which completes the geometric part of the argument once (55) has been certified.

Coefficient estimates

Write a lattice point as \(z=a(n+q/2+i\sqrt3 q/2)\), with \((n,q)\in\mathbb Z^2\). Then \[|z/a|^2=n^2+nq+q^2\geq\frac34\max(|n|,|q|)^2.\] There are \(8v\) integer pairs with \(\max(|n|,|q|)=v\). For every \(d=6m\geq6\) and every positive integer \(D\), the sum outside the square \(|n|,|q|\leq D\), before multiplication by \(a^{-d}\), is therefore at most \[ 8(4/3)^{d/2}\sum_{v>D}v^{1-d} \leq 8(4/3)^{d/2}\frac{D^{2-d}}{d-2}. \tag{56}\] The last inequality follows by comparison with the integral from \(D\) to infinity. The first shell has six points with squared norm \(1\) and two with squared norm \(3\). Moreover, \[\sum_{v\geq2}v^{1-d} \leq 2^{1-d}+\int_2^\infty x^{1-d}\,\,\mathrm dx =2^{1-d}\left(1+\frac2{d-2}\right) \leq3\,2^{1-d}.\] It follows that, for \(d=k_m\geq12\), \[ |c_m|\leq\frac{1}{d a^d} \left(6+2\,3^{-d/2}+24(4/3)^{d/2}2^{1-d}\right) \leq\frac7{d a^d}. \tag{57}\] Indeed, the bracket equals \(6+50/729<7\) at \(d=12\), and its two nonconstant terms decrease with \(d\). The finite coefficient enclosures may use (56) with \(D=256\) for \(m=1\) and \(D=48\) for \(2\leq m\leq38\), including the factor \(1/(k_m a^{k_m})\). Reflection symmetry makes each coefficient real.

Uniform derivative bounds for omitted modes

For \(k=6m\geq234\), define the rational number \[ w_k=10^6\left( 30k^2(13/20)^k+3\cdot10^6 k^3 2^{-k} +5600k^5(3/5)^k+7\cdot10^{11}k^2(289/500)^k \right). \tag{58}\] For \(0\leq j\leq3\) and \(0<r\leq l\), the construction satisfies \[ |u_m^{(j)}(r)|\leq w_k, \qquad |D_m^{(j)}(r)|\leq w_k. \tag{59}\] At spline knots these inequalities mean both one-sided derivatives. Here are explicit checks of the constants.

The polynomial cutoff \(S(y)=10y^3-15y^4+6y^5\) on \([0,1]\), extended constantly outside that interval, obeys \[|S|\leq1,\qquad |S'|\leq2,\qquad |S''|\leq6,\qquad |S'''|\leq60.\] The identities \(p^4=1/12\) and \(R_0^4=4/27\) give the convenient rational bounds \[.537<p<.538,\quad .620<R_0<.621,\quad .696<l<.697, \quad l/a<13/20,\quad R_0/a<289/500.\] All terminating decimals in this section denote rational numbers.

For \(U_m(r)=(-1)^m c_mr^k(1-S((r-R_0)/.1))\), Leibniz’ rule and (57) bound the third derivative by \[7k^2(l/a)^k \left(l^{-3}+\frac{60}{kl^2} +\frac{1800}{k^2l}+\frac{60000}{k^3}\right).\] Using \(l>.69\) and \(k\geq234\), the coefficient of \(k^2(13/20)^k\) is less than \(25.445<30\). The same calculation for orders \(0,1,2\) gives coefficients below \(0.000001\), \(0.000197\), and \(0.070572\), respectively. This proves the first term inside the parentheses in (58).

For \(p-c\leq r\leq p\), write \(q_m(r)=T_{0m}+T_{1m}(r-p)+T_{2m}(r-p)^2\), so that \(b_m(r)=S((r-p+c)/c)q_m(r)\). Since \(p/a=1/2\) and \(p>1/2\), \[|T_{0m}|\leq\frac7k2^{-k},\qquad |T_{1m}|\leq14k2^{-k},\qquad |T_{2m}|\leq\frac{14}{3}k^3 2^{-k}.\] The last inequality uses \(((k-2)(k-3)/6)-1/2\leq k^2/6\). Consequently \[\begin{align*} |q_m|&\leq(7/k+.84k+.0168k^3)2^{-k},\\ |q_m'|&\leq(14k+.56k^3)2^{-k},\\ |q_m''|&\leq(28/3)k^3 2^{-k}. \end{align*}\] In particular, \[|b_m'''|\leq\frac{60}{c^3}|q_m| +\frac{18}{c^2}|q_m'|+\frac6c|q_m''| <8500k^3 2^{-k}<3\cdot10^6 k^3 2^{-k}.\] The lower derivative orders satisfy the same final bound, and below \(p-c\) the function is zero.

On \([p,p+.01]\), put \(w=r^2-p^2\) and \(X=p+\sqrt w\). The polynomial formula for \(V_m\) gives, for \(0\leq\ell\leq3\), \[\left|\frac{\,\mathrm d^\ell V_m}{\,\mathrm dw^\ell}\right| \leq |c_m| k^{2\ell}X^{k-2\ell}.\] For \(j\geq\ell\), the coefficient comparison that proves this is \[j(j-1)\cdots(j-\ell+1)\binom{k}{2j} \leq k^{2\ell}\binom{k-2\ell}{2j-2\ell}.\] The bound \(p>.537\) implies \(.02p+.0001<.04p^2\), hence \(X\leq1.2p=(3/5)a\) on this interval. Also \(r<.55\) and \(X>1/2\). The chain rule gives \[V_m'''(r)=8r^3\frac{\,\mathrm d^3V_m}{\,\mathrm dw^3} +12r\frac{\,\mathrm d^2V_m}{\,\mathrm dw^2}.\] Thus its absolute value is at most \[7\left(8(.55)^3 2^6+ \frac{12(.55)2^4}{234^2}\right)k^5(3/5)^k <597k^5(3/5)^k.\] The lower derivatives satisfy smaller bounds. The common bound \(5600k^5(3/5)^k\) therefore applies through order three.

On \([p+.01,R_0]\), use \(z=p+i\sqrt{r^2-p^2}\), so that \(V_m=\Re(c_mz^k)\) and \(|z|=r\). Direct differentiation gives \[|z'|\leq7,\qquad |z''|\leq300,\qquad |z'''|\leq60000.\] For example, \(\sqrt{r^2-p^2}>.103\), \(p<.54\), and \(r<.63\) suffice in the formulas \(|z'|=r/\sqrt{r^2-p^2}\), \(|z''|=p^2/(r^2-p^2)^{3/2}\), and \(|z'''|=3p^2r/(r^2-p^2)^{5/2}\). Since \(1/r<2\), the third derivative of \(c_mz^k\) is bounded by \[7\left(8\cdot343+\frac{12\cdot7\cdot300}{234} +\frac{2\cdot60000}{234^2}\right) k^2(289/500)^k <20000 k^2(289/500)^k.\] All derivative orders through three are therefore bounded by \(3\cdot10^6 k^2(289/500)^k\) on this interval.

For the extension \(\mathcal HV_m=V_m(R_0)+P(r-R_0)V_m'(R_0)+Q(r-R_0)V_m''(R_0)\), direct bounds on the coefficients of \(P,Q\) show, for \(0\leq j\leq3\), \[|P^{(j)}|\leq2\cdot10^5,\qquad |Q^{(j)}|\leq4000.\] Thus the common bound for the extension is at most \[(1+200000+4000)\,3\cdot10^6 k^2(289/500)^k <7\cdot10^{11} k^2(289/500)^k.\] It remains valid on the constant part beyond \(R_0+h_s\).

Finally the cutoff \(\tau(r)=S((r-R_0-.001)/.075)\) satisfies \[|\tau^{(j)}|\leq1,27,1067,142223 \qquad(j=0,1,2,3).\] The maximal Leibniz factor through order three is \(1+3\cdot27+3\cdot1067+142223=145506<10^6\). Applying this to \(u_m=(1-\tau)U_m+\tau b_m\) and \(D_m=(1-\tau)(b_m-U_m)\) proves (59). The series and their indicated derivatives consequently converge absolutely and locally uniformly. For \(r<b\), the omitted \(u_m\) derivatives additionally have their original \(O(r^{k-3})\) behavior.

Norms and bounds for retained modes

For a function \(F\) with one-sided derivatives as needed, define \[|F|_{j,t}=\sum_{q=0}^j\frac{|F^{(q)}(t)|}{q!}, \qquad \|F\|_{j,I}=\sup_{t\in I}|F|_{j,t}.\] These norms are submultiplicative. For nonnegative weights \(\omega_i\), a weighted mode norm means \(\sup_{t\in I}\sum_i\omega_i|F_i|_{j,t}\): the sum is taken before the supremum.

The finite interval verification in Section 7 establishes the following bounds for the retained modes: \[ |\widetilde B_1|_{2,t}<6000, \qquad \sum_{i=0}^{38}(k_i+1)^3|\widetilde D_i|_{3,t}<2\cdot10^6, \qquad b\leq t\leq l. \tag{60}\] The finite coefficient enclosure also supplies \(|c_1|<1\). These assertions are not conclusions of the analytic estimates in this section. Conditional on them, the remaining estimates below are independent of the finite calculation.

On \([0,b]\), \[\widetilde B_1(t)=\pi t+ \sum_{m=1}^{38}(-1)^m k_mc_mt^{k_m-1}\] has norm \(|\widetilde B_1|_{2,t}<10\). For an explicit check, use \(t<b<2/7\), \(a>1\), \(\pi<22/7\), and \(|c_1|<1\) for the first mode. The modes \(k\geq12\) contribute at most \(\sum_{k=12,18,\ldots}7k^2(2/7)^{k-3}<.014\); the first term and its geometrically decreasing ratios prove this last rational inequality. Together with the linear term and the \(k=6\) term this is less than \(10\).

Summing the omitted modes

Set \[ \mathfrak a=\sum_{k=234,240,\ldots}(k+1)^4w_k <2\cdot10^{-21},\qquad W_*=4\sum_{k=234,240,\ldots}w_k\leq4\mathfrak a. \tag{61}\] Here is a rational verification of the numerical inequality. For each term \((k+1)^4k^sr^k\), with \[(s,r)=(2,13/20),(3,1/2),(5,3/5),(2,289/500),\] the ratio at the next allowed index is \[\left(\frac{k+7}{k+1}\right)^4 \left(\frac{k+6}{k}\right)^s r^6.\] It decreases with \(k\), and at \(k=234\) it is strictly below \(1/10\) for all four pairs. The first weighted summand obeys the exact rational comparison \[\begin{split} 235^4\,10^6\bigl(&30\cdot234^2(13/20)^{234} +3\cdot10^6\cdot234^3 2^{-234}\\ &+5600\cdot234^5(3/5)^{234} +7\cdot10^{11}\cdot234^2(289/500)^{234}\bigr) <1.56\cdot10^{-21}. \end{split}\] Multiplying by \(10/9\) gives a bound less than \(2\cdot10^{-21}\), as asserted in (61).

Since \(\Delta B_1=\sum_{m\geq39}u_m'\), the derivative estimates give \[|\Delta B_1|_{2,t}\leq W_*,\qquad 0\leq t\leq l,\] and, on \([b,l]\), \[\sum_{i\geq39}(k_i+1)^j|D_i|_{3,t}\leq4\mathfrak a, \qquad 0\leq j\leq4.\] Indeed the divided-factorial sums use factors \(5/2\) and \(8/3\), respectively, both less than \(4\).

Propagation to the sector functions

Recall the quantities entering the sector formula: \[\begin{align*} B_1(t)&=B'(t)+1/t,& S_1(t)&=\int_0^t B(r)r\,\,\mathrm dr,\\ S_2(t)&=\int_0^t(-2B_1(r)+rB_1(r)^2)\,\,\mathrm dr,& H_k(t)&=t^{-k}\int_0^t B_1(r)r^k\,\,\mathrm dr,\\ Q_i(t)&=tD_i'(t)-k_iD_i(t). \end{align*}\] Write \(N=2\cdot10^6\) temporarily. Direct differentiation of \(\Delta S_1\) gives \[|\Delta S_1|_{2,t}\leq3W_*.\] For example \(|\Delta B^{(j)}|\leq W_*/4\) through order three, so the norm is at most \([t^2/2+t+(1+t)/2]W_*/4\) for \(t\leq l<1\). For the other integral, use \[\Delta(-2B_1+tB_1^2) =\bigl(-2+t(2\widetilde B_1+\Delta B_1)\bigr)\Delta B_1.\] The second-order norm of its integral is at most twice the supremum of the first-order norm of the integrand, because \(l<1\). Submultiplicativity and \(|t|_{1,t}<2\) yield \[|\Delta S_2|_{2,t} \leq(4+8\cdot6000+4W_*)W_*<60000W_*.\]

The identities and elementary bounds \[H_k'=B_1-kH_k/t,\quad |H_k|\leq\frac{t\|B_1\|_{0,[0,t]}}{k+1},\quad |H_k'|\leq2\|B_1\|_{0,[0,t]},\quad |H_k''|\leq|B_1'|+\frac{2k}{t}\|B_1\|_{0,[0,t]}\] imply, for \(b\leq t\leq l\), \[ |H_k|_{2,t}\leq6(k+1)(6000+W_*),\qquad |\Delta H_k|_{2,t}\leq6(k+1)W_*. \tag{62}\] Indeed \(l/(k+1)+5/2+k/b<6(k+1)\). The bound for \(\Delta H_k\) follows by the same linear argument applied to \(\Delta B_1\).

The radial source and energy functions are \[\begin{align*} f_0(t)={}&-2\pi\left(S_1(t)+t^2\sum_i\frac{D_i(t)}{k_i+2}\right)\\ &-\frac12\left(\log t+S_2(t) +\sum_i(-2D_i+2k_iD_iH_{k_i})(t) +\sum_{i,j}\frac{k_ik_j}{k_i+k_j}D_i(t)D_j(t)\right),\\ e(t)={}&\frac12\sum_{i,j}\frac{Q_i(t)Q_j(t)}{k_i+k_j}. \end{align*}\] The preceding bounds give \[\begin{align*} |f_0-\widetilde f_0|_{2,t} &\leq2\pi(3W_*+16\mathfrak a)\\ &\quad+\frac12\bigl(60000W_*+8\mathfrak a +12(2\cdot10^6)W_* +48(6000+W_*)\mathfrak a\\ &\hspace{44mm}+8(2\cdot10^6)\mathfrak a +16\mathfrak a^2\bigr) <6\cdot10^7\mathfrak a. \tag{63}\end{align*}\] For completeness, \(|t^2|_{2,t}=(t+1)^2<4\) gives the \(16\mathfrak a\) source term. The sum of \(-2D_i\) contributes \(8\mathfrak a\). Split the mixed term as \[\sum_i k_i(D_iH_{k_i}-\widetilde D_i\widetilde H_{k_i}) =\sum_{i\leq38}k_i\widetilde D_i\Delta H_{k_i} +\sum_{i\geq39}k_iD_iH_{k_i}.\] The finite and omitted weighted norms give \(12NW_*\) and \(48(6000+W_*)\mathfrak a\) after including its factor \(2\). Bounding \(k_ik_j/(k_i+k_j)\) by \(k_ik_j\), the quadratic difference is at most \(8N\mathfrak a+16\mathfrak a^2\). These are precisely the terms in (63). Finally \(W_*\leq4\mathfrak a\), \(\pi<22/7\), and \(\mathfrak a<2\cdot10^{-21}\) make its total coefficient less than \(56,264,181<6\cdot10^7\).

Since \(k_i\geq6\) and \(t<1\), direct differentiation of \(Q_i\) shows \[|Q_i|_{2,t}\leq2(k_i+1)|D_i|_{3,t}.\] For example, its coefficients on \(|D_i|,|D_i'|,|D_i''|,|D_i'''|\) are bounded respectively by \(k_i\), \(t+k_i-1\), \(t+(k_i-2)/2\), and \(t/2\). The finite and omitted sums of these norms are bounded by \(2N\) and \(8\mathfrak a\). Dropping the denominator \(k_i+k_j\geq1\) in the energy difference therefore gives \[ |e-\widetilde e|_{2,t} \leq16(2\cdot10^6)\mathfrak a+32\mathfrak a^2. \tag{64}\] The norm contains half the absolute second derivative. Accounting for this factor, (63) and (64) imply absolute truncation errors less than \(10^{-11}\) for every member of \[ f_0,\quad e,\quad f_0',\quad e',\quad f_0'',\quad e'' \qquad (b\leq t\leq l). \tag{65}\] In fact the respective second derivative bounds are below \(2.4\cdot10^{-13}\) and \(1.3\cdot10^{-13}\).

The constant region

For \(r\geq l\), all \(D_i\) vanish and \(B\) is constant. The omitted constant satisfies \[|\Delta B(l)|\leq \sum_{k=234,240,\ldots}\frac7k(R_0/a)^k \leq\frac{(7/234)(289/500)^{234}}{1-(289/500)^6} <6.075\cdot10^{-58}<10^{-11}.\] Since \(B_1(t)=1/t\) there, the derivatives of \(\log t\) and \(S_2(t)\) cancel. The full infinite-mode functions consequently obey the exact formula \[ e(t)=0,\qquad f_0(t)=f_0(l)-\pi B(l)(t^2-l^2),\qquad t\geq l. \tag{66}\] Thus outside \([b,l]\) one propagates the enclosures of \(f_0(l)\) and \(B(l)\) through this exact formula; no uniform absolute error for arbitrarily large \(t\) is asserted.

Outward arithmetic and derivative intervals

The calculation follows the interval inclusion principle: every operation encloses its exact result, and analytic remainders are added explicitly. Revol and Rouillier (Revol and Rouillier 2005) discuss this principle and the dependence of enclosure width on expressions and subdivision. We specify the integer operations here, so the proof’s arithmetic assumptions and rounding rules are explicit.

Put \(q=2^{144}\). A scalar interval is represented by integers \(u\le v\), meaning \([u/q,v/q]\). Rational inputs are rounded down at the lower endpoint and up at the upper endpoint. Addition and negation are exact at this denominator. Multiplication takes the least and greatest of the four endpoint products, then divides by \(q\) with outward rounding. If \(0<u\le v\), inversion gives \[\left[\frac{\lfloor q^2/v\rfloor}{q}, \frac{\lceil q^2/u\rceil}{q}\right];\] negative intervals are treated by negation, and division by an interval containing zero is forbidden. Integer powers use repeated multiplication. For a positive interval, integer square roots give the enclosure \[\left[\frac{\lfloor\sqrt{uq}\rfloor}{q}, \frac{\lfloor\sqrt{vq}\rfloor+1}{q}\right].\] The extra unit in the upper endpoint is harmless when the upper square root is an integer.

For logarithms, first extract an integer power of two and reduce the argument to \([1/4,1]\). The arguments used below are positive table intervals or enclosures of fixed positive constants; their widths allow this reduction. With \(z=(v-1)/(v+1)\), evaluate \[ \log v=2\sum_{j=0}^{199}\frac{z^{2j+1}}{2j+1}+\mathcal R, \qquad |\mathcal R| \le\frac{2(3/5)^{401}}{401(1-(3/5)^2)}<q^{-1}. \tag{67}\] The same series at \(z=1/3\) encloses \(\log2\), used to restore the extracted power. All summation errors are already enclosed by outward arithmetic; the displayed remainder is added separately. Similarly, use \[ \pi=16\arctan(1/5)-4\arctan(1/239),\qquad \arctan(1/v)=\sum_{j=0}^{149} \frac{(-1)^jv^{-2j-1}}{2j+1}+\mathcal R_v. \tag{68}\] The alternating-series bound gives \(|\mathcal R_v|\le v^{-301}/301<q^{-1}\) for both values of \(v\). Adding \([-2q^{-1},2q^{-1}]\) to each arctangent enclosure before the linear combination is sufficient. These procedures also enclose \(p,R_0,x_0,l\) using their exact defining formulas.

To bound derivatives we use interval jets \[(F_0,F_1,F_2,F_3),\qquad F_j\supseteq F^{(j)}(t)/j!.\] Addition, multiplication, inversion, logarithms, and square roots are performed as formal power series through degree three. For example, the product coefficient is \(\sum_{r=0}^jF_rG_{j-r}\), and the inverse coefficients follow recursively from the product being one. A variable on a radial interval \(T\) is represented by \((T,1,0,0)\); the resulting coefficients enclose its derivatives at every point of \(T\) on the chosen spline piece. At a single point, the leading interval encloses that point and the same derivative entries are used. After differentiation only the coefficients supported by the available derivative orders are used. The calculation below never requires a fourth derivative of a spline: third derivatives of \(B,D_i\) suffice for the second derivatives of \(f_0,e\) and for the radial quadrature errors.

For a differential equation \(Y'=F(t,Y)\) with an enclosed value at the evaluation point, its jet coefficients are reconstructed in increasing order by \[ Y_{j+1}=\frac{[F(t,Y)]_j}{j+1},\qquad 0\le j\le2. \tag{69}\] Thus the enclosed value and the known differential equation supply the derivatives without numerical differentiation.

Spline polynomials are evaluated by interval Horner arithmetic. One may first translate a polynomial \(P(t)=\sum_{\ell=0}^Na_\ell t^\ell\) to the midpoint \(m\) of its leading interval: \[ P(m+h)=\sum_{j=0}^N \left(\sum_{\ell=j}^N\binom{\ell}{j}a_\ell m^{\ell-j}\right)h^j. \tag{70}\] This is an exact polynomial identity, also for jets. It reduces interval widths without changing the spline or introducing an approximation; direct Horner evaluation remains a valid enclosure.

Coefficients and the radial table

Retain the auxiliary mode \(0\) and modes \(1\le m\le38\). The two modes with exponent six remain separate in sums involving \(D_i\). Enclose each \(c_m\) by summing over \(|n|,|r|\le D\), where \(D=256\) for \(m=1\) and \(D=48\) otherwise, and then adding the symmetric error in (56), multiplied by \(a^{-k_m}/k_m\). The summands can be evaluated using only integers before that last factor: if \(k=k_m\), \(A=2n+r\), and \(N=n^2+nr+r^2\), then \[ \operatorname{Re}(n+r/2+i\sqrt3r/2)^{-k} =\frac{\displaystyle\sum_{j=0}^{k/2}(-1)^j \binom{k}{2j}A^{k-2j}(3r^2)^j}{2^kN^k}. \tag{71}\] The pair \((0,0)\) is omitted. The resulting enclosure also tests \(|c_1|<1\).

The radial breakpoints are \[ b,\quad p-.06,\quad p,\quad R_0,\quad R_0+.001, \quad R_0+.075,\quad l=R_0+.076. \tag{72}\] Each intervening interval is subdivided by affine interpolation between its exact endpoints. For the first interval use \(n=\lfloor8000(u-v)\rfloor+1\) subintervals, and for all later intervals use \(n=\lfloor24000(u-v)\rfloor+1\), where \(v,u\) are its endpoints. In execution an outward upper bound on the product is used to choose \(n\), which can only add subdivisions. All true node positions are enclosed by the corresponding interval arithmetic. The length of every actual subinterval is bounded above in the calculation; these upper bounds are used in all error estimates.

For the fixed table, write \([u_i/q,v_i/q]\) for the valid outward enclosure of its \(i\)th node, \(0\le i\le m\). The interval constructor enforces \(u_i\le v_i\), and the node arithmetic preserves this invariant. After appending the last node and before propagating any integral state, the program requires \[ v_i<u_{i+1}\qquad(0\le i<m); \tag{73}\] it throws an exception if a comparison fails. Thus, on a completed table construction, \(u_i\le v_i<u_{i+1}\le v_{i+1}\). Both endpoint arrays are strictly increasing, and the stored step interval is \[H_i=[(u_{i+1}-v_i)/q,(v_{i+1}-u_i)/q]\] with positive lower endpoint. This also justifies every division by \(H_i\) below. These facts concern the initialized table and its unchanged aligned node, point, and cell arrays in the fixed calculation.

On \(0<t\le p-.06\), the truncated functions have the monomial form \[\widetilde B(t)=-\log t+\pi t^2/2 +\sum_{m=1}^{38}a_mt^{k_m},\qquad \widetilde B_1(t)=\pi t+\sum_{m=1}^{38}k_ma_mt^{k_m-1}, \qquad a_m=(-1)^mc_m.\] For compact formulas introduce the finite list \[(e_0,v_0)=(2,\pi),\qquad (e_m,v_m)=(k_m,k_ma_m)\quad(1\le m\le38).\] This list indexes monomials of \(\widetilde B_1\), rather than the sector modes; its exponent \(e_0=2\) is unrelated to \(k_0=6\). The starting integrals, including their derivatives, are obtained from the exact identities \[ \begin{split} \widetilde S_1(t)&=\frac{t^2}{2}\left(\frac12-\log t\right) +\frac\pi8t^4+ \sum_{m=1}^{38}\frac{a_mt^{k_m+2}}{k_m+2},\\ \widetilde S_2(t)&=-2\sum_{i=0}^{38}\frac{v_it^{e_i}}{e_i} +\sum_{i,j=0}^{38}\frac{v_iv_jt^{e_i+e_j}}{e_i+e_j},\\ \widetilde H_k(t)&=\sum_{i=0}^{38}\frac{v_it^{e_i}}{e_i+k}. \end{split} \tag{74}\] In particular, the table starts propagation at the exact node \(p-.06\) with these values; preceding nodes and cells use (74) directly.

There is a useful exact simplification for the sector functions on this first interval. Since \(b_m(t)=0\), \(D_m(t)=-a_mt^{k_m}\) and \(D_0(t)=0\), substitution in (46) cancels every nonconstant mode and gives \(H(r,\varphi)=\gamma(r)\) for \(r\le t\). Equivalently, every \(Q_i\) vanishes. Direct radial integration gives \[ f_0(t)=(\pi t^2-\tfrac12)\log t-\frac{3\pi^2t^4}{8}, \qquad e(t)=0,\qquad 0<t\le p-.06. \tag{75}\] Indeed, \[\int_0^t\gamma(r)r\,\mathrm dr =-\frac{t^2}{2}\log t+\frac{t^2}{4}+\frac{\pi t^4}{8}, \qquad \int_0^t{}^{\rm ren}\gamma'(r)^2r\,\mathrm dr =\log t-\pi t^2+\frac{\pi^2t^4}{4}.\] Combining these with the source and energy factors in the sector formula proves (75). The identity holds for the full infinite series as well as for every finite truncation. It may therefore replace the mode sums when evaluating \(f_0,e\) and their first two derivatives on the first interval, producing tighter enclosures. The full finite-mode values (74) at \(p-.06\) are still required to start the continuation of \(S_1,S_2,H_k\).

For a subsequent radial step \([v,u]\), propagate the integrals by \[ \begin{split} \widetilde S_1(u)&=\widetilde S_1(v)+ \int_v^u r\widetilde B(r)\,\mathrm dr,\\ \widetilde S_2(u)&=\widetilde S_2(v)+ \int_v^u\bigl(-2\widetilde B_1(r) +r\widetilde B_1(r)^2\bigr)\,\mathrm dr,\\ \widetilde H_k(u)&=(v/u)^k\widetilde H_k(v) +\int_v^u(r/u)^k\widetilde B_1(r)\,\mathrm dr. \end{split} \tag{76}\] For each integrand \(Y\), the midpoint rule has the rigorous enclosure \[ \int_v^uY(r)\,\mathrm dr \in (u-v)Y((v+u)/2) +\left[-\frac{(u-v)^3}{24}\sup_{[v,u]}|Y''|, \frac{(u-v)^3}{24}\sup_{[v,u]}|Y''|\right]. \tag{77}\] Midpoint values use point jets; second-derivative bounds use a jet whose leading interval contains \([v,u]\). The formula applies to weak second derivatives bounded almost everywhere as well.

For values inside a step, use the integral equations rather than the midpoint approximation. If \(T=[v,u]\) and \(h=u-v\), then \[ \begin{split} \widetilde S_j(T)&\subseteq \widetilde S_j(v) +[0,h]\,\operatorname{range}_T Y_j, \qquad j=1,2,\\ \widetilde H_k(T)&\subseteq (v/T)^k\widetilde H_k(v) +[-h\|\widetilde B_1\|_{0,T}, h\|\widetilde B_1\|_{0,T}], \end{split} \tag{78}\] where \(Y_1=r\widetilde B\) and \(Y_2=-2\widetilde B_1+r\widetilde B_1^2\). For the last inclusion, \((r/t)^k\le1\) whenever \(v\le r\le t\le u\). Recover the positive jet orders from \[ \widetilde S_1'=t\widetilde B,\qquad \widetilde S_2'=-2\widetilde B_1+t\widetilde B_1^2,\qquad \widetilde H_k'=\widetilde B_1-k\widetilde H_k/t. \tag{79}\] Substitution into (48) produces both node enclosures and cell enclosures of the ordered six-tuple \[ (\widetilde f_0,\widetilde f_0',\widetilde f_0'', \widetilde e,\widetilde e',\widetilde e''). \tag{80}\] No numerical differentiation of table values is used.

Separately, subdivide each interval in (72) into steps of length at most \(.00022\). Direct interval jets for \(\widetilde B_1\) and \(\widetilde D_i\) test the two finite-core inequalities (60). Together with \(|c_1|<1\), these are the finite premises of the truncation analysis; their verification does not use any infinite-mode error bound. Only after these tests pass is (65) invoked.

The spline data are \(C^2\) and have bounded one-sided third derivatives. Accordingly \(f_0,e\) are \(C^1\) with locally bounded weak second derivatives on \((b,\infty)\). A radial cell uses its own spline formula at its endpoints; either one-sided second derivative is valid for its interior. When a query meets a knot, the hull includes both adjacent cells. Thus no jump of a second derivative is mistaken for global classical differentiability.

Queries for the infinite-mode functions

Let \(\epsilon=[-10^{-11},10^{-11}]\), and let \(T\) be a valid outward enclosure of the queried true radii. In the fixed production calls every such radius is at least \(b\), although \(T\) may extend below \(b\). For derivatives, take the componentwise hull of (80) over every radial cell meeting \(T\) and add \(\epsilon\) to each component. By (65), this encloses the full infinite-mode tuple \((f_0,f_0',f_0'',e,e',e'')\) at every queried true radius in \(T\cap[b,l]\). The exact continuation below covers queried true radii at least \(l\). Here is the cell-selection argument for the implementation. In integer endpoint units write the valid query as \([a,b']\), let \(j\) be the first index with \(v_j\ge a\), and let \(k\) be the first index with \(u_k>b'\), using \(m+1\) if either index is absent. The closed envelope of cell \(i\) is \([u_i,v_{i+1}]\). It meets \([a,b']\) exactly when \(v_{i+1}\ge a\) and \(u_i\le b'\). By the ordering proved above, these conditions are equivalent to \(i\ge j-1\) and \(i<k\). Therefore the binary searches on the upper and lower endpoint arrays select exactly \[\mathrm{im}=\max(0,j-1)\le i<\min(m,k)=\mathrm{ix}.\] The non-strict intersection inequalities include both cells at a shared endpoint. Every true cell meeting the query lies in its envelope and is therefore selected; any extra cells only widen the hull. This is a statement about the fixed initialized table, not a general-purpose lookup interface.

For a query requiring values alone, improve the enclosure by linear interpolation of the endpoint values. On a radial cell \([v,u]\) put \[X=[0,1]\cap\frac{T-v}{u-v},\qquad L_Y(T)=Y(v)+X\bigl(Y(u)-Y(v)\bigr), \quad Y=\widetilde f_0,\widetilde e.\] An empty intersection is discarded. The interpolation error is \[ |Y(t)-L_Y(t)|\le\frac{(u-v)^2}{8} \sup_{[v,u]}|Y''|. \tag{81}\] The cell enclosure supplies the derivative bound. Add this symmetric error and \(\epsilon\), then take the hull over all relevant cells. Shared coefficient uncertainties need not be independent: interval operations enclose all their possible values and hence also their correlated exact values. In particular, the positive stored interval \(H_i\) contains the true step \(u-v\), so outward evaluation of the quotient defining \(X\) contains the normalized coordinate of every true queried point in that cell. Intersecting with \([0,1]\) cannot discard such a point.

Queries that meet \([l,\infty)\) also include the exact continuation (66). Enclose its two constants by the last node value plus \(\epsilon\) and by \[ B(l)\in-\log R_0+\pi R_0^2/2+ \sum_{m=1}^{38}V_m(R_0)+\epsilon. \tag{82}\] Thus the continuation contributes the tuple \[ \bigl(f_0(l)-\pi B(l)(T^2-l^2),\, -2\pi B(l)T,\,-2\pi B(l),\,0,0,0\bigr). \tag{83}\] Using the entire query interval in this expression is conservative even if it extends below \(l\), since the relevant exterior subset is contained in it. This rule propagates the constant uncertainties through the exact formula and asserts no radius-independent error on the unbounded tail.

Enclosing the calibrated constants

We evaluate the exact integrals defining \(\lambda\) in (53). Put \(s_0=.01\) and \(t(s)=\sqrt{p^2+s^2}\). Taylor’s formula at \(p\) removes the apparent division by zero in the first integral. Specifically, \[ \begin{split} \int_0^{s_0}p\frac{f_0(t(s))-f_0(p)}{s^2}\,\mathrm ds \in{}& f_0'(p)\left(\frac{s_0}{2} -\frac{s_0^3}{24p^2} +\left[0,\frac{s_0^5}{80p^4}\right]\right)\\ &+\frac{s_0^3}{6} \left[\frac{p}{(\sqrt{p^2+s_0^2}+p)^2},\frac1{4p}\right] \operatorname{range}_{[p,t(s_0)]} f_0''. \end{split} \tag{84}\] Indeed, \(t-p=s^2/(t+p)\), while \[\frac{p}{t+p}=\frac12-\frac{s^2}{8p^2} +\left[0,\frac{s^4}{16p^4}\right].\] The quadratic Taylor remainder of \(f_0\) supplies the second line of (84). The derivative range comes from the radial cells with the infinite-mode error added.

The implementation obtains this range by a direct cell scan over \([p,t(s_0)]\), and obtains the derivative bounds below by a second direct scan over \([p,R_0]\). Both use the same closed-cell selection rule, with outward endpoint bounds, but neither adds the exterior formula. This is sufficient here: \(p^4=1/12\) gives \(p>1/2>b\), \(x_0=R_0/2>s_0\), and \[t(x_0)=\sqrt{p^2+(R_0/2)^2}=R_0<l.\] Hence both true scan domains are contained in the finite table \([b,l]\). A nonempty scan alone would not justify a query extending beyond that table.

On \([s_0,x_0]\) integrate \(p(f_0(t(s))-f_0(p))/s^2\) by a midpoint sum with \(45000\) equal steps. Integrate \(e(t(s))/p\) on \([0,x_0]\) by a separate midpoint sum with \(45000\) equal steps. Obtain integer upper bounds \[C\ge\sup_{[p,R_0]}|f_0'|,\quad C'\ge\sup_{[p,R_0]}|f_0''|,\quad E\ge\sup_{[p,R_0]}|e'|,\quad E'\ge\sup_{[p,R_0]}|e''|\] from the cell enclosures plus \(\epsilon\). The two integrands have second derivatives bounded respectively by \[ \frac{C'}p+\frac{8C}{s_0^2},\qquad \frac{E'}p+\frac E{p^2}. \tag{85}\] For the first bound, set \(V=f_0(t(s))-f_0(p)\); then \[|V|\le Cs^2/(2p),\qquad V''=f_0''(t)(s/t)^2+f_0'(t)p^2/t^3,\] and differentiate \(pV/s^2\) twice. The resulting terms are \(p(V''/s^2-4V'/s^3+6V/s^4)\), bounded by the first expression in (85). The second follows by twice differentiating \(e(t(s))/p\). For an interval of length \(L\) and equal step length \(h\), add the symmetric error \(Lh^2\sup|Y''|/24\) to each complete midpoint sum. Together with (84), the endpoint terms and the factors in (53), this encloses \(\lambda\). The same formula then encloses \(\mu\).

To enclose \(K=2\int_0^{x_0}G(p,s)\,\mathrm ds\), choose an integer \(n\) and equal steps of length \(h=x_0/n\le.0001\), and use their midpoint sum. Independently enclose \(G_{ss}(p,s)\) on a partition with step length at most \(.0002\). If \(h\) is the midpoint step length and \(L_K\) bounds this second derivative, the enclosure is \[ K\in2\left(\sum_jhG(p,s_j) +[-x_0h^2L_K/24,x_0h^2L_K/24]\right). \tag{86}\] The factor two multiplies both the sum and its error.

The first sign tests are \[ \begin{gathered} -6.731<\lambda<-6.729,\qquad K<0,\\ f_0(l)-\lambda l^2/2-\mu>0,\qquad -\pi B(l)-\lambda/2>0. \end{gathered} \tag{87}\] Each strict inequality is accepted only when the corresponding interval endpoint proves it. The bound on \(\lambda\) also serves as a check on the normalization. The last two signs and (66) give \[ f(t)=f(l)+\bigl(-\pi B(l)-\lambda/2\bigr)(t^2-l^2)>0, \qquad G(d,s)>0\quad\hbox{if }t=\sqrt{d^2+s^2}\ge l. \tag{88}\] Here \(e(t)=0\) and \(d/t^2>0\).

The compact rectangle and its derivative bounds

Because \(l<.697\) and \(\sqrt{b^2+.65^2}>l\), the signs (88) reduce the remaining work to \[ \mathcal D=[b,.697]\times[0,.65]. \tag{89}\] For \(d\ge.697\), \(G(d,s)>0\) for every \(s\ge0\) and \(g_d(x)\ge0\). For \(d\ge b\) and \(x\ge.65\), the primitive is increasing in \(x\). The negative part of \(G\) vanishes for \(s\ge.65\).

We obtain derivatives of \(G\) by differentiating its defining formula, using the full radial query enclosures. The following identities specify all needed factors. Put \[t=\sqrt{d^2+s^2},\quad A=d/t^2,\quad Z=-s^2/(dt^2)=A-1/d,\quad r_d=d/t,\quad r_s=s/t,\] and write \(f'=f_0'-\lambda t\), \(f''=f_0''-\lambda\). Then \[ \begin{aligned} A_d&=(1-2d^2/t^2)/t^2,&Z_d&=A_d+1/d^2,\\ A_s&=-2ds/t^4,&Z_s&=A_s,\\ A_{dd}&=2d(d^2-3s^2)/t^6,&Z_{dd}&=A_{dd}-2/d^3,\\ A_{ss}&=-2d/t^4+8ds^2/t^6,&Z_{ss}&=A_{ss}. \end{aligned} \tag{90}\] For \(v=d\) or \(v=s\) let \(r_v\) have the meaning above. The derivative formulas are \[ \begin{split} G_v={}&A_vf+Af'r_v+Z_ve+Ze'r_v,\\ G_{vv}={}&A_{vv}f+2A_vf'r_v +A\bigl(f''r_v^2+f'(1-r_v^2)/t\bigr)\\ &+Z_{vv}e+2Z_ve'r_v +Z\bigl(e''r_v^2+e'(1-r_v^2)/t\bigr). \end{split} \tag{91}\] These formulas hold on each smooth piece and almost everywhere across the radial knots. No mixed second derivative of \(G\) is needed.

Partition both sides of (89) into intervals of length at most \(.0008\). Interval substitution in (91) on every closed rectangle supplies global upper bounds \[ L_1\ge\|G_d\|_{\infty,\mathcal D},\quad L_2\ge\|G_s\|_{\infty,\mathcal D},\quad L_3\ge\|G_{dd}\|_{\infty,\mathcal D},\quad L_4\ge\|G_{ss}\|_{\infty,\mathcal D}. \tag{92}\] The value-only interpolation is not used to obtain these derivative bounds.

The primitive away from its stationary point

Let \[ \mathcal D_0=[.52528,.54930]\times[.29420,.32621]. \tag{93}\] Its interior contains \((p,x_0)\). To verify the primitive bound on the complement, partition the \(d\) interval at the two vertical sides of \(\mathcal D_0\), and partition the \(s\) interval at its two horizontal sides. Subdivide these pieces into steps with maximal lengths \[D_h\le.00019,\qquad S_h\le.00018.\] Every cell therefore lies either in \(\mathcal D_0\) or outside its interior. For each \(d\) node, form successive midpoint sums from \(0\) to every \(x\) node, using the actual step lengths: \[Q(d,x)=\sum_{[s_j,s_{j+1}]\subseteq[0,x]} (s_{j+1}-s_j) G\left(d,\frac{s_j+s_{j+1}}2\right).\] The sum at \(x=0\) is zero and is tested as well. The value queries in Section 7.3 may be used here. At every node outside the interior of \(\mathcal D_0\), enclose \(2Q(d,x)-K\) and let \(m_*\) be a lower bound for all these enclosures.

The quadrature error at grid nodes is bounded by \(E_0\), and the additional interpolation error between nodes by \(E_I\), where \[ E_0=\frac{2(.65)S_h^2L_4}{24},\qquad E_I=\frac{2\bigl(D_h^2(.65)L_3+S_h^2L_2\bigr)}8, \tag{94}\] Indeed, the bilinear interpolation estimate follows from \[\partial_d^2(2g_d(x)-K)=2\int_0^xG_{dd}(d,s)\,\mathrm ds, \qquad \partial_x^2(2g_d(x)-K)=2G_s(d,x).\] Apply the one-dimensional interpolation error in each variable successively; the resulting bounds add and require no mixed derivative. Thus the finite assertion \[ m_*-E_0-E_I>0 \tag{95}\] proves \(2g_d(x)>K\) throughout \(\mathcal D\setminus\operatorname{int}\mathcal D_0\). All the derivative estimates are weak derivative estimates, so the same conclusion holds when a quadrature segment or rectangle crosses a radial knot.

The negative part on the whole rectangle

At every \(d\) node, including those between the vertical sides of \(\mathcal D_0\), use the same \(s\) partition to enclose \[N(d)=\sum_j(s_{j+1}-s_j) \max\left\{0,-G\left(d,\frac{s_j+s_{j+1}}2\right)\right\}.\] Let \(N_*\) be an upper bound for all these sums. The function \(y\mapsto\max\{0,-y\}\) is \(1\)-Lipschitz. On a subinterval of length \(h\), its midpoint error is therefore at most \(L_2h^2/4\): integrate \(L_2|s-s_{\rm mid}|\). Summing the errors gives at most \(.65S_hL_2/4\). Every \(d\) lies within \(D_h/2\) of a node; changing \(d\) to that node changes the integral by at most \(.65D_hL_1/2\). Consequently the assertion \[ N_*+K+.65\left(\frac{S_hL_2}{4}+\frac{D_hL_1}{2}\right)<0 \tag{96}\] proves \[\int_0^{.65}\max\{0,-G(d,s)\}\,\mathrm ds<-K \qquad(b\le d\le.697).\] Together with (88), this gives the entire negative-part inequality in (55).

Convexity in the stationary rectangle

It remains to prove the primitive inequality on \(\mathcal D_0\). Partition its \(d\) interval into steps of length at most \(.00004\). Partition \([0,.32621]\) into steps of length at most \(.00008\), including \(.29420\) as a node. The accumulation in \(s\) must begin at zero even though the Hessian is tested only above \(.29420\).

For a \(d\) interval \(D\) and an \(s\) interval \(S_j\), use the full derivative queries to enclose \[F_1(D,S_j)\supseteq G_{dd}(D,S_j),\quad F_2(D,S_j)\supseteq G_d(D,S_j),\quad F_3(D,S_j)\supseteq G_s(D,S_j).\] For a cell \(D\times S_j\) meeting \(\mathcal D_0\), enclosures of the three Hessian entries of \((d,x)\mapsto g_d(x)\) are \[ \begin{split} A_j&=\sum_{i<j}|S_i|F_1(D,S_i) +[0,|S_j|]F_1(D,S_j),\\ B_j&=F_2(D,S_j),\qquad C_j=F_3(D,S_j). \end{split} \tag{97}\] The first line is direct interval integration, not midpoint quadrature: it bounds both the completed subintervals and the partial final subinterval. Test, on every such cell, that \[ \inf A_j>0,\qquad \inf C_j>0,\qquad \inf(A_jC_j-B_j^2)>0. \tag{98}\] Interval multiplication encloses the determinant even when its factors are correlated. These inequalities imply positive definiteness, and in particular positive semidefiniteness, of the Hessian almost everywhere on \(\mathcal D_0\).

The function \((d,x)\mapsto g_d(x)\) belongs to \(W^{2,\infty}\) on a neighbourhood of this rectangle. To justify the passage from Hessian bounds to convexity, mollify on each smaller interior rectangle: the mollified Hessian is positive semidefinite, so the mollified functions are convex. Local uniform convergence proves convexity in the interior, and continuity includes the boundary. By Lemma 16, \((p,x_0)\) is an interior stationary point. Hence \[ 2g_d(x)\ge2g_p(x_0)=K\qquad((d,x)\in\mathcal D_0). \tag{99}\] This step uses the exact stationary point and does not infer convexity or stationarity from sampled function values.

Finite assertions and the scalar conclusion

For clarity, the finite procedure has the following fixed parameters and obligations.

The spherical logarithmic-energy constant

An independent recurrence for the lattice coefficients

Fix integers \((n,q)\ne(0,0)\) and write \[A=2n+q,\qquad N=n^2+nq+q^2,\qquad \zeta=A+i\sqrt3q,\qquad w=\frac{\zeta}{2}.\] Thus \(w=n+q/2+i\sqrt3q/2\), \(N=|w|^2>0\), and \(|\zeta|^2=4N\). A point of the covolume-one triangular lattice is \(aw\), where \(a=\sqrt{2/\sqrt3}\). For \(m\geq0\), put \[r_m(n,q)=\operatorname{Re}(\zeta^{6m}).\] The two numbers \(u=\zeta^6\) and \(v=\overline{\zeta}^{6}\) have \[u+v=2r_1,\qquad uv=(4N)^6,\] and hence satisfy \(X^2-2r_1X+(4N)^6=0\). Multiplying this equation by \(u^{m-1}\) and \(v^{m-1}\), respectively, and adding gives the integer recurrence \[ \begin{split} r_0&=1,\\ r_1&=A^6-45A^4q^2+135A^2q^4-27q^6,\\ r_{m+1}&=2r_1r_m-(4N)^6r_{m-1}\qquad(m\geq1). \end{split} \tag{101}\] Here all the quantities on the right are integers. This proves the recurrence for every individual lattice point, without a numerical approximation to a complex power.

The binomial theorem identifies its numerator exactly. Since only even powers of \(i\sqrt3q\) have nonzero real part, \[ r_m(n,q)=\sum_{j=0}^{3m}(-1)^j\binom{6m}{2j} A^{6m-2j}(3q^2)^j. \tag{102}\] This is the numerator in (71), after writing \(q\) for its second integer coordinate. Indeed, \(w^{-6m}=\overline{w}^{6m}/N^{6m}\), and therefore \[ \operatorname{Re}(w^{-6m}) =\frac{r_m(n,q)}{2^{6m}N^{6m}} =\frac{r_m(n,q)}{(64N^6)^m}. \tag{103}\] Thus the recurrence and the binomial expression agree for every point and every mode, not only for a selected finite list of comparisons.

Alternative exact bounds for the omitted modes

For the allowed indices \(k=234,240,\ldots\), split the weighted summand \((k+1)^4w_k\) of (58) into four positive rationals: \[\begin{align*} t_1(k)&=10^6\,30(k+1)^4k^2(13/20)^k,\\ t_2(k)&=10^6(3\cdot10^6)(k+1)^4k^3(1/2)^k,\\ t_3(k)&=10^6\,5600(k+1)^4k^5(3/5)^k,\\ t_4(k)&=10^6(7\cdot10^{11})(k+1)^4k^2(289/500)^k. \end{align*}\] Then \[\mathfrak a=\sum_{j=0}^{\infty}\sum_{i=1}^4t_i(234+6j).\] This is the same \(\mathfrak a\) as in (61), distinct from the lattice scale \(a\). If \[(p_i,\beta_i)=(2,13/20),(3,1/2),(5,3/5),(2,289/500),\] in the order \(i=1,2,3,4\), the exact consecutive ratios are \[ \rho_i(k)=\frac{t_i(k+6)}{t_i(k)} =\left(\frac{k+7}{k+1}\right)^4 \left(\frac{k+6}{k}\right)^{p_i}\beta_i^6. \tag{106}\] Both rational factors decrease for \(k>0\), so each \(\rho_i(k)\) decreases. The four initial ratios, with no rounding, are \[\begin{array}{c|c} i&\rho_i(234)\\ \hline 1&(241/235)^4(40/39)^2(13/20)^6\\ 2&(241/235)^4(40/39)^3(1/2)^6\\ 3&(241/235)^4(40/39)^5(3/5)^6\\ 4&(241/235)^4(40/39)^2(289/500)^6. \end{array}\] Here is an exact simultaneous check that all four are less than \(1/10\). The integer comparisons \[\begin{align*} 9\cdot241^4=30\,360\,623\,049 &<30\,498\,006\,250=10\cdot235^4,\\ 7\cdot40^5=716\,800\,000 &<721\,793\,592=8\cdot39^5,\\ 800\cdot13^6=3\,861\,447\,200 &<4\,032\,000\,000=63\cdot20^6 \end{align*}\] give, since \(p_i\leq5\) and \(0<\beta_i\leq13/20\), \[ \rho_i(234) \leq(241/235)^4(40/39)^5(13/20)^6 <\frac{10}{9}\frac87(13/20)^6<\frac1{10} \quad(1\leq i\leq4). \tag{107}\]