    Density of Binary Disc Packings:Lower and Upper Bounds

We provide, for any r∈ (0,1), lower and upper bounds on the maximal density of a packing in the Euclidean plane of discs of radius 1 and r. The lower bounds are mostly folk, but the upper bounds improve the best previously known ones for any r∈[0.11,0.74]. For many values of r, this gives a fairly good idea of the exact maximum density. In particular, we get new intervals for r which does not allow any packing more dense that the hexagonal packing of equal discs.

Authors

07/15/2018

A new lower bound for classic online bin packing

We improve the lower bound on the asymptotic competitive ratio of any on...
01/01/2020

A new upper bound for spherical codes

We introduce a new linear programming method for bounding the maximum nu...
07/12/2021

Bounds for Multiple Packing and List-Decoding Error Exponents

We revisit the problem of high-dimensional multiple packing in Euclidean...
07/20/2021

Prior-Free Clock Auctions for Bidders with Interdependent Values

We study the problem of selling a good to a group of bidders with interd...
01/31/2019

On two-fold packings of radius-1 balls in Hamming graphs

A λ-fold r-packing in a Hamming metric space is a code C such that the r...
12/03/2017

Exact upper and lower bounds on the misclassification probability

Exact lower and upper bounds on the best possible misclassification prob...
06/30/2019

Upper bounds on the graph minor theorem

Lower bounds on the proof-theoretic strength of the graph minor theorem ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

A disc packing (or circle packing) is a set of interior-disjoint discs in the Euclidean plane. Its density is the proportion of the plane covered by the discs:

 δ:=limsupk→∞area of the square [−k,k]2 % covered by discsarea of the square [−k,k]2.

If all the discs have the same radius, it has been proven by Tóth [FT43] (see also [CW10]) that the density is at most

 δ1:=π2√3≈0.9069,

reached for the so-called hexagonal compact packing, where discs are centered on a triangular grid (of size twice the disc radius).

What about more sizes of discs? In particular, what about the maximal density of binary disc packings, that is, packings with discs of two different sizes? Up to scaling, one can always assume that the largest disc has radius and the smallest one has radius . We then denote by the maximal density of packings by discs of radius and . Can we find the exact value of for each ?

The first decisive step was done by Blind [Bli69], who proved that for where

 rB:= ⎷7tanπ7−6tanπ66tanπ6−5tanπ5≈0.74299.

In other words, binary disc packings cannot achieve higher densities than packings of equal discs when the disc sizes are too close. This yields the exact maximal density over the whole interval . Besides this, to the best of our knowledge, there are only specific ratios for which the exact value of is known. They are the ratios which allow triangulated binary packings, that is, packings with two sizes of discs whose contact graphs are triangulated (the vertices of the contact graph of a packing are the disc center, and the edges connect the centers of tangent discs). These ratios are algebraic numbers which have been characterized in [Ken06]. The maximal density for each of them is proven in [BF21] (see also [Hep00, Hep03, Ken04b] for the first cases).

If we cannot obtain an exact value for , can we find lower and upper bounds? In particular, for which do we have , that is, which ratios allow binary disc packings which achieve higher densities than packings of equal discs? This is a question that is of particular interest in materials science because the maximization of density seems to be a criterion for the formation of materials (due to attractive forces of the van der Waals type). A ratio such that thus suggests the possibility to obtain new materials by combining atoms or nanoparticles within this ratio. A striking example is given by the binary assemblies of nanoparticles whose synthesis is explained in [PDKM15]: they all correspond faithfully to one of the triangulated binary packings that have been proven to maximize the density.

To obtain a lower bound, it suffices to find a “good” packing. For example, if is small enough, namely , we can simply insert a small disc in each hole of a hexagonal compact packing of large discs to get . More interesting, we can get lower bounds over a whole interval of ratios by continuously modifying a “good” packing so that the density decreases as little as possible. A particularly efficient way to proceed is the so-called “flipping and flowing method”, first used by Tóth [FT64], see also [CP19, CG20]. In Section 2 we detail lower bounds obtained in this way. These lower bounds seem to be mostly “folk” (see, e.g., [FJFS20, Ken04a]), but we provide here a SageMath worksheet [Dev19] (the file lower_bounds.sage in supplementary materials) which give explicitly the transformations as well as the corresponding densities. They are depicted in Fig. 1 (green curve). As a corollary, we get:

Corollary 1

One has for any in where , , , and are algebraic number whose minimal polynomials are in the file lower_bounds.sage. Figure 1: Lower bound (green curve) and upper bound (red curve) on δ(r). The blue dots and lines indicate where lower and upper bounds coincide (the dots correspond to the 9 triangulated packings, the lines to the intervals Ii’s in Cor. 2). The best previous upper bound is indicated by the black dotted curve. See also zooms in Fig. 10, 11, 12 and 13.

Finding upper bounds is more challenging. For a given , it indeed amouts to proving that among the uncountably many packings of discs of size and , none has density larger than the claimed upper bound. A first milestone was the proof in [Flo60] that is less than the density inside a triangle with mutually tangent discs of size , and centered on its vertices,see Fig. 1, bottom center. This upper bound was later on enhanced for in [Bli69], by proving that is less than the density inside the union of a regular heptagon and a regular pentagon, with the heptagon (resp. pentagon) being circumscribed to a large disc (resp. small disc), see Fig. 1, bottom right (the above mentionned constant is the value of for which this bound reaches ). In Section 3, we explain how we obtain upper bounds by enhancing the computer-aided method used in [BF21] to prove the maximal density of the triangulated packings. In a nutshell, the interval of the possible disc size ratio is divided into sufficiently many small intervals on which a program (the C++ code is provided in the supplementary materials) searches by dichotomy an upper bound on the maximum density. We have tried to make this paper readable as independently of [BF21] as possible by recalling the main lines of the method used in [BF21] (Subsec. 3.1). It nevertheless relies on many technical details of [BF21] that are omitted here and is therefore not self-contained. The obtained upper bounds improve the previous ones for any but the values in this interval which allow a triangulated binary packing. They are depicted in Fig. 1 (red curve). As a corollary, we get:

Corollary 2

One has for any in

 [0.4445,0.4532]I1∪[0.4917,0.5145]I2∪[0.5666,0.6270]I3∪[0.6468,1)I4,

where the endpoints of the ’s are exact numerical values obtained by computer.

Combining the intervals given in Cor. 1 and 2 answers whether a given ratio allows a binary packing more dense than the hexagonal compact packing of equal discs in more than of the cases. We have deliberately called both these results “corollaries” to emphasize the fact that the main result of this paper are the lower and upper bounds depicted in Fig. 1, which unfortunately cannot be stated in the form of a classical theorem. The feasability of closing the gap between the lower and upper bounds, is discussed in Section 4.

2 Lower bounds by flipping and flowing

2.1 Principle

To date, all disc packings that have been proven to maximize density (the binary ones and the hexagonal compact packing of equal discs) are found to have a triangulated contact graph, i.e. maximum number of contacts between discs [FT43, BF21, Fer19, Bli69]. This backs up the rule of thumb that the more contact between discs in a packing, the denser it is. However, a triangulated contact graph is only possible for some very particular sizes of discs. This is where flipping and flowing comes into play. The principle is, starting from a particularly dense packing, to continuously modify the ratio of disc sizes while trying to keep as much contacts between discs as possible, in the hope of decreasing the density as little as possible. The term “flowing” refers to the continuous modification of disc sizes. The term “flipping” refers to the particular (but frequent) case where the transformation connects two packings whose graphs are triangulated and differ by one or more flips, i.e. a diagonal of a quadrilateral is replaced by the other diagonal.

2.2 The flows

Figures 27 describe flows between (or from) especially dense packings. The ratios ’s and ’s are algebraic numbers whose exact values can be found in the file lower_bounds.sage. These packings are all triangulated, except those for in Fig. 7. They are all periodic and described by their fundamental domains. The density is depicted (red curve), with the horizontal axe corresponding to the density of the hexagonal compact packing. This yields the lower bound depicted in Fig. 1.

Remark 1

One checks that when the flow on the right of and the one on the left of cross, for , the density is higher about than (it is important for Cor. 1). The bottom of the green valley between and in Fig. 1 is thus (very slightly) above the horizontal dotted line.

2.3 Interstitials packings

For , a small disc can fit in the holes of a hexagonal compact packing of large discs. To obtain the lower bound depicted in Fig. 1 for , we simply put as much as possible discs of a hexagonal compact packing of small discs inside the holes of a hexagonal compact packing of large discs (the center of each hole is the center of a small disc). This is not optimal (some improvements can be found in [UST04]). At least, this yields over .

2.4 Computation

To compute the density along a flow, we use the fact that the contacts between discs in the packings yield quadratic equations in the coordinates of disc centers and the ratio . This form a polynomial system which has to be of dimension at least to allow to vary and at most to have as much contacts as possible. Solving this system allows to compute the positions of the discs and the density as functions of the ratio . The complete calculations are provided in the file lower_bound.sage. The main tool is the simple function stick((x1,y1,r1),(x2,y2,r2),r3), which takes the coordinates x1,y1 of a disc of radius r1, the coordinates x2,y2 of a disc of radius r2, a real number r3 and returns the coordinates x3,y3 (if they exist) of the disc of radius r3 which is exteriorly tangent to and such that sees on the left of .

Consider, for example, the periodic binary packing whose fundamental domain is depicted on Fig. 8, left. Discs are numbered and we arbitrarily set the positions of the two first discs such that they are tangent:

d0=(0,0,1)
d1=(2,0,1)

We then stick one by one the three following discs:

d2=stick(d0,d1,r)
d3=stick(d2,d1,1)
d4=stick(d2,d3,1)

This allows to compute the density in the fundamental domain:

d=(2+2*r^2)/((d4+(d3-d1))*d1)

This yields for the density along the flow between and the expression

 π(r2+1)(r+1)416(r+2)√r+2r√r.

All the considered cases are similar, except two which are slightly more complicated. The first case is the periodic binary packing whose fundamental domain is depicted on Fig. 8, right. It corresponds, in the previous subsection, to the flow between and . In this case, it is not possible to describe the packing by sticking discs one by one to two previous discs. To get around this problem, we define a multivariate polynomial ring over whose variables are the coordinates of the disc centers, the disc ratio r and the density d:

K.<x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,r,d>=QQ[]

We then add the equations which describe the disc contacts or symmetries of the packing (each polynomial must be equal to zero):

eqs=[y1,
x1-2*x2,
x1+x2-(x3+x4+x5),
y1+y2-(y3+y4+y5),
x3^2+y3^2-(1+r)^2,
(x4-x3)^2+(y4-y3)^2-(r+r)^2,
(x4-x5)^2+(y4-y5)^2-(r+r)^2,
(x3-x5)^2+(y3-y5)^2-(r+r)^2,
(x4-x1)^2+(y4-y1)^2-(1+r)^2,
(x2-x5)^2+(y2-y5)^2-(1+r)^2,
(x1-2*x3)^2+(y3+y3)^2-(r+r)^2,
d*x1*y2-(1+6*r^2)]

We then compute the ideal defined by these equations and eliminate all the variables but r and d:

I=ideal(K,eqs)
J=I.elimination_ideal([x1,y1,x2,y2,x3,y3,x4,y4,x5,y5])
P=QQ[r,d](J.gens())

This yields a polynomial of degree in r and in d. Actually, we can remove a factor which would correspond to a density larger than . We get an irreducible polynomial of degree in d^2 whose coefficients are polynomials in r. This allows to get a closed-form expression for the density along the flow between and :

 π(6r2+1)√47r4+84r3+54r2+12r+3−(7r3+13r2+9r+3)√45r2−6r−3√6(r4+12r2+12r+3).

The second case corresponds, in the previous subsection, to the flow on the left of . It is similar, except that we eventually get a polynomial in degree in d^2: we cannot derive a closed-form expression, but an implicit plot is possible.

Remark 2

The “algebraic” method by ideal elimination can actually be used for all the flows. The more “pedestrian” method of adding the discs one by one leads however to much faster calculations in practice.

3 Upper bounds via localizing potentials

3.1 Checking an upper bound for a fixed ratio

Given a ratio and a candidate upper bound on the maximum density, we want to check whether or not. The strategy used is the one in [BF21], which resembles the one used by Hales to prove the Kepler conjecture, nicely exposed in [Lag02]. Here we will only sketch this strategy and omit all technical details that are not strictly necessary to understand what follows. We refer to [BF21] for the complete description (which is unfortunately necessary to fully understand the C++ code provided in the supplementary materials).

Given a packing, consider a triangulation of its disc center (we used the Fejes-Molnár triangulation [FTM58] for its suitable properties). The density inside the triangles, that is, the proportion of the triangle area covered by the discs of the packing, can greatly vary from one triangle to the other. In particular, it may be larger than the candidate upper bound on the maximum density for any packing. The idea is to prove that whenever a triangle is too dense, there are necessarily neighbor triangles whose densities are low enough so that the density drops below once averaged over all these triangles. More precisely, the density of each triangle is distributed among its three vertices, and it is proven that no matter how a vertex of the triangulation of a packing is surrounded, the average of the densities distributed to this vertex by the triangles to which it belongs is at most . Key parameters appear to be the so-called “base vertex potentials” , , , , and . Each specifies, in a triangle whose vertices are center of mutually tangent discs of radii , and , the amount of density distributed to the center of the disc of radius . Since the density of such a triangle is determined and since there are different such triangles (depending on the radii , , or of the discs centered on vertices), two of the base vertex potentials have to be chosen and the others follow. This choice has a strong influence on whether the maximum density can be proven to be less than . This choice turned out to be relatively simple to make in each of the 9 ratios considered in [BF21]. But here the ratio can be any, so it is necessary to automate this choice in a satisfactory way.

To explain how this choice is made, we need to mention two other parameters that appear in [BF21], namely and . The base vertex potential is indeed defined only in triangles with mutually tangent discs. For other triangles, grows linearly with the angle in with a proportionality factor . When the angle in increases, the triangle tends to be less dense and can “absorb” more density from its neighbors, which encourages taking large values of and . However, the total density that each triangle can absorb cannot is limited, which in turns limits the value that and can take. This can be formalized by a list of inequalities on , and two of the ’s (since the other ones follow as already mentioned). In the program upper_bound.cpp, we consider and . The inequalies are defined in the function add_vertex_positivity_constraints of the program parameters_xy.cpp. This yields a -dimensional parameter polytope (polyhedron P in the function find_delta in upper_bound.cpp) in which a point has to be chosen.

The following automatic choice was developed through trial and error, with the goal being to prove an upper bound as low as possible. It is implemented in the function set_xy_generic in the program upper_bound.cpp. First, it sets the ratio to be equal to the maximal possible value of divided by the maximal possible value of

. This amounts to intersect the parameter polytope with a hyperplane to get a new polytope. We then take either the barycenter of the vertices of this new polytope if

, or a vertex of it which minimizes otherwise. This defines the parameters and . We do not consider the values and it defines because our program deals with rational polytopes whereas can be not rational (actually, it can even be an interval as we shall later see). Instead, we keep only the numerical values of and and then proceed as in [BF21] to compute and (as well as the parameter , not mentioned here) which will decide exactly whether is a suitable upper bound or not.

3.2 Finding an upper bound for an interval of ratios

The previous subsection explained how to check a given upper bound for a given ratio. But the goal of this paper is to find an upper bound as good as possible for any ratio. This yields two problems:

1. we have no candidate upper bound on the maximal density;

2. we have a continuum of ratio to consider.

Assuming is fixed, the first problem can be fixed by dichotomy on the candidate density. Namely, we maintain two variables and , such that the proof that we have an upper bound on the exact maximal density succeeds for but fails for . We start with slightly less than and equal to the Florian upper bound. At each step, we check whether the proof works for and update and accordingly. We stop when is smaller than a fixed precision, namely . We finally output : it is an upper bound on the maximal density (and the best that we get by our method, up to the fixed precision). This is done in the function find_delta in upper_bound.cpp.

To fix the second problem, we can subdivide is many small intervals and rely on interval arithmetic to consider each small interval as a value for . The smaller the intervals, the better the precision, thus the better the upper bound on the maximal density. We thus want to subdivide is interval as small as possible, with the limiting factor being the computation time.

Actually, we found it more relevant to compute upper bound only for regularly spaced discrete values of with maximal precision. This indeed gives a fair idea of which bounds could be obtained because the maximal density is quite regular:

Proposition 1

For in , the maximal density satisfies

 |δ(y)−δ(x)||y−x|≤πy2√3.

The proof of this proposition is given in Appendix A. In particular, the maximal density is -Lipshitz over any interval . This makes the computation much lighter. For example, on our computer, subdividing the interval takes min. for subintervals of length and min. (h.) for subintervals of length . In comparison, discrete values of with a step of takes min. and yields for these discrete points an upper bound on the maximal density which is comparable to the one obtained with intervals of length , in around times less computation time. The regularity allows to extend this upper bound everywhere. Fig. 9 illustrates this point. Figure 9: The upper bound obtained with subintervals of length 0.0001 (resp. 0.00001) in light gray (resp. dark gray). The red line shows the upper bound obtained with discrete values of r, with the minima being reached for the discrete values of r. Since connecting these minima by segments is expected to give a fair idea of the best upper bound that can be obtained by our method, we only represented these minima on the figures Fig. 1, 10, 11, 12 and 13.

Table 1 gives the step and execution time on each connected component of the complement of (the intervals on which the maximal density has been proven to be ). A list of values can be found in the supplementary materials (file trace_upper_bound.txt), see also Fig. 10, 11 and 12 for a zoom of Fig. 1.

3.3 Checking δ(r)=δ1 over I1, I2, I3 and I4

The upper bound obtained in the previous subsection suggest over the intervals , , and defined in Corollary 2. To prove this, it is no longer enough to consider discrete values of r: we really have to use interval arithmetic. Instead of partitionning each in many equally small intervals, we use again the dichotomy to subdivide as little as possible. Namely, we bisect each while the precision does not suffice to conclude. But we no longer need to make a dichotomy on the candidate density since it is . Table 2 gives the total number of subintervals into which the program upper_bound_HCP.cpp eventually divided each interval as well as the execution time on our Laptop (i5-7300U CPU 2.60GHz).

4 Discussion

Figures 10, 11, 12 and 13 depicts lower and upper bound on the maximal density over each interval of interest. The upper bounds are depicted by red points and the lower bound by a green curve.

The black points indicate the smallest that we get if we only check that the parameter polytope is not empty (this corresponds to the first pass in the function find_delta in upper_bound.cpp). These points thus give a lower bound on the best upper bound that we can hope to obtain with our method, assuming that we know how to choose the parameters in an optimal way. This lower bound is however not necessarily achievable (this is in particular the case when the black points are under the green curve, i.e., under the lower bound proved on the maximum density): the parameter polytope could be not empty without containing any parameter that allows to prove the candidate bound. Indeed, the polytope is not calculated in an exact way, whereas the verification once the parameters are chosen is exact and it is used to establish the upper bound.

Nevertheless, for the ratios with a large difference between red and black points, the upper bound can probably be improved by choosing the parameters more finely. The upper bound turns out to be very sensitive to the choice of parameters in the polytope, which is therefore a delicate exercise (especially since it is done automatically here because of the very large number of different ratios considered).

The black points also show that the largest intervals over which one can expect our method to prove that the maximal density is are only slightly larger than the ’s in Corollary 2, namely:

 [0.4398,0.4644],[0.4862,0.5182],[0.5624,0.6285],[0.6455,1]. Figure 10: Between I3 and I4, that is, for r∈[0.6269,0.6469]. The upper bound (red points) is quite close to the lower bound (green curve). We conjecture that the lower bound is tight on this interval. For r1≈0.6375, the lower bound reaches its maximum and has been proven to be tight [BF21]. Figure 11: Between I2 and I3, that is, for r∈[0.5145,0.5666]. For r3≈0.5332 and r2≈0.5451, the lower bound reaches two local maxima and has been proven to be tight [BF21]. For r≤r2, the upper bound (red points) are quite close to the lower bound (green curve). We conjecture that the lower bound is tight on this interval. For r≥r2 the lower and upper bounds diverge. The difference between black and red dots is moreover relatively large between 0.552 and 0.558, suggesting (but not proving) that the choice of parameters could be optimized. We do not exclude the possibility that the lower bound is not optimal, i.e. that a better flow could be defined to the right of r2 (remind Fig. 6). Figure 12: Between I1 and I2, that is, for r∈[0.4532,0.4917], there is a sort of “mysterious island”. The lower bound is indeed δ1 over this whole interval: no packing more dense than the hexagonal compact packing is known. The difference between red and black points suggest that the choice of parameters could be optimized, but does not leave any hope of reducing the upper bound to δ1 over the whole interval. Our conjecture is that the lower bound is tight, that is, this mysterious island is an artefact of our method. This possible “artefact” is discussed in more detail in the text. Figure 13: On the left of I1, that is, r∈[0.106,0.445]. The lower bound has been proved to be tight on its local maxima at r8≈1547, r7≈0.2808, r6≈0.3292, r5≈0.3861 and r4≈0.4142 [BF21]. We conjecture that it still holds on a neighborhood of these ratios. We also conjecture that the lower bound is tight at its local maximum at rb≈0.3691. This is less clear around rc≈0.2168. Upper and lower bound are quite different in the valleys betweens peaks of the lower bound. This can be due to artefacts (this is our hypothesis - see text) or to unknown dense packings.

In Fig. 12 and 13, we used the term “artefact”. By this, we mean a peak in the upper bound which is though to be well above the exact maximal density. Such peaks appear for ratios such that the discs fit together particularly well around one disc. Indeed, the method developed in [BF21] is rather “local”: it distributes the densities between each disc and its close neighbors and bounds from above the resulting average density. However, these locally dense arrangements may not combine well on a more global scale, leading to packings in the whole plane that are actually much less dense than the obtained upper bound. Figure 14 illustrates the case . The smaller , the more often discs fit together particularly well around one disc. This explain why there are more and more peaks in the red or black curves in Fig. 1 or 13 when becomes small. To get around this problem, it will probably be necessary to modify the method to make it less local, i.e., to distribute the densities on a larger scale. This unfortunately makes the method even more complex. Figure 14: Discs of radius 1 and 0.48. The large disc can be surrounded by discs to form a pattern which is locally quite optimal in terms of density (three leftmost patterns). This is also the case, to a lesser extent, of the small discs (three rightmost patterns). However, it seems impossible to combine these patterns to form a dense packing in the plane. We can indeed start from one of these patterns, but the more we add discs, the less the local patterns can resemble those represented here.

Appendix A Proof of Prop. 1

Proof. Recall that the maximal density can be approached arbitrarily close by the density of a periodic packing. Fix .

First, consider a periodic packing of discs of radii and and density at least . Assume its fundamental domain has area and contains discs of radius and discs of radius . Hence

 pπy2+qπA≥δ(y)−ε.

By replacing each disc of radius by a smaller disc of radius with the same center, we get a packing of discs of radii and whose density is

 pπx2+qπA.

Since this density is, by definition, less or equal than , we have

 δ(y)−ε≤pπy2+qπA≤δ(x)+pπ(y2−x2)A=δ(x)+pπ(x+y)A(y−x).

In the initially considered packing, the fraction of covered by the discs of radius is at most , the maximal density of a packing by equal discs. Hence

 pπy2A≤π2√3.

This yields

 δ(y)−ε≤δ(x)+π(x+y)2y2√3(y−x)≤πy2√3(y−x).

Taking gives the first half of the claimed inequality.

Conversely, consider a periodic packing of discs of radii and and density at least . Assume its fundamental domain has area and contains discs of radius and discs of radius . Hence

 pπx2+qπA≥δ(x)−ε.

By scaling the whole packing by , we get a packing of discs of radii and whose fundamental domain has area . Then, by replacing each disc of radius by a smaller disc of radius with the same center, we get a packing of discs of radii and whose density is

 pπy2+qπA(y/x)2=pπx2A+qπx2Ay2.

Since this density is, by definition, less or equal than , we have

 δ(x)−ε≤pπx2+qπA=pπy2+qπA(y/x)2+qπA−qπx2Ay2≤δ(y)+qπ(x+y)Ay2(y−x).

In the initially considered packing, the fraction of covered by the discs of radius is at most , the maximal density of a packing by equal discs. Hence

 qπA≤π2√3.

This yields

 δ(x)−ε≤δ(y)+π(x+y)2y2√3(y−x)≤δ(y)+πy2√3(y−x).

Taking gives the second half of the claimed inequality.