Dmytro Taranovsky
Last modified: October 18, 2021

Asymptotic Optimal Sphericity

Abstract: We prove that maximum sphericity for polyhedra with $n$ faces is at least $1 - \frac{5 \sqrt{3} π}{27n} - O(n^{-3/2})$, and conjecture it to be $1 - \frac{5 \sqrt 3 π}{27n} - C_1 n^{-3/2} + (C_2±o(1)) n^{-5/3} \log^{2/3} n$. We discuss in detail the conjectured fractal-like behavior of nonhexagonal faces for maximum sphericity polyhedra, and similar behavior and bounds for other optimization problems on the sphere such as the Thomson problem.

Introduction

Sphericity is defined as $π^{1/3}(6V)^{2/3}/S$ given volume $V$ and surface area $S$. The simplicity of the formulation (maximum volume for $n$-faced polyhedra with given surface area) conceals the rich structure of such polyhedra as $n→∞$. Computing maximum sphericity for small $n$ is hard, even though by the decidability of the real closed field, it is an algebraic multiple $π^{1/3}$. However, asymptotics can bring a sense of clarity, and here apparently, a convergent large scale fractal structure of the polyhedra. The four asymptotic terms are:
  1  -  maximum sphericity is 1 (the sphere).
  $\frac{-5 \sqrt{3} π}{27n}$  -  the hexagonal limit, corresponding to same-sized regular hexagonal faces, which we (re)derive here.
  $-C_1 n^{-3/2}$  -  a fractal-like structure due to non-existence of a parallel incompressible vector field on a sphere (the main subject here).
  $(C_2±o(1)) n^{-5/3} \log^{2/3} n$  -  savings from using area deformations rather than defect lines at small scales (quantization costs are smaller).

We take the readers through the rich picture, both for sphericity and for other optimization problems on the sphere. The hexagonal limit was known previously, but the $O(n^{-3/2})$ term and the main content of the paper is new. $o,O,Θ,Ω,ω$ are from the Big O Notation.

General properties and small $n$

Maximum sphericity convex $n$-polyhedron (i.e. with $n$ faces) is known to have an inscribed sphere tangent to all faces at their centroids. Moreover, this applies to local optima (Schoen 1986 citing Steinitz 1928); for example, moving the tangent point towards the face center (even without moving the planes of the other faces) increases the sphericity. That maximum (and not just supremum) sphericity exists follows from compactness and monotonicity, as one can always increase sphericity by truncating the sharpest corner. Maximum sphericity $n$-polyhedron is conjectured to be simple (i.e. 3 edges incident on each vertex), and unique up to similarity transformations. Conjecturally, only $n=4$ or 5 has any triangles; only $n=5-11,13$ has any quadrilaterals; no faces have 8 or more edges; the first heptagon appears at $n=33$.

Per (Schoen 1986) (elaborated with additional text and illustrations at (Schoen, website)), for small $n$, maximum sphericity $n$-polyhedra tend to be somewhat irregular, but usually with some symmetry, and sometimes indicate a bigger symmetry that is broken. Most of the polyhedra are produced from the previous one by truncating a vertex, or an adjacent group of 2-4 vertices, and reoptimizing for sphericity. As of this writing, the maximum has only been proved for $n=4,5,6,12$. Only n=4, 5 (a trigonal prism), 6, 12, 32 (conjecturally; a truncated icosahedron) are vertex-transitive; and only 4, 6, 12 are face-transitive (and regular). The regular octahedron and icosahedron should maximize sphericity for their numbers of vertices (and even edges or faces+edges) but not faces.

Imposing a symmetry constraint tends to have a significant sphericity cost, except for appropriate $n$; see "The Roundest Polyhedra with Symmetry Constraints" (Lengyel, Gáspár, Tarnai 2017).

Highest sphericity among Archimedean solids is the snub dodecahedron, and among Archimedean and Catalan solids, the disdyakis triacontahedron (the dual of the truncated icosidodecahedron). Modifying the parameters should yield the highest sphericities among vertex transitive and face transitive polyhedra, respectively. Both polyhedra have the largest number of faces in their classes (at 92 and 120, respectively).

Polyhedra with $O(1)$ nonhexagonal faces

For some moderate $n$, Goldberg polyhedra might give the best sphericity; these use hexagons and 12 pentagons with rotational icosahedral symmetry. (Some definitions of Goldberg polyhedra also allow octahedral symmetry with 6 squares, and tetrahedral symmetry with 4 triangles; all vertices have degree 3.) Having hexagons and 12 pentagons is particularly efficient, and pentagons 'repel' each other (as each creates a local angle deficiency), hence icosahedral symmetry. For other moderate $n$, we might get 12 well-separated pentagons, but without the symmetry.

However, projecting a regular hexagonal grid on a sphere creates distortion, but for best sphericity, we want almost all faces to approximate regular hexagons of the same size.
- Approximately, local distortion is proportional to region width squared, with loss per unit area proportional to distortion squared.
- I conjecture that $n'$ nonhexagonal faces (and no holes) cannot get us within $o(\frac{1}{n'n})$ sphericity of the hexagon ideal (note: $Θ(\frac{1}{n'n})$  $(4≤n'<\sqrt{n})$ is possible by allowing nonadjacent pentagon-heptagon pairs below).

Actually, $O(1)$ hexagons can fill in a triangular face by using degree four vertices to add a handle (with tiny size to avoid wasting surface area), so let us assume that the polyhedron is simply connected.

Given a triangulation of a combinatorial sphere with equilateral triangles and $O(1)$ topological defects (vertices of degree other than 6), we can use a conformal mapping to get a spherical polyhedron with all but $O(1)$ of the triangle angles being within $ε$ of 60° (so all but $O(1)$ triangles are close to equilateral; the $O(1)$ depends on $ε$ and the number of defects). The dual will have all but $O(1)$ faces approximate regular hexagons. However, the triangles (or hexagons) have different sizes, and as we approach a singularity (i.e. defect), size goes $x→Θ(x^c)$ if angle is increased in $c$ times ($x$ (approaching 0) is size; $c$ can be below 1). Nonconformal mappings should give better sphericity, but due to distortion, we are still $Θ(1/n)$ below the hexagon ideal.

Computation of the hexagonal limit

For a polyhedron with an inscribed sphere radius $r$, $V=\frac{r}{3} S$ (and $\frac{r}{n} S$ in $n$ dimensions), where $S$ is the surface area. Thus, given $r=1$ and average distance from the sphere center to the surface $1+ε$ (with the averaging using the solid angle), we get sphericity $1 - ε + O(ε^2)$. In turn, $ε=d^2/2-O(d^4)$ for (root mean square) distance $d$ from the tangential point. Given a regular $m$-gon face tangent at its center and with $\frac{4π}{n}$ solid angle, side length $s=\sqrt{\frac{16π}{mn} \tan \frac{π}{m}}+O(n^{-1})$, and $ε=(\frac{1}{16}\cot^2 \frac{π}{m}+\frac{1}{48}) s^2+O(n^{-2}) = $$ (\cot \frac{π}{m} + {\frac 1 3}\tan \frac{π}{m}) \frac{π}{mn} + O(n^{-2})$. Thus, idealized (1-sphericity)*n approaches $\frac{5 \sqrt 3 π}{27}$ for $m=6$ (and $\frac{π}{3}$ for $m=4$, and $\frac{2\sqrt{3} π}{9}$ for $m=3$). For small $n$, some computations are in (Goldberg 1935).

A construction for $1 - \frac{5 \sqrt{3} π}{27n} - O(n^{-3/2})$ sphericity

We will use an approximately regular hexagonal grid (i.e. faces are hexagons) with $Θ(\sqrt n)$ topological defects; note that radius/typical_edge_length is $Θ(\sqrt n)$. A single adjacent pentagon-heptagon pair can delete a ray in a hexagonal grid (for visualization, see (Abrego 2020)); and intuitively, we need to delete $Ω(\sqrt n)$ rays to cover a large convex region of the sphere. A nonadjacent pair can delete more but with a quadratic area of $Θ(1)$ distortion, so we will not use them.

At unit grid size, deleting a ray requires $Ω(1/r)$ relative distortion at distance $r$, and thus has cumulative cost $Ω(\log r)$ (relative to local unit distortion) unless the geometry or other defects cancel the distortion earlier. This threatens to add a logarithmic factor to the $O(n^{-3/2})$ bound (if we use $Θ(\sqrt n)$ defects at $Ω(n^{-2} \log n)$ per-defect cost), which we will avoid by carefully grouping defects as follows.

Recursively partition the sphere into roughly $n^{1/3}$ pieces, each time cutting a piece along a geodesic (or otherwise with short length), and efficiently reducing maximum piece diameter. From top to bottom, given a piece, choose an approximately parallel unit vector field, and uniformly rotate it for low angle with the vector field of the parent piece. Fill the leaf pieces with an approximate regular hexagonal grid (with the right unit size) oriented parallel to the vector field. Then, recursively (bottom up), join sibling pieces, warping a narrow strip near the separating geodesic for alignment of the grids, with $O(\sqrt n d^3)$ topological defects per join, where $d$ is the diameter of the piece (for unit sphere radius). Note that the defects concentrate along the boundaries of the larger pieces, and the effective radius of distortion for each defect is limited by the distance between the defects, hence $Θ(n^{-2})$ average cost per defect.

An attempt at the matching upper bound

Conjecture: For every triangulation of a radius $R$ sphere into $n$ triangles, $Ω(\sqrt n)$ triangles have $Ω(1)$ relative difference from being an equilateral triangle of area $4πR^2/n$. (MathOverflow question: https://mathoverflow.net/questions/391508/covering-the-sphere-with-an-approximately-planar-grid )

Assuming the conjecture, the maximum sphericity convex polyhedron with $n$-faces must have $Ω(\sqrt n)$ faces with $Ω(1)$ relative deviation from being a regular hexagon of the right size, as otherwise its dual can be converted into an impossible triangulation. But each such face has $Ω(n^{-2})$ sphericity cost (see computation of the hexagonal limit above; for a given solid angle, the regular hexagon is most efficient, and if solid angles are different, the disparity leads to lost sphericity). Thus, we are $Ω(n^{-3/2})$ sphericity below the hexagon ideal.

To see that the conjecture should be true, note that approximately equilateral triangles have a locally unique triangulation topology (assuming no hinge vertices), allowing a grid-based local coordinate system and a global coordinate chart, with (under the coordinates) isometric connections. For a simply-connected region, the coordinate systems can be chosen to agree on all points of overlap. Furthermore, for a 'sturdy' convex region and similar triangle sizes, a single coordinate system gives correct distances within a factor of $1+ε$, which makes it approximately planar. However, a long thin strip can easily bend.
    Next, a nonsimply connected region can have topological defects. A topological defect is characterized by a dislocation and holonomy (which affects how dislocation propagates); orientation is preserved for us. If all defects have zero holonomy, then the dislocations are additive, and due to the intrinsic curvature (for example, circles on spheres have low perimeters for their areas), we need $Ω(\sqrt n)$ coordinate dislocation (assuming similar triangle sizes) to cover a large convex area of the sphere. Nonzero defect holonomy, being quantized, ordinarily gives many distorted (or wrong-sized) triangles, but I do not have a simple characterization to complete the proof.

Here is a generalization of the conjecture to $d>1$ dimensions: Given $n$ points on a $d$-sphere of area $n$, for at least $Ω(n^{1-1/d})$ points, the nearby space must have $Ω(1)$ difference (in terms of distances) from being isometric to a portion of the Cartesian grid. (The 'space' involves both the points and the empty space.)

The conjectured sphericity fractal

The conflict between Euclidean and spherical geometry should create a fractal-like behavior:
- at large scales, defect lines are more efficient than area deformations,
- denser defect lines are more efficient per defect,
- but we also need finer lines to reach all areas of the sphere.

We conjecture as follows.

Convergence: Optimal large scale behavior converges, giving us a picture of the topological defects in the optimal sphericity $n$-polyhedron. Let $T$ denote locations with high density of topological defects (in the limit).

Parameterization: $T$ is determined by the cost per unit length of joining planar hexagonal grids as a function of the angle of misalignment and the direction of the grid. These costs can be approximated computationally. A number of problems besides sphericity will have different costs, but with similar topological and qualitative behavior. For parameterized problems (such as $r^{-s}$ potential below), at some parameter values, we will have two solutions for $T$ (or more, if there are more parameters), but the set of such parameter values has measure zero, and should not include any algebraic values (unless one uses contrived parameterizations or goes outside the domain).

Basic structure: $T$ has 12 connected components, corresponding to equal area regions (60° holonomy). Each component is a tree. $T$ has zero area; no end of $T$ is in $T$. Topologically (with a possible exception for the degree of the root), each connected component of $T$ looks like a binary tree that branches (i.e. every branch branches) at every rational interval, with large denominators corresponding to small branches. Each path is differentiable except at all branch points.

Flow: We have conservation of flow — angle of misalignment between grids; corresponds to area (solid angle) of drainage. At each root (i.e. when the angle crosses 30°), the flow adds up to 60° and cancels out (though in a sense, the cancellation happens everywhere).

Symmetry: The number of components and distance minimization would suggest an approximate icosahedral symmetry, with 5 branches joined at each of the 12 roots (with each branch roughly bisecting an icosahedral face). Any symmetry is approximate as branches spontaneously break the symmetry (and repel each other) and do not have a symmetric effect on opposite sides. However, assuming 24° misalignment per unit length costs less than roughly $2 \cos \frac{π}{5}$ of 12° misalignment (which is somewhat probable), it becomes favorable to join two or four of the branches into pairs. As the savings increase, it becomes profitable to elongate the joined branches. It is unclear if there is an approximate symmetry (especially if the deviations are large), but one possibility is an approximate pyritohedral symmetry, with one direction of the deviation from regularity (on the pyritohedral continuum) if there are 4 branches at each root, and the other if there are 3 branches.

Branch angles and sizes: Branch angles depend on the parent flow $α$ and child flow $β$, locally minimizing cost of flow. Under isotropy, this corresponds to conservation of momentum (nonrelativistic), where speed equals cost per unit flow per unit length. For (differentiable) anisotropy, the conservation still holds if we add a sideways component to the flow momentum. Typically (in terms of diversion of flow) $\min(β,α-β)=Θ(α)$, in which case the angles are $Θ(\log^{-1/2} α^{-1})$; for other $β<α/2$, its angle is still $Θ(\min(1,\sqrt{\frac{\log(α/β)}{\log(1/α)}}))$. For every interval on $T$ with $O(1)$ relative reduction in flow, the distribution of side branch sizes has approximately power law behavior (possibly with a logarithmic or so correction, and without constraining the size of the largest branch).

Isotropy: For small flow $α$, we get approximate isotropy: Grid-induced relative cost anisotropy is $O(\log^{-1} α^{-1})$ per unit angle (deformations can use preferred alignments at scales below the inter-defect distance). Also, $T$ has small typical angles, and a constant-speed detour does not change your average angle. As an aside, if we manually inserted (see parameterization above) flow-independent smooth $Θ(1)$ anisotropy for small $α$: the first order term (change in cost per change in angle) would mostly cancel out (average direction being fixed), leaving rotation within $Θ(1)$ factor of random drift (which is still a bias for the preferred directions); the second order term would create local (not necessarily positive) correlations between branch directions that do not persist across scales (depending on sign, it may favor straight or curved branches).

The complement of $T$: The complement of $T$ is connected and contains a dodecahedral skeleton, but otherwise tree-like, with unique paths between points and the skeleton. Almost all points (all but zero area) lie on the ends. At almost all ends and on uncountably many points on each path, the compass direction of the path diverges with unbounded winding angles. I do not know whether path lengths are finite. A flow $α$ branch perpendicular to the boundary (and not too far from it) would push/perturb the boundary by a typical angle $Θ(\log^{-1} α^{-1})$ (based on savings from the branch versus extra length). This is not enough for infinite lengths, unless branches are nearly parallel to the boundary, in which case the perturbation angles are greater.

The end of branching: For finite $n$, branching of defect lines ends at typical flow $Θ(n^{-1/3} \log^{5/6} n)$ — length $Θ(n^{-1/6} \log^{2/3} n)$ and defect-free region width $Θ(n^{-1/6} \log^{1/6} n)$. At this scale, optimal branching cost for a region (roughly proportional to $\mathrm{area}⋅\mathrm{length}$ for $O(\log^{-1/2} n)$ angles) has the same order of magnitude as region deformation cost (roughly proportional to $\mathrm{area}⋅\mathrm{width}^4$). Due to these savings on branching, optimal sphericity is $1 - \frac{5 \sqrt 3 π}{27n} - C_1 n^{-3/2} + (C_2±o(1)) n^{-5/3} \log^{2/3} n$. By scale symmetry, there should be no terms between the $C_1$ and $C_2$.
Further terms: Further terms (hidden in the above $o(1)$) might bring only polylogarithmic extra precision, starting with a likely $O(n^{-5/3} \log^{1/6} n)$ term due to extra distances at branch ends. Quantization cost is apparently lower at $O(n^{-2})$ per branch join and $O(n^{-2} \log^{1/3} n)$ per branch end, plus quantization savings. In many problems, the quantization pseudorandom noise is (by approximate independence) polynomially lower than the quantization cost, which can be partially expressed here through various bounds on nonmonotonic or oscillatory behavior.

Some related problems

Unit circle packing on a sphere: Despite its superficial similarity (especially for small $n$), unit circle packing on a sphere (Tammes problem) likely has different asymptotic behavior due to apparent inability to smooth out fractional distances at $o(1)$ relative cost. We have $\frac{π\sqrt{3}}{6}-O(n^{-1/3})$ density, using regions of width $Θ(n^{-1/6}R)$ (at sphere radius $R$), with each region approximately hexagonally packed, and with higher scale structure (above $Θ(n^{-1/6}R)$) less important for the packing density. The slow convergence helps explain why it takes hundreds or thousands of circles to overcome the close-to-optimal packing density of just 12 circles; see (Clare and Kepert 1986, 1991).

Maximum volume polyhedra that can be cut from a sphere:
- If the number of vertices is fixed, this is quite similar to sphericity (even if the fractal shape is mildly different).
- If the number of faces is fixed, then (most faces being approximate regular hexagons now) $ε$-sized deformations have $O(ε)$ rather than $O(ε^2)$ relative cost, unless the deformations are conformal. Conformal mapping theory suggests that we can still get within $O(n^{-3/2})$ volume of the hexagonal limit (this does not generalize to higher dimensions). (And we can easily get $O(n^{-7/5})$ by using regions of width $Θ(n^{-1/10})$.)

Minimum length equal area partitions of the sphere: Using area-preserving deformations (see Thomson problem), we get within $Θ(1)$ total perimeter of the hexagonal limit (for unit sphere radius). This is presumably optimal.

Sphericity given numbers of vertices and edges:
* If the maximum sphericity polyhedron must instead have a given number of vertices, we get sphericity within $O(n^{-3/2})$ of the idealized triangle lattice. The resulting polyhedron need not have an inscribed sphere, but using a circumscribed sphere (though presumably suboptimal) works well enough for the bound.
* Approximately square lattice appears to be best when $n_V+kn_F$ is fixed for $≈2.8<k<≈24.5$ (sphericity depends more on the number of faces $n_F$ than the number of vertices $n_V$). For approximately square lattice, we cannot just freely place face planes (like for hexagons) or vertices (like for triangles). However, we can approximate uniform curvature, and though I did not verify it, there appears enough degrees of freedom to allow $O(ε)$ deformations at $O(ε^2)$ relative cost and get the $O(n^{-3/2})$ bound. The corresponding sphericity fractal should consist of 8 equal-area regions, perhaps with an approximate octahedral symmetry.

The Thomson problem

The Thomson problem asks for the minimum electric potential of $n$ classical electrons on the unit sphere, where the potential between each unordered pair of particles at distance $r$ is $r^{-s}$, with $s=1$. We will consider other $s$ as well, with the potential $-\log r$ for $s=0$ and $-r^{-s}$ for $s>0$ (i.e. repulsive field with power law strength).

For $s<-2$, the particles will be evenly divided between two antipodes. For $s=-2$ (corresponding to sum of squared distances), any arrangement with the center matching the center of the sphere will have the same potential (Copher 2020).

For $s>-2$, we expect the distribution to approximate the triangle grid (or hexagonal, if you look at Voronoi regions), with a fractal-like set of defects set analogous to the sphericity fractal.

If $s→∞$, the behavior at the largest scale is unclear (compare the Tammes problem), though it should still converge for fixed $s$, and with smaller scale behavior resembling the sphericity fractal. One possibility is chaotic behavior for scales above $s^{-Θ(1)}$ (though still deterministic and continuous in $s$, except for a measure 0 set of $s$).

Defining the lattice limit: Given a lattice or otherwise sufficiently regular distribution on the plane, we can compute the expected product of charge densities as a function of inter-charge distance $r>0$ (formalized as a distribution (expected value given a compact continuous distribution of $r$) to handle divergence for specific $r$). Then, given a sphere or another surface, the idealized lattice potential is obtained by integrating over (unordered) pairs of points on the surface using the geodesic distance to get the expected product of charge densities, and the actual distance to convert the product into the contribution for potential.
Notes:
* This also works in higher dimensions.
* If the desired number of points is fixed, we first implicitly rescale the lattice to have the same area density.
* If we want a nonuniform continuum limit charge distribution on a surface, we can first rescale the metric in the integral to make that distribution uniform (without changing actual or geodesic distances). This is especially useful for nonspherical $d$-dimensional surfaces with $s<d$, where the poppy-seed bagel theorem does not apply.

Cancelling long-distance interactions: For $s≤d$ (number of dimensions), local potentials are exceeded by long-distance potentials, but we can (mostly) cancel them out as follows. We will use sturdy ($O(1)$ perimeter^2/area) equal area regions with the electron in the center of the region, and add a uniform positive charge to make the total charge zero. This way, all regions have zero net charge and dipole moment, leaving only quadrupole-quadrupole interactions, whose potential falls as $r^{-4-s}$, and higher. For $s≥d$, we can avoid local divergence by zeroing out potential for inter-charge distance $<ε$; the choice of $ε$ does not affect the differences between potentials as along as all configurations have inter-electron distances above $ε$ and curvature is fixed (or $s<d+2$ and $ε→0$).

Modified idealized lattice potential: If for small $s$, the idealized lattice potential does not give the right cancellation, here is a modification (or if it agrees, a computation/approximation method) using cell-based computations. Essentially, given a remote cell, all points in the cell will be given the same weight. In two dimensions, the difference will be (at most) $O(1)$ for $s>0$, $O(\log n)$ for $s=0$ and $O(n^{-s/2})$ for $s<0$, which does not affect the below theorem for $s≥-1$ (or $s>d-3$ in $d$ dimensions). To avoid local divergence, no modification is made for $s≥d$.
    Consider a $d$-dimensional periodic charge distribution (with periodicity in $d$ dimensions), with each cell having zero net charge and dipole moment. Pick a cell center. For every center within distance $<π$ (assuming unit sphere radius) of the cell center, place the two cells on the sphere and at the same geodesic distance and relative orientation between the two centers as they were on the plane/space, and compute their interaction potential (assigning half to the current cell). The choice of the projection does not affect the asymptotic analysis as long as areas are preserved (or charge densities modified accordingly) and we have low distortion near the above geodesic (which is also used for transporting orientations; if we do not get better cancellation elsewhere, Lambert cylindrical equal-area projection with the geodesic is reasonable). Add up the potentials using weighting $(\sin θ \, / \, θ)^{d-1}$ for inter-center geodesic distance $θ$ (this compensates for the spherical geometry), then include the internal potential for the initial cell, and multiply by sphere area over cell area to get the ideal spherical potential. A (perhaps more natural) alternative that agrees to sufficient precision for us is to use all pairs of points (with distance $<π$ and weighting as above) instead of region centers, using the geodesic between the points to place their cell(s) on the sphere (the normalization factor is now sphere area over cubed cell area provided that the first point is restricted to one cell).
    Nonzero net charge is handled by adding uniform compensating charge to the grid, and then removing its contribution, which depends only on the total charge, to the sphere potential. Nonzero dipole moment can be handled by translating the cell boundaries by a fractional amount. For the triangle grid, each cell is the hexagon of the dual grid.

Theorem: For the Thomson problem with $s>-3/2$ (and two dimensions), there is a configuration within $O(n^{(1+s)/2})$ potential of the modified triangle grid ideal. For $s=-3/2$, we get $O(n^{-1/4} \log n)$, and for $-3/2>s>-2$, we get $O(n^{s/6})$.
Notes:
* For $s≤-3/2$, optimality of the above depends on whether there are additional cancellations of long-distance potential changes beyond those in the proof here.
* Just like for sphericity, the difference from the modified triangle grid ideal is expected to be $C_1 n^{(s+1)/2} - (C_2±o(1)) n^{s/2+1/3} \log^{2/3} n$ ($C_1$ and $C_2$ depend on $s$) for $s>-1$. For $-1>s>-3/2$, we should get $C_1 n^{(s+1)/2}±O(n^{s/6})$, and for $s=-1$, $C_1±o(n^{-1/6} \log^2 n)$.
Proof sketch:
    We refine the sphericity construction to use equal-area deformations (which correspond to incompressible flow). Normalize the grid interval at $Θ(1)$ by rescaling distance, and then normalize the force between two charges to be $Θ(1)$ at distance $Θ(1)$ (i.e. at typical distance between adjacent charges). At distance $r$ from a given nearby defect, the displacement is $O(1)$, with $O(1/r)$ gradient (deformation) and $O(1/r^2)$ second order correction. Furthermore, we do not need deformations at scales above the inter-defect distances.
    Now, given two charges at distance $r_0$ between each other and distance $Θ(r)$ from the defect, the uncompensated change in quadrupole-to-quadrupole (and higher) potential is $O(r_0^{-s-3}/r^2)$: This corresponds to $r_0^{-s-4}$ base for the potential, $r$-fold reduction due to small deformation size, and $r/r_0$ reduction due to cancellation of local first order effects (below). The contribution between charges at disparate distance scales from the defect is smaller. Thus, the cumulative cost is $O(\log r)$ if $s>-1$, $O(\log^2 r)$ if $s=-1$, and $O(r^{-1-s})$ otherwise, which suffices.
    In more detail, consider two cells at distances $r_1$ and $r_2$ ($r_1≤r_2$) from the defect, at distance $r_0$ from each other.
Case 1: $r_0$ is $Ω(r_1)$ (and thus $Ω(r_2)$). The change in quadrupole moments is $O(r_1^{-1})$ and $O(r_2^{-1})$, and the change in distance is $O(1)$. Thus, the change in potential is $O(r_1^{-1} r_2^{-4-s})$, which is small enough.
Case 2: $r_0$ is much smaller than $r_1$. We get sufficient cancellation of the linear component because an area-preserving linear transformation equals a product of shear transformations (of comparable magnitude) along the symmetry (reflection) axes of the lattice. Specifically, the interaction potential between two cells centered (before deformation) at 0 and $C$ equals the integral over points $a$ and $b$ (within range) of the potential of the six point system: $-a,0,a,C-b,C,C+b$ with charge 1 at $0$ and $C$ and charge -1/2 at $±a$ and $±b$ (the arithmetic on points is done before the deformation; $|C|=r_0$). The change due to deformation can be decomposed into:
- $O(1)$ shear linear transformations, each of magnitude $O(r_1^{-1})$ (at distance 1), and where for each, all but $O(r_0/r_1)$ fraction of the six-point systems are matched with a system that exactly compensates the length changes
- translation of the $(-b,C,b)$ triple by $O((r_0/r_1)^2)$ (dependent on the triple)
- movement of points by $O(r_1^{-2}d^2)$ where $d=O(1)$ is the distance to the relevant center point.
Each of these has a small enough impact on the potential.
    Also, regarding the total modified idealized lattice potential, we rely on small distance changes compared to what the ideal assumes (smaller scales as above; for all scales, $O(1)$ grid intervals effective change due to area preservation), and use isotropy of the quadrupole moment for regular $n$-gons to deal with large-scale orientation changes.

Higher dimensions

Similar constructions exist for $d$-dimensional spheres and polyhedra approximating them. For all $d$, we can get within $O(n^{1-1/d})$ defects of any fixed $d$-dimensional lattice, with smoothing of deformations as above. This amounts to $1-C n^{-2/d}-O(n^{-3/d})$ sphericity ($C$ depends on the lattice). If the lattice is optimal, the optimal choice of defects should converge to a $d-1$ dimensional fractal.

Typically (across problems), at misalignment angle $α$, the relative cost per defect is still $Θ(\log (1/α))$ (where $α$ is small and the cost at large angles is normalized at $Θ(1)$). And $O((1/α)^{1/d-ε})$ relative cost suffices to get within the cost of having $O(n^{1-1/d})$ unit-cost defects.

For $d=8$, the $E_8$ lattice is generally optimal, and for $d=24$, the Leech lattice is. For other $d>2$, especially not too small $d$, it is unclear if there is an optimal lattice (which may depend on the optimization criteria). If not, then the divergence from the Euclidean $d$-dimensional limit might be less, as instead of $d-1$-dimensional defect streams, we might tweak the arrangements to avoid defects (less likely, the divergence can be more than in the lattice case due to difficulty of smoothing deformations).

For $r^{-s}$ potential, the quadrupole-quadrupole potentials exceed local connections for $d≥s+4$, perhaps giving a new behavior (unless we get a higher cancellation). For $d<s+3$, using our smoothing of deformations, we are within $O(n^{(d-1+s)/d})$ potential of the lattice limit/ideal, at least if the lattice has enough reflection symmetries to get the cancellation used in the proof of the theorem. For $d≥s+3$, the quadrupole-based divergence from the nonmodified lattice limit might be larger. It is plausible that we have enough cancellations (and costs if we do otherwise) to give the fractal-like optimal behavior all the way for $s>-2$.

References

(Abrego 2020) Abrego, Andoni. "Particles and fields from topological defects in spacetime mesh." https://community.wolfram.com/groups/-/m/t/2030734?p_p_auth=e4lQ49tu . Accessed: Oct 6, 2021.
(Clare and Kepert 1986) Clare B.W. and Kepert D.L. "The closest packing of equal circles on a sphere." Proc. R. Soc. Lond. A405329–344 (1986). http://doi.org/10.1098/rspa.1986.0056.
(Clare and Kepert 1991) Clare B.W. and Kepert D.L. "The optimal packing of circles on a sphere." J Math Chem 6, 325–349 (1991). https://doi.org/10.1007/BF01192589.
(Copher 2020) Copher, Jessica. "Sums of Squared Distances between Points on a Unit n-Sphere." (2020) arXiv:2001.02745.
(Goldberg 1935) Goldberg, Micahel. "The isoperimetric problem for polyhedra." Tôhoku Math. J. 1935, 40, 226–236.
(Lengyel, Gáspár, Tarnai 2017) Lengyel, András, Zsolt Gáspár and Tibor Tarnai. "The Roundest Polyhedra with Symmetry Constraints." Symmetry 9 (2017): 41.
. https://doi.org/10.3390/sym9030041.
(Schoen 1986) Schoen, Alan. "A defect-correction algorithm for minimizing the volume of a simple polyhedron which circumscribes a sphere." Proceedings of the second annual symposium on Computational geometry, August 1986, pp 159–168. https://doi.org/10.1145/10515.10533.
(Schoen, website) Schoen, Alan. "Roundest Polyhedra." https://schoengeometry.com/a_poly.html . Accessed: Oct 6, 2021.