Optimal Random Packing of Spheres and Extremal Effective Conductivity

A close relation between the optimal packing of spheres in Rd and minimal energy E (effective conductivity) of composites with ideally conducting spherical inclusions is established. The location of inclusions of the optimal-design problem yields the optimal packing of inclusions. The geometrical-packing and physical-conductivity problems are stated in a periodic toroidal d-dimensional space with an arbitrarily fixed number n of nonoverlapping spheres per periodicity cell. Energy E depends on Voronoi tessellation (Delaunay graph) associated with the centers of spheres ak (k=1,2,…,n). All Delaunay graphs are divided into classes of isomorphic periodic graphs. For any fixed n, the number of such classes is finite. Energy E is estimated in the framework of structural approximations and reduced to the study of an elementary function of n variables. The minimum of E over locations of spheres is attained at the optimal packing within a fixed class of graphs. The optimal-packing location is unique within a fixed class up to translations and can be found from linear algebraic equations. Such an approach is useful for random optimal packing where an initial location of balls is randomly chosen; hence, a class of graphs is fixed and can dynamically change following prescribed packing rules. A finite algorithm for any fixed n is constructed to determine the optimal random packing of spheres in Rd.


Introduction
Geometry and physics are deeply interconnected in modern science. The present paper is devoted to a recently discovered relation between the optimal packing of spheres and the effective conductivity of composites. The keyword here is "optimal", and it concerns the deterministic and random packing of inclusions in R d . Studies of packing advance our knowledge about the effective properties of densely packed composites, porous and granular media, metallic glasses, and emulsions [1][2][3].
This study focuses on nonoverlapping equal spherical inclusions. The spatially symmetric arrays of disks, i.e., the hexagonal (triangular) lattice in R 2 and the FCC or HCP lattice in R 3 , naturally arise in a deterministic statement. Random optimal packing is usually related to random-walk algorithms and their implementation (see [1][2][3][4] and references therein). In the present paper, we develop the theory of optimal random packing. The derived optimal structures can be considered to be structures having hidden random symmetry. An optimal structure achieves the local density maximum, which cannot be improved by any local perturbations within the considered class of packing. Such a local maximum translates into the physical result that effective conductivity does not exceed the conductivity of FCC and HCP lattices, which achieve the highest packing density in R 3 .
The packing problem refers to geometrical optimization problems. One of the most popular problems concerns the highest packing density of equal spheres in R d . Various geometrical methods were developed to study this problem [5,6]. The complete solution for d = 2 was given by Toth [7], and for d = 3 by Hales and Ferguson [6].
It was established in [8] that the solution to the physical problem of the optimal effective conductivity for d = 2 implies a solution to the geometrical packing problem. Effective conductivity can be defined as the energy density per a representative cell, and is determined by conductivities of constituents and by their location [9]. An optimal design problem in the theory of composites can be stated as the optimization of effective properties when the concentrations of phases are fixed, and the distribution of phases must be established [10]. Other variants of optimal-design problems with external loading were discussed in [11][12][13] and works cited therein.
The present paper concerns the optimal effective conductivity of composites with spherical inclusions. The deterministic problem can be stated in the following way. Given infinitely many equal spherical perfect conductors D i (i = 1, 2, . . .), it is necessary to locate them in the host medium of finite conductivity in such a way that the homogenized medium is macroscopically isotropic, and its effective conductivity attains the minimal value. This statement requires explanations concerning the infinite number of inclusions. The proper homogenization problem must be stated in a periodicity cell Q 0 displayed in Figure 1 with a finite number of inclusions D i (i = 1, 2, . . . , n). Of course, one may take limit n → ∞ with simultaneously decreasing radius r, but the homogenization problem must be periodic [9]. More precisely, d linearly independent vectors determine a periodicity cell in R d .
In particular, the concentration of inclusions is a given parameter in the periodic problem. Different locations of inclusions yield different values of effective conductivity. The minimal value of effective conductivity is achieved for a location of inclusions, which solves the deterministic optimal-design problem. This minimum is global since it is found over all the locations of inclusions. In the present paper, we investigate local minima that describe optimal random locations. Moreover, the physical optimal-design problem is shown to solve the geometrical optimal-packing problem in deterministic and random statements.
Recent results in structural approximations [14,15] showed that densely packed composites including random media [16] can be investigated by means of the functional associated to discrete energy. It has the following structure: where the minimum is taken over values t j prescribed to inclusion D j . Value g kj expresses the main term of the interparticle flux between neighbor inclusions D j and D k when distance δ jk between them tends towards zero. The discrete consideration of inclusions in Functional (1) yields a graph Γ with vertices at the centers of inclusions connected by edges with prescribed flux g kj . This relates minimization Problem (1) to the discrete Laplacian on a graph, with the eigenvalues and other objects of graph theory developed in [17][18][19][20]. Functional (1) was applied to the many-bodies capacity problem to explain Tamm's screening effect [21,22].
The most complete results on applications of Functional (1) were obtained in the twodimensional (2D) case when value g kj had the square-root singularity in δ jk (see Keller's formula (14) below). This immediately restricts the set of graphs to graphs generated by the Voronoi diagram. Moreover, the theory of homogenization requires the consideration of periodic problems to determine effective conductivity [9]. This also concerns the structural approximation that should generally be applied to periodic discrete structures. Therefore, periodic Functional (1) must be considered on toroidal manifolds. This restricts the direct applications of general graph theory [17][18][19][20] that is valid for nonperiodic graphs. Some results on 2D structural approximations were established in the papers outlined below. Random doubly periodic graphs were studied and applied to the homogenization problem in [16]. The minimization of periodic Functional (1) with respect to its geometry was applied to pattern formation in [23,24]. The 2D diffusive flux in reaction-diffusion systems described by nonlinear PDE was approximated by a discrete problem for Functional (1). In particular, the general problem of plane-pattern formation based on the Turing mechanism was reduced to the geometric problem of random packing.
In higher-dimensional spaces, g kj does not always have a singularity as δ jk → 0; see capacity of a cone-plane pair for d = 3 in ( [14], Section 3.7.1). We see below that the flux between two spheres in R d for d > 3 is always regular for the linear conductivity on gradient.
This paper is devoted to extension of previous results [16,23,24] to R d , and is partially presented in arXiv:1412.7527. Section 2 presents structural approximation theory in R d . Section 3 is devoted to the construction of discrete energy in R d for the p-Laplacian. The general algorithm to minimize discrete energy (1) and energy estimations are derived in Section 4. Results are applied to the sphere-packing problem in deterministic and random statements in Section 5. Section 6 contains the discussion and concluding remarks.

Structural Approximation in R d
Let ν j (1, 2, . . . , d) be the fundamental translation vectors in space R d (d ≥ 2), forming a lattice Q = {∑ d j=1 m j ν j : m j ∈ Z}. Fundamental parallelotope Q 0 is defined by its 2 d vertices 1 2 ∑ d j=1 (±ν j ). Two points a, b ∈ R d are identified if their difference a − b = ∑ d j=1 m j ν j belongs to lattice Q. Hence, toroidal topology is introduced on Q 0 with glued opposite faces. In case R 2 , fundamental parallelogram Q 0 becomes the classical flat torus. Distance a − b between two points a, b ∈ Q 0 is introduced as where the modulus means the Euclidean distance between points a and b in R d . Consider n nonoverlapping balls D k = {x ∈ R d : |x − a k | < r} of radius r with centers a k in the cell Q 0 (see Figure 1). Let D 0 denote parallelotope Q 0 without balls |x − a k | ≤ r (k = 1, 2, . . . , n). Following [14,15] and works cited therein, we now present the main results on the nonlinear homogenization of periodic composites governed by the p-Laplacian equation (p ≥ 2) ∇ · |∇u| p−2 ∇u = 0, x ∈ D 0 . (3) Scalar function −u(x) is called the potential (sometimes the potential is taken with the plus sign), and vector function J = |∇u| p−2 ∇u is called the flux. p-Laplace Equation (3) can be smoothly continued into perforated domain D 0 + ∑ d j=1 ∑ m j ∈Z m j ν j . Take a unit vector ξ = (ξ 1 , ξ 2 , . . . , ξ d ) ∈ R d . The external potential is determined by linear function u 0 (x) = ξ · x up to an arbitrary additive constant. Vector Ω m 1 ,...,m d = ∑ d j=1 m j ξ j ν j is introduced with m j ∈ Z. It is assumed that the potential is quasiperiodic: The flux is periodic: The potential satisfies boundary conditions where t k are undetermined constants. The total normal flux through each sphere vanishes: ∂D k J(x) · n ds = 0, k = 1, 2, . . . , n, where n denotes the outward unit normal vector to sphere ∂D k . Problem (4)-(7) describes the field in the periodic composite when inclusions D k are occupied by a perfect conductor, and host conductivity is governed by Equation (3). In electrostatics, E = ∇u denotes the electric field, J = |∇u| p−2 ∇u the electric current density. Energy passing through the cell per unit volume (effective conductivity) is calculated by Formula (6.1.6) from [15] where |Q 0 | stands for volume of Q 0 . It is assumed that centers a k (k = 1, 2, . . . , n) are distributed in Q 0 in such a way that the corresponding composite is isotropic in the macroscale, i.e., the effective conductivity of the composite is expressed by a scalar λ and does not depend on unit external flux ξ.
Energy λ can be found as the minimum of the functional [15]: where space V consists of the quasiperiodic functions from the Sobolev space W 1,p (Q 0 ): Here, quasiperiodicity means that Conditions (4) and (5) are fulfilled for v. Though the definition of space V depends on external flux ξ, the energy does not because of macroscopic isotropy.
Consider a periodic graph Γ in cell Q 0 with the vertices at a k (k = 1, 2, . . . , n) and the edges corresponding to the necks between neighbor inclusions. Neighbors are defined as balls that share a common edge of the Voronoi tessellation of Q 0 with respect to their centers in toroidal topology. Designation a k ∼ a j is used for two neighboring vertices. For each fixed a k , set J k of indices is introduced for neighbor vertices. Their total number N k = #J k is called the degree of vertex a k . The Voronoi tessellation in a bounded domain is precisely described in [14,15]. The evident modification to periodic structures can be taken over. We use the term "Delaunay graph" for the constructed graph Γ following [16]. The corresponding triangulation slightly differs from commonly used Delaunay triangulation in degenerate cases. For example, consider a square and its four vertices. Traditional Delaunay triangulation has four sides of the square and one of the diagonals. In our approach, the Delaunay graph has only four sides (see Example 2 in Section 4).
The discrete network model is based on the justification that the flux is concentrated in the necks between closely spaced inclusions (see [15] for nonlinear Equation (3) and [14,15] in linear case p = 2). First, we discuss the linear case when p = 2. Let δ kj denote the gap between balls D k and D j , The computation of relative interparticle flux g kj between D k and D j (transport coefficient in terms of [14] and capacity in [25]) relies on Keller's formulas [25] g kj = g where in 3D g and in 2D Flux (13) or (14) is assigned to every edge of graph Γ.
In the nonlinear case for p > 2, we have the following formulas due to [15] in 3D: (15) and in 2D kj depends only on the geometry of the problem, i.e., on radius r and gap δ kj given by (11). [15]. Formulas (15) and (16) were taken from [15] with slight corrections (see Equations (34) and (36) and a remark below them).
In order to formulate the main asymptotic result of [15], the maximal length of edges δ = max k max j∈J k δ jk of the graph Γ is introduced. Following [14,15] we consider class D of macroscopically isotropic composites with densely packed inclusions. The term "densely packed inclusions" means that any such a location of balls has a percolation chain for δ = 0. More precisely, consider a set of nonoverlapping balls in Q 0 with δ ≥ 0 endowed with the toroidal topology. A change in δ means that the centers of balls are fixed, but their radii change and remain equal. Let for δ = 0 exist a chain of touching balls connecting the opposite faces of parallelotope Q 0 . Such a chain is called a percolation chain. (In percolation theory, a percolation chain of balls is usually defined as a chain of overlapping open balls). Macroscopic isotropy implies that each pair of the opposite faces of Q 0 possesses a percolation chain for δ = 0 (see Figure 2). Any macroscopically isotropic location not belonging to class D of densely packed inclusions can be replaced by an element of D having higher concentration. This can be performed by parallel translations of nontouching groups of balls to make them touch.
The double sum over the edges of Γ is introduced.
The main term of discrete energy has the following form: where g (0) kj is given by (15) or by (16).

Theorem 1 ([15]
). Continuous energy (8) and discrete energy (18) for p ≥ 2 in spaces R d (d = 2, 3) tend to +∞ as δ → +0. Moreover, λ can be approximated by σ for a sufficiently small δ: Theorem 1 justifies the theory of structural approximation [14,15] on the basis of a mesoscopic discretization when edges and vertices of graph Γ correspond to inclusions. This theorem was applied to the systematic study of densely packed deterministic composites [26] with given locations of inclusions, and to random two-dimensional composites [16]. Delaunay graph Γ determines summation (17) in energy (18).

Discrete Network in R d
Though the theory of structural approximation [14][15][16]26] was developed and justified in R d for d = 2, 3, its main results hold for d > 3 under the condition that interparticle flux g kj tends to infinity as δ kj → 0. The latter can be established by a formal repetition of the arguments from [14,15] in R d . Here, we partially fill this gap and calculate the main term g Following Keller [25] and explanations by Kolpakov ([14], pp. 18-25), we estimate interparticle flux g jk between two balls D j and D k with the centers located at points (0, 0, . . . , 0, ± a 2 ) of radius r < a 2 in R d (d = 2, 3, . . .). Let potential u(x) be equal to constants ±C on the spheres ( The spheres can be approximated by the paraboloids near the gap: where δ jk = a − 2r. The potential can be approximated up to O(1) by the function (see (Section 6.4.12) from [15] for 2D and 3D cases).
We now proceed to investigate the general case. Integral (24) in the spherical coordi- This can be analogously calculated to the volume of a d-dimensional ball. First, let us calculate where the Γ function is used. The integral in R is calculated by applying formula where hypergeometric function 2 F 1 is used. Therefore, Let p be a natural number greater than 2. Hereafter, in order to obtain a singular behavior for nearly touching balls, we assume that The following asymptotic formula takes place for Z → ∞: where q > d−1 2 . Formula (32) was derived by using of Mathematica . Applying (32) to (30), we obtain the main asymptotic term of g jk as δ jk → 0: Let d be an odd number. Then, (33) becomes Let d be an even number. Using formula we rewrite (33) in the form Formulas (33), (34), and (36) are slightly different with the corresponding formulas from [15] (see for instance (Section 6.2.13)) by multipliers. This can be related to the error being taken to power p in [15] (e.g., Section 6.4.15 from [15] for 2D and 3D cases) instead of the correct p − 1 from (24). Let us introduce vector t = (t 1 , t 2 , . . . , t n ) ∈ R n and matrix a = (a 1 , a 2 , . . . , a n ) ∈ R n × R d . The main term of interparticle flux g (0) jk can be written as the function of one variable a k − a j ≡ x ≥ 2r: Constant c depends on r, p, and d, and can be explicitly written using (34) and (36).
Consider the functional associated with energy and the main term of the discrete energy Theorem 1 can be extended to general space R d as follows.
Theorem 2. Consider class D of densely packed balls in R d described in Figure 2. Continuous energy (8) (see also (9)) and discrete energy (39) for p ∈ N satisfying (31) tend to +∞ as δ := max kj δ kj → +0. Moreover, λ can be approximated by σ for a sufficiently small δ: Proof of the theorem repeats the proof of Theorem 1 from [15] by its extension to periodicity cell Q 0 in space R d .

Theorem 3.
Consider class D of densely packed balls in R d . Let σ = σ(a) attain the global minimum at a location a * (δ) for sufficiently small δ. Then, optimal packing is attained at a * (0).
Let optimal packing be attained at another location a * for which concentration φ * > φ * . Location a * contains a percolation chain. Take such a radius r 0 < r for which the concentration is reduced to φ * . Then, all balls in this location a * with radius r 0 are separated from each other; hence, corresponding conductivity σ(a * ) is a finite number. However, minimal conductivity σ(a * (δ)) tends to infinity as δ → 0. This leads to a contradiction.
The theorem is proved.

Extremal Energy
Let all periodic Delaunay graphs Γ with n vertices per cell be divided into the equivalence classes of isomorphic graphs. A class of graphs are denoted as G. We are looking for the global minimum of Functional (38) in t and a in the torus topology: in a fixed class G. Minimum (41) exists since continuous function f (x) decreases and 0 ≤ f ( a k − a j ) ≤ +∞ for all a k − a j ≥ 2r. Function f (x) as a convex function for 2r ≤ x < +∞ satisfies Jensen's inequality where sum of positive numbers q i is equal to unity. Equality holds if and only if all x i are equal.
The following finite algorithm can be constructed to determine the optimal random packing for any fixed number of spheres n in periodicity cell Q 0 . The number of graph classes G is finite since the number of vertices n is fixed, and the number of edges for any vertex is bounded. Discrete energy E(t, a) algebraically depends on a finite number of variables t 1 , t 2 , . . . , t n and a 1 , a 2 , . . . , a n given on the closed finite sets. An exception arises in space R 3 when f (x) is expressed through a logarithmic function by (13). Function f (x) depends on radius r = min k =j a k − a j . Therefore, Minimization (41) in a fixed class G is realized for an elementary function of the finite number of variables. In the considered algorithm, the radius r is unknown. Optimal packing is obtained for maximal concentration c(G) = nπr 2 . The densest packing is equal to the minimum taking over a finite number of classes min G c(G).
The remaining part of the section is devoted to estimation of Minimum (41). The derived relations simplify the algorithm for a class of regular graphs in which each vertex has the same number of neighbors.
Let finite sum ∑ k,j in (41) be arranged in such a way that x i = a k − a j and Hölder's inequality states that for non-negative a i and b i This implies that Function f (x) decreases; hence, (43) and (45) give The minimum of the right-hand part of (46) on a = (a 1 , a 2 , . . . , a n ) is independently achieved on t k for max a g(a), where

Lemma 1. For any fixed class of graphs G, any local maximizer of g(a) is the global maximizer that fulfils system of linear algebraic equations
where s k can take values 0, ±1, · · · , ±d in accordance with class G. System (48) always has a unique solution a = (a 1 , a 2 , . . . , a n ) up to an arbitrary additive constant vector.
Proof. It follows from Definition (2) and the properties of the Voronoi tessellation that for some s kj which can take the values 0, ±1. The extremal points of (47) can be found from system of equations ∇ k g(a) = 0, k = 1, 2, . . . , n, where a k = (x (k) Parallelotope Q 0 endowed with toroidal topology is a compact manifold without boundary; hence, all extremal points of (47) satisfy this system. Equations (50) can be written in equivalent form (48), where The sum of all equations (48) gives an identity; hence, they are linearly dependent. Moreover, if a = (a 1 , a 2 , . . . , a n ) is a solution of (48), then (a 1 + c, a 2 + c, . . . , a n + c) is also a solution of (48) for any c ∈ R d . Consider the homogeneous system corresponding to (48): where a k belongs to cell Q 0 . System (51) can be decoupled by coordinates onto independent systems: Each l-th system (52) has only constant solutions: This follows from the consideration of the quadratic form: It is symmetric and positive semidefinite. Therefore, quadratic form (47) has a global minimum attained at a linear set of R n . All local minima of (47) coincide with the global minimum. The quadratic form attains the global minimum at Set (53). Therefore, all solutions of (52), and hence of (51), are only constants. Then, System (48) has only one condition of solvability for the right-hand part, which is fulfilled. Therefore, System (48) always has a unique solution up to an arbitrary additive constant.
The lemma is proved.
The l-th coordinates x (k) l of point a k (l = 1, 2, . . . , d and k = 1, 2, . . . , n) are located on real axis R (coincidence is permitted). For definiteness, let the n-th point to be fixed as x (n) l = 0 (l = 1, 2, . . . , d). Then, real numbers x (k) l satisfy uniquely solvable system where ν ( ) l denotes the l-th coordinate of basic vector ν . Solutions of (55) can be presented as linear combination where x (k ) l satisfy system All computations can be separately performed for l = 1, 2, . . . , d. As a result, we obtain expressions for optimal a k in class G The next numerical problem consists in determining of such a basis {ν } d =1 that possesses d percolation chains connecting the opposite faces of parallelotope Q 0 . This is a combinatorial problem that is not discussed in the present paper.
Lemma 1 gives energy bound min t E(t, a), with a satisfying (48) for the packing problem in fixed class G. This bound yields the exact global minimum when Inequality (46) becomes an equality. Let some of the differences |t k − t j | vanish, and the rest of the differences |t k − t j | are considered to be equal. If, for equal differences |t k − t j |, differences a k − a j are also equal, we have the desired global minimum. Such a case takes place, for instance, for graphs with the same degree of vertices (kissing number for packing) N = N k for all vertices. Example 1. Consider class of graphs A 2 in R 2 , with n = m 2 (m ∈ N) vertices containing hexagonal lattice A 2 . This class is determined by kissing number N k = 6 for all disks. The first step in the proof of optimality of the hexagonal lattice is proof of equality N k = 6. Here, we refer to [7]. Hence, if we restrict ourselves by class A 2 , we do not arrive to a set of graphs containing optimal packing.
The second step consists in the direct check that the regular hexagonal lattice satisfies (48). Then, the simple external flux is applied, and the corresponding |t k − t j | and a k − a j are calculated. Consider the regular hexagonal lattice generated by vectors The radius of disks equals r = 1 2 . Periodicity cell Q 0 is determined by vectors (59). It is convenient to consider A 2 as a set of perpendicular layers to axis x 2 . The number of layer q ∈ Z is introduced, where the q-th layer consists of disks of which the centers have the x 2 coordinate equal to √ 3 2 q (see Figure 3). Let the external flux be directed along axis x 2 and determined by external potential u 0 (x) = √ 3 2 x 2 . Then, the continuous and discrete potentials take value q on all disks of the q-th layer. Therefore, difference |t k − t j | is equal to zero if a k and a j belong to the same layer, and equals unity if a k and a j belong to neighboring layers. Difference a k − a j takes the same value of 2r = 1 for touching disks from neighboring layers. Inequality (46) becomes an equality in this case. It follows from Lemma 1 that global Minimum (41) on class of graphs A 2 is attained at the regular hexagonal graph. Optimal location a * satisfies (48), thereby depending on basic vectors ν 1 and ν 2 . If we take other basic vectors, different from (59), we may arrive at location a * , which is not macroscopically isotropic.
This example gives an alternative "physical" proof of the optimal 2D packing attained for the hexagonal array, but the first step of the proof contains equation N k = 6, established by geometrical arguments.

Example 2.
Consider class of graphs Z 2 in R 2 , with n = m 2 (m ∈ N) vertices containing square lattice Z 2 . This class is determined by contact number N k = 4 for all disks.
Consider the regular square lattice generated by vectors ν 1 = (m, 0) and ν 2 = (0, m). Periodicity cell Q 0 is determined by these basic vectors. In accordance with structural approximation theory, the corresponding Delaunay graph consists of edges parallel to axes x 1 and x 2 . Lattice R 2 consists of layers perpendicular to axis x 1 . The q-th layer consists of disks of which the centers have the x 2 -coordinate be equal to q. The external flux is determined by external potential −u 0 (x) = −x 2 . Then, the continuous and discrete potentials take value q on all disks of the q-th layer. Lemma 1 implies that global Minimum (41) on class of graphs Z 2 is attained at the regular square graph.
Regular square lattice Z 2 gives optimal packing in class Z 2 . Class Z 2 is geometrically unstable, since any perturbation of a vertex yields flipping, and changes the structure of Z 2 . Example 3. Consider class of graphs A 3 in R 3 with n = m 3 (m ∈ N) vertices containing regular face-centered cubic lattice A 3 . Since this lattice consists of layers, we can consider the external flux to be directed perpendicular to the layers. Then, the potential is constant in each layer, and differences a k − a j = 2r are constant, where neighboring points a k , a j belong to the neighboring layers. This implies that lattice A 3 reaches optimal packing in class A 3 . Optimal location a * in class A 3 depends on the basic vectors prescribed to FCC lattice A 3 since a * satisfies System (48), of which the right-hand side depends on the basic vectors. The class of graphs containing the lattice corresponding to the hexagonal close-packing (HCP) determines another optimal location.
The class of graphs containing A 3 with arbitrary fixed basic vectors determines its optimal location a * through System (48). This location a * depends on the basic vectors and must determine a macroscopically isotropic structure.

Optimization in Various Classes of Graphs
We now proceed to discuss the described algorithm in general form at the beginning of Section 4. Section 5.1 contains remarks concerning packing problems in the deterministic statement. Section 5.2 contains an example of the optimal packing of balls in the random statement. This example illustrates how the algorithm can be adjusted to random optimal packing.

Deterministic Problem
In order to define the set of all optimal packing, it is sufficient to describe all classes G and to construct the optimal graph for every G by Minimization (41). This problem is referred to as a finite combinatorial optimization problem in fixed space R d when the number of balls per cell n is finite. The set of graph classes can be reduced in the deterministic statement by considering only such a class G of which the optimal graph has equal edges. Unfortunately, this property of equal edges is not known at the beginning. The optimization procedure is performed in a class G, and afterwards to check the lengths of its edges. The considered numerical examples in R 2 suggest that only regular arrays (square, hexagonal, etc.) have such a property. The suggestion is formulated below in the form of a hypothesis.
Let us fix class of graphs G. Delaunay graph Γ ⊂ G of which the vertices satisfy (48) is called optimal in considered class G. It follows from the theory of structural approximation [14,15] that the energy (41) of every optimal Delaunay graph Γ must only contain singular terms in δ. In limit case δ → 0, we call edges of the Delaunay graph having length 2r the principal edges. Adding or deleting longer edges to or from graph Γ can change the energy only up to a value of order O(1). Subgraph Γ 0 of Γ containing only the principal edges is called the principal graph. It is also convenient to consider class S in which the number of vertices n is represented in form n = n d 1 , where n 1 is a natural number. This condition can be always fulfilled by the addition of d adjacent cells along a fundamental vector. For definiteness, let n 1 = 2M + 1 be an odd number. Below, designations Q 0 , n, etc. are used for the modified class of graphs S. Consider set of periodically located lattice points Let any point of L be connected to other κ d closest points of L. Such a regular graph is denoted by Γ d . This Γ d graph determines class G d .

Hypothesis 1.
Any principal graph Γ 0 can be extended to a regular graph Γ d where each vertex has the same degree κ d .

Hypothesis 2.
Any graph Γ 0 is a subgraph of a graph Γ ∈ G d .
The main point of the hypothesis is based on the eventual possibility to replace one graph by another by addition and deleting edges greater than δ, since such a replacement does not change the energy up to O(1).
If the hypothesis is true, layered lattices (60) solve the optimal packing problem in the corresponding R d for which the kissing number is known.

Stochastic Problem
Here, we demonstrate by an example how the algorithm can be adjusted to optimal random packing.
Graphs (A , E ) and (A , E ) are isomorphic. Therefore, graph (A , E ) is one of the optimal graphs that can be obtained as the optimal result of a random packing.   It is interesting to compare our optimal structures with the known structures. As an example, consider the patterns in Figure 6 simulated in [27]. First, a boundary-value problem in the square for the Turing system was numerically solved, and the concentration distribution was found as a function U(x, y). Gray spots in Figure 6 show the hills of levels U(x, y) = constant. In the present paper, these gray spots are processed by the following method. First, we find the gravity centers of the spots and consider them to be vertices of Delaunay graph. Next, we consider the inner points of this graph and calculate new set of points a k by Lemma 3. These points are the centers of the blue disks displayed in Figure 6. The edges of the optimal Delaunay graph are shown. Vertices of degrees 5, 6, and 7 are visible. This is the optimal-packing graph corresponding to the considered class of graphs.

Discussion
In this paper, we show that the solution to the physical problem of minimal energy (effective conductivity) implies a solution to the geometrical problem of optimal packing. A general algorithm is proposed to determine the densest packing in different statements. Lemma 1 yields bounds for the classical optimal-packing problem. For graphs with the same degree of vertices, these bounds are exact, hence solving the optimal-packing problems within the considered classes. Therefore, we can "algorithmically" solve the optimal-packing problem in any fixed class of graphs G. Let us discuss these classes of graphs. Though the number of spheres per cell n is arbitrary and can tend to infinity, the periodicity cell is fixed at the beginning; hence, optimal location a * depends on the basic vectors. It also follows from the observation that the right-hand part of System (48) contains the basic vectors. For a sufficiently large n, the majority of equations of (48) are homogeneous, and only "boundary" equations are inhomogeneous. Their number is of order √ n. It is interesting to investigate the asymptotic dependence of a * on the basic elements as n → ∞. This shows, for instance, the effectiveness of the hexagonal-type packing in the cubic cell in R d . We do not know a priori that solution a * of System (48) forms a macroscopically isotropic structure for a given cell Q 0 and for a fixed structure of G. However, a simple necessary condition of isotropy can be checked: location a * must contain d percolation chains connecting the opposite face of parallelotope Q 0 (see Figure 2).
A straightforward verification that class G yields optimal packing can be proposed for a fixed n. Below, the analogy with the laminated lattices ( [5], Chapter 6) is used. First, for clarity, consider a 1D case. Then, there exists exactly one class of graph Z, including regular 1D lattice Z. A graph from Z consists of sequentially located edges along the axis. It follows from System (48) that each point a k lies in the middle of segment (a k−1 , a k+1 ) that immediately implies points a k form lattice Z. The above scheme can be extended to R d . Ultimately, one must perform such a procedure for all classes of graphs, and take the best packing in the deterministic statement. Though optimal packing can be constructed in every class of isomorphic graphs, their comparison to select the best packing seems to be a hard combinatorial problem. Vertices in the considered class can be considered to be a perturbation of the rigid laminated lattice. Therefore, Theorem 2 from [5] (p. 164) can be extended to the classes of graphs containing the corresponding laminated lattices.
The constructed algorithm can be effectively applied to packing in classes G which can be considered to be a randomly chosen structure. Corresponding concentration c(G) of balls attains the maximal value in class G on optimal graph Γ. Random packing can be considered to be a dynamical process of finding ball locations by obeying some stochastic and deterministic rules. Such a process terminates if the reached packing concentration cannot increase in the framework of the considered algorithm. All these terminated locations belong to the set of optimal locations within the classes of isomorphic graphs.