Robert Botet
email: robert.botet@universite-paris-saclay.fr

Sylvie Kwok

Bernard Cabane

Filling space with polydisperse spheres in a non-Apollonian way

Completing a Swiss-Cheese theorem of Lieb and Lebowitz, we prove that any population of spheres with power-law radius distribution ∝ 1/r d f +1 can completely fill 3D Euclidean space if the exponent is such that 2.8 ≤ d f < 3. This sufficient condition extends considerably the known part of the ensemble of spacefilling populations of polydisperse spheres. The self-similar spatial arrangement of the polydisperse spheres related to the theorem is discussed using a numerical example with d f = 2.875. By calculating the small-angle scattering structure factor of the resulting packing, we found it to present several crystalline peaks indicating some regularity. This is significantly different from the featureless structure factor of an Apollonian packing which represents total disorder. We thereby argue that the Lieb and Lebowitz algorithm for filling space with spheres is fundamentally different from Apollonian constructions.

Introduction

Dense vitreous materials are much sought after because of their unique technological applications. Indeed, they manifest remarkable mechanical properties such as high resilience and high formability [1]; they often possess amazing optical and electrical properties [START_REF] Kobwittaya | Materials Chemistry of Ceramics[END_REF]; their thermal properties are anomalous [START_REF]Mechanical and Thermal Properties of Ceramics[END_REF] etc. Two classic examples of dense vitreous materials are 1) vitrified ceramics [4] and 2) ultra-high performance concrete [START_REF] Schmidt | Ultra High Performance Concrete (UHPC) -Proceedings of the International Symposium on Ultra High Performance Concrete[END_REF]. A main issue for the practical preparation of such materials is that the disordered local structure requires a lot of energy (thermal and mechanical) to obtain and control [6]. Therefore, any avenue to lower energy consumption is of paramount interest.

An efficient way to get around the problem is to pack together a large population of small particles to form a paste as dense as possible before melting or sintering the paste to eliminate tiny voids and obtain the final vitreous material [7]. Intuitively, it is clear that such an approach requires a very wide initial size distribution since one has to use smaller and smaller pieces to occupy the likewise increasingly smaller voids remaining in the packed system during the filling process. Hence the crucial question: is there an optimal particle size distribution that can give the most compact paste by simple mixing?

We address here this question when the particles are all spherical in shape, because spherical particles are either natural in some cases (e.g. extremely dense emulsions of spherical droplets may result from elastic energy minimization [8]), or provide good approximation for small pebbles or convex grains. In spite of Ulam's packing conjecture [START_REF] Gardner | The Colossal Book of Mathematics: Classic Puzzles. Paradoxes, and Problems[END_REF] stating that spheres are one of the worst possible cases of all the convex bodies for occupying space densely [START_REF] Conway | [END_REF]11], the number of space-filling populations of spheres is nonetheless expected to be very limited and thus relatively easy to identify. From the current state of the art, we know that polydisperse sphere populations with power-law distributions, n(r) ∝ 1/r d f +1 , of the sphere radii r, can geometrically fill space entirely without overlaps when the exponent d f has specific values: d f = 2.47 (numerical result) [12], d f = 2.73 (numerical result) [13], any value of d f such that 2.47 < d f < 2.88 (numerical range) [14] and an infinite number of discrete values of d f between 2.9885 and 3 (exact result) [15].

Interestingly, the last result, engendered by a little-known theorem of Lieb and Lebowitz (LL) [15], is often neglected in modern discussions pertaining to dense spherepackings. A brief literature search in the domain reveals a multitude of publications relating to Apollonian packing algorithms, but few mentions of the LL algorithm (likely due to its description having been casually tucked away by the authors within a dense, mathematical, work). In the present work, we seek to revive general interest in the LL theorem, by extending it to demonstrate that any population of spheres with 2.8 ≤ d f < 3 (exact range) can be packed as densely as desired without overlaps. In fact, several naturally-occurring distributions are empirically known to be within these bounds, e.g. (see [16] for extensive data) gravel (d f = 2.82) or glacial tills (d f = 2.88). We compare LL and Apollonian packings in the 3D space, deducing that the extended LL theorem authorises space-filling populations of spheres previously inaccessible by Apollonian constructions.

Our initial motivation was to compare packings resulting from the LL algorithm and packings found in our own experiments of extremely dense emulsions of spherical droplets [8,[START_REF] Kwok | [END_REF]. We then discuss in the present work such sphere arrangements in the only 3D Euclidean space. Similar work could also be interesting to do in the 2D space, though comparison with real experiments might be more limited in that case.

The mathematical context

The enduring popularity of Apollonian constructions has a historical root. From the mathematical point of view, Leibniz [START_REF] Leibniz | The Leibniz -Des Bosses correspondence[END_REF] was the first to propose an algorithm to cover the 2D Euclidean space with an infinite number of disks of decreasing radii, based on a geometric construction introduced by Apollonius de Perga that builds a disk tangent to three existing disks. Much later, Mandelbrot [START_REF] Mandelbrot | The Fractal Geometry of Nature[END_REF] conjectured that one could use a similar algorithm to fill the 3D Euclidean space with spheres arranged in a self-similar way. A rigorous self-similar packing algorithm based on inversive geometry was proposed in [14]. All Apollonian algorithms are based on a maximal geometric condition and is generically called 'osculatory packing' [12] or 'Apollonian packing' [START_REF] Farr | [END_REF]. In 3D, it is based on the iterative addition of the largest interstitial sphere compatible (that is without overlaps) with the biggest void remaining in the system [21]. The algorithm can be used to build regular and disordered packings alike, attaining the final volume fraction ϕ = 1 [22] (in the mathematical sense that the interior of the set complementary to the sphere packing is of null measure).

An osculatory packing definitely different from Apollonian has been suggested by Manna in [23]. In the Apollonian algorithm, the n th sphere is the largest sphere consistent with the empty pore space left by the n -1 spheres, while, in the Manna's algorithm, the centre of the n th sphere is chosen randomly in the empty pore space of the n -1 spheres and its radius is the largest possible in that place. Both algorithms are in the group of the 'osculatory algorithms' since the n th sphere is tangent to at least one of the n -1 spheres, but the Apollonian algorithm is known to require significantly less number of small spheres than the Manna's case.

The self-similar assumption

In the following, we consider the problem of totally filling a unit cube with a polydisperse population of perfect spheres.

By firstly packing spheres all of the same radius, one can obtain a system of maximum volume fraction ϕ = π/3 √ 2 ≃ 0.74 if the arrangement is regular (facecentered cubic lattice), or ϕ ≃ 0.64 if the packing is disordered (random close packing) [24]. The next step is to fill with smaller spheres the voids remaining in the system. Using again monodisperse spheres of radius smaller than the radius of the largest-void insphere, one can increase the volume fraction of the system. But because of their complicated shapes, one cannot totally fill the voids in such a way: when this second step is completed, many smaller voids still remain. These voids can be partially filled in the same manner; by iterating the process using increasingly smaller spheres, all voids are eventually filled completely when the sphere radius tends to 0. Note that such schematic description of this algorithm does not require any length scale to be defined, and one can then expect the resulting packing to be self-similar. Therefore, we are now looking for a radius distribution of the spheres, n(r), in the form of a scale-free power law [25,26]:

n(r) ∝ 1 r d f +1 . (1)
The radius distribution is written here as a continuous distribution of the radii, r, but discrete distribution can be used as well (and will be used later on). The exponent, d f , of the power law is a positive constant equal to the fractal dimension of the matter-void interface [27,28], that is, asymptotically, the union of the surfaces of all the spheres forming the packing [23].

The fractal dimensions of the Apollonian packing and related sphere packings

The value, d Apo f , of the fractal dimension of the 3D Apollonian packing has been variously and precisely estimated from intensive numerical works [12,21]:

d Apo f = 2.4739... (2)
It would appear that the value ( 2) is characteristic of an osculatory packing process as it is the same whether the initial packing is regular or random. Slightly different values of d f have also been found for non-osculatory, packing-limited growth algorithms. [13,29].

Besides the commonly discussed Apollonian packing, Lieb and Lebowitz [15] also described filling space with populations of spheres whose radii are power-law distributed as long as the corresponding value of d f is larger than 2.9885. In the following sections, we examine how the original LL theorem may be extended down to d f = 2.8, by gradually approaching a spherical geometry from known space-filling polyhedra (known as plesiohedra).

Wesler's theorem

A theorem by Wesler [START_REF] Wesler | Proc. Amer. Math. Soc[END_REF] states that the total surface of any sphere-packing filling the unit cube is infinite. Then, applying the theorem to the radius-distribution (1), the value of d f must be in the range:

2 ≤ d f < 3 . ( 3
)
This condition was tentatively improved upon by Aste [26], who argued that among all the space-filling polydisperse sphere packings, the Apollonian packing is the one with the lowest possible value of d f (i.e. with the smallest proportion of tiny spheres). Then, a necessary condition stronger than (3) for spheres filling space could be: 2.4739 ≤ d f < 3; up to now, this argument remains a conjecture.

In the following, we demonstrate and discuss sufficient conditions of the form (3) on the value of d f for a population of spheres with radius distribution (1) to be spacefilling. As the argument uses space-filling polyhedra, we start first in setting out a short review of these objects.

Space-filling polyhedra

Definitions of two geometric quantities related to convex bodies

Let us call ω(r) a finite convex body, in which r is the radius of the insphere, that is the largest sphere that can fit inside ω(r) [START_REF] Sahy | [END_REF]. We introduce two characteristic numbers related to ω(r):

• the asphericity number [11] γ c is the ratio between the radius of the circumscribed sphere (that is the smallest sphere containing ω(r)) and the radius of the insphere. The index c , in "γ c ", is to remind that this number is related to the circumscribing sphere

• the number a v comes from the volume of ω(r), written: a v 4πr 3 /3. The index v , in "a v ", reminds one that this number is related to the volume of the convex body These characteristic numbers are practical for measuring the distances in terms of r for the convex body ω(r).

Note that if ω(r) is the sphere of radius r, these quantities are respectively:

γ c = 1, a v = 1.

Introducing the plesiohedra

Perfect tessellation of space with identical polyhedra has been a long-pondered problem. An anecdote attributed to Aristotle (4 th century BC) is well known: he proclaimed without proof in his book De Caelo that space could be filled by identical tetrahedra ("among solids [it is agreed that ] only two [fill the place that contain them ], the regular tetrahedron and the cube "). This conjecture was disproved in the 15 th century by Regiomontanus in a lost manuscript [START_REF] Hilbert | Geometry and the Imagination[END_REF][START_REF] Simplicius | On Aristole: On the Heavens[END_REF]. Even back then, it was already recognised that drawing up a complete list of space-filling polyhedra was not as straightforward a problem as one would have thought.

Many different convex polyhedra have been found to fill 3D Euclidean space by tiling identical copies of them. They are generically named plesiohedra [START_REF] Grünbaum | [END_REF]. In the following we shall restrict the discussion to three convex plesiohedra selected among the five parallelohedra, because they are simple and rather "close" to the sphere. We list these selected plesiohedra in decreasing order of γ c (smaller γ c corresponds to more sphere-like polyhedra):

• the cube is the only Platonic solid to be a plesiohedron. Its characteristics are:

γ c = √ 3 ≃ 1.
73 and a v = 6/π ≃ 1.90. The Cartesian coordinates of the 8 vertices of the cube ω(r) centred in the origin (0, 0, 0) are: (± r, ± r, ± r). The cube is the Wigner-Seitz cell of the regular simple cubic lattice.

• the rhombic dodecahedron is an Archimedean solid that has:

γ c = √ 2 ≃ 1.41 and a v = √ 18/π ≃ 1.35.
The Cartesian coordinates of the 14 vertices of the rhombic dodecahedron ω(r) centred in the origin are:

(± r/ √ 2, ± r/ √ 2, ± r/ √ 2), (± r √ 2, 0, 0), (0, ± r √ 2, 0) and (0, 0, ± r √ 2)
. The rhombic dodecahedron is the Wigner-Seitz cell of the regular fcc lattice.

• the Kelvin's truncated octahedron is another space-filling Archimedean solid that has:

γ c = √ 5/3 ≃ 1.

General sufficient space-filling conditions

Let us consider the ensemble whose elements are all the populations of perfect spheres with discrete distribution of their radii, r k , such that:

{n k , r k } k=0,1,2,... means n k spheres (n k ∈ N ⋆ ) of radius r k , (4) k → r k is a decreasing sequence of positive real numbers , n k r 3 k ≤ 1 √ 32 for any k = 0, 1, 2, . . . ( 5
) ∞ ∑ k=0 n k r 3 k = 3 4π (space-filling condition of the packing) , ( 6
) ∞ ∑ k=0 n k r 2 k = ∞ (infinite surface of the packing) . ( 7
)
For every population in the ensemble above, the index k is called the rank. The condition [START_REF] Schmidt | Ultra High Performance Concrete (UHPC) -Proceedings of the International Symposium on Ultra High Performance Concrete[END_REF] expresses that the ensemble of monodisperse spheres of rank k cannot constitute a volume fraction larger than the maximum value π/3 √ 2.

Though necessary for spheres with discrete radius distribution, ( 5)-( 7) do not constitute sufficient conditions to represent a space-filling population of spheres. To do so, one has to find a geometrical algorithm able to build iteratively the sphere packing into the unit cube, adding n k spheres of radius r k without overlap at each rank k. This way is analogous to invoking a kind of Maxwell's demon able to put each sphere at its correct place in space [35].

Sufficient condition using the LL algorithm

The following algorithm -named "LL algorithm" since it was introduced by Lieb and Lebowitz in [15] -can be used to build perfect packing of spheres complying with ( 5)- (7). We describe below the algorithm starting from the 3D cube of edge length 1, though any initial 3D compact domain can, in principle, be used (in [15], the unit ball was used instead). In the case of the unit cube, the unit length is defined as the cube edge, i.e. all particle radii in the formulae below will be expressed on the scale of the cube-edge length.

Running through the LL algorithm

At rank 0 of the sphere packing, n 0 spheres of radius r 0 are placed without overlap inside the empty cube of edge length 1 (periodic boundary conditions are used). Placement of these spheres can be regular or random, and they act as seeds for the packing. Such placement is possible because of [START_REF] Schmidt | Ultra High Performance Concrete (UHPC) -Proceedings of the International Symposium on Ultra High Performance Concrete[END_REF]. When the rank 0 is completed, the volume V 0 of the domain Ω 0 left empty in the cube is:

V 0 = 1 - 4π 3 n 0 r 3 0 < 1 . ( 8
)
In the LL algorithm, a loop of iteration builds the system of rank j from the system completed at rank j -1, that is when all the n k spheres of radius r k , k = 0, 1, • • • , j -1, have been positioned inside the cube by the Maxwell's demon. The main loop is as follows:

• loop "j -1 → j" :

-step 1 : part of the empty domain Ω j-1 is covered with a dense regular array of ν j non-overlapping identical plesiohedra ω(r j ), each of them totally included inside Ω j-1 (that is none of these polyhedra intersects the boundaries of Ω j-1 ). This step can always be completed as long as the value of of r j is so small that ν j ≥ 1.

-step 2 : n j plesiohedra are chosen amongst the ν j , and replaced by their inspheres.

The selection of the plesiohedra is random. The value of n j is given by (4). The ν j -n j remaining polyhedra are left empty. Clearly, this step can be achieved provided n j ≤ ν j .

After the loop has been completed, it results in the void Ω j-1 of volume

V j-1 = 1 - 4π 3 j-1 ∑ k=0 n k r 3 k < 1 , (9)
in the system. The loop above starts from j = 1, and is continued for all integer positive values of j. The radii r k being all > 0, the sequence of void volumes V j-1 , as written in [START_REF] Gardner | The Colossal Book of Mathematics: Classic Puzzles. Paradoxes, and Problems[END_REF], is decreasing and its limiting value vanishes for j → ∞ because of the complete space filling condition (6). A sketch of the iterative process as a Swiss-Cheese construction, is shown on Fig. 2 in the 2D case.

It is clear from the construction shown in Fig. 2 that the LL algorithm is a variant of the Sierpiński algorithm used to build 3D Menger cubes [START_REF] Mandelbrot | The Fractal Geometry of Nature[END_REF]. Let us recall that, similarly to the asymptotic complete sphere-packing, Menger cubes have asymptotically zero volume (corresponding to ( 6)) and infinite surface (corresponding to (7)). Noticeable differences are: 1) open spheres (open disks in 2D space) are removed at each rank, instead of open cubes for the Menger sponge ; 2) LL algorithm is basically random (though it is not mandatory) whereas Sierpiński algorithm is usually regular.

The only condition (by induction) required for the LL algorithm to be continued to infinity is that n j ≤ ν j at each rank j. This is an explicit sufficient condition, but hard to use because the intricate shape of void domains complicates the task of calculating the number ν j of plesiohedra included in the void domain Ω j-1 . To overcome this difficulty, one could search analytically for the lower-bound number of plesiohedra, depending on the actual shape of Ω j-1 . If such a lower-bound can be found (named R j ), the less stringent sufficient condition n j < R j can be used instead for each value of j = 1, 2, • • • .

Determining R j

Considering the system at rank j -1, the domain Ω j-1 is the union of ν j complete polyhedra ω(r j ) ('complete polyhedron' means here: 'totally inside Ω j-1 ') and of a domain δΩ j , of volume δV j , which is too thin to accept any complete polyhedron ω(r j ):

V j-1 = ν j v j + δV j . ( 10
)
By construction, the domain δΩ j cannot include any sphere of diameter 2γ c r j (that is the diameter of the sphere circumscribing a polyhedron ω(r j )), then δV j is smaller than the sum of the volumes of the outer shells of width 2γ c r j around each sphere of rank k ≤ j -1: The inequality (11) is the key ingredient that allows one to write explicitly a sufficient condition for the population (4) of spheres be space-filling. This inequality is more precise than the corresponding inequality written by Lieb and Lebowitz in their work (the expression (3.4) in [15]) and it results in a much better estimate of the lower bound R j . Using ( 6), ( 9) and ( 11), the following lower bound for ν j is obtained:

δV j ≤ 4π 3 j-1 ∑ k=0 n k [ (r k + 2γ c r j ) 3 -r 3 k ] . ( 11
)
ν j ≥ 3 4πα v 1 r 3 j - 1 α v j-1 ∑ k=0 n k ( r k r j + 2γ c ) 3 . ( 12
)
It results in the sufficient condition (named: (CS 1)):

the LL algorithm builds a space-filling packing of spheres with distribution ( 4)-( 7) if the inequality:

3 4π 1 n j r 3 j - j-1 ∑ k=0 n k n j ( r k r j + 2γ c ) 3 ≥ α v ( 13
)
holds true for every value of j = 1, 2, • • • In ( 13), the parameters α v and γ c depend only on the plesiohedron used in the LL algorithm.

Sufficient condition for a generic fractal sphere population

In this section we shall derive the sufficient condition on the exponent d f , for the discrete distribution of spheres ( 4)-( 7) with:

n k = n 0 n k ; r k = r 0 n k/d f , k = 0, 1, 2, • • • , ∞ (14)
such that these spheres are space-filling using the LL algorithm. This discrete distribution writes n k ∝ 1/r

d f k
similar to the continuous distribution (1) written: n(r) dr = (1/r d f ) dr/r. Indeed, from (14), the index k is proportional to the logarithm of the corresponding radius, hence the term dr/r may be put aside in the continuous distribution.

There are four constant parameters in (14), namely: n 0 and n are two integer numbers such that n 0 ≥ 1 and n ≥ 2, and two real parameters r 0 and d f such that 0 < r 0 < 1/(2 5/6 n 1/3 0 ) and 2 ≤ d f < 3, to be consistent with ( 5) and (7). We have also the rough but useful inequality: γ c r 0 < 1, which comes from r 0 < 1/2 (because the edge length of the cube is 1), and γ c < 2.

We note first that, using the distribution ( 14), the space-filling condition (6) writes as the equality:

4π 3 n 0 r 3 0 = 1 -n 1-3/d f , ( 15
)
from which the volume fraction of the system at rank j ≥ 0, namely:

ϕ j = ∑ k=0,••• ,j n k 4πr 3 k /3
, is found to depend only on the two parameters n and d f :

ϕ j = 1 - 1 n (j+1)(3/d f -1) . ( 16
)
Since n ≥ 2, ϕ j tends clearly to 1 when j → ∞. This formula can also be written:

1 -ϕ j = ( r j+1 r 0 ) 3-d f , (17)
a popular formula for the porosity of fractal sphere-packing [36] (in theoretical physics) and of concrete [START_REF]Caquot Le role des matériaux inertes dans le béton in Société des Ingénieurs Civils de France pps[END_REF] (in applied physics).

The seed-particle radius r 0

From the equality (15), the condition ( 5) gives an upper bound of the parameter n:

2 ≤ n < 2 2d f /(3-d f ) , ( 18
)
where we used for the sake of simplicity the inequality: 1 -π/ √ 18 > 1/4. Also from (15), one finds the value of r 0 as a function of the three other parameters n 0 , n and d f :

r 0 = ( 3 4πn 0 ( 1 -n 1-3/d f ) ) 1/3 . ( 19
)
4.2. The condition R j /n j ≥ 1

Using ( 14) and ( 19), the sufficient condition ( 13) is written:

f 1 + 6γ c n j(1-2/d f ) g 2 + 12γ 2 c n j(1-1/d f ) g 1 + 8γ 3 c n j g 0 ≥ α v , ( 20
)
in which the four functions f 1 , g 0 , g 1 , g 2 , independent of the rank j, are given by:

f 1 ≡ 1 1 -n 1-3/d f - ( 6γ c n 1-2/d f -1 + 12γ 2 c n 1-1/d f -1 + 8γ 3 c n -1 ) , ( 21
)
g m ≡ 1 n 1-m/d f -1 , ( 22
)
and m = 0, 1, 2. Since all the functions g m are > 0, the relation ( 20) is fulfilled if the weaker inequality (independent on the rank j): f 1 > α v is realized.

Summarizing the sufficient condition

The corresponding sufficient condition (named: (CS 2)) is written as:

the LL algorithm builds a perfect space-filling packing of spheres with the distribution (14) if there is at least one integer number n ∈ [2, 2 2d f /(3-d f ) ) (the notation [a, b) denotes the interval a ≤ n < b), such that:

f 1 ≡ 1 1 -n 1-3/d f - ( 6γ c n 1-2/d f -1 + 12γ 2 c n 1-1/d f -1 + 8γ 3 c n -1 ) > α v . ( 23
)
Once the inequality ( 23) is fulfilled for a value of n, we may select any integer number n 0 such that:

n 0 > 6 π ( 1 - 1 n 3/d f -1 ) , ( 24
)
which comes from the fact that r 0 < 1/2 (no sphere can be larger than the edge length of the box) and relation (15). At last, the value of r 0 is given by ( 19) as a function of d f , n and n 0 .

5. Explicit sufficient conditions for various space-filling polyhedra approaching the derivation for spheres

The cube

The cube has parameters:

γ c = √ 3 ; α v = 6 π , ( 25
)
hence the expression of the function f 1 :

f 1 ≡ 1 1 -n 1-3/d f - ( 6 √ 3 n 1-2/d f -1 + 36 n 1-1/d f -1 + 24 √ 3 n -1 ) . ( 26
)
The function f 1 is an increasing function of d f for any fixed value of n > 2. The inequation f 1 ≥ 6/π, with f 1 given in (26) admits solutions in integer numbers n ∈ [ 2, 2 2d f /(3-d f ) ) for any value of d f in the range:

2.855 ≤ d f < 3 , ( 27
)
The evolution of the f 1 as a function of n is plotted in Fig. 3 for several relevant values of d f . The sufficient condition (27) is much wider than the set of discrete values d f = 2 + log p/ log(p + 1) (with p any integer ≥ 26) derived in [15] and which are all included in the range [2.9885, 3).

The rhombic dodecahedron

This polyhedron has parameters:

γ c = √ 2 ; α v = 8 . ( 28
)
It is closer to the sphere in the sense that the value of γ c is smaller than for the cube. The same analysis as before leads to the sufficient condition:

2.847 ≤ d f < 3 , ( 29
)
that is a slightly wider range than for the cube case.

The truncated octahedron

This polyhedron is more interesting than the previous one since its parameters:

γ c = √ 5/3 ; α v = 8/π √ 3 (30)
show that its shape is much closer to the shape of a sphere (the ratio between the radius of the circumscribed sphere and the one of the inscribed sphere is only 1.29). In this case, the same analysis as before leads to the sufficient condition:

2.800 ≤ d f < 3 , (31)
The stricter sufficient condition (CS 1) allows the replacement of the lower threshold 2.800 by 2.798 in (29).

In short, using the sufficient condition (CS 2) with, successively, the cube, the rhombic dodecahedron and the truncated octahedron as the plesiohedra required in the LL algorithm, we showed the following theorem:

• for any value of d f in the range 2.8 ≤ d f < 3, it is possible to pack iteratively

n k ∝ 1/r d f spheres of radii in geometric progression r k ∝ ρ k (with 0 < ρ < 2 -1/3
) in the unit 3D cubic box. The packing is complete in the sense that the packing density equals 1 for k → ∞.

A note about the best possible sufficient condition

The only parameters of the sufficient condition (23) are γ c and α v of the plesiohedron used in the LL algorithm. The more sphere-like is the plesiohedron, and the closer to 1 are these two parameters, no matter how complicated the shape of the polyhedron. Even aperiodic space tiling [START_REF] Grünbaum | Tilings and patterns Chap[END_REF] can be used as well. So far, we have considered only three plesiohedra, with the truncated octahedron having a shape relatively close to the sphere (this is the "Kelvin's polyhedron" [START_REF] Weaire | The Kelvin Problem: Foam Structures of Minimal Surface Area London[END_REF]).

At this juncture, one wonders how much more the range of fractal dimensions given by ( 31) could be expanded, if all possible plesiohedra (known or yet to be discovered) were considered in this gradual geometrical approach towards perfect sphericity. A partial answer may be gleaned from the remark that the parameters γ c > 1 and α v > 1 for any plesiohedron, with the limiting values = 1 for both parameters (these limit values are the respective values for the perfect sphere). Let us then suppose that we know a sequence of different plesiohedra with decreasing values of γ c approaching 1. Such sequence could begin with: /cube (γ c = 1.73)/ rhombic dodecahedron (γ c = 1.35)/ truncated octahedron (γ c = 1.29)/ • • • . Provided such a sequence exists, the smallest value of d f compatible with a space-filling sphere population (14) under the LL algorithm, is given by the sufficient condition (23) written with γ c = α v = 1:

1 1 -n 1-3/d f - ( 6 n 1-2/d f -1 + 12 n 1-1/d f -1 + 8 n -1 ) > 1 , (32)
We first solve numerically the corresponding equation in d f , with the variable n > 1:

1 1 -n 1-3/d f - ( 6 n 1-2/d f -1 + 12 n 1-1/d f -1 + 8 n -1 ) = 1 , ( 33
)
For any value of n, there is only one solution of [START_REF] Simplicius | On Aristole: On the Heavens[END_REF], which is shown on the Fig. 4. Moreover, when n ≫ 1, the two first terms of (33) are dominant, and the equation leads to:

n 2-5/d f ≃ 6 (34)
that is: d f ≃ 5/(2 -log 6/ log n) which tends to 5/2 when n → ∞.

There is no acceptable solution corresponding to d f ≤ 5/2. The more general sufficient condition (CS 1) leads to the same conclusion since the left-hand term of (13), calculated for d f = 5/2, grows from a negative value for n = 2, to the value 1 -(as 1 -5/n 1/5 when n → ∞) for any value of j. We conclude that, using the algorithm LL, the range of values of d f sufficient for the packing of spheres ( 14) be space-filling is always included in the interval:

5/2 < d f < 3 .
An interesting remark here is that the Apollonian fractal dimension d Apo f ≃ 2.474 ± 0.001 being smaller than 2.5, the Apollonian sphere-packing cannot be a population (14) packed using the LL algorithm, whatever the plesiohedron we use in this algorithm. In Section 6, we will examine more closely the difference between packings resulting from LL and from Apollonian constructions.

Conceivable improvements of the theorem

As discussed in Section 5.4, the theorem derived in Section 5.3 can be straightforwardly improved using space-filling polyhedra 'rounder' than the truncated octahedron adhering to the sufficient condition (23). These polyhedra remain to be found, and their values of α v and γ c can then be calculated to see how the range of acceptable d f is affected.

A more complicated improvement might result from finding a better lower bound R j of the coefficients ν j , following the reasoning detailed in Section 3.3. This would require an analysis of the shapes of the voids left in the intermediary sphere packings, more carefully than the simple calculation presented in 3.3. This approach could probably be effective, but requires much work.

Another improvement consistent with the LL algorithm would include tessellation of space with two (or more) different sorts of polyhedra [START_REF] Conway | [END_REF]. This would change the formula appearing in the sufficient condition (CS 1), involving then two (or more) values of γ c and α v . That way, two (or more) sphere radii in fixed ratios between them would have to be considered at each rank of the LL algorithm iteration.

A numerical example

In this Section, we propose using the LL algorithm to build numerically the first stages of a sphere packing in cube of edge-length L = 1 with periodic boundary conditions. The fractal dimension d f = 2.875 > 2.8 is selected, and we use the truncated octahedron in the algorithm.

The inequation (23) solves into:

74 ≤ n ≤ 2.2 10 11 (35)
Selecting the value n = 74 and the number n 0 = 2 of initial spheres, we give in

• • • • • • • • • • • Table 1
. Parameters of a possible numerical simulation of sphere packing using the LL algorithm with truncated octahedron as the plesiohedron. The sphere population is chosen with fractal dimension d f = 2.875, n = 74 and n 0 = 2. When the rank j of the iteration increases, the total number, N s , of spheres increases exponentially, while the volume fraction ϕ j increases very smoothly. Numerical simulation of such sphere packing is hardly conceivable beyond the rank j = 4 using this method.

We performed a numerical simulation up to j = 3 (number total of spheres: N s = 821 550) using the LL algorithm starting from 2 randomly-placed spheres of radius r 0 = 0.2731 in the unit cube. At every rank k = 1, 2, 3, the selection of the n k truncated octahedra among the ν k , is performed randomly. Visualization of a small part of the packing showing the particles of ranks 0, 1, 2 and 3 is shown on Fig. 5, and the successive regular lattices, on which the truncated octahedra are lying, are quite visible on it.

A similar conclusion about the spurious evidence of local lattices of plesiohedra in the LL algorithm case can be drawn from the analysis of the sphere-packing in the reciprocal space. The reciprocal space is the natural framework to analyze the Small-Angle X-ray Scattering data from systems of particles. In this space, the natural coordinate frame is represented by the scattering wave vector coordinate q = (2π/λ)(u sca -u inc ), in which λ is the wavelength of the scattered X-ray, and u inc , u sca are the unit vectors in the incident and scattered directions respectively. The static structure factor [41]:

S(q) = ∑ i,j f i f j e iq(r i -r j ) / ∑ i f 2 i (36)
of an ensemble of spheres (with the centre of the sphere labelled i is located in r i , and its form factor is f i ) is a measured quantity which includes most of the information about the spatial arrangement of the spheres in the system. The structure factor depends only on the modulus q = 4π/λ sin(θ/2), in which θ is the scattering angle, when spherical symmetry applies.

In the case of the LL algorithm packing of spheres, the structure factor function is almost constant, with a number of strong fine peaks reminiscent of Bragg peaks coming 5, built using the LL algorithm up to level 3. The restricted Braggs peaks are due to the limited regular packing sets inherent to the LL algorithm. The variable q is counted in units determined by the choice of the unit length (the edge of the cube). For comparison, the dotted red curve is the structure of an Apollonian packing, as discussed in [START_REF] Kwok | [END_REF], and the blue curve quite close to value 1 is for packing obtained using the Manna's algorithm, as discussed in [13]. In these three examples, the systems were made of about 800 000 spheres, and angle-averaging was performed with 25 000 random Euler angles. No Bragg peak is seen either in the Apollonian packing or in the Manna packing structure factor.

from the BCC symmetry of the tessellation of space by truncated octahedra (see Fig. 6). The LL structure factor is quite different from the structure factor of an Apollonian packing (dotted red curve) and from the structure factor of a Manna packing (blue curve), both shown on the Fig. 6 for comparison. This result confirms our conjecture that the osculatory algorithms are different in essence from the LL algorithm.

Conclusion

In the present work, we extended a theorem from Lieb and Lebowitz to derive a wide range of values of the exponent d f for which a population of spheres with continuous radius-distribution ∼ 1/r d f +1 is space-filling. The method can probably be improved and hints are given in the text. Although the basis of the LL algorithm does not allow for the inclusion of all possible space-filling populations of spheres -as evidenced by the comparison with Apollonian packing -the extended LL theorem nevertheless grants access to a much larger set of sphere populations, previously unknown or unexplored. Physical properties of extremely dense sphere packings, as generated by the LL algorithm or by the Apollonian osculatory algorithm, remain to be discovered. However, unusual properties can be expected if one refers to the case of the Menger sponges since both sorts of systems share basic features (e.g. extreme density, extreme surface, self-similar spatial distribution of matter, etc.). Indeed, Menger sponges are known for the spectacular -and still unexplained -localization of electromagnetic waves [42] or strongly anomalous thermal diffusion [43]. But unlike the Menger cubes -which are artificial structures -, extremely dense sphere packings appear naturally as high internal-phase-ratio emulsions [START_REF] Kwok | [END_REF]8]. For studies and applications involving physical properties of extreme sphere-packing, the role of the exponent d f is probably of utmost importance and we may now begin to understand how the value of d f is related to the building process.

3 )

 3 29 and a v = 8/π √ 3 ≃ 1.47. The Cartesian coordinates of the 24 vertices of the truncated octahedron ω(r) centred in the origin are: (0, ± r/ √ 3, ± 2 r/ √ and the six permutations of the three coordinates. The truncated octahedron is the Wigner-Seitz cell of the regular bcc lattice.

Figure 1 .

 1 Figure 1. Sketch of the three plesiohedra used as examples in this work. From left to right: the cube, the rhombic dodecahedron and the truncated octahedron. The respective circumscribing spheres are shown to realise how 'round' they are.

Figure 2 .

 2 Figure 2. The 2D LL algorithm in the unit square, described as a Swiss-Cheese model. At rank k = 0, an initial tessellation with big squares is used and the initial holes are the open disks inscribed in n 0 squares. These initial circular holes are kept at rank k = 1 . Then, n 1 smaller squares of edge length r 1 are created to tile the plane, wherein open disks are randomly inscribed so long as they do not intersect with any other holes. At rank k = 2, smaller tiling squares of edge length r 2 are used in which newer open disks are now inscribed without overlaps with all existing holes, etc. In this manner, the iteration is performed such that the plane at rank k is punched by a number of non-overlapping open disks, and, for k → ∞, the unit square is full except the complementary set of the union of all the open holes. This complementary set is of null measure if the condition (6) holds.

Figure 3 .

 3 Figure 3. Shape of the function f 1 as given in (26) in the case of the cube as the plesiohedron used in the LL algorithm. Three values of the fractal dimension d f are used to exemplify the threshold of (23). The maximum of the function f 1 increases regularly with the value of d f . For d f > 2.855, the values of f 1 are larger than α v = 6/π ≃ 1.91 (the horizontal line) on some interval of n, making true the sufficient condition for the sphere population (14) be space-filling.

Figure 5 .

 5 Figure 5. Cut through the rank 3 sphere packing obtained using the LL algorithm with truncated octahedron as the plesiohedron. The parameters are: d f = 2.875, n = 74, n 0 = 2. Colors correspond to the sizes of the spheres (that is: r 0 = 0.273 (blue), r 1 = 0.061 (green), r 2 = 0.014 (yellow), r 3 0.003 (red)). The size of that part of the packing is 10% of the unit cube. Since it is a 2D projection of a 3D arrangement of spheres, this image is not as informative as the 2D model shown on Figure 2. However, one can see clearly the underlying regular truncated-octahedra honeycomb used by the LL algorithm, despite the random selection of the sphere locations amongst the available sites at each iteration.

Figure 6 .

 6 Figure 6. Structure factor (black curve) of the ensemble of spheres shown on Fig.5, built using the LL algorithm up to level 3. The restricted Braggs peaks are due to the limited regular packing sets inherent to the LL algorithm. The variable q is counted in units determined by the choice of the unit length (the edge of the cube). For comparison, the dotted red curve is the structure of an Apollonian packing, as discussed in[START_REF] Kwok | [END_REF], and the blue curve quite close to value 1 is for packing obtained using the Manna's algorithm, as discussed in[13]. In these three examples, the systems were made of about 800 000 spheres, and angle-averaging was performed with 25 000 random Euler angles. No Bragg peak is seen either in the Apollonian packing or in the Manna packing structure factor.

  Figure 4. Solution (black curve) in d f of the equation (33) versus the variable n. The red dashed line is the special value d f = 5/2. The asymptotic solution for n → ∞ is d f ≃ 5/(2 -log 6/ log n) → 5/2. figure shows that, for given value of d f , the inequation (32) has solution in n only if d f > 5/2.

	d f
	2.8
	2.6
	2.4
	2.2
	10 000	100 000 000	n

Table 1

 1 the characteristics of the numerical simulation compatible with the parameters deduced from the sufficient condition(23).

	j	N s	r j	ϕ j
	0	2	0.2731 0.171
	1	150	0.0611 0.312
	2	11 102	0.0137 0.430
	3	821 550	0.003	0.527
	4	6.0 10 7 6.8 10 -4 0.608
	5	4.5 10 9 1.5 10 -4 0.675
	6	3.3 10 11 3.4 10 -5 0.730
	• • •	• • •	• • •	• • •
	15 2.2 10 28 4.8 10 -11 0.950
	•


Etienne Cantin

Olivier Lopez

Nicolas Cahuzac

Yuhao Liu

Mathieu Manceau

Benoît Darquié

Christian Chardonnet

Anne Amy-Klein

Nicolas Quintin

Martin Rabault

Grégoire Coget

Peter Rosenbusch

Vincent Ménoret

Nicolas Quintin Renater

Mads Tønnes

Benjamin Pointard

Maxime Mazouth

Hector Alvarez-Martinez

Caroline B Lim

Philip Tuckey

Michel Abgrall

Rodolphe Le Targat

Paul- Eric Pottie

REFIMEVE Fiber Network for Time and Frequency Dissemination and Applications

Keywords: optical link, frequency dissemination, fiber network, high-precision molecular spectroscopy

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

REFIMEVE is a national research infrastructure composed of optical fiber links to disseminate time and frequency reference signals from the French National Metrological Institute, LNE-SYRTE, to more than 30 research labs and research infrastructures all over France, over the optical fiber backbone of RENATER, the French National Research and Education Network. Consolidation and expansion of the network, together with its extension to time transfer and new users were recently funded and are currently being implemented. We have the objective to provide the ultra-stable signals with an uptime better than 90% and thus we have developed effective remote control and supervision of all links [START_REF] Cantin | An accurate and robust metrological network for coherent optical frequency dissemination[END_REF]. Last year the network was extended up to 2x3500 km and it has currently connections to UK, Germany, Italy and Switzerland, as shown in Figure 1.

II. NETWORK EXTENSION

The REFIMEVE signal source is an ultra-stable laser whose drift is actively compensated using a continuous comparison to a H-Maser with an operational optical frequency comb. The fiber network uses essentially the academic network of RENATER, thanks to Dense Wavelength Division Multiplexing with data traffic. It is mainly equipped with Repeater Laser Stations, which repeat the ultra-stable signal from one fiber link to another [START_REF] Guillou-Camargo | First industrial-grade coherent fiber link for optical frequency standard dissemination[END_REF], and Multibranch Laser Stations acting as hubs [START_REF] Cantin | An accurate and robust metrological network for coherent optical frequency dissemination[END_REF]. We have also implemented some industrial-grade Extraction Stations for local use in urban areas [START_REF] Bercy | In-line extraction of an ultrastable frequency signal over an optical fiber link[END_REF]. An ensemble of low noise phase-and-frequency recorders, local oscillators and GPS receivers enables remote measurements and clock comparisons at border crossings.

All equipment is remotely controlled and for every link from locations A to B, we implement a second parallel link from B to A in order to assess and trace the proper dissemination of the signal and supervise the network in real time. The link's residual stability is evaluated from the phase Figure 1: schematic of the REFIMEVE network comparison of the signal injected in a link and the signal after dissemination from A to B and back to A. Figure 2 displays the residual frequency fluctuations of such an end-to-end phase comparison for four REFIMEVE links during 5 months. The link's length ranges from 2x340 to 2x900 km. The uptime of the links varies between 71 and 85%. The residual frequency fluctuations are typically 10 -17 or below for an integration time of 1000 s, enabling precise clocks comparison within less than one hour. Figure 3 displays the corresponding end-to-end stabilities, calculated with the modified Allan deviation. The short-term instability is typically below 10 -15 at 1-s averaging time and the long-term instability at a level of 10 -19 or below. We also checked that the bias of the transmitted frequency is compatible with zero within an uncertainty of respectively 4×10 -20 , 6×10 -20, 1×10 -19 and 2×10 -19 for the four links (a) to (d) considered in Figures 2 and3.

Due to the missing data in long-term measurements, the estimation of the link stabilities is not straightforward. We either concatenated the data or replaced the missing data with simulated data of the same statistical properties. This last approach enables us to mitigate the effects of missing data as shown in [START_REF] Tønnes | Coherent fiber links operated for years: effect of missing data[END_REF]. We applied both approaches to a very long data set of 5 years of operation of the Paris-Strasbourg link, with 44% uptime, as displayed on Fig. 4. The graph shows that the stability obtained when missing data are replaced with simulated data (violet) is fully compatible with the stability calculated from a data set without any missing data at shortterm (orange), and with concatenated data when the minimal duration of a data subset exceeds the statistical coherence time. Without a proper set of filtering parameters, a concatenated approach can lead to erroneous evaluation of the relative frequency stability and introduce frequency biases (green curve). Furthermore, with this processing method we were able to assess the link residual instability at 1-year averaging time and the link uncertainty at a level of a few 10 -21 .

III. RECENT DEVELOPMENT

We have recently extended the REFIMEVE network both in Paris region and on a national scale, to the south-west part of France and to CERN in Switzerland, which is a major European interconnection point for academic networks and GÉANT.

A schematic of the current network in Paris region is shown on Fig. 5. It disseminates the optical signal and a White Rabbit time signal to all the main universities in north, center and south of Paris in a star network configuration. All the links depart from LNE-SYRTE where a multibranch station enables us to coherently disseminate the optical reference to all the users [START_REF] Cantin | An accurate and robust metrological network for coherent optical frequency dissemination[END_REF]. Extraction stations are used at Sorbonne Université and Université Paris-Saclay to drop the signal to the user labs [START_REF] Bercy | In-line extraction of an ultrastable frequency signal over an optical fiber link[END_REF].

The link to CERN is a branching link of length 200 km along the main link between Paris and Italy, as displayed on Fig. 6. Its stability is limited by the residual noise of the previous segments of the link, which totalize a length of 720 km. However its stability reaches the 10 -18 scale in less than 2 hours, and a few 10 -19 at 1-day averaging time. This performance could be improved in the future with a better filtering of the accumulated noise at the branching node.

IV. CURRENT APPLICATIONS

The ultra-stable frequency signal disseminated with REFIMEVE has been used for more than 10 years for precision measurements in the Paris regional area, mainly at LKB and LPL [START_REF] Santagata | High-precision and wide tuneable mid-IR spectrometer with traceability to primary frequency standards[END_REF]. It is disseminated to a few laboratories in Besançon, Lille, Grenoble [START_REF] Votava | Comb coherencetransfer and cavity ring-down saturation spectroscopy around 1.65 µm: kHz-accurate frequencies of transitions in the 2ν3 band of 12 CH4[END_REF], Marseille and Paris, where current applications of REFIMEVE concern photonics, laser stabilisation or control, atomic and molecular spectroscopy. Current developments are very active in precision molecular spectroscopy in the mid-infrared spectral range around 10 µm, with molecules of interest for atmospheric and astrophysics studies [START_REF] Santagata | High-precision and wide tuneable mid-IR spectrometer with traceability to primary frequency standards[END_REF][START_REF] Votava | Comb coherencetransfer and cavity ring-down saturation spectroscopy around 1.65 µm: kHz-accurate frequencies of transitions in the 2ν3 band of 12 CH4[END_REF] or for test of fundamental physics [START_REF] Bielsa | HCOOH high-resolution in the 9.18 µm region[END_REF]. At LPL, molecules are probed with Quantum Cascade Lasers stabilised to the REFIMEVE optical reference, through an optical frequency comb. Fig. 7 shows a saturated absorption spectrum of the line P(E,co,0,2,33) of methanol, a molecule involved in numerous atmospheric reactions, and present in the interstellar medium. The spectrum is recorded in transmission of a Fabry-Perot cavity, with a resolution around 25 kHz limited by the pressure broadening. A Lorentzian fit enables us to determine the line center with an uncertainty around 600 Hz.

V. CONCLUSIONS AND PERSPECTIVES

The REFIMEVE network is mature and enables more and more labs to benefit from an ultra-stable and accurate frequency signal for their precision experiments. Moreover, thanks to its European connection, it gives the opportunity to compare the best atomic clocks in Europe, as performed with Italy during 4 months in 2021-2022 [START_REF] Clivati | Coherent optical fiber link across Italy and France[END_REF] or through a two-months measurement campaign within Europe in 2022. It is also very interesting for applications to chronometric geodesy, using these clocks comparison campaigns [START_REF] Lisdat | A clock network for geodesy and fundamental science[END_REF] or with the planned development of a mobile platform which can be placed at the output end of any links providing transportable comb and clock for clock comparisons.

REFIMEVE contributes to the roadmap towards a pan-European fiber network, as actively prepared through the project CLONETS-DS. The development of such networks is also in progress in other EU countries and all around the world, for instance in China and Japan [START_REF] Zang | Cascaded transfer of optical frequency with a relay station over a 224 km deployed fiber link[END_REF][START_REF] Akatsuka | Optical frequency distribution using laser repeater stations with planar lightwave circuits[END_REF]. Interconnection and combination with free-space techniques, such as the wellestablished radio-communication Two-Way Satellite Time and Frequency Transfer (TWSTFT) and more-ground breaking techniques such as optical TWSTFT, is a challenge to be faced to enable world-wide multi-technique optical timescale comparisons.

FUNDING

The authors have received financial support from Agence Nationale de la Recherche ANR T-REFIMEVE (ANR-21-ESRE-0029), REFIMEVE+ (ANR-11-EQPX-0039), Labex

Figure 2 :

 2 Figure 2: Relative frequency fluctuations (1000s/point) with time of 4 REFIMEVE links, during 5 months a) Paris-Lille-Paris (2x340 km); b) Paris-Strasbourg-Paris (2x650 km); c) Paris-Lyon-Modane-Lyon-Paris (2x900 km) ; d) Lyon-Marseille-Lyon (2x440 km).

Figure 3 :

 3 Figure 3: End-to-end stability (modified Allan deviation) for the 4 REFIMEVE links of Fig. 2.

Figure 4 :

 4 Figure 4: End-to-end stability (modified Allan deviation) for the Paris-Strasbourg-Paris link along 5 years of data. Missing data are taken into account either by concatenating the data or by replacing missing data with simulated data having the same statistical properties.

Figure 6 :

 6 Figure 6: End-to-end stability (modified Allan deviation) of the Grenoble-CERN-Grenoble link of 2x200 km.

Figure 5 :

 5 Figure 5: schematic of the REFIMEVE regional network in Paris.

Figure 7 :

 7 Figure 7: saturated absorption spectrum of the line P(E,co,0,2,33) of methanol at a pressure of 0,05 Pa, around 29 THz (10 µm), with 1s/point. The signal is detected in transmission of a Fabry-Perot cavity with a Quantum Cascade Laser stabilised to the REFIMEVE signal.

  (ANR-10-LABX-48-01), LNE, Domaine d'intérêt majeur SIRTEQ ATH-2019 (ONSEPA), and the CNRS through the 80|Prime program and through the Imperial College London -CNRS PhD joint programme.