Locally optimal 2-periodic sphere packings

TL;DR

This paper extends Voronoi's method to enumerate all locally optimal 2-periodic sphere packings in various dimensions, showing that in dimensions 3, 4, and 5, no such packing exceeds the density of optimal lattices.

Contribution

It generalizes Voronoi's enumeration algorithm to 2-periodic packings and implements it for dimensions up to 6, providing new insights into local optimality.

Findings

No 2-periodic packing exceeds the density of optimal lattices in dimensions 3, 4, and 5.

Partial enumeration performed in dimension 6.

Method successfully enumerates all locally optimal 2-periodic packings when finitely many exist.

Abstract

The sphere packing problem is an old puzzle. We consider packings with m spheres in the unit cell (m-periodic packings). For the case m = 1 (lattice packings), Voronoi proved there are finitely many inequivalent local optima and presented an algorithm to enumerate them, and this computation has been implemented in up to d = 8 dimensions. We generalize Voronoi's method to m > 1 and present a procedure to enumerate all locally optimal 2-periodic sphere packings in any dimension, provided there are finitely many. We implement this computation in d = 3, 4, and 5 and show that no 2-periodic packing surpasses the density of the optimal lattices in these dimensions. A partial enumeration is performed in d = 6.

Figures

Tables

Table 1. Table 1: Summary of enumeration results. For d = 3 , 4 , 5 𝑑 3 4 5 d=3,4,5 we list the number of vertices of the Ryshkov-like polyhedron, the number of algebraically extreme 2-periodic sets, the highest number density achieved for packing radius 1 1 1 , and the number of algebraically extreme sets that achieve this density. In parentheses, we indicate how many of the corresponding 2-periodic sets are not also lattices.

vertices of $ℛ\left(4\right)$	4 (2)	10 (6)	34 (25)
algebraically extreme 2-periodic sets	3 (1)	7 (3)	29 (20)
highest density	$1/\left(2\sqrt{2}\right)$	$1/8$	$1/\left(8\sqrt{2}\right)$
multiplicity of highest density	3 (1)	2 (0)	5 (2)

Table 2. Table 2: Some invariants of the 29 algebraically extreme 2-periodic sets in d = 5 𝑑 5 d=5 dimensions, including density and kissing number. The L/N column indicates whether this is the representation of a lattice as a 2-periodic set ( L ) or a truly nonlattice arrangement ( N ). Here α 1 = 1 3 ​ ( tan − 1 ⁡ ( 3 ​ 762 ) − π ) subscript 𝛼 1 1 3 superscript 1 3 762 𝜋 \alpha_{1}=\tfrac{1}{3}(\tan^{-1}(3\sqrt{762})-\pi) , α 2 = 1 3 ​ ( tan − 1 ⁡ ( 9 ​ 47 / 17 ) − π ) subscript 𝛼 2 1 3 superscript 1 9 47 17 𝜋 \alpha_{2}=\tfrac{1}{3}(\tan^{-1}(9\sqrt{47}/17)-\pi) , α 3 = 1 3 ​ ( tan − 1 ⁡ ( 486 ​ 1077 / 5867 ) − π ) subscript 𝛼 3 1 3 superscript 1 486 1077 5867 𝜋 \alpha_{3}=\tfrac{1}{3}(\tan^{-1}(486\sqrt{1077}/5867)-\pi) , and x 3 = 8189832078 + 22432320 ​ 661 ​ cos ⁡ α 3 − 35353733064 ​ cos ⁡ 2 ​ α 3 − 3050582422 ​ cos ⁡ 4 ​ α 3 + 6990736 ​ 661 ​ cos ⁡ 5 ​ α 3 subscript 𝑥 3 8189832078 22432320 661 subscript 𝛼 3 35353733064 2 subscript 𝛼 3 3050582422 4 subscript 𝛼 3 6990736 661 5 subscript 𝛼 3 x_{3}=8189832078+22432320\sqrt{661}\cos{\alpha_{3}}-35353733064\cos{2\alpha_{3}}-3050582422\cos{4\alpha_{3}}+6990736\sqrt{661}\cos{5\alpha_{3}} .

#	$\delta$	$\delta$	L/N	$\kappa$
1	0.0884	$1/\left(8\sqrt{2}\right)$	L	40
2	0.0884	$1/\left(8\sqrt{2}\right)$	L	40
3	0.0884	$1/\left(8\sqrt{2}\right)$	L	40
4	0.0884	$1/\left(8\sqrt{2}\right)$	N	40
5	0.0884	$1/\left(8\sqrt{2}\right)$	N	40
6	0.0795	$\sqrt{5-2\sqrt{6}}/4$	N	35
7	0.0795	$\sqrt{5-2\sqrt{6}}/4$	N	35
8	0.0795	$\sqrt{5-2\sqrt{6}}/4$	N	34
9	0.0786	$1/\left(9\sqrt{2}\right)$	L	30
10	0.0786	$1/\left(9\sqrt{2}\right)$	L	30
11	0.0786	$1/\left(9\sqrt{2}\right)$	L	30
12	0.0786	$1/\left(9\sqrt{2}\right)$	N	30
13	0.0786	$1/\left(9\sqrt{2}\right)$	N	30
14	0.0771	$\sqrt{419+1011\sqrt{33}}/1024$	N	30
15	0.0765	$-\sqrt{2}+\left(2\sqrt{5}\right)/3$	N	29
16	0.0765	$-\sqrt{2}+\left(2\sqrt{5}\right)/3$	N	29
17	0.0758	$3/\left(2\sqrt{126+4\sqrt{19}\cos {\alpha }_1-266\cos 2{\alpha }_1}\right)$	N	30
18	0.0750	$3/40$	N	30
19	0.0748	$1/\left(36\sqrt{-58+26\sqrt{5}}\right)$	N	26
20	0.0748	$\left(3\sqrt{3}\right)/\left(\left(\sqrt{\left(3-4\cos {\alpha }_2-8\cos 2{\alpha }_2\right)}\right)\left(4\left(5+4\cos {\alpha }_2\right)\right)\right)$	N	28
21	0.0738	$59049/\left(4\sqrt{x_3}\right)$	N	25
22	0.0737	$5/\left(48\sqrt{2}\right)$	N	30
23	0.0722	$1/\left(8\sqrt{3}\right)$	L	30
24	0.0722	$1/\left(8\sqrt{3}\right)$	L	30
25	0.0722	$1/\left(8\sqrt{3}\right)$	L	30
26	0.0722	$1/\left(8\sqrt{3}\right)$	N	30
27	0.0722	$1/\left(8\sqrt{3}\right)$	N	30
28	0.0722	$1/\left(8\sqrt{3}\right)$	N	30
29	0.0722	$1/\left(8\sqrt{3}\right)$	N

Equations

$$\delta (\Xi )=\frac{m}{\det \Lambda }\text{,}$$

$$\rho (\Lambda )=\frac{1}{2}\underset{\mathbf{l},{\mathbf{l}}^{\prime }\in \Lambda }{\min }\Vert \mathbf{l}-{\mathbf{l}}^{\prime }\Vert =\frac{1}{2}\underset{\mathbf{l}\in \Lambda }{\min }\Vert \mathbf{l}\Vert \text{.}$$

$$Q(\mathbf{n})=⟨\mathbf{n}{\mathbf{n}}^T,Q⟩\ge \lambda \text{ for all }\mathbf{n}\in {\mathbb{Z}}^d$$

$$J(n_1,\dots ,n_d,l_1,\dots ,l_{m-1})=\Vert n_1{\mathbf{a}}_1+\dots +n_d{\mathbf{a}}_d+l_1{\mathbf{b}}_1+\dots +l_{m-1}{\mathbf{b}}_{m-1}{\Vert }^2\text{.}$$

$$J=\left(\begin{array}[]{cc}Q&R^{T}\\ R&S\end{array}\right)\text{,}$$

$$\rho (\Xi )$$

$$={\left(\underset{\mathbf{l},\in \Lambda ,0\le i,j\le m-1}{\min }\Vert \mathbf{l}+{\mathbf{b}}_i-{\mathbf{b}}_j{\Vert }^2\right)}^{1/2}$$

$$=\frac{1}{2}(\min J{)}^{1/2}\text{,}$$

$$J\in \mathcal{R}(\lambda )\text{ iff }J(\mathbf{k})=⟨\mathbf{k}{\mathbf{k}}^T,J⟩\ge \lambda \text{ for all }\mathbf{k}\in \mathcal{M}\setminus \{0\}\text{.}$$

$$J(\mathbf{k})$$

$$=(\mathbf{n}+\mathbf{q}{)}^{T}Q(\mathbf{n}+\mathbf{q})+{\mathbf{l}}^{T}S\mathbf{l}-{\mathbf{q}}^{T}Q\mathbf{q}\text{,}$$

$$T=\left(\begin{array}[]{cc}U&V\\ 0&W\end{array}\right)\text{.}$$

$$W=\left(\begin{array}[]{cccc}0&1&0&\cdots\\ 0&0&1&\cdots\\ \vdots&\vdots&\vdots&\ddots\end{array}\right)\Pi\left(\begin{array}[]{ccc}-1&-1&\cdots\\ 1&0&\cdots\\ 0&1&\cdots\\ \vdots&\vdots&\ddots\end{array}\right)\text{,}$$

$$\text{minimize }\mathrm{obj}(J)=\det Q(J)\text{ subject to }J\in {\mathcal{R}}_0(4{\rho }_0^2)\text{,}$$

$$T_J=J+{\left(\mathcal{S}[N(J)]\right)}^{\perp }=\{J^{\prime }:⟨J^{\prime }-J,\mathbf{u}{\mathbf{u}}^T⟩=0\text{ for all }\mathbf{u}\in N(J)\}\text{,}$$

$$C_J=\{J^{\prime }:J^{\prime }(\mathbf{k})\ge \lambda \text{ for all }\mathbf{k}\in \mathcal{M}\text{ such that }J(\mathbf{k})=\lambda \}\text{.}$$

$$G_{J}=\left(\begin{array}[]{cc}Q^{-1}&0\\ 0&0\end{array}\right)\text{.}$$

$$J^{\prime}(t)=\left(\begin{array}[]{cc}Q&R+t\tilde{R}\\ R^{T}+\tilde{R}^{T}&S+t\tilde{S}+t^{2}\tilde{R}^{T}Q^{-1}\tilde{R}\end{array}\right)\text{,}$$

$$J^{\prime}(t)=\left(\begin{array}[]{cc}AA^{T}&(1-t)A\mathbf{b}(t)\\ \mathbf{b}(t)^{T}A^{T}&\|\mathbf{b}(t)\|^{2}\end{array}\right)\text{,}$$

$$\mathrm{obj}(J^{\prime })\text{,}$$

$$g_i(J^{\prime })\triangleq 4-⟨{\mathbf{k}}_i{\mathbf{k}}_i^T,J^{\prime }⟩\le 0\text{ for }i=1,\dots ,s\text{,}$$

$$h(J^{\prime })\triangleq {\lambda }_{\text{min}}(J^{\prime })=0\text{.}$$

$$C_0=\{\tilde{J}:$$

$$⟨\nabla g_i,\tilde{J}⟩\ge 0\text{ for all }i\in I^0,$$

$$⟨\nabla h,\tilde{J}⟩=0\}\text{,}$$

$$J^{\prime }([\mathbf{v}](t))=t^2{c}^2{\mathbf{u}}^T{Q}^{\prime }\mathbf{u}-tc{\mathbf{u}}^T{Q}^{\prime }{\mathbf{e}}^T+\mathbf{e}{Q}^{\prime }\mathbf{e}-2t{c}^2+{\mathbf{l}}^{T}S\mathbf{l}\text{.}$$

$$J^{\prime }(U(\mathbf{n}),\mathbf{l})=({\mathbf{n}}^{\prime }{)}^T{U}^T{Q}^{\prime }U\mathbf{n}+2({\mathbf{n}}^{\prime }{)}^T{U}^T(R^{\prime }{)}^T\mathbf{l}+2({\mathbf{n}}^{\prime \prime }{)}^T{U}^T(R^{\prime }{)}^T\mathbf{l}+{\mathbf{l}}^T{S}^{\prime }\mathbf{l}\text{.}$$

$$({\mathbf{k}}_{\sigma (i)}^{\prime }{)}^T{J}_2{\mathbf{k}}_{\sigma (j)}^{\prime }={\mathbf{k}}_i^T{T}^T{J}_{1}T{\mathbf{k}}_j\text{ for all }1\le i,j\le d+m-1\text{.}$$

$$J_{3,1}=2\left(\begin{array}[]{cccc}2&0&0&-1\\ 0&2&0&1\\ 0&0&4&2\\ -1&1&2&2\end{array}\right)\quad J_{3,2}=2\left(\begin{array}[]{cccc}2&-1&1&0\\ -1&2&0&1\\ 1&0&6&3\\ 0&1&3&2\end{array}\right)$$

$$J_{3,3}=\frac{2}{3}\left(\begin{array}[]{cccc}6&3&0&0\\ 3&6&0&3\\ 0&0&16&8\\ 0&3&8&6\end{array}\right)$$

$$J_{3,1v}=2\left(\begin{array}[]{cccc}2&0&0&1\\ 0&2&0&-1\\ 0&0&2&1\\ 1&-1&1&2\end{array}\right)\quad J_{3,2v}=2\left(\begin{array}[]{cccc}2&1&0&0\\ 1&2&0&1\\ 0&0&2&1\\ 0&1&1&2\end{array}\right)$$

$$J_{3,3v}=2\left(\begin{array}[]{cccc}2&1&-1&0\\ 1&2&0&1\\ -1&0&2&1\\ 0&1&1&2\end{array}\right)\quad J_{3,4v}=2\left(\begin{array}[]{cccc}2&-1&-1&0\\ -1&2&1&1\\ -1&1&2&0\\ 0&1&0&2\end{array}\right)\text{.}$$

$$J_{4,1}=2\left(\begin{array}[]{ccccc}2&0&0&0&1\\ 0&2&0&0&-1\\ 0&0&2&0&-1\\ 0&0&0&2&-1\\ 1&-1&-1&-1&2\end{array}\right)\quad J_{4,2}=2\left(\begin{array}[]{ccccc}2&-1&-1&0&-1\\ -1&2&1&0&1\\ -1&1&2&0&0\\ 0&0&0&4&2\\ -1&1&0&2&2\end{array}\right)\text{.}$$

$$J_{4,3}=2\left(\begin{array}[]{ccccc}2&-1&-1&-1&-1\\ -1&2&1&1&0\\ -1&1&2&1&1\\ -1&1&1&2+2\tau&-\tau\\ -1&0&1&-\tau&2\end{array}\right)\text{,}$$

Full text

††journal: Journal of LaTeX Templates

1 Introduction

The sphere packing problem asks for the highest possible density achieved by an arrangement of nonoverlapping spheres in a Euclidean space of $d$ dimensions. The exact solutions are known in $d=2$ [1], $d=3$ [2], $d=8$ [3] and $d=24$ [4]. This natural geometric problem is useful as a model of material systems and their phase transitions, and, even in nonphysical dimensions,it is related to fundamental questions about crystallization and the glass transition [5].The sphere packing problem also arises in the problem of designing an optimal error-correcting code for a continuous noisy communication channel [6]. Perhaps even more than the obvious applications of the sphere packing problem, contributing to its importance is the unexpected wealth of remarkable structures that appear as possible solutions and merit study in their own right [7, 8].

There does not appear to be a systematic solution or construction that achieves the optimal packing in every dimension, and every dimension seems to have its own quirks, unearthing new surprises [6]. In some dimensions, such as $d=8$ and $d=24$, the solution is unique and given by exceptionally symmetric lattices. In others, such as $d=3$, there is a lattice that achieves the highest density, but this density can also be achieved by other periodic packings with larger fundamental unit cells and even by aperiodic packings. For $d=10$ it seems that the highest density, achieved by a periodic arrangement with $40$ spheres per unit cell, cannot be achieved by a lattice. In any given dimension, the densest packing known is periodic, but it is not known whether there is some dimension where the densest packing is not periodic. In fact, frighteningly little is known in general as we go up in dimensions. It is possible that in all dimensions the densest packing is achieved by a periodic packing with a universally bounded number of spheres in the unit cell, that in all dimensions it is achieved by a periodic packing, but only with unboundedly many spheres in the unit cell, or that in some dimension it is not achieved by a periodic packing at all. In any case, periodic packings with arbitrarily many spheres in the unit cell can approximate the optimal density in any dimension to arbitrary precision.

Therefore, it is reasonable to ask for the optimal density achieved by a periodic packing in $d$ dimensions with $m$ spheres per unit cell, which we denote $\phi_{d,m}$. In the limit $m\to\infty$, this density approaches (and possibly equals for some finite $m$) the optimal packing density $\phi_{d}$. Apart from the dimensions where $\phi_{d}$ is known, there are no cases with $m>1$ where $\phi_{d,m}$ is known. In this paper, we describe a method for calculating $\phi_{d,2}$, and use it to obtain $\phi_{4,2}=\phi_{4,1}$ and $\phi_{5,2}=\phi_{5,1}$. Our method enumerates all inequivalent local optima, and therefore does not terminate if there are infinitely many. However, we conjecture that there are finitely many inequivalent local optima in each dimension, and therefore that our method reduces the problem of calculating $\phi_{d,2}$ to a finite calculation. Our method could also be generalized to larger $m$, but becomes more complicated.

For the case $m=1$, corresponding to lattices, Voronoi gave an algorithm to enumerate all the locally optimal solutions to the problem. In geometrical terms, Voronoi’s algorithm uses the fact that the space of lattices up to isometry can be parameterized (redundantly) by positive definite quadratic forms. The subset in the linear space of quadratic forms that corresponds to lattices with no two points less than a certain distance apart is a locally finite polyhedron (a closed convex set, every intersection of which with a compact polytope is a compact polytope; see A). The additional fact that the density of lattice points is a quasiconvex function implies that local maxima can only occur at the vertices of this polyhedron. Because there are only finitely many vertices that correspond to distinct lattices, the local optima of the lattice sphere packing problem can be fully enumerated. This interpretation of Voronoi’s algorithm was suggested by Ryshkov [9] and so the polyhedron is known as the Ryshkov polyhedron. The lattices corresponding to vertices of the polyhedron are known as perfect lattices, but not every perfect lattice is locally optimal, or an extreme lattice. This algorithm has been implemented and executed to determine $\phi_{d,1}$ for $d\leq 8$, as well as to enumerate all perfect lattices in these dimensions [10]. However, the exploding number of perfect lattices as $d$ increases make the execution of this algorithm currently unfeasible even for $d=9$ [11].

Schürmann transported some of the concepts from Voronoi’s theory to the case of periodic packings and was able to show that any extreme lattice, when viewed as a periodic packing, is still locally optimal [12, 13]. However the nonlinearities that arise in the case $m>1$ have made it difficult to use this version of the theory to provide an analog of Voronoi’s algorithm. In this paper we transport the concepts in a slightly different way than Schürmann, permitting us to reconstruct many of the elements of Voronoi’s algorithm in the case of periodic packings with fixed $m$. These elements behave nicely enough for $m=2$ to allow us to provide an algorithm for the enumeration of all locally optimal packings, and we execute this algorithm for $d\leq 5$.

2 Theoretical preparation

2.1 Periodic point sets

For a set of points, $\Xi\subset\mathbb{R}^{d}$, the packing radius is the largest radius, such that if balls were centered at each of the points, they would not overlap. Namely, the packing radius is half the infimum distance between any two points, $\rho(\Xi)=\tfrac{1}{2}\mathrm{inf}_{x,y\in\Xi}\|x-y\|$. A set of points is periodic if there are $d$ linearly independent vectors $\mathbf{a}_{1},\ldots,\mathbf{a}_{d}$, such that the set is invariant under translation by any of these vectors. The group generated by $\mathbf{a}_{1},\ldots,\mathbf{a}_{d}$ is a lattice, $\Lambda=\{\sum_{i=1}^{d}n_{i}\mathbf{a}_{i}:n_{1},\ldots,n_{d}\in\mathbb{Z}\}$. If all translations that fix the set are in $\Lambda$, we say that $\Lambda$ is a primitive lattice for the periodic set. The primitive lattice is unique. A set that is periodic under translation by a lattice $\Lambda$ and has $m$ orbits under translations by $\Lambda$ called $m$-periodic. We will only deal with sets of positive packing radius, and therefore $m$ is necessarily finite.

The number density of a periodic set $\Xi$ is given by

[TABLE]

where $\det\Lambda$ is the volume of a parallelotope generated by the generators of $\Lambda$, $\{\sum_{i=1}^{d}x_{i}\mathbf{a}_{i}:0\leq x_{1},\ldots,x_{d}\leq 1\}$. This parallelotope and its translates under $\Lambda$ tile $\mathbb{R}^{d}$, and its volume is the same as the volume of any other fundamental cell of $\Lambda$, that is, a polytope whose $\Lambda$-translates tile $\mathbb{R}^{d}$. Therefore, this is a property of $\Lambda$, and not of the particular basis. Also, since $\Lambda$ itself is not uniquely determined from $\Xi$, but can be any sublattice of the primitive lattice of $\Xi$, it is important to note that the formula for $\delta(\Xi)$ above is independent of the choice of $\Lambda$. The largest density that can be achieved by a packing of equal-sized balls centered at the points of $\Xi$ is $\phi(\Xi)=V_{d}\rho(\Xi)^{d}\delta(\Xi)$, where $V_{d}$ is the volume of a unit ball in $\mathbb{R}^{d}$. In any family of point sets closed under homothety, such as periodic sets or $m$-periodic sets, maximizing the packing density $\phi$ is equivalent to maximizing the number density under the constraint $\rho(\Xi)\geq\rho_{0}$ for some fixed $\rho_{0}$.

2.2 Voronoi’s algorithm

A $1$-periodic packing is a translate of a lattice. The packing radius of a lattice is given by

[TABLE]

The packing radius of a lattice and the volume of its fundamental cell are both invariant under rotations of the lattice. Therefore, to find the densest lattice packing, we only need to consider lattices up to rotation. Consider a lattice $\Lambda=A\mathbb{Z}^{d}$ generated by $\mathbf{a}_{1},\ldots,\mathbf{a}_{d}$. The quadratic form $Q\colon\mathbb{Z}^{d}\to\mathbb{R}$, $Q(n_{1},\ldots,n_{d})=\|n_{1}\mathbf{a}_{1}+\ldots+n_{d}\mathbf{a}_{d}\|^{2}$ determines $\Lambda$ up to rotations. However, since a lattice can have different generating vectors, different quadratic forms can correspond to the same lattice. $Q$ and $Q^{\prime}$ correspond to the same lattice if and only if $Q=Q\circ U$, where $U\in GL_{d}(\mathbb{Z})$. This is precisely the group of linear maps that map $\mathbb{Z}^{d}$ to itself. The packing radius of a lattice $\Lambda$ corresponding to the quadratic form $Q$ is $\rho(\lambda)=\tfrac{1}{2}(\min Q)^{1/2}$, where $\min Q$ is the minimum of $Q$ over nonzero vectors. The determinant of $\Lambda$ is given by $(\det Q)^{1/2}$.

The linear space of quadratic forms can be identified with the linear space of symmetric matrices $\mathcal{S}^{d}$ using the standard basis of $\mathbb{Z}^{d}$, so that $Q(\mathbf{n})=\mathbf{n}^{T}Q\mathbf{n}$. And the natural inner product in this space is given by $\langle Q,Q^{\prime}\rangle=\operatorname{tr}QQ^{\prime}$. The condition $\min Q\geq\lambda$ can be written as the intersection of infinitely many linear inequalities:

[TABLE]

It can be easily checked that for $\lambda>0$, Eq. (3) implies that $Q$ is a positive definite matrix, and a lattice corresponding to it can be recovered, e.g.  by Cholesky decomposition.

The set $\{Q:\min Q\geq\lambda\}$ is a locally finite polyhedron, that is, a closed convex set such that every intersection of it with a compact polytope is a compact polytope. Loosely speaking, though it is not a polyhedron (the intersection of finitely many halfspaces), it behaves locally like one. We were not able to find a satisfactory reference on locally finite polyhedra, and so we prove some basic results about them in A. In the remainder of the text, we use the less cumbersome polyhedron to refer to any locally finite polyhedron. Ryshkov reinterpreted Voronoi’s algorithm in terms of the polyhedron $\{Q:\min Q\geq\lambda\}$, which is therefore known as the Ryshkov polyhedron [9]. Though the Ryshkov polyhedron has infinitely many vertices, there are only finitely many orbits under the action of $GL_{d}(\mathbb{Z})$. Therefore, if we start at some vertex, enumerate the vertices with which it shares an edge, and repeat the process for any enumerated vertex not of the same orbit as a previous vertex, we eventually have representatives of all the orbits of vertices in a connected component of the 1-skeleton of the Ryshkov polyhedron. Since the 1-skeleton of a locally finite polyhedron is connected, these are all the orbits of vertices.

We wish to find the minimum of $\det Q$ in this polyhedron. Because the determinant is a quasiconcave function, its local minima in the polyhedron can occur only on the vertices, and to find the global minimum it is enough to compare its value at all the vertices.

2.3 The Ryshkov-like Polyhedron

To parameterize $m$-periodic point sets, we can use a similar scheme. An $m$-periodic point set is the union of $m$ translates of a lattice, $\Xi=\bigcup_{i=0}^{m-1}(\Lambda+\mathbf{b}_{i})$, where $\Lambda=A\mathbb{Z}^{d}$. Without loss of generality, we can take $\mathbf{b}_{0}=0$. Now, consider the set $\mathcal{M}\subseteq\mathbb{Z}^{d+m-1}$ given by $\mathcal{M}=\mathbb{Z}^{d}\times(E-E)$, where $E=\{0,\mathbf{e}_{1},\ldots,\mathbf{e}_{m-1}\}$ is the standard basis of $\mathbb{Z}^{m-1}$ plus the zero vector. We define the following function $J:\mathcal{M}\to\mathbb{R}$:

[TABLE]

This function can be extended uniquely to a quadratic form over $\mathbb{R}^{d+m-1}$, represented by a symmetric matrix $J\in\mathcal{S}^{d+m-1}$, namely

[TABLE]

where $Q_{ij}=\mathbf{a}_{i}\cdot\mathbf{a}_{j}$, $R_{ij}=\mathbf{b}_{i}\cdot\mathbf{a}_{j}$, and $S_{ij}=\mathbf{b}_{i}\cdot\mathbf{b}_{j}$. The quadratic form, and therefore also the function $J:\mathcal{M}\to\mathbb{R}$, determines $\Xi$ up to rotation.

The packing radius of $\Xi$ is given by

[TABLE]

where $\min J$ is the minimum of $J$ over $\mathcal{M}\setminus\{0\}$.

Let $\mathcal{R}(\lambda)=\{J:\min J\geq\lambda\}$. This set is defined as the intersection of linear inequalities:

[TABLE]

We want to show that $\mathcal{R}(\lambda)$, like the Ryshkov polyhedron, is a locally finite polyhedron.

Lemma 1**.**

Let $Q:\mathbb{R}^{d}\to\mathbb{R}$ be a positive definite quadratic form, satisfying $Q(\mathbf{n})\geq\lambda$ for all $\mathbf{n}\in\mathbb{Z}^{d}$ and $\mathrm{tr}\,Q\leq C$. Then any vector $\mathbf{x}\in\mathbb{R}^{d}$ satisfying $Q(\mathbf{x})\leq 1$ also satisfies $\|\mathbf{x}\|\leq M$, where $M$ depends on $d$, $\lambda$, and $C$.

The proof is given in Theorem 3.1 of [10].

Theorem 1**.**

Let $W=\{J\in\mathcal{R}(\lambda):|J_{ij}|\leq L\}$. There are only finitely many $\mathbf{k}\in\mathcal{M}$ such that for some $J\in W$, $J(\mathbf{k})\leq\lambda$ .

Proof.

Let $\mathbf{k}=(\mathbf{n},\mathbf{l})$. There are only finitely many choices for $\mathbf{l}$. So, we need only show that for any fixed choice of $\mathbf{l}$, there are finitely many $\mathbf{n}\in\mathbb{Z}^{d}$ such that $J(\mathbf{k})\leq\lambda$ for some $J\in W$. We have

[TABLE]

where $\mathbf{q}=Q^{-1}R^{T}\mathbf{l}$. Therefore, if $J(\mathbf{k})\leq\lambda$ then $(\mathbf{n}+\mathbf{q})^{T}Q(\mathbf{n}+\mathbf{q})\leq\lambda+\mathbf{q}^{T}Q\mathbf{q}-\mathbf{l}^{T}S\mathbf{l}\leq\lambda+\mathbf{q}^{T}Q\mathbf{q}=\lambda+\mathbf{l}^{T}RQ^{-1}R^{T}\mathbf{l}$.

Since $\mathbf{l}^{T}RR^{T}\mathbf{l}\leq 4dL^{2}$ and $\mathrm{tr}Q\leq dL$, we have from the lemma that $\mathbf{l}^{T}RQ^{-1}R^{T}\mathbf{l}\leq 4dL^{2}M^{2}$. Therefore, if $J(\mathbf{k})\leq\lambda$ then $(\mathbf{n}+\mathbf{q})^{T}Q(\mathbf{n}+\mathbf{q})\leq\lambda+4dL^{2}M^{2}$. Again, from the lemma, we have that $(\mathbf{n}+\mathbf{q})^{T}(\mathbf{n}+\mathbf{q})\leq\lambda M^{2}+4dL^{2}M^{4}$. Therefore, there can only be finitely many choices of $\mathbf{n}$ for each choice of $\mathbf{l}$ such that $J(\mathbf{k})\leq\lambda$ for some $J\in W$. ∎

So, any intersection of $\mathcal{R}(\lambda)$ with a compact polytope is a compact polytope and $\mathcal{R}(\lambda)$ is a locally finite polyhedron. We call it the Ryshkov-like polyhedron. Moreover, since each vertex of $\mathcal{R}(\lambda)$ is the unique solution of a linear system of equations with integer coefficients, the vertices of $\mathcal{R}(\lambda)$ have rational coordinates.

2.4 Symmetries of the Ryshkov-like polyhedron

The symmetries of $\mathcal{R}(\lambda)$ are tightly linked with the symmetries of the set $\mathcal{M}\subset\mathbb{Z}^{d+m-1}$. Namely, if $T:\mathbb{R}^{d+m-1}\to\mathbb{R}^{d+m-1}$ is a linear map such that $T(\mathcal{M})=\mathcal{M}$, then $J\mapsto J\circ T$ maps $\mathcal{R}(\lambda)$ to itself.

We decompose $T$ into blocks with a top-left block of size $d\times d$. It is easy to see that the bottom-left block must be zero, or else there is always some $\mathbf{k}\in\mathcal{M}$ such that $T(\mathbf{k})=(\mathbf{n},\mathbf{l})$, where $\mathbf{l}\not\in(E-E)$. Therefore we write

[TABLE]

As a consequence, we also have that $U$ and $W$ must be invertible as maps $\mathbb{Z}^{d}\to\mathbb{Z}^{d}$ and $(E-E)\to(E-E)$. Therefore, $U\in GL_{d}(\mathbb{Z})$ and $W$ is a permutation of $E$, namely

[TABLE]

where $\Pi$ is a $m\times m$ permutation matrix.

It is easy to check that whenever $U\in GL_{n}(d)$, $W$ is a permutation of $E$, and $V$ is an arbitrary $d\times(m-1)$ integer matrix, then $T$ is a symmetry of $\mathcal{M}$. Let us call this group of symmetries $\Gamma$, and use its members to act on elements of $\mathcal{M}$ by $\mathbf{k}\mapsto T\mathbf{k}$ or elements of $\mathcal{R}(\lambda)$ by $J\mapsto J\circ T$. We conjecture that for any $\lambda>0$, $\mathcal{R}(\lambda)$ has finitely many faces (in particular vertices) up to the action of $\Gamma$, but we do not have a proof, except in the cases where we have enumerated the vertices (see Sec. 4), and found, by the fact of the algorithm halting, that there were finitely many vertices.

2.5 The rank constraint

Equation (5) described how to obtain a quadratic form $J\in\mathcal{R}(\lambda)$ from an $m$-periodic packing $\Xi$ of packing radius $\rho(\Xi)\geq\tfrac{1}{2}\lambda^{1/2}$. However, the reverse operation is not always possible. Clearly, $\mathrm{rank}\,J=d$ is a necessary condition. In fact, it is also sufficient, since $J\in\mathcal{R}(\lambda)$ for $\lambda>0$ implies that $Q$ is positive definite, and therefore $\mathbf{a}_{1},\ldots,\mathbf{a}_{d}$ can be recovered through, e.g., Cholesky decomposition, and $\mathbf{b}_{1},\ldots,\mathbf{b}_{m-1}$ by solving $R_{ij}=\mathbf{b}_{i}\cdot\mathbf{a}_{j}$.

Let $\mathcal{R}_{0}(\lambda)=\{J\in\mathcal{R}(\lambda):\mathrm{rank}\,J=d\}$. Since $Q$ is positive definite within $\mathcal{R}(\lambda)$, we can write this condition as the vanishing of the Schur complement $S-RQ^{-1}R^{T}=0$. For $m=2$, we also have two simpler equivalent forms for this constraint: $\det J=0$ or $\lambda_{\text{min}}(J)=0$, where $\lambda_{\text{min}}$ denotes the smallest eigenvalue.

2.6 The density objective

We wish to find the maximum of $\delta(\Xi)$ among $m$-periodic sets $\Xi$ of packing radius at least $\rho_{0}$. This is equivalent to the minimization problem

[TABLE]

where $Q(J)$ is the top-left block of $J$ (see (5)). Since the problem is equivalent for different values of $\rho_{0}$, we will often use for convenience the explicit value $\rho_{0}=1$, which corresponds to $\lambda=4\rho_{0}^{2}=4$. While $\mathrm{obj}(J)$, like the objective in Voronoi’s algorithm, is quasiconcave, the nonlinearity of the rank constraints does not allow for a straightforward characterization of the local minima.

A vivid illustration of a new type of local minima that can arise is a family of 9-dimensional 2-periodic packing arrangements known collectively as the fluid diamond packing [10, p. 41]. The lattice $D_{9}=\{\mathbf{n}\in\mathbb{Z}^{n}:n_{1}+n_{2}+\ldots+n_{9}\in 2\mathbb{Z}\}$ has a packing radius $\rho=1/\sqrt{2}$ and is an extreme lattice, so that any nearby lattice has a smaller packing radius or a larger determinant. The deep holes of a lattice are the points that maximize the distance to the nearest lattice point. For $D_{9}$, the deep holes are at $D_{9}+\mathbf{b}_{0}$, where $\mathbf{b}_{0}=(\tfrac{1}{2},\tfrac{1}{2},\ldots,\tfrac{1}{2})$, and are at a distance of $3/2>2\rho$ from the lattice. So, $D_{9}^{+}(\mathbf{b})=D_{9}\cup(D_{9}+\mathbf{b})$ has the same packing radius, but double the density when $\mathbf{b}$ is sufficiently close to $\mathbf{b}_{0}$. As a 2-periodic point set, this is clearly a local maximum of the density since a nearby 2-periodic point must be of the form $\Lambda\cup(\Lambda+\mathbf{b}^{\prime})$, where $\Lambda$ is nearby $D_{9}$, and therefore has a larger determinant or a packing radius that is no larger. Thus we have a 9-dimensional family 2-periodic point sets that are all locally optimal. In fact, the fluid diamond packing achieves the highest known packing density for spheres in 9 dimensions. A similar situation can be constructed with any extreme lattice whose covering radius (the distance of the deep holes from the lattice) is more than twice its packing radius. Another example, therefore can be constructed with the lattice $A_{8}$, but this example is less often discussed, presumably because other packings achieve higher density in that dimension.

This vivid example should discourage us from attempting to transport the methods of the Voronoi algorithm, where a critical result was that all local optima lie on vertices and are therefore isolated. However, we will show that, at least in the case $m=2$, the packing constructed this way are the only possible problem we may encounter. Namely, for $m=2$, a local minimum of $\mathrm{obj}(J)$ in $\mathcal{R}_{0}(\lambda)$ is either an isolated minimum lying on an edge of $\mathcal{R}(\lambda)$ or is the union of two translates of an extreme lattice whose covering radius is more than twice the packing radius.

To start, we define a sufficient condition for local optimality by linearizing all of the constraints, which we call algebraic extremeness, and prove that it is indeed a sufficient condition. Let $J\in\mathcal{R}_{0}(\lambda)$. The set $\{J^{\prime}:\mathrm{rank}\,J^{\prime}=d\}$ is a $[\tfrac{1}{2}(d+m)(d+m-1)-\tfrac{1}{2}m(m-1)]$-dimensional manifold in the neighborhood of $J$, and its tangent space is given by

[TABLE]

where $N(J)$ is the null space of $J$, and $\mathcal{S}[N(J)]$ is the space of quadratic forms over $N(J)$. The linear inequality constraints, encoded in the polyhedron $\mathcal{R}(\lambda)$, give rise in the neighborhood of $J$ to the cone

[TABLE]

Finally, the gradient of the objective $\mathrm{obj}(J)$ is proportional to

[TABLE]

Definition 1**.**

A configuration $J\in\mathcal{R}_{0}(\lambda)$ is called algebraically extreme if $\langle J^{\prime}-J,G_{J}\rangle>0$ for all $J^{\prime}\in C_{J}\cap T_{J}$ except $J^{\prime}=J$.

Theorem 2**.**

Let $J\in\mathcal{R}_{0}(\lambda)$ be algebraically extreme. Then

$J$* is an isolated minimum of $\mathrm{obj}(J^{\prime})$ (11) among $J^{\prime}\in\mathcal{R}_{0}(\lambda)$.* 2. 2.

$J$* lies on a $\tfrac{1}{2}m(m+1)$-dimensional face of $\mathcal{R}(\lambda)$.*

Proof.

The first claim follows as the usual sufficient optimality criterion for a smooth optimization problem with smooth equality and inequality constraints [14].

To prove the second claim, note that if $J$ does not lie on a $\tfrac{1}{2}m(m+1)$-dimensional face of $\mathcal{R}(\lambda)$, then $C_{J}$ contains an affine space passing through $J$ of dimension $\tfrac{1}{2}m(m+1)+1$. Since $T_{J}$ is an affine space of codimension $\tfrac{1}{2}m(m+1)$ passing through $J$, then $C_{J}\cap T_{J}$ must include a line, and it is not possible for $\langle J^{\prime}-J,G_{J}\rangle$ to be positive everywhere on this line, even with $J$ excepted. Therefore, $J$ is not algebraically extreme. ∎

Definition 2**.**

A configuration $J\in\mathcal{R}_{0}(\lambda)$ is called a fluid packing if it sits on a nonconstant continuous curve $J^{\prime}\colon[0,1]\to\mathcal{R}_{0}(\lambda)$ of the form

[TABLE]

such that $J^{\prime}(0)=J$ and $J^{\prime}(t)$ is a local minimum of $\mathrm{obj}$ over $\mathcal{R}_{0}(\lambda)$ for all $0\leq t\leq 1$.

One example of a fluid packing is $D_{9}^{+}(\mathbf{b}_{0})$, which is part of the 9-dimensional fluid diamond packing described above. Specifically, let $A$ be a matrix whose rows form a basis for $D_{9}$, let $\mathbf{b}_{1}$ be a point such that $\|\mathbf{b}_{0}-\mathbf{b}_{1}\|<3/2-\sqrt{2}$, and let

[TABLE]

where $\mathbf{b}(t)=(1-t)\mathbf{b}_{0}+t\mathbf{b}_{1}$.

When $m=2$, every fluid packing is the union of two extreme lattices. The continuous curve in Definition 2 when $m=2$ must be of the form (13). Therefore, $J^{\prime}(t)(\mathbf{k})$ is a nonconstant quadratic function of $t$ for every $\mathbf{k}=(\mathbf{n},\pm 1)$ and there is some $0<t^{\prime}<1$ such that $J^{\prime}(t)(\mathbf{k})>\lambda$ for every such $\mathbf{k}$. Because $J^{\prime}(t^{\prime})$ is a local minimum of $\mathrm{obj}$ over $\mathcal{R}_{0}(\lambda)$, $Q$ is a local minimum of $\det Q^{\prime\prime}$ over $\{Q^{\prime\prime}:\mathbf{n}^{T}Q^{\prime\prime}\mathbf{n}\geq\lambda\text{ for all }\mathbf{n}\in\mathbf{Z}^{d}\}$.

Theorem 3**.**

Let $m=2$, and let $J\in\mathcal{R}_{0}(\lambda)$. If $J$ is a local minimum of $\mathrm{obj}(J^{\prime})$ among $J^{\prime}\in\mathcal{R}_{0}(\lambda)$ then one of the following is true:

$J$* is algebraically extreme.* 2. 2.

$J$* is a fluid packing.*

Proof.

Without loss of generality, let $\lambda=4$. Let $\{\mathbf{k}_{1},\ldots,\mathbf{k}_{s}\}=\{\mathbf{k}:J(\mathbf{k})=4\}=J^{-1}(4)$. The quadratic form $J$ is a local optimum of the problem

[TABLE]

As a consequence, it satisfies the Karush-Kuhn-Tucker condition. Namely, there exist $u_{\mathbf{k}}\geq 0$ and $v$ such that $\nabla\mathrm{obj}+\sum_{i=1}^{s}u_{i}\nabla g_{i}+v\nabla h=0$. We consider the derivative with respect to $J^{\prime}_{d+1,d+1}$, which we denote $\partial_{d+1,d+1}$. Since $\partial_{d+1,d+1}\mathrm{obj}=0$, $\partial_{d+1,d+1}g_{i}\leq 0$, and $\partial_{d+1,d+1}h>0$, we have that $v\geq 0$. The Lagrangian corresponding to these KKT coefficients is $L=\mathrm{obj}+\sum u_{i}g_{i}+vh$. Let $C_{0}$ be the marginal cone

[TABLE]

where $I^{+}=\{i:u_{i}>0\}$ and $I^{0}=\{i:u_{i}=0\}$. Suppose that $J$ is not algebraically extreme, then the marginal cone is nontrivial. Consider the Hessian of the Lagrangian $H=\mathrm{hess}\,L=\mathrm{hess}\,\mathrm{obj}+v\mathrm{hess}\,h$, where, e.g., $\langle\tilde{J},(\mathrm{hess}\,\mathrm{obj})\tilde{J}\rangle=(d^{2}/dt^{2})\mathrm{obj}(J+t\tilde{J})$. A necessary condition for $J$ to be a local optimum of (14) is that $\langle\tilde{J},H\tilde{J}\rangle$ is nonnegative for all $\tilde{J}\in C_{0}$ [14, p. 216]. However, both $\mathrm{hess}\,\mathrm{obj}$ and $\mathrm{hess}\,h$ are negative semidefinite (for $\mathrm{obj}$, this is the quasiconcavity of the determinant function; for $h$, this is a well-known result of second-order eigenvalue perturbation theory). It follows that $C_{0}$ must lie in the null-space of $\mathrm{hess}\,\mathrm{obj}$, that is, the $Q$ component of any vector in the marginal cone is zero.

Consider a nonzero element $\tilde{J}\in C_{0}$, with blocks $\tilde{Q}=0$, $\tilde{R}$, and $\tilde{S}$, and let $J^{\prime}(t)$ be a one-parameter family such that $Q^{\prime}(t)=Q$, $R^{\prime}(t)=R+t\tilde{R}$, and $S^{\prime}(t)=S+t\tilde{S}+t^{2}\tilde{R}^{T}Q^{-1}\tilde{R}$. By the Schur complement condition, we see that $\mathrm{rank}\,J^{\prime}(t)=d$ for all $t$. Also, $g_{i}(J^{\prime}(t))=4-\mathbf{k}_{i}^{T}J^{\prime}(t)\mathbf{k}_{i}=\langle\nabla g_{i},\tilde{J}\rangle-t^{2}\mathbf{l}^{T}\tilde{R}^{T}Q^{-1}\tilde{R}\mathbf{l}$ is nonpositive for $t>0$ by the positive definiteness of $Q^{-1}$. Therefore, for some $\epsilon$, $J^{\prime}(t)\in\mathcal{R}_{0}(4)$ for all $0\leq t<\epsilon$. Since $\mathrm{obj}(J^{\prime}(t))=\mathrm{obj}(J)$ and $J$ is a local minimum of $\mathrm{obj}$ over $\mathcal{R}_{0}(4)$, so must $J^{\prime}(t)$ for all sufficiently small $t$. Therefore $J$ is a fluid packing. ∎

We restricted the statement of Theorem 3 to the case $m=2$, where its analysis is more tractable thanks to the simplified form of the rank constraint (see Section 2.5). However, we conjecture that it is also true for $m>2$, but the analysis of optimality conditions becomes more complicated. Theorems 2 and 3 suggest a straightforward generalization of Voronoi’s algorithm to 2-periodic sets, which we elaborate in the next section.

3 The generalized Voronoi algorithm

3.1 Outline

Our algorithm seeks to enumerate all the locally optimal 2-periodic sets in $d$ dimensions. As we proved in Theorem 3, those are either fluid packings or algebraically extreme packings. The problem of enumerating the fluid packings is rather straightforward for $m=2$, since the two component lattices must themselves be extreme.

The enumeration of the algebraically extreme 2-periodic sets of packing radius $1$ consists of two steps that can be conceptually thought of as occurring one after the other, but in our implementation are actually interleaved. The first step is to enumerate all the vertices and edges of $\mathcal{R}(4)$ up to equivalence under the action of $\Gamma$. The second step is to take each edge, represented in the form $\{J+tJ^{\prime}:0\leq t\leq 1\}$ or $\{J+tJ^{\prime}:t\geq 0\}$, solve for all $t$ such $J+tJ^{\prime}\in\mathcal{R}_{0}(4)$, and check whether $J+tJ^{\prime}$ is algebraically extreme. We elaborate on these steps below.

3.2 Enumerating connected components of fluid packings

As noted in Section 2.6, a 2-periodic fluid packing is the union of two translates of an extreme lattice. The relative translation must be a vector that is at a distance from the nearest lattice point of at least twice the packing radius. The local maxima of the distance from the nearest lattice point are called the holes of a lattice, and are well-studied objects in the theory of lattices [6]. Namely, they are the circumcenters of the Delaunay cells of the lattice (when those lie inside the Delaunay cell). Every fluid packing is in the same connected component as the fluid packing obtained by taking the union of an extreme lattice and a translation of it by a hole of that lattice. The extreme lattices of a given dimension can be enumerated using Voronoi’s algorithm, and their holes can be determined, for example, by finding the vertices of the Voronoi cell of the lattice [6]. This enumeration might be redundant if two holes of different lattice orbits are connected by a path whose distance to the lattice does not go below twice the packing radius, but it is guaranteed to give us at least one representative of every connected component of fluid packings.

3.3 The shortest vector problem and related problems

Throughout the algorithm, we need to solve a problem analogous to the problem known as the shortest vector problem (SVP) in the context of lattices [15]. In our context, given a quadratic form $J:\mathbb{R}^{d+m-1}\to\mathbb{R}$ with positive definite upper-left $d\times d$ block $Q$, we wish to find its minimum over nonzero vectors in $\mathcal{M}=\mathbb{Z}^{d}\times(E-E)$. We also, in some cases, want to enumerated the vectors attaining this minimum, or more generally enumerate all vectors that attain a value below some threshold.

For $\mathbf{k}=(\mathbf{n},\mathbf{l})\in\mathcal{M}$ with fixed $\mathbf{l}$, $J(\mathbf{k})$ is an inhomogeneous quadratic form of $\mathbf{n}$ (see Eq. (8)). Finding the minimum of an inhomogeneous quadratic form over the integer vectors is a problem known as the closest vector problem (CVP) for lattices and is closely related to the SVP [15]. Since there are only finitely many possible values of $\mathbf{l}$, the problem reduces to a finite number of instances of the CVP problem for the underlying lattice [10].

As in the case of lattices, it might be expected that reduction of the quadratic form (using elements of $\Gamma$) would significantly decrease the average time to compute the SVP. To reduce the form, we first reduce the underlying lattice (the $Q$ block) by applying an operation $T$ of the form (9), where $V=0$ and $W=1$. We then perform size reduction of the translation vectors, by applying the appropriate operation $T$, where $U=1$ and $W=1$, so that $Q^{-1}R^{T}$ has entries in $[-\tfrac{1}{2},\tfrac{1}{2}]$.

3.4 Enumeration of vertices

The algorithm to enumerate the vertices of $\mathcal{R}(4)$ is similar to the one accomplishing the analogous task in Voronoi’s algorithm. We start with a known vertex of $\mathcal{R}(4)$, denoted $J_{1}$. We compute the extreme rays of its cone $C_{J_{1}}$ [16]. For each such ray, $\{J_{1}+tJ^{\prime}:t\geq 0\}$ there are two possibilities: either it lies entirely in $\mathcal{R}(4)$ or there is some $t>0$ such that $J_{1}+tJ^{\prime}$ is another vertex of $\mathcal{R}(4)$, which we say is contiguous to $J_{1}$. The first possibility does not exist in the original algorithm for lattices, but does occur in our case, and we must check for this possibility. We discuss this problem in Section 3.5. If it is determined that the edge is not unbounded, we need to determine $t^{*}=\max\{t:J_{1}+tJ^{\prime}\in\mathcal{R}(4)\}$. This can be done by using bisection to narrow down the possible range of values for $t^{*}$ the point where we can identify the inequalities that become violated for $t>t^{*}$. We can then solve for $t^{*}$ by setting them to equalities. This would be a variation of the Algorithm $2$, $3.1.5$ from Ref. [10]. As $t^{*}$ is known to be rational, we find that in the low dimensional cases we tackle in this work, a practical and efficient method is to check whether $J_{1}+(a/b)J^{\prime}$ is a vertex of $\mathcal{R}(4)$ for increasing natural values of $a$ and $b$. However, as the number of dimensions increases we expect this method to stop being practical, as it does in the case of the lattice Voronoi theory.

For each contiguous vertex, we check if it is equivalent to $J_{1}$ (more on this step in Section 3.6), and if not, we add it to a queue of vertices to be processed and to our partial enumeration of vertices. At each subsequent step of the algorithm, we remove a vertex from the queue, compute its contiguous vertices, check them for equivalence against all vertices in our partial enumeration, and add the ones that are not equivalent to previously enumerated vertices to the queue and to the partial enumeration. The enumeration is complete when the queue is empty.

As we mentioned in Section 2.4, we conjecture that $\mathcal{R}(4)$ has finitely many inequivalent vertices for all $d$, and therefore that the enumeration always halts. However, even if there are infinitely many vertex orbits, and the enumeration is performed breadth-first, then any orbit will eventually be enumerated by the algorithm.

The vertices of $\mathcal{R}(4)$ and the facets of their cones are all rational since they are determined by equations of the form $\mathbf{k}^{T}J\mathbf{k}=4$ with integer coefficients. Therefore, all these calculations can be performed using exact arithmetic.

One way to obtain a starting vertex to initialize the algorithm is as follows: consider $A_{d}$ (or $D_{d}$, or any extreme, nonfluid $d$-dimensional lattice). It can be represented as 2-periodic set by taking a sublattice of index 2 instead of the primitive lattice. This 2-periodic set is necessarily algebraically extreme [13], and so it lies on an edge of $\mathcal{R}(4)$, which must terminate at a vertex in at least one direction.

3.5 Detection of unbounded edges

Given a vertex $J$ of $\mathcal{R}(4)$ and an extreme ray of $C_{J}$, of the form $\{J+tJ^{\prime}:t\geq 0\}$, we want to determine if the ray lies entirely in $\mathcal{R}(4)$. This is the case if and only if $J^{\prime}(\mathbf{k})\geq 0$ for all $\mathbf{k}\in\mathcal{M}$. So, we have a problem very similar to the SVP above, except that we may not assume that $Q^{\prime}$ is positive definite. If $Q^{\prime}$ is not even positive semidefinite, then there is some $\mathbf{k}=(\mathbf{n},0)$ such that $J^{\prime}(\mathbf{k})=Q^{\prime}(\mathbf{n})<0$, and the ray is not unbounded. So the only remaining case is when $Q^{\prime}$ is positive semidefinite and has nontrivial null space, $N(Q^{\prime})$.

We break this remaining case into two cases. First, consider the case that there exist $\mathbf{l}\in(E-E)$ and $\mathbf{u}\in N(Q^{\prime})$ with $c=\mathbf{u}^{T}(R^{\prime})^{T}\mathbf{l}>0$. Let $\mathbf{v}(t)=(-tc\mathbf{u},\mathbf{l})$, let $[\mathbf{v}](t)\in\mathcal{M}$ be the closest integer vector to $\mathbf{v}(t)$, and let $(\mathbf{e},0)=\mathbf{v}(t)-[\mathbf{v}](t)$ be the remainder. Then

[TABLE]

All terms except $-2tc^{2}$ are either zero or bounded as $t\to\infty$ and so $J^{\prime}([\mathbf{v}](t))<0$ for large enough $t$ and the ray is not unbounded.

The final case is that $\mathbf{u}^{T}(R^{\prime})^{T}\mathbf{l}=0$ for all $\mathbf{l}\in(E-E)$ and $\mathbf{u}\in N(Q^{\prime})$. In this case, the problem reduces to the first case, where $Q^{\prime}$ is positive definite, albeit in a smaller dimension: Since $Q^{\prime}$ is a rational matrix, there is a unimodular transformation $U\in GL_{d}(\mathbb{Z})$ such that $Q^{\prime}\circ U$ has $\mathrm{span}\,\{\mathbf{e}_{r+1},\ldots,\mathbf{e}_{d}\}$ as its null space and is positive definite on $\mathrm{span}\,\{\mathbf{e}_{1},\ldots,\mathbf{e}_{r}\}$. We can find this transformation by computing the Hermite normal form or the Smith normal form of $Q^{\prime}$ scaled to an integer matrix. Let $\mathbf{n}=\mathbf{n}^{\prime}+\mathbf{n}^{\prime\prime}$ with $\mathbf{n}^{\prime\prime}\in N(Q^{\prime}\circ U)$ and $\mathbf{n}^{\prime}\in N(Q^{\prime}\circ U)^{\perp}$. Then

[TABLE]

Since, by assumption of this case, $(\mathbf{n}^{\prime\prime})^{T}U^{T}(R^{\prime})^{T}\mathbf{l}=0$, the value of $J^{\prime}(U(\mathbf{n}),\mathbf{l})$ depends only on $\mathbf{n}^{\prime}$ and $\mathbf{l}$, and we may find its minimum using our SVP method.

3.6 Equivalence checking

Given two quadratic forms $J_{1}$ and $J_{2}$, we wish to determine if $J_{2}=J_{1}\circ T=T^{T}J_{1}T$ for some $T\in\Gamma$. Let us denote this relation as $J_{1}\sim J_{2}$. This problem may apply to vertices of $\mathcal{R}(4)$, as part of the algorithm for enumerating vertices, but we may also apply it to any pair of quadratic forms in $\mathcal{R}(4)$ that are not necessarily vertices. Let us first prove some useful results.

Definition 3**.**

A set $M\subset\mathcal{M}$ is perfect if $J(\mathbf{k})=J^{\prime}(\mathbf{k})$ for all $\mathbf{k}\in M$ implies $J=J^{\prime}$.

A set $M$ is perfect if and only if the set $\{\mathbf{k}\mathbf{k}^{T}:\mathbf{k}\in M\}$ spans the space of symmetric $(d+m-1)\times(d+m-1)$ matrices. Denote by $J^{-1}(A)\subset\mathcal{M}$ the set of vectors in $\mathcal{M}$ attaining values in the set $A$. If $A=\{4\}$ and $J$ is a vertex of $\mathcal{R}(4)$, then $J^{-1}(A)$ is perfect. A direct consequence of the definition of a perfect set is the following lemma:

Lemma 2**.**

Let $A$ be a set of real values, and let $M_{1}=J_{1}^{-1}(A)$ and $M_{2}=J_{2}^{-1}(A)$ be the set of vectors in $\mathcal{M}$ that achieve these values. If $M_{1}$ and $M_{2}$ are perfect, then the following are equivalent:

$J_{2}=J_{1}\circ T$* for some $T\in\Gamma$.* 2. 2.

$T(M_{2})=M_{1}$, and $J_{2}|_{M_{2}}=(J_{1}\circ T)|_{M_{2}}$ for some $T\in\Gamma$.

In particular a $T$ that satisfies one condition also satisfies the other.

Proof.

First, suppose $J_{2}=J_{1}\circ T$, and let $\mathbf{k}\in M_{2}$. Since $J_{2}(\mathbf{k})=J_{1}(T\mathbf{k})\in A$, we have that $T\mathbf{k}\in M_{1}$. Similarly, if $\mathbf{k}\in M_{1}$, then $T^{-1}\mathbf{k}\in M_{2}$. So $T(M_{2})=M_{1}$, and clearly $J_{2}|_{M_{2}}=(J_{1}\circ T)_{M_{2}}$ follows a fortiori from the unrestricted equality.

The other direction follows immediately from the definition of a perfect set. ∎

Therefore, a simple algorithm to check for equivalence is as follows: first, construct a set $A$ such that $M_{1}=J_{1}^{-1}(A)$ is perfect. When $J_{1}$ is a vertex of $\mathcal{R}(4)$, the set $A=\{4\}$ suffices. Otherwise, we find the smallest $a>4$ such that $A=[4,a]$ suffices. Next, compute $M_{2}=J_{2}^{-1}(A)$. If $M_{2}$ is not perfect, $|M_{1}|\neq|M_{2}|$, or $J_{2}$ does not take values in $A$ over $M_{2}$ with the same frequency as $J_{1}$ does over $M_{1}$, then $J_{1}\not\sim J_{2}$. We give labels $1,\ldots,s$ to the elements of $M_{1}=\{\mathbf{k}_{1},\ldots,\mathbf{k}_{s}\}$ and $M_{2}=\{\mathbf{k}^{\prime}_{1},\ldots,\mathbf{k}^{\prime}_{s}\}$, such that $\mathbf{k}_{1},\ldots,\mathbf{k}_{d+m-1}$ are linearly independent (there must be such a linearly independent subset for $M_{1}$ to be perfect). We now try to construct an injective map $\sigma:\{1,\ldots,d+m-1\}\to\{1,\ldots,s\}$ such that

[TABLE]

We can do this by a backtracking search, constructing $\sigma$ on $1,\ldots,n<d+m-1$ for increasing $n$, and backtracking when no possible assignment of $\sigma(n+1)$ satisfies (18) for $1\leq i,j\leq n+1$. For each such complete map produced by the backtracking search, there is a unique linear map $T^{\prime}$ such that $T^{\prime}\mathbf{k}_{i}=\mathbf{k}^{\prime}_{\sigma(i)}$ for $1\leq i\leq d+m-1$. If $T^{\prime}\in\Gamma$ and $T^{\prime}(M_{1})=M_{2}$, we are done and $J_{1}\sim J_{2}$. Otherwise, we continue with the backtracking search. If the backtracking search concludes without finding any equivalence, then $J_{1}\not\sim J_{2}$.

3.7 Enumeration of algebraically extreme forms

Given an edge of $\mathcal{R}(4)$ of the form $J+tJ^{\prime}$, where $0\leq t\leq 1$ or $t\geq 0$, we wish to identify the points of this edge that lie in $\mathcal{R}_{0}(4)$. As we are limiting ourselves to the case $m=2$, we simply need to solve $\det(J+tJ^{\prime})=0$ for $t$, which is a univariate polynomial equation. This equation can be solved for $t$ as a generalized eigenvalue problem $J\mathbf{x}=-tJ^{\prime}\mathbf{x}$. It is possible that the entire edge lies in $\mathcal{R}_{0}(4)$, as happens for one edge in our enumeration for $d=5$. In that case, either $Q$ is constant over the edge, and any locally optimal packing on the edge is necessarily a fluid packing, or $Q$ is not constant, and therefore neither is $\mathrm{obj}$. So, only points that minimize $\mathrm{obj}$ over the edge — necessarily endpoints by quasiconcavity — can be locally optimal packings.

For each possible candidate configuration (each root of $\det(J+tJ^{\prime})$ if it is not identically zero, or the endpoints if it is), we test to see if it is algebraically extreme. If $J$ is algebraically extreme, we can certify that it is by considering the dual problem: $J$ is algebraically extreme if and only if the cone $\{\sum_{\mathbf{k}\in J^{-1}(4)}\eta_{\mathbf{k}}\mathbf{k}\mathbf{k}^{T}+\mathbf{x}\mathbf{x}^{T}:\eta_{\mathbf{k}}>0,\mathbf{x}\in N(J)\}$ is full-dimensional and includes $G_{J}$. If $J$ is not algebraically extreme, we can certify that by the direct problem of finding $J^{\prime}$ that gives a counterexample for the definition: $J^{\prime}\neq J$, $\langle J^{\prime}-J,G_{J}\rangle\leq 0$, and $J^{\prime}\in C_{J}\cap T_{J}$. Since $t$ comes from the root of a univariate polynomial, it is not, in general, rational.

In our calculation, we used floating point representation for the candidate algebraically extreme configurations and for the implementation of the test for whether they are algebraically extreme.

4 Numerical results

We use an implementation of the algorithm described in the previous section to fully enumerate the vertices of the Ryshkov-like polyhedron $\mathcal{R}(4)$ and the algebraically extreme $m$-periodic packings for $m=2$ and $d=3,4,5$. The results of the enumeration are summarized in Table 1. The extreme lattices in these dimensions are the classical root lattices $A_{d}$, $D_{d}$, and $E_{5}$, whose covering radii are known [6], and do not exceed twice the packing radius. Since the existence of a 2-periodic fluid packing requires the existence of an extreme lattice with a hole that is at a distance of more than twice the packing radius from the nearest lattice point (see Section 3.2), there are no fluid packings for $d=3,4,5$. Our enumeration of the algebraically extreme 2-periodic algebraically extreme packings exhausts all locally optimal 2-periodic packings.

Our attempted enumeration for $d=6$ did not terminate after over a month of running. The main bottleneck appears to be the enumeration of extreme rays of $C_{J}$ for vertices $J$ with large number of minimal vectors $J^{-1}(4)$ (see Sec. 4.4). This is similar to the main bottleneck in the enumeration of perfect lattices in $d=8$, where only by exploiting the symmetries of these high-kissing-number lattices, was the full enumeration made tractable [17]. We are hopeful that a similar approach could be used for $m=2$ to make full enumeration in higher dimensions than $d=5$ tractable, but we do not attempt to implement it in this work.

We use a mix of floating-point arithmetic and exact arithmetic in various parts of the algorithm. Therefore, we do not claim the results in this section as rigorous results, but rather as numerical results. We relied on floating-point arithmetic for the shortest vector problem, giving sufficient tolerance so that no candidate shortest vectors are rejected and then using exact arithmetic to verify that they are in fact all of the same length. Determining the rays of vertex cones and finding contiguous vertices along those rays was done with exact rational arithmetic. The final step of finding candidate optima along the edges of the polyhedron and checking whether these candidates are algebraically extreme was done using floating-point arithmetic. After obtaining the set of algebraic extreme packings, we were able to use a computer algebra system to obtain exact algebraic expressions for their coordinates.

4.1 d=3

The enumeration in $d=3$ gave 4 inequivalent vertices and 3 inequivalent algebraically extreme 2-periodic sets. All the algebraically extreme 2-periodic sets in 3 dimensions of packing radius $1$ have the same density $\delta=1/(2\sqrt{2})\approx 0.354$, and they are represented by the following quadratic forms:

[TABLE]

[TABLE]

The forms $J_{3,1}$ and $J_{3,2}$ are two inequivalent representations of the fcc lattice as a 2-periodic set. $J_{3,1}$ is a stacking of square layers, whereas $J_{3,2}$ is a stacking of triangular layers (see Figure 1). The forms are not equivalent under the action of $\Gamma$ because the corresponding two sublattices of the fcc lattices are not equivalent under the symmetries of the fcc lattice. The form $J_{3,3}$ represents the hexagonal-close-packing 2-periodic arrangement, which is not a lattice. Note that a form $J$ is the representation of a lattice as a 2-periodic set if and only if $Q^{-1}R^{T}\in(\tfrac{1}{2}\mathbb{Z})^{d}$.

All the algebraically extreme 2-periodic sets have the same kissing number, 12, but as quadratic forms they have different number of minimal vectors: $|J_{3,1}^{-1}(4)|=20$, $|J_{3,2}^{-1}(4)|=|J_{3,3}^{-1}(4)|=18$. This difference occurs because contacts between spheres in the same orbit contribute the same minimal vector as the analogous contact in a different orbit. Denote by $|\cdot|_{*}$ a counting measure that gives weight $m$ to vectors of the form $\mathbf{k}=(\mathbf{n},0)$. Then the average kissing number is $\kappa=|J^{-1}(\mathrm{min}\,J)|_{*}/m$.

The vertices of the Ryshkov-like polyhedron are also interesting to consider. They are

[TABLE]

[TABLE]

They have full rank, so they do not correspond to 2-periodic sets in 3 dimensions. However, because they are positive definite, they correspond to 4-dimensional sets that are periodic (with two orbits) under a 3-dimensional lattice. When they are projected to the space spanned by the lattice, they can be interpreted as binary packings of non-additive spheres, where the two sphere species have equal self-radius, but smaller radius when interacting with each other. The first is a simple cubic lattice of one species with its body-center holes filled with spheres of the other species. The second is a hexagonal lattice with one of the two inequivalent triangular prism-shaped holes in each unit cell filled. The third and fourth are the fcc lattice with its octahedral or (one of its) tetrahedral holes filled, giving, respectively, a simple cubic lattice with alternating species (the NaCl crystal structure) and the diamond crystal structure (see Figure 2).

Finally, we point out as an example of an unbounded edge, the edge connecting $J_{3,3}$ and $J_{3,2v}$.

4.2 d=4

In 4 dimensions, we obtain 7 inequivalent algebraically extreme 2-periodic sets. Six of those lie in the relative interior of edges of the Ryshkov-like polyhedron, and one is a vertex. Two algebraically extreme sets achieve the maximum density, $\delta=1/8=0.125$, and both are representations of the lattice $D_{4}$ as a 2-periodic set. They are represented by the forms,

[TABLE]

The form $J_{4,1}$ is a vertex of the Ryshkov-like polyhedron. Both have the kissing number of $D_{4}$, $\kappa=24$, as the average kissing number. The next-highest density achieved is $\delta=2/\sqrt{144+64\sqrt{5}}\approx 0.118$, which is achieved by a single algebraically extreme 2-periodic set:

[TABLE]

where $\tau=(1+\sqrt{5})/2$ is the golden ratio. Its average kissing number is $\kappa=22$. Finally, there are four algebraically extreme 2-periodic sets that achieve the same density as the $A_{4}$ lattice, $\delta=1/(4\sqrt{5})\approx 0.112$. Two of them are in fact representations of the $A_{4}$ lattice, but the other two are not also lattices. All four have the same kissing number as the lattice, $\kappa=20$. Including $J_{4,1}$, there are 10 inequivalent vertices of the Ryshkov-like polyhedron.

4.3 d=5

In 5 dimensions, we obtain 29 inequivalent algebraically extreme 2-periodic sets (see Table 4.3 for a tabulation of their invariants). Five inequivalent 2-periodic sets achieve the maximum density $\delta=1/(8\sqrt{2})\approx 0.0884$, three of which are representations of the lattice $D_{5}$ as a 2-periodic set, and two are nonlattices that achieve the same density. The ones that represent the lattice $D_{5}$ are given by the forms,

[TABLE]

[TABLE]

[TABLE]

where $J_{5,1}$ is also a vertex of the Ryshkov-like polyhedron (the only vertex that is also in $\mathcal{R}_{0}$). The two that are not lattices but achieve the same density are

[TABLE]

[TABLE]

All have the same kissing number $\kappa=40$.

The next highest density, $\delta=(\sqrt{5-2\sqrt{6}})/4\approx 0.0795$, is achieved by three algebraically extreme 2-periodic sets, all of which are not also lattices:

[TABLE]

[TABLE]

[TABLE]

The first two of these have a kissing number of $\kappa=35$, and the third has $\kappa=34$. The third largest density is that achieved by one of the three extreme lattices in 5 dimensions, $\delta=1/(9\sqrt{2})\approx 0.0786$, and is achieved by five 2-periodic sets, three of which are representations of the lattice, and all have the same kissing number as the lattice. The third extreme lattice, $A_{5}$, has the lowest density of all extreme lattices, $\delta=1/(8\sqrt{3})\approx 0.722$, and this density is achieved by seven 2-periodic sets, three of which represent the lattice, and all having the same kissing number as the lattice. This is also the lowest density achieve by any algebraically extreme 2-periodic set, although there is a form of $\mathcal{R}_{0}(4)$ on an edge of $\mathcal{R}(4)$ that has lower density, but fails to be algebraically extreme.

Including $J_{5,1}$, the Ryshkov-like polyhedron has 34 vertices. Two of the vertices are not positive semidefinite, and this is the lowest dimension where such vertices occur. When interpreted as nonadditive binary sphere packings, these packings would have a nonself-radius larger than the self-radius. In this sense they are similar to the fluid packings (the distance between orbits is larger than the distance within each orbit), but the underlying lattices are not extreme. Another interesting phenomenon that first occurs in 5 dimensions is an edge of $\mathcal{R}(4)$ that is completely contained in $\mathcal{R}_{0}(4)$. Such an edge is the one connecting $J_{5,1}$ to $J_{5,1}\circ T$, where

[TABLE]

as well as all the equivalent edges. Along the edge $J_{5,1}+t(J_{5,1}\circ T-J_{5,1})$, the objective determinant is $\mathrm{obj}(t)=512(1+t-t^{2})$, so any internal point cannot be a locally optimal 2-periodic set.