Loop corrections for hard spheres in Hamming space
Abolfazl Ramezanpour Saman Moghimi-Araghi



 


aramezanpour@gmail.com
Department of Physics, College of Science, Shiraz University, Shiraz 71454, Iran
Medical Systems Biophysics and Bioengineering, Leiden Academic Centre for Drug Research, Faculty of Science, Leiden University, Leiden, The Netherlands	

samanimi@sharif.edu
Physics Department, Sharif University of Technology, P.O. Box 11155-9161 Tehran, Iran




§ ABSTRACT 
We begin with an exact expression for the entropy of a system of hard spheres within the Hamming space. This entropy relies on probability marginals, which are determined by an extended set of Belief Propagation (BP) equations. The BP probability marginals are functions of auxiliary variables which are introduced to model the effects of loopy interactions on a tree-structured interaction graph. We explore various reasonable and approximate probability distributions, ensuring they align with the exact solutions of the BP equations. Our approach is based on an ansatz of (in)homogeneous cavity marginals respecting the permutation symmetry of the problem. Through thorough analysis, we aim to minimize errors in the BP equations. Our findings support the conjecture that the maximum packing density asymptotically conforms to the lower bound proposed by Gilbert and Varshamov, further validated by the solution of the loopy BP equations.




Loop corrections for hard spheres in Hamming space
    Saman Moghimi-Araghi
    September 9, 2024
==================================================





§ INTRODUCTION

Finding the maximum packing density of hard spheres is a challenging problem in physics and coding theory, especially in high dimensions <cit.>. For example, exact solutions are only known for very small and specific dimensions like 1, 2, 3, 8, and 24<cit.>. Currently, there's a significant gap between the best lower and upper bounds on the maximum packing density in very high dimensions <cit.>. One approach to tackle this problem is to start with a crystalline structure made of elementary cells with a flexible basis. These cells form a periodic density distribution of the spheres, determined by the shape of the elementary cell and the arrangement of the spheres in the basis. By adjusting these parameters, we can find an upper bound for the maximum packing density that satisfies certain equations for the spatial density distribution.

On the other hand, the problem can be understood as mapping to a physical system of interacting particles, wherein the relevant packing configurations, weighted by the Boltzmann factor, are examined as the number density of hard spheres increases. Mean-field approximations are commonly employed to construct a phase diagram of the system, particularly in very high (infinite) dimensions where this approximation is expected to perform well if handled with care <cit.>.

In the following, we will focus on the packing problem within the binary Hamming space. At zero temperature, this translates into a constraint satisfaction problem where any two spheres cannot overlap. The Bethe approximation (cavity method) has been utilized to estimate the maximum packing density in the Hamming space. Nevertheless, the degree to which these results accurately reflect the behaviour as the dimensionality approaches infinity in this fully-connected model with numerous interconnected interactions remains uncertain.

In this study we try a different approach by first introducing some auxiliary variables to map the original problem to an extended problem with a tree interaction graph where the Bethe approximation is exact <cit.>. Then we write an exact expression for the entropy in terms of the (unique) solution to the Bethe, or Belief Propagation (BP), equations <cit.>. To solve these equations and estimate the system's entropy, we make reasonable and manageable approximations.

The structure of the paper is as follows. We begin with a more precise statement of the problem and a summary of the results.  In Sec. <ref> we present the mapping to a tree interaction graph and write the associated Bethe equations and the entropy.   
Section <ref> is devoted to finding approximate BP solutions, starting with a naive homogeneous (liquid) solution and ending with inhomogeneous solutions which respect the permutation symmetry of the problem. The concluding remarks and some details of the results are given in Sec. <ref> and Appendices <ref>-<ref>, respectively. 





 §.§ The problem statement

Consider the binary Hamming space in n dimensions. In other words, an n-dimensional discrete space where each dimension can take either zero or one. We want to place N hard spheres of diameter d on the points of this lattice in way that the spheres do not overlap. However, if they have no overlap, they are considered completely non-interacting. Placing a hard sphere at any point σ⃗_i∈ (0,1)^n produces a forbidden or occupied space which is the set of discrete points within the ball of radius d at which the other hard spheres cannot be located. This volume is given by

    V_d(σ⃗_i)={σ⃗_j:D_ij<d},

as the distance between any two points in Hamming space is equal to 

    D_ij=∑_a=1^n(σ_i^a-σ_j^a)^2.

And if the coordinates of any two points are not identical in exactly q dimensions, their distance from each other is equal to D_ij=q. Thus, the number of inaccessible points is equal to 

    V_d=|V_d(σ⃗_i)|=∑_l=0^d-1C(l:n),

where C(l:n)=n!/(l!(n-l)!). Note that this number differs from the volume of a spherical shell with radius d.


Let's write the partition function of this system. Since the spheres are considered non-interacting, the partition function of the system is just the number of possible configurations for the packing. Any arbitrary configuration of these N spheres in the Hamming space is represented by 

    σ⃗={σ⃗_i ∈ (0,1)^n: i=1,⋯,N}.

and thus, the partition function is given by 

    Z_n(N,d)=∑_σ⃗∏_i<j𝕀(D_ij≥ d)=e^S_n(N,d).

where the indicator function 𝕀(C)=1 if constraint C is satisfied, otherwise 𝕀(C)=0.
As mentioned, the partition function represents the number of possible configurations for packing. Naturally, in a space of a certain size, if we increase the number of spheres, we reach a point where no valid configuration can be found and the spheres would surely have overlaps. We call the maximum number of spheres for which a valid configuration is found N_max. In other words, we must have 

    Z_n(N_max,d)> 0,1cm Z_n(N_max+1,d)= 0.

This quantity and other properties of the system depend on the space dimension n and the diameter considered for the spheres d. Increasing n enlarges the space, while increasing d reduces the number of states. We are interested in the case where d, n →∞ while their ratio remains constant and equal to δ.


Within the Bethe approximation and assuming the replica symmetry <cit.>, the entropy S_n(N,d)=ln Z_n(N,d) is estimated by

    S_n^LBP(N,d)=∑_iΔ S_i-∑_i<jΔ S_ij.

The entropy contributions of the nodes (spheres) Δ S_i and edges (interactions) Δ S_ij depend on the solutions to the Bethe equations

    η_i→ j(σ⃗_i) ∝∏_k i,j(∑_σ⃗_k𝕀(D_ik≥ d)η_k→ i(σ⃗_k)).

Here the cavity marginal η_i→ j(σ⃗_i) is the probability of finding sphere i in position σ⃗_i in the absence interaction with sphere j. Given a fixed point of the above equation, one computes the entropy contributions

    e^Δ S_i   = ∑_σ⃗_i∏_j i(∑_σ⃗_j𝕀(D_ij≥ d)η_j→ i(σ⃗_j)),
    
    e^Δ S_ij   = ∑_σ⃗_i,σ⃗_j𝕀(D_ij≥ d)η_i→ j(σ⃗_i)η_j→ i(σ⃗_j).

 
For a tree interaction graph, the Bethe equations with the unique solution η_i→ j(σ⃗_i)=1/2^n provide the exact entropy:

    S_n^T(N,d)=Nln(2^n)+(N-1)ln (1-V_d/2^n),

which is greater than zero for any d<n, irrespective of the tree structure.

The entropy obtained by the loopy Belief Propagation (LBP) equations with a homogeneous or liquid solution which respects the translation symmetry of the problem (η_i→ j(σ⃗_i)=1/2^n) is:

    S_n^LBP(N,d)=Nln(2^n)+N(N-1)/2ln (1-V_d/2^n).

The maximum N here is given by

    N_max^LBP=1-2ln 2^n/ln (1-v_d)→ (2ln 2)n/v_d,

where v_d=V_d/2^n. 

Note that the Gilbert-Varshamov (GV) lower bound states that

    N_max≥2^n/V_d=1/v_d=N_max^GV.

It seems that even other solutions to the loopy Bethe equations do not result in an exponentially larger maximum number of spheres than the GV lower bound <cit.>.
     



 §.§ Summary of the results

The main results of this paper are listed below:



  * The interaction graph of spheres is partitioned into a tree interaction graph and a set of induced loopy interactions. The effect of loopy interactions is represented by messages passing along the edges of the tree graph. We formulate the BP equations with these auxiliary variables in an extended space to obtain an exact expression for the entropy of the packing problem in the Hamming space.


  * A naïve ansatz for the BP probability marginals reproduces the entropy, previously obtained by the loopy Bethe equations (assuming replica symmetry or one step of replica symmetry breaking). This asymptotically coincides with the GV lower bound for the maximum number of spheres.


  * We observe that the aforementioned naïve ansatz asymptotically satisfies (“on average”) the extended BP equations. Similar results are obtained with other reasonable candidates for homogeneous (liquid) BP marginals, which are expected to be closer to the exact solution of the extended BP equations.


  * Numerically, we observe that an initially inhomogeneous solution to the BP equations approaches the homogeneous one as the number of spheres increases, starting from two neighboring spheres localized in the Hamming space. However, this does not conclusively prove that there is no packing density that asymptotically exceeds the GV lower bound.






§ AN EXACT TREE REPRESENTATION

Consider a connected tree graph T of N nodes and the associated local branches or cavity trees T_i→ j. More precisely, T_i→ j is the subgraph that is obtained after removing edge (ij) with root node i. This can recursively be defined as follows

    T_i→ j=i ∪_k∈∂ i∖ j T_k→ i.

Here ∂ i denotes the set of neighbors of node i in the graph. Figure <ref> displays such a tree interaction graph.  
The size of such cavity tree is denoted by

    N_i→ j=|T_i→ j|.






For a given configuration σ⃗, let us define the cavity messages h_i→ j to represent the space occupied by all the spheres in the cavity tree T_i→ j∖ i, that is

    h_i→ j=∪_k∈ T_i→ j∖ i V_d(σ⃗_k).

Or, in terms of the other incoming messages 

    h_i→ j=∪_k∈∂ i ∖ j( h_k→ i∪ V_d(σ⃗_k) ).

Equivalently, the messages h_i→ j can be considered as the positions of all the spheres in T_i→ j except sphere i.
These message are defined to have access to the positions of all the spheres locally at each node of the tree graph T.


Now with the extended set of variables σ⃗ and 𝐡={h_i→ j: i=1,…,N, j∈∂ i}, we rewrite the partition function of the sphere packing problem,

    Z_n(N,d)=∑_σ⃗∏_(ij)∈ T𝕀(D_ij≥ d) 
    ×(∑_𝐡∏_i ∏_j∈∂ i𝕀(V_d(σ⃗_i)∩ h_j→ i=∅)𝕀(h_i→ j=∪_k∈∂ i∖ j h_k→ i∪ V_d(σ⃗_k))).

The constraints are to ensure that the messages h_i→ j satisfy the necessary equations and at the same time any two spheres have a distance larger than or equal to d. Note that there is only one solution to the constraints on the messages h_i→ j for a given configuration σ⃗ of the spheres.    

The Bethe equations for the above problem are recursive equations for the cavity probability distributions μ_i→ j(σ⃗_⃗i⃗,h_i→ j:σ⃗_⃗j⃗,h_j→ i) along the directed edges of the tree graph. Here μ_i→ j(σ⃗_⃗i⃗,h_i→ j:σ⃗_⃗j⃗,h_j→ i) is the probability of having (σ⃗_⃗i⃗,h_i→ j) as the position of sphere i and the message h_i→ j conditioned on the values of (σ⃗_⃗j⃗,h_j→ i). In a tree graph, the incoming variables are independent of each other, therefore <cit.>, 

    μ_i→ j(σ⃗_⃗i⃗,h_i→ j:σ⃗_⃗j⃗,h_j→ i) ∝𝕀(σ⃗_⃗i⃗,h_j→ i)
    ×∑_{σ⃗_k,h_k→ i:k∈∂ i∖ j}𝕀(h_i→ j)∏_k∈∂ i∖ j[𝕀(D_ik) 𝕀(σ⃗_⃗i⃗,h_k→ i)𝕀(h_i→ k)μ_k→ i(σ⃗_⃗k⃗,h_k→ i:σ⃗_⃗i⃗,h_i→ k)],

where for brevity's sake, we defined

    𝕀(σ⃗_⃗i⃗,h_k→ i)    =𝕀(V_d(σ⃗_i)∩ h_k→ i=∅),
    𝕀(D_ik)    =𝕀(D_ik≥ d),
    𝕀(h_i→ k)    =𝕀(h_i→ k=∪_j∈∂ i∖ k h_j→ i∪ V_d(σ⃗_j)).

In the following, the messages h_i→ j are called internal messages and the cavity probabilities μ_i→ j are called the external BP messages. In practice, the BP equations are solved by iteration in a random sequential way starting with initial BP messages μ_i→ j. Here the tree structure of the interaction graph insures that there a unique solution to the above equations.   


It is straightforward to start from the partition function Z_n(N,d) of the tree interaction graph T and relate the free entropy S_n(N,d) to the BP cavity marginals <cit.>. This can be written in term of the variables (nodes) and interactions (edges) contributions to the entropy 

    S_n(N,d)=∑_i Δ S_i-∑_(ij)∈ TΔ S_ij,

where

    e^Δ S_i=∑_σ⃗_i∑_{σ⃗_j,h_j→ i:j∈∂ i}∏_j∈∂ i[𝕀(D_ij) 𝕀(σ⃗_⃗i⃗,h_j→ i)𝕀(h_i→ j)μ_j→ i(σ⃗_⃗j⃗,h_j→ i:σ⃗_⃗i⃗,h_i→ j)],

and

    e^Δ S_ij= ∑_σ⃗_i,h_i→ j,σ⃗_j,h_j→ i𝕀(D_ij) μ_i→ j(σ⃗_⃗i⃗,h_i→ j:σ⃗_⃗j⃗,h_j→ i)μ_j→ i(σ⃗_⃗j⃗,h_j→ i:σ⃗_⃗i⃗,h_i→ j).



In the following we shall work with a tree structure T which is represented by a chain of interacting spheres.
This allows us to simplify the BP equations and obtain simpler expressions for the BP cavity marginals and the entropy.      

  



 §.§ Chain representation

Here we assume that the spheres are arranged in a chain from i=1,…,N. Besides the σ⃗_i we need the messages h_i passing from left to right. Thus on each edge (i,i+1) we have three constraints: D_i,i+1≥ d, h_i+1=h_i∪  V_d(σ⃗_i), and  h_i ∩ V_d(σ⃗_i+1)=∅. Let the indicator function 𝕀_i,i+1 represent these constraints. Then, the left-to-right BP messages μ_i(σ⃗_i,h_i) and right-to-left BP messages ν_i(σ⃗_i,h_i) are given by the following BP equation,

    μ_i(σ⃗_i,h_i) ∝∑_σ⃗_i-1,h_i-1𝕀_i-1,iμ_i-1(σ⃗_i-1,h_i-1),
    ν_i(σ⃗_i,h_i) ∝∑_σ⃗_i+1,h_i+1𝕀_i,i+1ν_i+1(σ⃗_i+1,h_i+1).

The boundary messages are

    μ_1(σ⃗,h_1)    = 1/2^nδ_h,∅,
    ν_N(σ⃗,h_N)    = 1/2^n1/2^n(N-1).

Here we are considering h_i as the set of positions of spheres j=1,…,i-1.
In this way, the right-to-left messages are uniform distributions. Also away from the boundary

    ν_i(σ⃗,h) = 1/2^n1/2^n(i-1).



The entropy within the chain representation reads

    S_n(N,d)=∑_i=1^N Δ S_i-∑_i=1^N-1Δ S_i,i+1,

where now the node contribution is

    e^Δ S_i=∑_σ⃗_i,h_i(∑_σ⃗_i-1,h_i-1𝕀_i-1,iμ_i-1(σ⃗_i-1,h_i-1))(∑_σ⃗_i+1,h_i+1𝕀_i,i+1ν_i+1(σ⃗_i+1,h_i+1)),

and the edge contribution  is given by

    e^Δ S_i,i+1= ∑_σ⃗_i,h_i,σ⃗_i+1,h_i+1𝕀_i,i+1μ_i(σ⃗_i,h_i)ν_i+1(σ⃗_i+1,h_i+1).



In the following we shall work with the chain representation.




§ APPROXIMATE SOLUTIONS OF THE BETHE EQUATIONS





 §.§ A naive approximation of the BP messages

Let us write the BP probability marginals μ_i(σ⃗_i,h_i) in term of probability of position of sphere i, p_i(σ⃗_i) and a conditional probability q_i(h_i:σ⃗_i), 

    μ_i(σ⃗_i,h_i)= p_i(σ⃗_i) q_i(h_i:σ⃗_i),

with h_i={σ⃗_j=1,…,i-1}.
Then, by symmetry we take a uniform distribution for σ⃗_i,

    p_i(σ⃗_i)=1/2^n.

We also assume that the probability distribution of the other spheres is factorized and uniform

    q_i(h_i:σ⃗_i)=∏_j=1^i-1𝕀(D_i,j≥ d)/2^n-V_d.


In this way, for the node contributions to the entropy from Eq. <ref> we obtain

    e^Δ S_i=1/2^n+ni∑_l_1,l_2,l_12=d^n C(l_12:n)
    Ω(l_1,l_2:l_12)
    ×(1-2V_d-O_01(l_1)-O_02(l_2)-O_12(l_12)+2O_012(l_1,l_2,l_12)/2^n-V_d)^i-2.

The other terms in the entropy come from the interactions along the chain and are obtained from Eq. <ref>

    e^Δ S_i,i+1= 1/2^n+ni∑_l=d^n C(l:n)
    (1-V_d-O_01(l)/2^n-V_d)^i-1.

Here O_ij(l) is the overlap of two spheres of radius d at distance l. O_ijk(l_i,l_j,l_ij) is the overlap of three spheres of radius d when (i,j) have distance l_ij and the other sphere is at distances l_i,l_j form i and j, respectively.
The function Ω(l_i,l_j:l_ij) is the number of possible points for the third sphere given the distances. See Fig. <ref> for a schematic representation of these quantities.






Let us take the limit n,d→∞ with δ=d/n finite. All distances are scaled with n, for instance r_ij=l_ij/n. Then v_d=V_d/2^n=e^-n[ln(2)-H(δ)], with H(δ)=-δln(δ)-(1-δ) ln(1-δ) for the binary entropy, is an exponentially small quantity. In Appendix <ref> we obtain the asymptotic entropy contributions of the nodes and edges in the total entropy,  

    Δ S_i    ≃ nN(3/2ln 2-ṽ_d),
    Δ S_i,i+1   ≃ nN1/2(ln 2-ṽ_d),
    
    S^naive   =Δ S_i-Δ S_i,i+1≃ nN(ln 2-1/2ṽ_d).

where we defined ṽ_d=Nv_d/n. Therefore, the entropy vanishes at

    N_max^naive=(2ln 2)n/v_d,

at the same point that the loopy BP entropy S_n^LBP(N,d) in Eq. <ref> goes to zero.
 
 


 §.§ Beyond the naive approximation

Note that the probability distribution of the spheres is invariant under the permutation of the spheres. Therefore, by de Finetti theorem <cit.> the probabilities q_i(h_i:σ⃗_i) can well be approximated by a convex combination of factorized distributions,   

    q_i(h_i:σ⃗_i)=∑_αc_i(α)∏_j=1^i-1f_i,j^α(σ⃗_j:σ⃗_i),

where

    ∑_αc_i(α)    =1,
    ∑_σ⃗_jf_i,j^α(σ⃗_j:σ⃗_i)    =1.

Recall that h_i is equivalent to the set of positions {σ⃗_j=1,…,i-1}. 
Let us for simplicity consider only one of the factorized terms. Then the BP equations are

    p_i(σ⃗_i)∏_j=1^i-1f_ij(σ⃗_j:σ⃗_i)
    ∝∑_σ⃗_i-1,h_i-1𝕀_i-1,i p_i-1(σ⃗_i-1)∏_j=1^i-2f_i-1,j(σ⃗_j:σ⃗_i-1).

It means that 

    p_i(σ⃗_i)f_i,i-1(σ⃗_i-1:σ⃗_i)∏_j=1^i-2f_i,j(σ⃗_j:σ⃗_i)
    
    = p_i-1(σ⃗_i-1)/z_i𝕀(D_i,i-1≥ d)∏_j=1^i-2𝕀(D_i,j≥ d)f_i-1,j(σ⃗_j:σ⃗_i-1),

where the normalization factor z_i is

    z_i= ∑_σ⃗_i,h_i p_i-1(σ⃗_i-1)𝕀(D_i,i-1≥ d)∏_j=1^i-2𝕀(D_i,j≥ d)f_i-1,j(σ⃗_j:σ⃗_i-1).

Taking p_i(σ⃗_i)=p_i-1(σ⃗_i-1)=1/2^n, the above equation suggests the following solution

    f_i,j(σ⃗_j:σ⃗_i)=𝕀(D_i,j≥ d)/2^n-V_d,

for j=1,…,i-2, and

    f_i,i-1(σ⃗_i-1:σ⃗_i)=𝕀(D_i,i-1≥ d)g_i,i-1(σ⃗_i-1:σ⃗_i),

for j=i-1. Given the {f_i-1,j:j=1,…,i-2}, the function g_i,i-1(σ⃗_i-1:σ⃗_i) should be chosen according to the following constraints,

    ∑_σ⃗_i-1g_i,i-1(σ⃗_i-1:σ⃗_i)𝕀(D_i,i-1≥ d)    =1,
    
    g_i,i-1(σ⃗_i-1:σ⃗_i)
       =(2^n-V_d)^i-2/z_i∏_j=1^i-2f_i-1,j(σ⃗_j:σ⃗_i-1).

The equations can not be true for any configuration. To overcome the problem of dependence on σ⃗_j, we can simply average over these variables using a uniform measure, leading to

    g_i,i-1(σ⃗_i-1:σ⃗_i)
    = (2^n-V_d)^i-2/z_i⟨∏_j=1^i-2f_i-1,j(σ⃗_j:σ⃗_i-1)⟩
    
    =(1-v_d)^i-2/z_i∑_σ⃗_i-2 g_i-1,i-2(σ⃗_i-2:σ⃗_i-1)𝕀(D_i-1,i-2≥ d)
    
    =(1-v_d)^i-2/z_i.

In this way z_i=2^n(1-v_d)^i-1, thus

    f_i,i-1(σ⃗_i-1:σ⃗_i)=𝕀(D_i,i-1≥ d)/2^n-V_d.

It means that all the f_i,j for j=1,…,i-1 are given by the expression we used in the naive approximation of the messages. Moreover, the f_i,j are the same function, as expected from the permutation symmetry of the problem.

Note that given the above f_i,j, according to Eq. <ref>,

    z_i=1/(2^n-V_d)^i-2∑_l_12=d^nC(l_12:n)
    (∑_l_1,l_2=d^n
    Ω(l_1,l_2:l_12))^i-2.

For n→∞, after ignoring the exponentially small overlaps,

    z_i≃1/(2^n-V_d)^i-2(2^n-V_d)(2^n-2V_d)^i-2=2^n(1-v_d)(1-V_d/2^n-V_d)^i-2≃ 2^n(1-v_d)^i-1,

which is consistent with the expression we obtained for this quantity after Eq. <ref>. 





 §.§ Minimizing the BP errors

In this subsection, we look for a better solution for f_i,j(σ⃗_j:σ⃗_i)=𝕀(D_i,j≥ d)g_i,j(σ⃗_j:σ⃗_i) by minimizing the distance

    ℒ[{g_i,j=1,…,i-1}]=∑_{σ⃗_j=1,…,i-1:D_i,j≥ d}( ∏_j=1^i-1g_i,j-1/z_i∏_j=1^i-2𝕀(D_i-1,j≥ d)g_i-1,j)^2,

constrained with

    ∑_σ⃗_j:D_i,j≥ dg_i,j=1.

The stationary equations for the g_i,j(σ⃗_j:σ⃗_i) result in

    λ_i,j=g_i,j(σ⃗_j:σ⃗_i)∏_k=1:k≠ j^i-1(∑_σ⃗_k:D_i,k≥ dg_i,k^2)
    
    -1/z_i∑_σ⃗_i-1:D_i,i-1,D_i-1,j≥ dg_i,i-1g_i-1,j∏_k=1:k≠ j^i-2(∑_σ⃗_k:D_i,k,D_i-1,k≥ dg_i,kg_i-1,k),

for j=1,…,i-2, and

    λ_i,i-1=g_i,i-1(σ⃗_i-1:σ⃗_i)∏_k=1^i-2(∑_σ⃗_k:D_i,k≥ dg_i,k^2)-1/z_i∏_k=1^i-2(∑_σ⃗_k:D_i,k≥ d,D_i-1,k≥ dg_i-1,kg_i,k).

The Lagrange multipliers λ_ij are to ensure the normalization constraints.
Let us assume that the g_i,j depend only on the Hamming distances D_i,j.
By symmetry we also assume that g_i,j=g_i is the same for all j=1,…,i-1 and consider only the later equation for g_i,i-1. In this way, from Eq. <ref> we get

    g_i(l)=λ_i+1/z_i(∑_l_1,l_2=d^n Ω(l_1,l_2:l)g_i-1(l_1)g_i(l_2))^i-2/(∑_l_1=d^nC(l_1:n)g_i(l_1)^2)^i-2,


    z_i=∑_l=d^nC(l:n)
    (∑_l_1,l_2=d^nΩ(l_1,l_2:l)g_i-1(l_1) )^i-2,

and for the Lagrange multiplier

    λ_i=1/2^n-V_d(∑_l_1=d^nC(l_1:n)g_i(l_1)^2)^i-2
    
    -1/z_i(2^n-V_d)∑_l=d^nC(l:n)(∑_l_1,l_2=d^n Ω(l_1,l_2:l)g_i-1(l_1)g_i(l_2))^i-2.

The above equations can be rewritten in a compact form as  

    g_i(l)=1/2^n-V_d+(X_i(l)-⟨ X_i(l)⟩),

where we defined

    X_i(l)=1/z_i( 
    ∑_l_1,l_2=d^n Ω(l_1,l_2:l)g_i-1(l_1)g_i(l_2)/∑_l_1=d^nC(l_1:n)g_i(l_1)^2)^i-2,

and

    ⟨ X_i(l)⟩=1/2^n-V_d∑_l=d^nC(l:n)X_i(l).


Given the g_i-1(l), the equations can be solved by iteration for g_i(l), starting from g_2(l)=1/(2^n-V_d).
We see that the trial solution g_i(l)=1/(2^n-V_d), gives 

    ⟨ X_i(l)⟩   =1/2^n-V_d,

and

    X_i(l)=(1-v_d1-o_ij(l)/1-v_d)^i-2/∑_l=d^nC(l:n)(1-v_d1-o_ij(l)/1-v_d)^i-2,

where v_d=V_d/2^n and o_ij=O_ij/V_d are exponentially small quantities. Approximating (1-o_ij(l))/(1-v_d)≈ 1 is therefore consistent with the solution g_i(l)=1/(2^n-V_d) obtained in the previous subsection. 








 §.§ Considering the permutation symmetry

Let us consider the general case where the conditional probability marginal in the BP messages is a convex combination of factorized forms

    q_i(h_i:σ⃗_i)=∑_αc_i(α)∏_j=1^i-1( 𝕀(D_i,j≥ d)g_i,j^α(σ⃗_j:σ⃗_i) ).

Then we try to minimize the consistency error of the above trail solutions in the BP equations

    ℒ[{c_i^α,g_i,j=1,…,i-1^α}]
    
    =∑_{σ⃗_j=1,…,i-1:D_i,j≥ d}( ∑_αc_i(α)∏_j=1^i-1g_i,j^α-1/z_i∑_αc_i-1(α)∏_j=1^i-2𝕀(D_i-1,j≥ d)g_i-1,j^α)^2,

with respect to the variables c_i(α), and g_i,j^α.

The stationary equations for the unknown variables lead to the following expressions

    λ_i=∑_α'c_i(α')∏_j=1^i-1(∑_σ⃗_j:D_i,j≥ dg_i,j^αg_i,j^α')
    
    -1/z_i∑_α'c_i-1(α')∑_σ⃗_i-1:D_i,i-1≥ dg_i,i-1^α∏_j=1^i-2(∑_σ⃗_j:D_i,j,D_i-1,j≥ dg_i,j^αg_i-1,j^α'),

and

    λ_ij^α=c_i(α)∑_α'c_i(α')g_i,j^α'∏_k=1:k≠ j^i-1(∑_σ⃗_k:D_i,k≥ dg_i,k^αg_i,k^α')
    
    -1/z_ic_i(α)∑_α'c_i-1(α')∑_σ⃗_i-1:D_i,i-1≥ dg_i,i-1^αg_i-1,j^α'∏_k=1:k≠ j^i-2(∑_σ⃗_k:D_i,k,D_i-1,k≥ dg_i,k^αg_i-1,k^α'),

if j=1,…,i-2, and for j=i-1

    λ_i,i-1^α=c_i(α)∑_α'c_i(α')g_i,i-1^α'∏_k=1^i-2(∑_σ⃗_k:D_i,k≥ dg_i,k^αg_i,k^α')
    
    -1/z_ic_i(α)∑_α'c_i-1(α')∏_k=1^i-2(∑_σ⃗_k:D_i,k,D_i-1,k≥ dg_i,k^αg_i-1,k^α').

The Lagrange multipliers λ_i and λ_ij^α ensure that 

    ∑_αc_i(α)    =1,
    ∑_σ⃗_j:D_i,j≥ dg_i,j^α(σ⃗_j:σ⃗_i)    =1.





To be specific, we take 

    ∑_αc_i(α)[…]=∑_l_j=1,…,i-1=d^nc_i(l_1,…,l_i-1)1/(i-1)!∑_P[…].

This means that the i-1 spheres are distributed in distances l_1,…,l_i-1, and to ensure the permutation symmetry we sum over all the (i-1)! permutations with equal weights. More precisely, now the the conditional part of the BP messages read as follows

    q_i(h_i:σ⃗_i)=∑_l_j=1,…,i-1≥ dc_i(l_1,…,l_i-1)1/(i-1)!∑_P∏_j=1^i-1g_i,j^l_Pj(σ⃗_j:σ⃗_i),

where as before h_i={σ⃗_1:i-1}, and j'=Pj shows the effect of permutation P on j.

By symmetry, the coefficients c_i(l_1,…,l_i-1)=c_i({l_1:i-1}) depend on the number of spheres at distance l.
As before, we assume that the g_i,j^l_Pj depend on Hamming distances D_i,j and consider a complete and normalized representation

    g_i,j^l_Pj(σ⃗_j:σ⃗_i)=1/w(l_Pj)δ_l_ij,l_Pj.

This fixes the values of the g_i,j variables in the cavity marginals.  
We also defined w(l)=C(l:n) for the number of points at distance l from an arbitrary point in the Hamming space.
In this way, the stationary equations are given by

    λ_i=1/(i-1)!(i-1)!∑_{l_1:i-1'}∑_P,P'c_i({l_1:i-1'})∏_j=1^i-1(δ_l_Pj,l_P'j'/w(l_Pj))
    
    -1/z_i1/(i-1)!(i-2)!∑_{l_1:i-2'}∑_P,P'c_i-1({l_1:i-2'})∑_σ⃗_i-1:D_i,i-1≥ dg_i,i-1^l_P(i-1)∏_j=1^i-2(δ_l_Pj,l_P'j'/w(l_Pj)),

which can be simplified to

    λ_i=1/∏_j=1^i-1w(l_j)c_i({l_1:i-1})
    -1/z_i1/(i-1)∑_k=1^i-11/∏_j kw(l_j)c_i-1({l_1:i-1}∖ l_k).

Or, in terms of the c_i variables we have   

    c_i({l_1:i-1})=λ_i∏_j=1^i-1w(l_j)
    +1/z_i1/(i-1)∑_k=1^i-1w(l_k)c_i-1({l_1:i-1}∖ l_k).

On the other hand, from the normalization constraints 

    ∑_l_1,…,l_i-1=d^n c_i({l_1:i-1})    =1,
    1/2^n∑_σ⃗_i,{σ⃗_1:i-1}q_i({σ⃗_1:i-1}:σ⃗_i)    =1,

one obtains

    λ_i=1/(2^n-V_d)^i-1(1-2^n-V_d/z_i),

and

    z_i=∑_l_1,…,l_i-1=d^n c_i-1({l_1:i-2})w(l_i-1)∏_j=1^i-2(∑_l'=d^nΩ(l',l_j:l_i-1)/w(l_j) ).

Finally, after some straightforward algebra, we obtain a recursive equation for the c_i in terms of the previously defined quantities  

    c_i({l_1:i-1})=∏_j=1^i-1(w(l_j)/2^n-V_d)
    
    +∑_k=1^i-1w(l_k)c_i-1({l_1:i-1}∖ l_k)/(i-1)-(2^n-V_d)∏_j=1^i-1(w(l_j)/(2^n-V_d))/∑_l_1',…,l_i-1'=d^n c_i-1({l_1:i-2'})w(l_i-1')∏_j=1^i-2(∑_l'=d^nΩ(l',l_j':l_i-1')/w(l_j') ).

A reasonable initial condition for the above equation is

    c_2(l_1)=w(l_1)/2^n-V_d.

This makes the second term for computing c_3(l_1,l_2) zero and results in

    c_3(l_1,l_2)=w(l_1)/2^n-V_dw(l_2)/2^n-V_d.

The same happens for other coefficients and we obtain

    c_i({l_1:i-1})=∏_j=1^i-1(w(l_j)/2^n-V_d),

which is equivalent to a uniform distribution of the sphere positions {σ⃗_1:i-1}. 


From the above solutions we obtain the node contributions to the entropy Eq. (<ref>),

    e^Δ S_i=1/2^n(i+2)∑_l_1,…,l_i-2=d^n c_i-1({l_1:i-2})1/(i-2)!∑_P
    ∑_σ⃗_i-1,σ⃗_i,σ⃗_i+1:D_i,i-1,D_i,i+1,D_i-1,i+1≥ d∏_j=1^i-2(∑_σ⃗_j:D_j,i-1,D_j,i,D_j,i+1≥ d g_i-1,j^l_Pj(σ⃗_j:σ⃗_i-1)),

which simplifies to

    e^Δ S_i=1/2^n(i+1)∑_l_1,l_2,l_3=d^n C(l_1:n)Ω(l_2,l_3:l_1)(1/2^n-V_d∑_l_i-1,l_i,l_i+1=d^n Y(l_i-1,l_i,l_i+1:l_1,l_2,l_3))^i-2.

The edge contributions are obtained from Eq. (<ref>), 

    e^Δ S_i,i+1= 1/2^n(i+2)∑_l_1,…,l_i-1=d^n c_i({l_1:i-1})1/(i-1)!∑_P
    ∑_σ⃗_i,σ⃗_i+1:D_i,i+1≥ d∏_j=1^i-1(∑_σ⃗_j:D_j,i,D_j,i+1≥ d g_i,j^l_Pj(σ⃗_j:σ⃗_i)),

or, after simplification

    e^Δ S_i,i+1=1/2^n(i+1)∑_l=d^n C(l:n)(1/2^n-V_d∑_l_i,l_i+1=d^n Ω(l_i,l_i+1:l))^i-1.

As expected, these are identical expressions to those obtained for the entropy using the naive approximation of the BP messages.





 §.§ Inhomogeneous solutions

Let us break the translation symmetry and look for more general solutions to the BP equations μ_i(σ⃗_i,h_i)= p_i(σ⃗_i) q_i(h_i:σ⃗_i) where p_i(σ⃗_i) is not necessarily 1/2^n. 
Assuming again that the conditional part is given by

    q_i(h_i:σ⃗_i)=∑_αc_i(α)∏_j=1^i-1( 𝕀(D_i,j≥ d)g_i,j^α(σ⃗_j:σ⃗_i) ),

we minimize the resulted error in the BP equations

    ℒ[p_i(σ⃗_i),{c_i^α,g_i,j=1,…,i-1^α}]=∑_σ⃗_i,{σ⃗_j=1,…,i-1:D_i,j≥ d}
    ( p_i(σ⃗_i)∑_αc_i(α)∏_j=1^i-1g_i,j^α-1/z_ip_i-1(σ⃗_i-1)∑_αc_i-1(α)∏_j=1^i-2𝕀(D_i-1,j≥ d)g_i-1,j^α)^2.

The main equations and details of computations for this general case are given in Appendix <ref>. 


In the following, we continue with the last expression for the q_i(h_i:σ⃗_i) in the previous subsection which respects the permutation symmetry, that is

    q_i(h_i:σ⃗_i)=∑_l_j=1,…,i-1≥ dc_i(l_1,…,l_i-1)1/(i-1)!∑_P∏_j=1^i-1g_i,j^l_Pj(σ⃗_j:σ⃗_i).

As before g_i,j^l_Pj(σ⃗_j:σ⃗_i)=1/w(l_Pj)δ_l_ij,l_Pj is fixed and we have only two stationary equations for the p_i(σ⃗_i) and c_i({l_1:i-1}). Derivation of the error function Eq. (<ref>) with respect to the p_i(σ⃗_i) in presence of the normalization constraints leads to

    p_i(σ⃗_i)=γ_i'
    +1/∑_{l_1:i-1≥ d}c_i^2({l_1:i-1})/∏_j=1^i-1w(l_j)∑_{l_1:i-1≥ d}c_i({l_1:i-1})/∏_j=1^i-1w(l_j)
    ×1/(i-1)!∑_P1/z_i∑_{l_1:i-2'≥ d}c_i-1({l_1:i-2'}) ∑_σ⃗_i-1:D_i,i-1=l_P(i-1)p_i-1(σ⃗_i-1)
    ∏_j=1^i-2( Ω(l_Pj,l_j':l_P(i-1))/w(l_j')).

The stationary equations for the c_i read as follows 

    c_i({l_1:i-1})=λ_i'∏_j=1^i-1w(l_j)
    +1/∑_σ⃗_ip_i^2(σ⃗_i)∑_σ⃗_ip_i(σ⃗_i)
    ×1/(i-1)!∑_P1/z_i∑_{l_1:i-2'≥ d}c_i-1({l_1:i-2'})∑_σ⃗_i-1:D_i,i-1=l_P(i-1)p_i-1(σ⃗_i-1)∏_j=1^i-2( Ω(l_Pj,l_j':l_P(i-1))/w(l_j')).

The normalization factor z_i is 

    z_i=∑_{l_1:i-2'≥ d}^n c_i-1({l_1:i-2'})∑_σ⃗_i,σ⃗_i-1:l_i,i-1≥ dp_i-1(σ⃗_i-1)∏_j=1^i-2(∑_l=d^nΩ(l,l_j':l_i,i-1)/w(l_j')),

and the coefficients γ_i' and λ_i' are Lagrange multipliers.

Note that the first terms of the two equations are expected from a homogeneous solution and in both equations we have a term like

    Q(σ⃗_i,{l_1:i-1})≡1/(i-1)!∑_P1/z_i∑_{l_1:i-2'≥ d}c_i-1({l_1:i-2'})
    
    ×∑_σ⃗_i-1:D_i,i-1=l_P(i-1)p_i-1(σ⃗_i-1)∏_j=1^i-2( Ω(l_Pj,l_j':l_P(i-1))/w(l_j')).

Normalization conditions give the Lagrange multipliers

    γ_i'    =1/2^n(1-
    1/∑_{l_1:i-1≥ d}c_i^2({l_1:i-1})/∏_j=1^i-1w(l_j)∑_{l_1:i-1≥ d}c_i({l_1:i-1})/∏_j=1^i-1w(l_j)∑_σ⃗_i Q(σ⃗_i,{l_1:i-1}) ),
    λ_i'    =1/(2^n-V_d)^i-1(1-1/∑_σ⃗_ip_i^2(σ⃗_i)∑_σ⃗_ip_i(σ⃗_i)∑_{l_1:i-1≥ d}Q(σ⃗_i,{l_1:i-1}) ).

We rewrite the above equations for the p_i and c_i in a more compact form   

    p_i(σ⃗_i)    =1/2^n+Δ p_i(σ⃗_i),
    
    c_i({l_1:i-1})    =∏_j=1^i-1(w(l_j)/2^n-V_d)+Δ c_i({l_1:i-1}),

with deviations Δ p_i,Δ c_i from the liquid solution defined by

    Δ p_i(σ⃗_i)≡1/∑_{l_1:i-1≥ d}c_i^2({l_1:i-1})/∏_j=1^i-1w(l_j)
    ×∑_{l_1:i-1≥ d}c_i({l_1:i-1})/∏_j=1^i-1w(l_j)(Q(σ⃗_i,{l_1:i-1})-1/2^n∑_σ⃗_i'Q(σ⃗_i',{l_1:i-1})),

and

    Δ c_i({l_1:i-1})≡1/∑_σ⃗_ip_i^2(σ⃗_i)∑_σ⃗_ip_i(σ⃗_i)
    
    ×(Q(σ⃗_i,{l_1:i-1})-∏_j=1^i-1(w(l_j)/2^n-V_d)∑_{l_1:i-1'≥ d}Q(σ⃗_i,{l_1:i-1'})).




Let us rewrite the equations for the p_i and c_i as 

    [p_i(σ⃗_i)-γ_i']( ∑_{l_1:i-1≥ d}c_i^2({l_1:i-1})/∏_j=1^i-1w(l_j))   =∑_{l_1:i-1≥ d}c_i({l_1:i-1})/∏_j=1^i-1w(l_j)Q(σ⃗_i,{l_j}),
    
    [c_i({l_1:i-1})-λ_i'∏_j=1^i-1w(l_j)](∑_σ⃗_ip_i^2(σ⃗_i))   =∑_σ⃗_ip_i(σ⃗_i)Q(σ⃗_i,{l_1:i-1}).

Now, we multiply the first equation by p_i(σ⃗_i) and sum over σ⃗_i. Using the normalization condition and the second equation we get

    [∑_σ⃗_ip_i^2(σ⃗_i)-γ_i']( ∑_{l_1:i-1≥ d}c_i^2({l_1:i-1})/∏_j=1^i-1w(l_j))
    
    =(∑_σ⃗_ip_i^2(σ⃗_i))∑_{l_1:i-1≥ d}c_i({l_1:i-1})/∏_j=1^i-1w(l_j)[c_i({l_1:i-1})-λ_i'∏_j=1^i-1w(l_j)].

Simplifying the equation results in

    γ_i'( ∑_{l_1:i-1≥ d}c_i^2({l_1:i-1})/∏_j=1^i-1w(l_j))
    =λ_i'(∑_σ⃗_ip_i^2(σ⃗_i)),

where we also used normalization condition ∑_{l_1:i-1≥ d}c_i({l_1:i-1})=1. 


For the liquid case

    p_i(σ⃗_i)    =1/2^n,
     
    c_i({l_1:i-1})    =∏_j=1^i-1(w(l_j)/(2^n-V_d)),

a consistent solution to the equations is given by

    Q(σ⃗_i,{l_1:i-1}=(i-1)/2^n∏_j=1^i-1(w(l_j)/2^n-V_d),

with Δ p_i=Δ c_i=0. In addition, for Eqs. <ref>, <ref>, and <ref> we obtain

    γ_i'    =1/2^n(1-∑_σ⃗_i,{l_1:i-1≥ d}Q(σ⃗_i,{l_1:i-1})),
    λ_i'    =1/(2^n-V_d)^i-1(1-∑_σ⃗_i,{l_1:i-1≥ d}Q(σ⃗_i,{l_1:i-1})),
    λ_i'/2^n   =γ_i'/(2^n-V_d)^i-1,

which are satisfied by the liquid solution. 



On the other side, we may consider a frozen solution, where the probability distributions are concentrated on a single configuration, for instance, 

    p_i(σ⃗_i)    =δ_σ⃗_i,σ⃗_i^*,
     
    c_i({l_1:i-1})    =δ_{l_1:i-1},{l_1:i-1^*}.

In this case, from Eqs. <ref>, <ref>, and <ref> we find

    γ_i'    =1/2^n(1-∑_σ⃗_iQ(σ⃗_i,{l_1:i-1^*})),
    λ_i'    =1/(2^n-V_d)^i-1(1-∑_{l_1:i-1≥ d}Q(σ⃗_i^*,{l_1:i-1})),
    λ_i'    =γ_i'/∏_j=1^i-1w(l_j^*).

This results in a consistency equation which should be satisfied by the (σ⃗_i^*,{l_1:i-1^*}) given the previous assignments for j=1,⋯,i-1,

    1/2^n(1-∑_σ⃗_iQ(σ⃗_i,{l_1:i-1^*}))=
    ∏_j=1^i-1w(l_j^*)/(2^n-V_d)^i-1(1-∑_{l_1:i-1≥ d}Q(σ⃗_i^*,{l_1:i-1})).

Starting from an initial condition one can try to satisfy the above equation by iteration to find a configuration (σ⃗_i^*,{l_1:i-1^*}) for sphere i. An initial condition starting from i=2 could be 

    p_2(σ⃗_2)=δ_σ⃗_2,0⃗,
    
    c_2(l_1)=δ_l_1,d.

Let us assume that the probabilities c_i({l_j}) are concentrated on a single configuration {l_j^*}. Then we compute the p_i(σ⃗_i) from Eq. <ref> using the above initial condition. The deviation Δ P=∑_σ⃗|p_i(σ⃗)-1/2^n| is then maximized over all possible configurations of {l_j^*}. 
Figure <ref> shows that the above process converges quickly to a uniform solution p_i(σ⃗_i)=1/2^n as the number of spheres N increases. The minimum of deviation is zero for all the points which are reported in the figure.   








§ CONCLUSION

We used an exact representation of the entropy of a system of hard spheres in the Hamming space to investigate the asymptotic behavior of the maximum packing density of the spheres. The method is based on a decomposition of the interaction graph into a tree structure and the induced loopy interactions, which their effects are taken into account by passing some (internal) messages along the tree interaction graph. The solutions to the Bethe equations for reasonable approximations of the BP cavity marginals are asymptotically consistent with the Gilbert-Varshamov lower bound for the packing density, but we can not exclude the possibility of solutions of higher densities.       

We considered an ansatz of BP marginals which are represented by a linear superposition of factorized probability distributions which respect the permutation symmetry of the problem. It would be interesting to try other classes of tractable cavity marginals as trial BP messages which are constrained by satisfying the exact BP equations. For instance, another possibility is working with simpler internal messages h_i→ j instead of approximating the BP cavity marginals. Given an ansatz of reasonable BP probability marginals, one can also try to maximize the Bethe entropy to find an upper bound within a specific class of solutions to the BP equations.                  



We would like to thank Parsa Rangriz and Francesco Zamponi for helpful discussions and useful comments. This work was performed using the Academic Leiden Interdisciplinary Cluster Environment (ALICE), the compute resources provided by Leiden University.





§ THE NAIVE APPROXIMATION OF THE BP MESSAGES: ASYMPTOTIC ENTROPY

The BP probability marginals μ_i(σ⃗_i,h_i) in general depend on the position probability of sphere i, p_i(σ⃗_i) and a conditional probability q_i(h_i:σ⃗_i), 

    μ_i(σ⃗_i,h_i)= p_i(σ⃗_i) q_i(h_i:σ⃗_i).

The h_i={σ⃗_j=1,…,i-1} represent the positions of all spheres j=1,⋯,i-1.
By symmetry we take a uniform distribution for σ⃗_i,

    p_i(σ⃗_i)=1/2^n.

We also assume that the probability distribution of the other spheres is factorized and uniform

    q_i(h_i:σ⃗_i)=∏_j=1^i-1𝕀(D_i,j≥ d)/2^n-V_d.


In this way, for the nodes' contribution to the entropy we have

    e^Δ S_i=1/2^n+ni∑_l_1,l_2,l_12=d^n C(l_12:n)
    Ω(l_1,l_2:l_12)
    ×(1-2V_d-O_01(l_1)-O_02(l_2)-O_12(l_12)+2O_012(l_1,l_2,l_12)/2^n-V_d)^i-2.

The edges' contribution in the entropy are

    e^Δ S_i,i+1= 1/2^n+ni∑_l=d^n C(l:n)
    (1-V_d-O_01(l)/2^n-V_d)^i-1.

Here O_ij(l) is the overlap of two spheres of radius d at distance l. And O_ijk(l_i,l_j,l_ij) is the overlap of three spheres of radius d when (i,j) have distance l_ij and the other sphere is at distances l_i,l_j form (i,j).
The function Ω(l_i,l_j:l_ij) is the number of possible points for the third sphere given the distances,

    Ω(l_i,l_j:l_ij)= C(l_i-l_j+l_ij/2:l_ij)C(l_i+l_j-l_ij/2:n-l_ij).

Thus for the overlaps we get

    O_ij(l_ij)=∑_l_i,l_j=0^dΩ(l_i,l_j:l_ij).

Let us define Y(l_i,l_j,l_k:l_ij,l_ik,l_jk) as the number of points at distances l_i,l_j,l_k from spheres i,j,k with distances l_ij,l_ik,l_jk,

    Y(l_i,l_j,l_k:l_ij,l_ik,l_jk)
    
    =∑_x_00=0^n C(x_11:x)C(x_01:l_ik-x)C(x_00:n-l_ik-l_jk+x)C(x_10:l_jk-x)

with

    2x    =l_ik+l_jk-l_ij, 
    
    2x_01   =l_j+l_k-l_jk-2x_00, 
    
    2x_10   =l_i+l_k-l_ik-2x_00, 
    
    2x_11   =l_ik+l_jk-l_i-l_j+2x_00.

Then

    O_ijk(l_i,l_j,l_ij)=∑_l_i,l_j,l_k=0^dY(l_i,l_j,l_k:l_ij,l_ik,l_jk).



Let us take the limit n,d→∞ with δ=d/n finite. All distances are scaled with n, for instance r_ij=l_ij/n. Then v_d=V_d/2^n=e^-n[ln(2)-H(δ)] with H(δ)=-δln(δ)-(1-δ) ln(1-δ) is an exponentially small quantity. The node's contribution in the entropy is given by  

    Δ S_node=∑_i=1^NΔ S_i=≃ nN∫_0^1 dx[H^*(x)+H_1^*(x)+H_2^*(x)
    
    -xṽ_d(2-o_01^*(x)-o_02^*(x)-o_12^*(x)+2o_012^*(x))]
    -nNln(2)/2,

after approximating ∑_i ≈ N∫_0^1 dx.
The overlaps are scaled o=O/V_d and ṽ_d=Nv_d/n. The star means the above quantities are computed at r^*(x),r_1^*(x),r_2^*(x) which

    r^*(x),r_1^*(x),r_2^*(x)=
    max_r,r_1,r_2 ∈ (δ,1)[ H(r)+H_1(r_1,r_2,r)+H_2(r_1,r_2,r)
    
    -xṽ_d(2-o_01(r_1)-o_02(r_2)-o_12(r)+2o_012(r_1,r_2,r))].

Here

    H_1(r_1,r_2,r)=rln r-r_1-r_2+r/2lnr_1-r_2+r/2- r_2-r_1+r/2lnr_2-r_1+r/2,

and

    H_2(r_1,r_2,r)=(1-r)ln(1-r)
    
    -r_1+r_2-r/2lnr_1+r_2-r/2- (1-r_1+r_2+r/2)ln(1-r_1+r_2+r/2).





The overlaps are exponentially small here. Consider for instance the maximum overlap o_ij(r) at r=δ,

    o_ij(δ)=O_ij/V_d≃ e^n[H_1(r',r',δ)+H_2(r',r',δ)-H(δ)].

By symmetry, the maximum is for r_1=r_2=r'. A plot of the exponent H_1(r',r',δ)+H_2(r',r',δ)-H(δ) as a function of r' in Fig. <ref> shows that this quantity is always negative for δ<1/2.

Since the overlaps are exponentially small, they can be ignored to get

    Δ S_node≃ nN(H(1/2)+H_1(1/2,1/2,1/2)+H_2(1/2,1/2,1/2)-1/2(2ṽ_d+ln(2)) ),

 

The link's contribution in the entropy is given by  

    Δ S_link=∑_i=1^N-1Δ S_i,i+1≃ nN(∫_0^1 dx[H^*(x)-xṽ_d(1-o_ij^*(x))]-ln(2)/2).

The star means the above quantities are computed at r^*(x) which

    r^*(x)=max_r∈ (δ,1)( H(r)-xṽ_d(1-o_ij(r))).

The exponentially small o_ij can be ignored to get

    Δ S_link≃ nN(H(1/2)-1/2(ṽ_d+ln(2)) ),



Putting all together the entropy reads

    Δ S_node≃ nN(3/2ln 2-ṽ_d),
    Δ S_link≃ nN1/2(ln 2-ṽ_d),
    
    S^naive=Δ S_node-Δ S_link≃ nN(ln 2-1/2ṽ_d).

Therefore, the entropy vanishes at

    N_max^naive=(2ln 2)n/v_d,

at the same point that the Bethe entropy in Eq. <ref> goes to zero.


 


§ THE BP EQUATIONS FOR INHOMOGENEOUS SOLUTIONS

Let us break the translation symmetry and look for more general solutions to the BP equations μ_i(σ⃗_i,h_i)= p_i(σ⃗_i) q_i(h_i:σ⃗_i) where p_i(σ⃗_i) is not necessarily 1/2^n. 
Assuming that

    q_i(h_i:σ⃗_i)=∑_αc_i(α)∏_j=1^i-1( 𝕀(D_i,j≥ d)g_i,j^α(σ⃗_j:σ⃗_i) ),

we minimize the expected error in the BP equations

    ℒ[p_i(σ⃗_i),{c_i^α,g_i,j=1,…,i-1^α}]=∑_σ⃗_i,{σ⃗_j=1,…,i-1:D_i,j≥ d}
    ( p_i(σ⃗_i)∑_αc_i(α)∏_j=1^i-1g_i,j^α-1/z_ip_i-1(σ⃗_i-1)∑_αc_i-1(α)∏_j=1^i-2𝕀(D_i-1,j≥ d)g_i-1,j^α)^2,

The above function should be minimized with respect to the variables p_i(σ⃗_i),{c_i^α,g_i,j=1,…,i-1^α} conditioned to the following normalization constraints:

    ∑_σ⃗_ip_i(σ⃗_i)    =1,
    ∑_αc_i(α)    =1,
    ∑_σ⃗_j:D_i,j≥ dg_i,j^α(σ⃗_j:σ⃗_i)    =1.



The stationary equations with respect to the variables p_i(σ⃗_i),{c_i^α,g_i,j=1,…,i-1^α} read as follows 

    γ_i=p_i(σ⃗_i)∑_α,α'c_i(α)c_i(α')∏_j=1^i-1(∑_σ⃗_j:D_i,j≥ dg_i,j^αg_i,j^α')
    
    -1/z_i∑_α,α'c_i(α)c_i-1(α')∑_σ⃗_i-1:D_i,i-1≥ dp_i-1(σ⃗_i-1)g_i,i-1^α∏_j=1^i-2(∑_σ⃗_j:D_i,j,D_i-1,j≥ dg_i,j^αg_i-1,j^α'),



    λ_i=∑_σ⃗_ip_i^2(σ⃗_i)∑_α'c_i(α')∏_j=1^i-1(∑_σ⃗_j:D_i,j≥ dg_i,j^αg_i,j^α')
    
    -1/z_i∑_σ⃗_ip_i(σ⃗_i)∑_α'c_i-1(α')∑_σ⃗_i-1:D_i,i-1≥ dp_i-1(σ⃗_i-1)g_i,i-1^α∏_j=1^i-2(∑_σ⃗_j:D_i,j,D_i-1,j≥ dg_i,j^αg_i-1,j^α'),

and

    λ_ij^α=p_i^2(σ⃗_i)c_i(α)∑_α'c_i(α')g_i,j^α'∏_k=1:k≠ j^i-1(∑_σ⃗_k:D_i,k≥ dg_i,k^αg_i,k^α')
    
    -1/z_ip_i(σ⃗_i)c_i(α)∑_α'c_i-1(α')∑_σ⃗_i-1:D_i,i-1≥ dp_i-1(σ⃗_i-1)g_i,i-1^αg_i-1,j^α'∏_k=1:k≠ j^i-2(∑_σ⃗_k:D_i,k,D_i-1,k≥ dg_i,k^αg_i-1,k^α'),

if j=1,…,i-2, and for j=i-1

    λ_i,i-1^α=p_i^2(σ⃗_i)c_i(α)∑_α'c_i(α')g_i,i-1^α'∏_k=1^i-2(∑_σ⃗_k:D_i,k≥ dg_i,k^αg_i,k^α')
    
    -1/z_ip_i(σ⃗_i)p_i-1(σ⃗_i-1)c_i(α)∑_α'c_i-1(α')∏_k=1^i-2(∑_σ⃗_k:D_i,k,D_i-1,k≥ dg_i,k^αg_i-1,k^α').

The Lagrange multipliers γ_i, λ_i and λ_ij^α are to ensure the normalization constraints. 





prsty

CS-book-1999 Conway JH, Sloane NJA. Sphere packings, lattices and groups, third edition, Grundlehren der Mathematischen Wissenschaften 290, Springer-Verlag, New York, 1999.

Z-book-2008 Zong C. Sphere packings. Springer Science & Business Media; 2008 Jan 20.

PZ-revmp-2010 Parisi G, Zamponi F. Mean-field theory of hard sphere glasses and jamming. Reviews of Modern Physics. 2010 Mar 16;82(1):789.

TS-revmp-2010 Torquato S, Stillinger FH. Jammed hard-particle packings: From Kepler to Bernal and beyond. Reviews of modern physics. 2010 Sep 15;82(3):2633.




F-mathz-1940 Fejes L. Über einen geometrischen Satz. Mathematische Zeitschrift. 1940 Dec;46(1):83-5.

H-anmath-2005 Hales TC. A proof of the Kepler conjecture. Annals of mathematics. 2005 Nov 1:1065-185.

V-anmath-2017 Viazovska MS. The sphere packing problem in dimension 8. Annals of mathematics. 2017 May 1:991-1015.

CKMRV-2017 Cohn H, Kumar A, Miller S, Radchenko D, Viazovska M. The sphere packing problem in dimension 24. Annals of Mathematics. 2017 May 1;185(3):1017-33.




MRRW-ieee-1977 McEliece R, Rodemich E, Rumsey H, Welch L. New upper bounds on the rate of a code via the Delsarte-MacWilliams inequalities. IEEE transactions on Information Theory. 1977 Mar;23(2):157-66.

B-mathres-1992 Ball K. A lower bound for the optimal density of lattice packings. International Mathematics Research Notices. 1992 Jan 1;1992(10):217-21.

DL-ieee-1998 Delsarte P, Levenshtein VI. Association schemes and coding theory. IEEE Transactions on Information Theory. 1998 Oct;44(6):2477-504.

S-jcomb-2001 Samorodnitsky A. On the optimum of Delsarte's linear program. Journal of Combinatorial Theory, Series A. 2001 Nov 1;96(2):261-87.

JV-ieee-2004 Jiang T, Vardy A. Asymptotic improvement of the Gilbert-Varshamov bound on the size of binary codes. IEEE Transactions on Information Theory. 2004 Jul 26;50(8):1655-64.

SST-jmat-2008 Scardicchio A, Stillinger FH, Torquato S. Estimates of the optimal density of sphere packings in high dimensions. Journal of Mathematical Physics. 2008 Apr 1;49(4).

C-geo-2002 Cohn H. New upper bounds on sphere packings II. Geometry & Topology. 2002 Jun 25;6(1):329-53.

KLV-mathres-2004 Krivelevich M, Litsyn S, Vardy A. A lower bound on the density of sphere packings via graph theory. International Mathematics Research Notices. 2004;2004(43):2271-9.

ACHT-jhep-2020 Afkhami-Jeddi N, Cohn H, Hartman T, de Laat D, Tajdini A. High-dimensional sphere packing and the modular bootstrap. Journal of High Energy Physics. 2020 Dec;2020(12):1-45.

CT-mathc-2022 Cohn H, Triantafillou N. Dual linear programming bounds for sphere packing via modular forms. Mathematics of Computation. 2022 Jan;91(333):491-508.





PS-jstat-1999 Procacci A, Scoppola B. Statistical mechanics approach to coding theory. Journal of statistical physics. 1999 Aug;96:907-12.

PS-pre-2000 Parisi G, Slanina F. Toy model for the mean-field theory of hard-sphere liquids. Physical Review E. 2000 Nov 1;62(5):6554.

PZ-jstat-2006 Parisi G, Zamponi F. On the high density behavior of Hamming codes with fixed minimum distance. Journal of statistical physics. 2006 Jun;123(6):1145-67.







CKPUZ-anrevcm-2017 Charbonneau P, Kurchan J, Parisi G, Urbani P, Zamponi F. Glass and jamming transitions: From exact results to finite-dimensional descriptions. Annual Review of Condensed Matter Physics. 2017 Mar 31;8:265-88.





RZ-pre-2012 Ramezanpour A, Zecchina R. Cavity approach to sphere packing in Hamming space. Physical Review E. 2012 Feb 6;85(2):021106.

R-pre-2013 Ramezanpour A. Computing loop corrections by message passing. Physical Review E. 2013 Jun 28;87(6):060103.





MP-physc-2001 Mezard M, Parisi G. The Bethe lattice spin glass revisited. The European Physical Journal B-Condensed Matter and Complex Systems. 2001 Mar;20:217-33.

KFL-ieee-2001 Kschischang FR, Frey BJ, Loeliger HA. Factor graphs and the sum-product algorithm. IEEE Transactions on information theory. 2001 Feb;47(2):498-519.

YFW-artint-2003 Yedidia JS, Freeman WT, Weiss Y. Understanding belief propagation and its generalizations. Exploring artificial intelligence in the new millennium. 2003 Jan 8;8(236-239):0018-9448.

MM-book-2009 Mezard M, Montanari A. Information, physics, and computation. Oxford University Press; 2009 Jan 22.

baldassi-polito-2009 Baldassi, C., Braunstein, A., Ramezanpour, A., Zecchina, R. (2015). Statistical Physics and Network Optimization Problems. In: Fagnani, F., Fosson, S., Ravazzi, C. (eds) Mathematical Foundations of Complex Networked Information Systems. Lecture Notes in Mathematics(), vol 2141. Springer, Cham.


KLC-pra-2013 Kraus CV, Lewenstein M, Cirac JI. Ground states of fermionic lattice Hamiltonians with
permutation symmetry. Physical Review A. 2013 Aug 28;88(2):022335.































