Also at ]Advanced Materials and Liquid Crystal Institute, Kent State
University, Kent, OH, USA

# Leaky Cell Model of Hard Spheres

Thomas G. Fai<EMAIL_ADDRESS>Department of Mathematics and Volen Center
for Complex Systems, Brandeis University, Waltham, MA, USA Jamie M. Taylor
Basque Center for Applied Mathematics (BCAM), Bilbao, Bizkaia, Spain  Epifanio
G. Virga  Dipartimento di Matematica, Università di Pavia, Pavia, Italy
Xiaoyu Zheng  Peter Palffy-Muhoray  [  Department of Mathematical Sciences,
Kent State University, Kent, OH, USA

###### Abstract

We study packings of hard spheres on lattices. The partition function, and
therefore the pressure, may be written solely in terms of the accessible free
volume, i.e. the volume of space that a sphere can explore without touching
another sphere. We compute these free volumes using a leaky cell model, in
which the accessible space accounts for the possibility that spheres may
escape from the local cage of lattice neighbors. We describe how elementary
geometry may be used to calculate the free volume exactly for this leaky cell
model in two- and three-dimensional lattice packings and compare the results
to the well-known Carnahan-Starling and Percus-Yevick liquid models. We
provide formulas for the free volumes of various lattices and use the common
tangent construction to identify several phase transitions between them in the
leaky cell regime, indicating the possibility of coexistence in crystalline
materials.

††preprint: AIP/123-QED

##  I Introduction

The thermodynamical properties of hard spheres, i.e. objects whose interior
points cannot overlap with those of any other hard object, have been a subject
of great theoretical interest given the conceptual simplicity and rich
phenomena exhibited by these systems Torquato and Stillinger (2010). Systems
of hard spheres are one of the few systems that are analytically tractable,
allowing for exact and/or approximate calculations of physical properties as
reviewed in (2; 3; 4). Despite the relative simplicity of hard spheres, they
exhibit non-trivial behaviors such as phase transitions and allow for
experimentally testable predictionsSantos (2016).

This subject has experienced a resurgence because of new computational tools
for many-body problems. Computer simulations, first used to generate random
packings of hard disks in two dimensions by the method of random sequential
addition Hoover, Hoover, and Hanson (1979) and to calculate freezing/melting
transitions by dynamic Monte Carlo simulation Hoover and Ree (1967), have been
used to study packings of circles and ellipses in two dimensions Vieillard-
Baron (1972) and spheres, ellipsoids, and other shapes in three dimensions
Chaikin  _et al._ (2006); Ni  _et al._ (2012), in tandem with sophisticated
computational methods to analyze packings based on techniques such as
computational topology Carlsson  _et al._ (2012) and discrete geometry Sastry
_et al._ (1997, 1998); Kapfer  _et al._ (2012); Maiti and Sastry (2014); Chen
and Holmes-Cerfon (2017). Boltzmann generators based on deep learning have
recently been developed to sample equilibrium states of confined systems of
hard disks Noé  _et al._ (2019).

This resurgence of interest and new data has highlighted the need for
predictive theories grounded in statistical mechanical first principles. A key
quantity that appears throughout these works is the free volume, i.e. the
volume of space available to the center of mass of each particle in the
packing. The partition function may be approximated in terms of the free
volume accessible to each sphere, which only requires a purely local
calculation in terms of a cage of neighboring spheres. In the special case of
lattice packings with a prescribed crystal structure, only a single unit cell
needs to be considered, and this observation led to the development of a cell
theory (CT) of hard sphere liquids Eyring and Hirschfelder (1937);
Hirschfelder, Stevenson, and Eyring (1937); Lennard-Jones and Devonshire
(1937, 1938).

This cell theory (alternatively called the Lennard-Jones-Devonshire cell
model) allows for various analytically tractable calculations. One can derive
exact expressions for the free volume in the case that the spheres are
arranged in a latticeBuehler  _et al._ (1951) and use these formulas in the
partition function. In this work, we further explore the cell theory by
computing expressions for the free volume for several crystal structures in
one, two, and three dimensions. Whereas the cell theory has thus far primarily
been used to model solids in the high-density limit, here we introduce a leaky
cell model which represents materials that retain their crystal structure at
low density. We use the resulting expressions to predict phase transitions
between lattices and interpret these predictions in the context of numerical
and experimental observations.

The extension of cell theories to low densities was already considered in
Buehler et al. Buehler  _et al._ (1951), in which they average between _hard
center_ and _soft center_ cell models to capture the possibility that spheres
may wander out of their local cages at sufficiently low densities. However, in
all the scenarios they considered, the center of the wandering sphere is
always restricted to a local Wigner-Seitz (or Voronoi) cell. In (22), spheres
are taken to be “bounded either by the collision spheres of its nearest
neighbors or else by mathematical partitions which bisect the distance between
lattice points, whichever barrier appears first.” It is natural then to
consider what happens when these mathematical partitions are absent.

It has been shown that lattice structures lose mechanical stability at low
density Woodcock (1997); Warshavsky, Ford, and Monson (2018). However, the
“ghost lattice” obtained at low density is nevertheless a useful theoretical
instrument which allows us to predict behaviors such as phase transitions
within a fully consistent cell theory. Here, we further develop this concept
by pushing the cell theory into the low and intermediate-density regime in
which spheres may leak out of their local cages.

In the leaky cell model proposed here, we allow for the possibility that
spheres are able to explore an extended region outside of the Wigner-Seitz
cell. We identify an intermediate regime in which spheres escape the Wigner-
Seitz cell—i.e. they may be brought outside the cell by a continuous
translation without making contact with any other sphere—and we find the
critical packing fractions at which the leaky cell model begins to deviate
from the classical one. We consider face-centered cubic (FCC), body-centered
cubic (BCC), and simple cubic (SC) lattices, and we find phase transitions
between these lattices in their corresponding leaky regimes.

To characterize the behavior of lattice packings in this leaky, moderate
density regime, we perform exact calculations of the free volume. Given its
analytic tractability, the leaky cell model extends lattice models outside
their typical range of application and yields predictions for crystalline
materials at moderate densities. The formulas for free volume that we compute
have interesting features such as discontinuities in the compressibility
factor observed for certain lattices. By allowing for a direct comparison of
the free volume of several lattice packings, this work elucidates the
thermodynamical importance of this quantity and provides a foundation for the
study of packings with local or global crystal structure.

We further consider an apparent limitation of the cell theory. Whereas using
the cell theory to model the high-density (solid) regime and empirical
equations of state to model the low-density (liquid) regime leads to good
agreement with observations from experiment and Monte Carlo simulation, our
results indicate that equations of state based purely on cell theory are
unable to accurately predict the freezing/melting transitions in a hard sphere
gas. We rectify this by introducing a new quantity, called the quasi-communal
entropy, which characterizes the difference between the hard sphere gas and
cell theory in the low-density regime.

This article is structured as follows. First, we describe the formulation of
the leaky cell model and its application to several three-dimensional (3D)
lattices. (The two-dimensional (2D) and one-dimensional (1D) cases are also
described in the appendices.) Next, we explain how to derive exact formulas
for the corresponding free volume measures based on a straightforward
geometric argument. Given these formulas, it is straightforward to compute
various physical quantities such as the free energy and compressibility
factor. We use the common tangent construction to identify packing fractions
at which different lattices may coexist with one another or with liquids
described by the empirical Percus-Yevick and/or Carnahan-Starling equations of
state.

##  II Leaky cell model

We calculate free volumes exactly in regular lattice arrangements such as
cubic lattices and hexagonal close-packed (HCP) lattices. The number of
spheres and total volume are fixed. The locations of all but one sphere are
fixed on the lattice. The free volume is defined to be the volume of the
region accessible to the center of mass of a sphere that wanders continuously
from its lattice site without making contact with any other sphere Hoover,
Hoover, and Hanson (1979).

###  II.1 Compressibility factor from free volume

In a hard sphere gas, the equation of state may be determined purely from the
free volume, as we now derive from first principles following closely the
reasoning of previous works Rice (1944); Kirkwood (1950); Wood (1952); Speedy
and Reiss (1991). Once we have a formula $\mathcal{F}=\mathcal{F}(v)$ for the
free volume $\mathcal{F}$ as a function of the Voronoi cell volume $v$, we
construct the partition function $\mathcal{Z}$ in terms of $\mathcal{F}$.

For $N$ hard spheres with center of mass positions $\mathbf{x}_{i}$ and
momenta $\mathbf{p}_{i}$ for $i=1,\dots,N$, the partition function in the case
of single occupancy of each lattice cell is given byEyring and Hirschfelder
(1937)

$\displaystyle\mathcal{Z}$
$\displaystyle=\frac{1}{h^{3N}}\idotsint\mathrm{d}\mathbf{x}_{1}\dots\mathrm{d}\mathbf{x}_{N}\idotsint\exp{\left(\frac{-\|\mathbf{p}_{i}\|^{2}}{2mk_{B}T}\right)}\mathrm{d}\mathbf{p}_{1}\dots\mathrm{d}\mathbf{p}_{N}$
$\displaystyle=\left(2\pi
mk_{B}T/h^{2}\right)^{3N/2}\idotsint\mathrm{d}\mathbf{x}_{1}\dots\mathrm{d}\mathbf{x}_{N}$
$\displaystyle=\left(2\pi mk_{B}T/h^{2}\right)^{3N/2}\mathcal{F}^{N},$ (1)

where the free volume $\mathcal{F}$ depends on the lattice configuration and
the excluded volumes of neighboring spheres, $h$ is Planck’s constant, $k_{B}$
is Boltzmann’s constant, $T$ is the temperature, and $m$ is the mass of each
sphere. (Note that the usual factor of $1/N!$ does not appear explicitly in
the denominator of (1). This is because there is another factor of $N!$ in the
numerator to account for all permutations of spheres among lattice sites, and
these two factors exactly cancel.) Therefore

$\log\mathcal{Z}=N\log\mathcal{F}-3N\log\Lambda,$ (2)

where we have introduced the de Broglie wavelength ${\Lambda=h/\sqrt{2\pi
mk_{B}T}}$. The Helmholtz free energy $A$ is defined in terms of $\mathcal{Z}$
via:

$A:=-k_{B}T\log\mathcal{Z}.$ (3)

Substituting the expression (2) for $\mathcal{Z}$ into this equation yields

$A=-k_{B}T\left(N\log\mathcal{F}-3N\log\Lambda\right).$ (4)

The pressure $p$ is defined by

$p:=-\frac{\partial A}{\partial V},$ (5)

i.e. the pressure is the negative gradient of the free energy with respect to
the total volume (at constant temperature). The second term in (4) is
independent of volume, so that differentiating yields

$p=k_{B}T\frac{\partial\log\mathcal{Z}}{\partial
V}=\frac{Nk_{B}T}{\mathcal{F}}\frac{\partial\mathcal{F}}{\partial V}.$ (6)

The Voronoi cell volume $v$ and total volume $V$ are related on the lattice by
$V=Nv$, so that the pressure may be rewritten in terms of local quantities as

$p=\frac{k_{B}T}{{\mathcal{F}}}\frac{\partial\mathcal{F}}{\partial v}.$ (7)

Finally, the compressibility factor $Z$ is defined via

$Z:=\frac{pV}{Nk_{B}T},$ (8)

and can therefore be expressed in terms of the free volume by

$Z=\frac{v}{\mathcal{F}}\frac{\partial\mathcal{F}}{\partial v}.$ (9)

The above formula, which also appears in Buehler et al. Buehler  _et al._
(1951), provides the foundation for our calculations. We use it to compute the
compressibility factor $Z$ directly from the free volume $\mathcal{F}$ and the
number density $\rho=v^{-1}$.

##  III Results

We next describe how the free volumes of lattices of hard spheres may be
computed analytically for the leaky cell model. We will denote the packing
fraction by $\eta$ throughout the remainder of the article. Letting $v_{0}$
denote the volume of each sphere in the system, $\eta=v_{0}/v$.

###  III.1 Exact formulas on lattices

The free volume in cell theory is defined in terms of the region accessible to
the center of mass of a sphere holding the position of every other sphere
fixed. It is the volume of the region accessible by continuous translation
without making contact with another sphere–this is a connected component in
the set of points located a distance of at least $2R$ from any other sphere,
where $R$ is the sphere radius.

Given that lattice structures are typically associated with solids, the cell
model has classically been applied in the high-density limit, in which case a
sphere cannot escape from the unit cell without making contact with a
neighboring sphere, which is forbidden in the hard sphere model. Therefore, to
compute the free volume in the high-density limit it is sufficient to consider
the unit cell only. In the leaky cell model explored here, we extend the cell
theory to lower packing fractions in which spheres are able to escape from the
unit cell.

Figure 1:  Development of the leaky cell model. Each polyhedron represents the
cage of neighbors for a chosen sphere in the lattice, with neighboring spheres
making up the vertices of the polyhedron. The left column illustrates the cage
at high packing densities near the close-packing fraction $\eta_{\text{cp}}$,
whereas the right column illustrates the extended cage of spheres between the
leaky packing fraction $\eta_{\text{l}}$ and the percolation limit
$\eta_{\text{p}}$. Cell structures were visualized using VESTA Momma and Izumi
(2011).

####  III.1.1 FCC/HCP lattice

We first describe packings of spheres in the FCC lattice, which involves
planar layers of hexagonally-arranged spheres. Note that, as pertains to the
free volume, the HCP and FCC lattices are equivalent; although we discuss only
the FCC lattice, the resulting formulas apply to the HCP lattice as well. (In
the HCP lattice, the layers are stacked in an ABAB pattern, whereas, in the
FCC lattice they are stacked in an ABCABC pattern. It follows that the FCC
lattice and HCP lattice are equivalent up to rotating one half of the unit
cell, so that the resulting free volume formulas are the same.)

As we discuss in detail later in this paper, for packings denser than the
leaky packing fraction $\eta^{\text{FCC}}_{\text{l}}$, the formulas we derive
are equivalent to those given previously by (22). They differ for packings
less dense than $\eta^{\text{FCC}}_{\text{l}}$ because in the leaky cell model
we must also account for the accessible volume outside the Wigner-Seitz cell.

The void space surrounding a sphere in the lattice can be partitioned into
$N_{t}=8$ tetrahedra and $N_{o}=6$ octahedra. As the packing fraction
increases, spheres becomes caged first by 18 neighboring spheres forming a
polyhedron with triangular octahedral faces and subsequently by 12 neighboring
spheres forming a polyhedron with both triangular and square faces (the
octahedral midplanes). See Figure 1(a)–(b).

Figure 2: Comparison between the free volumes in the leaky cell model and the
classical cell theoryBuehler  _et al._ (1951). The asymptotic approximation
corresponding to the curves labeled “approx CT” is given in Appendix B, the
black and red diamonds indicate the percolation and leaky packing fractions,
respectively, of Table 1, and the points labeled Monte Carlo are obtained by
sampling points uniformly within the cage and rejecting any that are inside of
a neighboring exclusion sphere. (a) FCC/HCP lattice, (b) BCC lattice, (inset)
Zoom-in showing that the free volume of the leaky cell model is bounded above
by the classical cell theory between the percolation and leaky cell
transitions, and (c) SC lattice.

As mentioned above, there is a transition in the free volume as the spheres
lose the ability to escape through the square faces at octahedral midplanes.
This leaky cell transition occurs when the edge length $a$ satisfies
$a=2\sqrt{2}R$, which corresponds to the packing fraction
$\eta^{\text{FCC}}_{\text{l}}\approx 0.26$. Interestingly, this is a smooth
transition; there is no discontinuity in the free volume or its first
derivative.

The percolation transition $\eta^{\text{FCC}}_{\text{p}}$ takes place at a
smaller packing fraction at which the spheres cannot escape through triangular
faces. Assuming the spheres have radius $R$, this occurs when $a=2\sqrt{3}R$,
and at this point spheres become locally caged and unable to move freely
throughout the domain. The packing fraction at which this occurs is
$\eta^{\text{FCC}}_{\text{p}}=\Omega/(3\sqrt{3/2})\approx 0.15$, where
$\Omega=\cos^{-1}(23/27)$ is the solid angle subtended by each vertex of a
regular tetrahedron.

The free volume is defined as a piecewise function between these transitions,
with the formula given in Eq. (31) of Appendix A. As detailed there, the
formula we use involves the three-dimensional double intersection volume
$V_{2I}(a)$ of intersecting spheres, as well as the intersection volume of
three and four spheres. In order to calculate these intersection volumes, we
use general formulas for volumes of higher-order intersections as given in
(30; 31).

In Fig. 2(a) we show that at high density our results agree with those of
(22). The disagreement at low densities is because of different formulations
are used for the free volume: in (22) the free volume is taken to be a subset
of the dodecahedral Voronoi cell associated to a given sphere, whereas in our
case we include all accessible volume in the lattice including the voids that
emerge when $a\geq 2\sqrt{2}R$ so that spheres can escape through square faces
of the cuboctahedron associated with the FCC lattice but remain caged by
triangular faces.

####  III.1.2 BCC lattice

There are notable differences in the leaky cell model for the body-centered
cubic (BCC) lattice. Initially, spheres are caged by 8 spheres at the vertices
of a cube (Figure 1(c)–(d)). At the leaky cell fraction
$\eta^{\text{BCC}}_{\text{l}}\approx 0.37$, spheres become caged by 14
spheres, which include the 6 spheres located across cubic faces in addition to
the 8 spheres at cubic vertices. Below the leaky cell fraction, the void space
may be tesselated by 6 octahedra. The free volume is computed by subtracting
the exclusion spheres of the 14 neighbors from the volumes of these octahedra,
accounting for all intersections (going up to quintuple intersections in this
case). See Appendix A for the explicit free volume formulas.

Unlike the FCC and SC lattices, in which spheres can escape from their local
cages at the leaky packing fraction, spheres in a BCC lattice are unable to
escape from their cubic cage at the leaky cell fraction. This is because the
neighboring spheres in the lattice positioned on the opposite sides of the
faces of the cubic cage also exclude volume within the cube. Consequently, the
free volume for the leaky cell model is actually less than the free volume
predicted by the cell theory (Fig. 2(b)).

Figure 3:  Properties of the leaky cell model for FCC/HCP, BCC, and SC
lattices. (a) Free volume $\mathcal{F}$ in units of $R^{3}$ as a function of
the packing density $\eta$, (b) free energy density $f$ vs. the packing
density $\eta$, (b) compressibility factor $Z$ vs. the packing density $\eta$.

####  III.1.3 SC lattice

We next consider the simple cubic lattice. In addition to the percolation
threshold at $\eta^{\text{SC}}_{\text{p}}\approx 0.19$, at which diagonal
neighbors on square faces begin to intersect and the center sphere becomes
caged by 26 neighbors, the SC lattice has a transition at
$\eta^{\text{SC}}_{\text{l}}\approx 0.29$ beyond which the sphere can no
longer escape through triangular faces (Figure 1(e)–(f)), so that opposite
corners of cubes are lost and the number of neighbors decreases to 6.

In Appendix C, we use analogous geometrical arguments on lattices in 2D to
compute free volumes and the resulting compressibility factors.

###  III.2 Comparison between different lattices

The critical packing fractions of the 3D lattice models are summarized in
Table 1.

Table 1: Percolation ($\eta_{\text{p}}$), leaky ($\eta_{\text{l}}$), and close-packed ($\eta_{\text{cp}}$) packing fractions in the FCC, BCC, and SC lattices. 3D Lattice | $\eta_{\text{p}}$ | $\eta_{\text{l}}$ | $\eta_{\text{cp}}$
---|---|---|---
FCC | $\frac{\pi}{9\sqrt{6}}\approx 0.14$ | $\frac{\pi}{12}\approx 0.26$ | $\frac{\pi}{3\sqrt{2}}\approx 0.74$
BCC | $\frac{9\pi}{128\sqrt{2}}\approx 0.16$ | $\frac{\pi}{6\sqrt{2}}\approx 0.37$ | $\frac{\pi\sqrt{3}}{8}\approx 0.68$
SC | $\frac{\pi}{12\sqrt{2}}\approx 0.19$ | $\frac{2\pi}{9\sqrt{6}}\approx 0.29$ | $\frac{\pi}{6}\approx 0.52$

As mentioned previously, above the leaky packing fraction the leaky cell model
reduces to the classical cell theory. To illustrate the difference between the
models below the leaky packing fraction, in Fig. 2 we compare the free volumes
(in units of $v_{0}$) obtained by either the leaky cell model or by classical
cell theory. Interestingly, the differences between the leaky and classical
models depend strongly on the lattice. In particular, there are qualitiative
differences between the FCC and BCC lattices; in the FCC lattice, the leaky
cell model has a free volume that is bounded below by the classical cell
theory, whereas in the BCC lattice it is bounded above by the classical cell
theory for $\eta>\eta^{\text{BCC}}_{\text{p}}$. This is a direct consequence
of the lattice geometry. In the FCC lattice, the wandering sphere leaks out of
the unit cell, whereas in the BCC lattice the exclusion spheres of neighbors
outside the unit cell leaks in.

In Fig. 3, we plot the free volumes (in units of $R^{3}$), free energy
densities (in $k_{B}T$ units), and compressibility factors obtained for the
three lattice models described above. Fig. 3(c) includes a comparison of the
compressibility factor $Z$ from these lattice free volumes—computed according
to Eq. (9)—to standard liquid models such as Percus-Yevick and Carnahan-
Starling for which

$\displaystyle Z^{c}_{PY}(\eta)$
$\displaystyle=\frac{1+\eta+\eta^{2}}{(1-\eta)^{3}}$ (10) $\displaystyle
Z_{CS}(\eta)$ $\displaystyle=\frac{1+\eta+\eta^{2}-\eta^{3}}{(1-\eta)^{3}},$
(11)

where $Z^{c}_{PY}$ and $Z_{CS}$ denote the Percus-Yevick and Carnahan-Starling
compressibility factors, respectively. Note that the Percus-Yevick
approximation does not lead to a unique compressibility factor; the
superscript in $Z^{c}_{PY}$ denotes that this compressibility factor is
obtained through the compressibility route (as opposed to the energy, virial,
or chemical-potential routes Santos (2016)). For convenience, in what follows
we will refer to this expression obtained through the compressibility route
simply as the Percus-Yevick approximation. Note the jump discontinuity in $Z$
for the simple cubic lattice that arises because $\mathcal{F}(v)$ is not
continuously differentiable at $\eta^{\text{SC}}_{\text{l}}$.

Figure 4: Phase transitions in the leaky cell model. Top row: common tangent
construction, bottom row: Maxwell construction via the chemical potential
$\mu$, in which the shaded areas represent equal areas. (a) We recover the
known PY-FCC freezing/melting transition. There is an analogous PY-BCC
freezing transition at $\eta\approx 0.55$ (not shown). (b) SC-FCC phase
transition at $\eta\approx 0.22$, (c) SC-BCC phase transition at $\eta\approx
0.22$, (d) SC-SC phase transition at $\eta\approx 0.3$. A polynomial
interpolant is used to smooth out the discontinuity in the free energy density
that occurs at $\eta_{\text{l}}^{\text{SC}}$.

###  III.3 Coexistence

We apply the common tangent construction to the free energy densities obtained
by the leaky cell model in order to determine the existence of phase
transitions. We first recover the theoretical freezing/melting transition
between the liquid phase and FCC lattice at packing fractions of approximately
0.47 & 0.53 respectively (with 0.47 corresponding to the liquid packing
fraction and 0.53 corresponding to the solid packing fraction) Hoover and Ree
(1968). Note that the empirical values of the freezing and melting packing
fractions are actually closer to 0.49 and 0.54 Hoover and Ree (1968). As shown
in Fig. 4(a) and Table 2,

Table 2: Phase transitions in the leaky cell model. Phase transition | Corresponding packing fractions
---|---
PY–FCC | 0.47–0.53
PY–BCC | 0.54–0.57
SC–FCC | 0.19–0.25
SC–BCC | 0.20–0.25
SC–SC (quasi)111We regard this as a quasi-phase transition, since only a near common tangent is observed. | 0.28–0.32
PY–SC | none
BCC–FCC | none

we are able to recover this result. Since this phase transition occurs above
the leaky packing fraction of $\eta^{\text{FCC}}_{\text{l}}\approx 0.26$, the
leaky cell model simply recovers the freezing/melting transition predicted by
the classical cell theory. In addition to the common tangent construction, we
also verify our result using the equivalent equal area Maxwell construction
Binder  _et al._ (2012). See Fig. 4. For the BCC lattice, we find a PY-BCC
freezing/melting transition at packing fractions of 0.54 & 0.57 analogous to
the known PY-FCC transition. As discussed next, using the appropriate
calibration to high-density is critical for obtaining this phase transition.

We find SC–FCC and SC–BCC inter-lattice phase transitions that occur at
packing fractions of approximately 0.19 & 0.24, respectively, which is between
the leaky packing fraction and the percolation limit (Figs. 4(b) and 4(c)).
This indicates that different lattice structures may coexist in the leaky
regime. Further, we find an SC–SC quasi intra-lattice phase transition around
packing fractions of 0.29 & 0.33, which is again within the leaky regime. Upon
close inspection of the free energy density for the SC lattice, one finds that
there is a near common tangent between these packing fractions, but that there
is only near-tangency, i.e. the tangent misses the curve by approximately
0.5%. We have verified this phase transition by constructing the convex hull
of the free energy, an approach that generalizes to non-differentiable
functions such as the present one.

Moreover, we can rule out the existence of phase transitions between other
lattice structures and between Percus-Yevick or Carnahan-Starling liquids and
the SC lattice. This is because the associated free energy densities do not
intersect, i.e. they are bounded above or below by one another so there can be
no common tangent (Fig. 3(b)).

Figure 5:  (a) There is an ambiguity in the absolute scale of the free energy
densities of the leaky or classical cell theories depending on the model
assumptions Eyring and Hirschfelder (1937). However, with the exception of the
BCC model (see inset (b)) the existence of a phase transition does not depend
on whether the high-density or low-density calibration is used.

###  III.4 Calibration

The free energy density is calculated in terms of the exact free volume
through Eq. (4), which invokes the single occupancy assumption and therefore
can only expected to be accurate in the high-density regime. If one replaces
the single-occupancy assumption by allowing movement from cell to cell and
positing that the total free volume is simultaneously accessible to all
spheres as it would be in a low-density gas, the partition function is
increased by a factor of $e$ per molecule Eyring and Hirschfelder (1937). The
resulting expressions for the free energy density differ only up to an
additive constant, and because this additive constant is the same across all
lattice models (provided they are treated consistently) it does not affect the
predicted phase transitions. However, the choice of additive constant does
affect the phase transitions when lattice models are compared to other free
energy densities on an absolute scale, e.g. the Percus-Yevick or Carnahan-
Starling liquids. From a mathematical standpoint, this is because the common
tangent construction that defines these phase transitions involves the
absolute free energy density and not only its derivatives.

More generally, the free energy is often obtained by integrating the pressure
from a reference value Ponce and Renon (1976); Mirzaeinia, Feyzi, and
Hashemianzadeh (2017):

$\displaystyle\beta\frac{A}{N}-\beta\frac{A^{\text{IG}}}{N}=\int_{0}^{\eta}\frac{Z-1}{\eta}\mathrm{d}\eta,$
(12)

in which $A^{\text{IG}}$, the free energy of an ideal gas, satisfies
$A^{\text{IG}}=-1+3\log{\frac{\Lambda}{2R}}+\log{\frac{6\eta}{\pi}}$. These
two approaches for normalizing the free energy, i.e. whether one directly uses
the absolute free energy given by (4) or integrates the compressibility factor
from a low-density reference according to (12), are essentially a question of
whether one calibrates the free energy using either the low-density or high-
density reference value, an issue also considered by previous authors Alder,
Hoover, and Young (1968).

We find that the phase transitions may depend sensitively on this calibration,
and in some cases this calibration affects whether or not a transition exists
at all. Whereas the PY-FCC phase transition does not depend on this particular
assumption—i.e. regardless of whether the free energy resulting from Eq. (4)
is used or whether a factor of $e$ is added to the free energy density, the
same phase transitions exist—the actual location of the phase transition
varies by as much as 10%. Moreover, the result is robust to the choice of the
empirical model for the liquid state (Fig. 5).

In contrast, the freezing/melting transition in the BCC lattice is highly
sensitive to the choice of empirical model. If the free energy is calibrated
to high-density and the Carnahan-Starling equation of state is used, no
freezing/melting transition is predicted, whereas using the low-density
calibration or the Percus-Yevick liquid model leads to the existence of a
freezing transition at packing fractions of approximately 0.55–0.65 (Fig.
5(b)).

Figure 5 illustrates the range of free energy densities obtained through the
leaky cell model depending on whether one includes the additive constant. It
shows that this ambiguity does not affect whether a freezing/melting
transition exists between the FCC lattice and the Percus-Yevick or Carnahan-
Starling liquids. Further, it does not affect any intra- or inter-lattice
phase transitions, since the additive constant is applied uniformly to all
models and therefore the common tangents do not change. To summarize, we find
that some phase transitions e.g. liquid-FCC are robust in that they exist for
any combination of the liquid theory (PY or CS) and lattice free energy
calibration constant (1 or $e$), whereas others e.g. liquid-BCC exist for
certain pairs and not for others. These findings are consistent with previous
observations of the robustness of the phase transition to the underlying
approximations Colot, Baus, and Xu (1986).

###  III.5 Cell theory at low density

Next, to further explore the importance of calibration and to test whether the
choice of calibration may affect phase transitions between cell models, we
compare versions of the leaky cell model calibrated to low and high density to
analogous versions of the standard cell theory (Fig. 6). In the figure legend,
those expressions which use the low-density calibration are labeled as “low-
density”; otherwise the high-density calibration is used.

Figure 6:  Cell theory at low density. First row contains 1D results, second
row contains 2D results, and third row contains 3D results. (a)
Compressibility factor, (b) free energy density, and (c) quasi-communal
entropy (the latter two in units of $k_{B}T$).

In one-dimension, we find that the compressibility factor predicted by the
leaky cell model is exact, and calibrating the free energy to low-density
yields a free energy indistinguishable from the exact solution. In two and
three dimensions, the leaky cell model is no longer exact. Neither the low-
density or high-density calibrations are able to match the empirical equations
of state over an appreciable range of packing fractions.

This agreement can be improved somewhat by considering weighted averages of
the standard and leaky cell models—weighting the leaky model by $\frac{1}{2}$
and $\frac{1}{3}$ in the 2D and 3D cases, respectively—in a manner similar to
the weighting between the hard center and soft center models used in Buehler
et al. Buehler  _et al._ (1951). However, it is evident that no combination of
the cell models provides a strong agreement with the empirical model (unless
these weights are varied with $\eta$ and considered as fitting parameters, in
which case the cell models lose their straightforward interpretations). This
reveals a limitation of the cell theory. No combination of the standard and
leaky cell theories can on its own predict the freezing/melting transition in
the hard sphere gas. Some additional information is required to capture the
behavior at low densities.

###  III.6 Quasi-communal entropy

To address this limitation, inspired by the notion of communal entropy used to
capture the error incurred by using the single occupancy assumption at low
density Rice (1938); Hoover and Ree (1968), we augment the free volume by
incorporating a quasi-communal entropy $s=s(\eta)$ via:

$\mathcal{F}^{\prime}(\eta)=\mathcal{F}^{CT}(\eta)\exp{\left(s(\eta)\right)},$
(13)

where $\mathcal{F}^{CT}$ is the free volume in the standard cell theory.

We then fit $s$ to the empirical liquid equation of state using the boundary
condition $s(0)=1$ in the low-density limit to recover the additional factor
of $e$ per molecule Eyring and Hirschfelder (1937). The quasi-communal entropy
$s$ reports on the discrepancy between the standard cell theory and the
empirical models, which generally agree with the first several virial
coefficients and are therefore a reliable model at low-density. It differs
from the notion of communal entropy explored in the literatureKirkwood (1950);
Hoover and Ree (1968) in that it is an aggregate measure of the additional
entropy in the hard-sphere gas rather than the entropy incurred solely from
the single-occupancy assumption.

In the one-dimensional case, the quasi-communal entropy is identical to the
communal entropy reported in Hoover and AlderHoover and Alder (1966), which
used systematic calculations to isolate the effects of the assumptions in cell
theory noted by Kirkwood Kirkwood (1950). Interestingly, the quasi-communal
entropy found in two and three dimensions is a non-monotonic function. As
mentioned previously, in addition to single-occupancy the cell theory also
assumes independence of neighboring lattice sites, which may be more realistic
in the low-density regime Kirkwood (1950). Presumably, the non-monotonicity
observed in the quasi-communal entropy corresponds in Kirkwood’s description
of the errors that arise from increased correlations between neighboring
cells, which has the effect of increasing the effective free volumes.

To better understand the quasi-communal entropy and how it compares with the
actual communal entropy, we revisit Tonks’ gas in one dimension, for which we
can compute the free volume exactly. We first review the formulation of Tonks’
gas model, following closely Tonks’ original paperTonks (1936). Consider a
segment of length $\ell$, which hosts $N$ rods of diameter $\sigma$, which
glide on the segment without going past each other. We call

$\lambda:=\frac{\ell}{N}$ (14)

the _specific length_ , and

$\eta:=\frac{\sigma N}{\ell}=\frac{\sigma}{\lambda}$ (15)

the packing fraction as above. In the thermodynamic limit, $N\gg 1$, the
entropy $S_{\mathrm{T}}$ of the system is

$S_{\mathrm{T}}=k_{B}TN\left[1+\ln\lambda(1-\eta)\right].$ (16)

Letting the force $f$ be defined as

$f:=T\frac{\partial S_{\mathrm{T}}}{\partial\ell},$ (17)

we arrive at the following expression for the compressibility factor

$Z:=\frac{f\ell}{Nk_{B}T}=\frac{1}{1-\eta}.$ (18)

Next, we consider Tonks’ gas in the context of cell theory, closely following
the presentation of RiceRice (1944). In this formulation, each rod is assumed
to be caged in a cell of length $\lambda$ (there are exactly $N$ cells, each
confining a single rod). It is a simple matter to derive the entropy
$S_{\mathrm{R}}$ of the system in this case,

$S_{\mathrm{R}}=k_{B}TN\ln\lambda(1-\eta),$ (19)

so that the _communal_ entropy $S_{\mathrm{c}}$ is simply

$S_{\mathrm{c}}:=S_{\mathrm{T}}-S_{\mathrm{R}}=k_{B}TN$ (20)

and the compressibility factor $Z$ is still delivered by (18) (because the
communal entropy is constant).

Next, we consider a leaky version of the cell theory above. This was first
done by Hoover and Alder Hoover and Alder (1966), in which the authors
considered a first version of leaky cells for Tonks’ gas. They allowed the
center of mass of each rod (and not just the whole rod) to be confined in a
cell and found that, as expected, the communal entropy in this case is smaller
than in the case of full confinement. Moreover, the compressibility factor is
affected, as the communal entropy is no longer constant.

Here, we take a more general view and consider the case in which the center of
mass of each rod can get as far as $(\alpha-\frac{1}{2})\sigma$ (either way)
from the end-points of the cell that bounds it. Here $0\leqq\alpha\leqq 1$ is
a fixed _leak parameter_. Formally, the coordinates $x_{j}$ of the centers of
mass of the rods obey the following inequalities,

$(j-1)\lambda-\left(\alpha-\frac{1}{2}\right)\sigma\leqq x_{j}\leqq
j\lambda+\left(\alpha-\frac{1}{2}\right)\sigma,\quad 2\leqq j\leqq N-1.$ (21)

The case $\alpha=0$ corresponds to the model of RiceRice (1944), whereas the
case $\alpha=\frac{1}{2}$ corresponds to the model of Hoover and AlderHoover
and Alder (1966). The limiting value $\alpha=1$ is the one where each rod is
almost allowed to leave its cell, as it can remain in contact with it at only
one point.

Figure 7: (a) Communal entropy vs. packing fraction in the interval (22) in
which (24) is valid, plotted for several values of the leak parameter; the red
curve corresponds to $\alpha=\frac{1}{2}$. The brown straight line is the
graph of the linear function in (27). (b) Corresponding graphs of the
compressibility factor $Z$ given by (25); the red curve corresponds to
$\alpha=\frac{1}{2}$.

As shown in detail in Appendix D, it turns out that, at least in the interval

$0\leqq\eta\leqq\frac{1}{2\alpha+1},$ (22)

the entropy $S$ of the system can be computed exactly (in the thermodynamic
limit). It is given by

$S=k_{B}TN\ln\frac{1}{2}\lambda\left(1+(2\alpha-1)\eta+\sqrt{[1+(2\alpha-1)\eta]^{2}-8\alpha^{2}\eta^{2}}\right),$
(23)

so that the communal entropy becomes

$S_{\mathrm{c}}=k_{B}TN\ln\frac{2(1-\eta)\mathrm{e}}{1+(2\alpha-1)\eta+\sqrt{[1+(2\alpha-1)\eta]^{2}-8\alpha^{2}\eta^{2}}}.$
(24)

Correspondingly, the compressibility factor reads as

$Z=\frac{1}{\sqrt{1-2\eta+4\alpha\eta+(1-4\alpha-4\alpha^{2})\eta^{2}}}.$ (25)

Both these functions are plotted in Fig. 7 for different values of $\alpha$,
the red graphs corresponding to $\alpha=\frac{1}{2}$. A simple computation
shows that at the upper limit of validity for (24), that is, at
$\eta=1/\left(2\alpha+1\right)$,

$\frac{S_{\mathrm{c}}}{k_{B}TN}=1+\ln\left(\frac{2}{2+\sqrt{2}}\right)\approx
0.47.$ (26)

This is precisely the lower bound of the window in which
$S_{\mathrm{c}}/k_{B}TN$ is plotted in Fig. 7(a).

Approximating the communal entropy for $\alpha=\frac{1}{2}$ by the linear
function

$S_{\mathrm{c}}\approx k_{B}TN(1-\eta),$ (27)

as suggested by the comparison between red and brown graphs in Fig. 7(a), one
obtains the following approximate form for the compressibility factor,

$Z_{\mathrm{a}}=\frac{1}{\sqrt{1-2\eta^{2}}}+\eta,$ (28)

which is plotted (in brown) in Fig. 8 against the exact compressibility factor
from (18) (in blue).

Figure 8: Comparison between the approximate compressibility factor
$Z_{\mathrm{a}}$ in (28) and the exact compressibility factor given by (18)
for $\alpha=\frac{1}{2}$.

The calculation above makes clear some of the ways in which the quasi-communal
entropy differs from the bona fide communal entropy. In particular, the true
communal entropy is a monotonic function bounded between 0 and 1 for any value
of the leak parameter. The fact that the quasi-communal entropy does not
satisfy these criteria reveals that it cannot be interpreted as the error
incurred by the single-occupancy assumption. Rather, it is properly
interpretated as a net error that includes not only the single-occupancy
assumption but also errors due to other approximations e.g. that the free
volumes of neighboring cells are uncorrelated.

##  IV Discussion

In this work, we extend the classical cell theory into an intermediate leaky
regime in which the local neighborhood of spheres expands, with consequences
on the resulting free volume. The relevant leaky packing fractions are
$\eta_{l}\approx 0.26$ for the FCC lattice, $\eta_{l}\approx 0.29$ for the SC
lattice, and $\eta_{l}\approx 0.37$ for the BCC lattice. Below this leaky cell
fraction and above the corresponding percolation thresholds, spheres remain
locally caged by their neighbors but the cage involves additional neighbors.
Interestingly, the effect of the leaky cell model is not uniform for all three
of these lattices. In the FCC and SC lattices, it leads to additional free
volume, whereas in the BCC lattice it leads to a decrease. This is a direct
consequence of the geometry of the lattice, as we show explicitly through
exact calculations of the free volume through the leaky regime.

In the leaky cell model explored here, we allow for an intermediate regime in
which spheres may escape their Voronoi cells while they remain sterically
trapped within an expanded cage. By expressing the free energy in terms of the
leaky cell volume, which may be calculated exactly, we apply the leaky cell
model to various lattice arrangments, and predict phase transitions between
lattices in the leaky regime.

We have identified two interlattice phase transitions within this intermediate
leaky regime. There are SC-FCC and SC-BCC phase transitions around
$\eta=0.22$. In addition, there is a quasi intralattice phase transition
around $\eta=0.3$ within the SC lattice, which arises because of a
discontinuity in the free energy density. These phase transitions indicate the
possibility that different lattice structures may coexist within crystalline
materials. Such cubic-cubic phase transitions have previously been observed
experimentally in fullerene crystals David  _et al._ (1992); Wochner  _et al._
(1997), and the leaky model may be an apt description of such systems.
However, it has been shown previouslySalsburg and Wood (1962) that not all
lattice configurations are stable, and that stability to shear requires
certain topological constraints on the number of contacts to be satisfied.
Therefore, physical forces beyond those of pure hard sphere models may be
required to realize some of the lattices discussed in this work.

An early critique of the cell theory was that the mean-field assumptions
underlying the model, such as the single occupancy of cells and the
independence of neighboring free volumes, were not always clearly specified
Kirkwood (1950); Wood (1952); Hoover and Alder (1966). In particular,
depending on the assumptions made, the resulting free energy may differ by an
additive constant. In this work we consider the effect of such translations in
the free energy in the context of the leaky cell model. We have found that
some of the phase transitions found are sensitive to this calibration,
particularly the freezing/melting transition in the BCC lattice, in which the
existence of a phase transition is sensitive to the choice of liquid model and
calibration. On the other hand, we find that predicting the known
freezing/melting transition in a hard sphere gas requires information beyond
the cell theory, such as an empirical liquid model (e.g. Carnahan-Starling):
regardless of whether it is calibrated to low or high density or some simple
combination, when the liquid model is based purely on a cell theory we find no
freezing transition reminiscent of the hard sphere gas.

Finally, we comment on possible extensions of this model. Although hard sphere
models are unable to capture the behavior of systems such as liquid crystals
in which the particles are rod-shaped, this work represents a step along the
path to understanding how cell theories based on free volume may be used in
non-spherical systems, as has recently been done in the context of predicting
the glass transition Gujrati (2021). Although it may not be possible to obtain
exact formulas in these more complicated settings, phase transitions could be
identified by using numerical methods to approximate the free energy, as has
been done previously using Monte Carlo simulation Barroso and Ferreira (2002);
Schilling and Schmid (2009). The formulas presented here are strictly valid on
lattices with periodic boundary conditions. In the case of a finite domain
with walls, it may be possible to calculate finite size corrections for the
leaky cell model as has been done for other lattice structures Polson  _et
al._ (2000).

###### Acknowledgements.

We thank Zachary M. Geballe for reading an early draft of this manuscript and
providing helpful feedback. T.G.F. acknowledges support from National Science
Foundation grant DMS-1913093. P.P-M. acknowledges support from the Office of
Naval Research (ONR N00014-18-1-2624). J.M.T. has been partially supported by
the Basque Government through the BERC 2018-2021 program; and by Spanish
Ministry of Economy and Competitiveness MINECO through BCAM Severo Ochoa
excellence accreditation SEV-2017-0718 and through project MTM2017-82184-R
funded by (AEI/FEDER, UE) and acronym “DESFLU". Finally, we wish to
acknowledge the support of the Institute for Mathematics and its Applications
(IMA), where this work was initiated during the 2018 program on “Multiscale
Mathematics and Computing in Science and Engineering”.

## Data Availability

The code used to perform the calculations presented here is available on
GitHub: https://github.com/thomasgfai/LeakyCellModel.

##  Appendix A  Exact free volumes on 3D lattices

To calculate exact free volumes accessible to spheres, we use the following
geometric idea. From the volume of the void space made up of a central sphere
and its neighbors, we first subtract the volumes of the exclusion spheres of
radius $2R$ with their corresponding solid angles. If the exclusion spheres
overlap, we must correct for this double-counting by adding the volumes of
double intersections between exclusion spheres. If the double intersections
overlap, we must subtract a triple intersection, and so on. Because there are
general formulas for the volumes of spherical double, triple, and quadruple
intersections Gibson and Scheraga (1987a, b), this method leads to an exact
formula for the free volume that can in be applied to any lattice.

###  A.1 FCC/HCP lattice

The FCC lattice involves planar layers of hexagonally-arranged spheres. The
surrounding void space can be partitioned into $N_{t}=8$ tetrahedra and
$N_{o}=6$ octahedra. As the packing fraction increases, a centrally-located
sphere becomes caged first by triangular faces (with 18 neighboring spheres)
and subsequently by octahedral midplanes (with 12 neighboring spheres). See
Fig. 1(a)–(b). Note that the cage in the left column of Fig. 1 is always a
subset of the corresponding cage in the right column. For example, for the FCC
lattice illustrated in the first row of Fig. 1, the cage in the right column
is obtained by adding outward-pointing right square pyramids to the cage in
the left column.

We express the free volume in terms of the edge length $a$. It is
straightforward to write these in terms of the Voronoi cell volume, as the
edge length $a$ is related to the Voronoi cell volume $v$ through the relation

$v=\frac{N_{t}}{4}V_{t}+\frac{N_{o}}{6}V_{o},$ (29)

where the numerical factors $4$ and $6$ come from the number of vertices of a
tetrahedron and octahedron and where $V_{t}=a^{3}/\left(6\sqrt{2}\right)$ and
$V_{o}=\sqrt{2}a^{3}/3$ are their respective volumes. It follows that

$v(a)=\frac{a^{3}}{\sqrt{2}}.$ (30)

We use the notation $V_{s}:=\frac{4}{3}\pi(2R)^{3}$ to denote the volume of a
sphere of radius $2R$. The free area $\mathcal{F}(v)$ on the FCC lattice
satisfies

$\mathcal{F}(v)=\begin{cases}0,\qquad&a<2R\\\ \\\\[-14.45377pt]
N_{t}\left(V_{t}-3V_{s}\frac{\Omega_{t}}{4\pi}+3V_{2I}(a)\frac{\alpha_{t}}{2\pi}\right)\qquad&\\\
\quad+N_{o}\left(V_{o}/2-2V_{s}\frac{\Omega_{o}}{4\pi}+2V_{2I}(a)\frac{\alpha_{o}}{2\pi}\right.\\\
\left.\quad+V_{2I}(\sqrt{2}a)-2V_{3I}^{\text{rt}}(a)+\frac{1}{2}V_{4I}(a)\right)\\\
\quad-4V_{3I}^{\text{eq}}(a),\qquad&2R<a<2\sqrt{2}R\\\ \\\\[-14.45377pt]
N_{t}\left(V_{t}-3V_{s}\frac{\Omega_{t}}{4\pi}+3V_{2I}(a)\frac{\alpha_{t}}{2\pi}\right)\qquad&\\\
\quad+N_{o}\left(V_{o}-5V_{s}\frac{\Omega_{o}}{4\pi}+8V_{2I}(a)\frac{\alpha_{o}}{2\pi}\right)\\\
\quad-16V_{3I}^{\text{eq}}(a),\qquad&2R<a<2\sqrt{3}R\\\ \\\\[-14.45377pt]
N_{t}\left(V_{t}-3V_{s}\frac{\Omega_{t}}{4\pi}+3V_{2I}(a)\frac{\alpha_{t}}{2\pi}\right)\qquad&\\\
\quad+N_{o}\left(V_{o}-5V_{s}\frac{\Omega_{o}}{4\pi}+8V_{2I}(a)\frac{\alpha_{o}}{2\pi}\right),\qquad&2\sqrt{3}R<a<4R\\\
\\\\[-14.45377pt]
N_{t}\left(V_{t}-3V_{s}\frac{\Omega_{t}}{4\pi}\right)\qquad&\\\
\quad+N_{o}\left(V_{o}-5V_{s}\frac{\Omega_{o}}{4\pi}\right),\qquad&4R<a,\end{cases}$
(31)

where $\Omega_{t}=\cos^{-1}(23/27)$, $\Omega_{o}=4\sin^{-1}(1/3)$ are the
solid angles subtended by the vertices of regular tetrahedra and octahedra and
$\alpha_{t}=\cos^{-1}(1/3)$ and $\alpha_{o}=\cos^{-1}(-1/3)$ are the
corresponding dihedral angles. Note that beyond the percolation threshold, at
which $a>2\sqrt{3}R$, the sphere can escape through tetrahedral faces.
However, we do not count the volume outside of the cage so that free volume
remains an intensive quantity and the system size $L$ does not enter the
formula.

In addition to the double intersection volume $V_{2I}(a)$, in Eq. (31) we have
used the following notation: $V_{3I}^{\text{eq}}(a)$ is the intersection
volume of three spheres at the vertices of an equilateral triangle of edge-
length $a$, $V_{3I}^{\text{rt}}(a)$ is the intersection volume of three
spheres along a right triangle at vertices of a square of edge-length $a$, and
finally $V_{4I}(a)$ is the quadruple intersection volume of four spheres at
the vertices of a square of edge-length $a$.

The volume of intersection between two spheres of radius $2R$ separated by a
distance $a$ is given by

$V_{2I}(a)=\frac{\pi}{12}\left(4(2R)+a\right)\left(2(2R)-a\right)^{2}.$ (32)

General formulas for volumes of higher-order intersections are given in Gibson
and Scheraga Gibson and Scheraga (1987a, b). In particular,

$\displaystyle V_{3I}^{\text{eq}}(a)$
$\displaystyle=\frac{a^{2}w^{\text{eq}}}{6}-\frac{3a}{2}\left(2(2R)^{2}-\frac{a^{2}}{6}\right)\tan^{-1}\left(\frac{2w^{\text{eq}}}{a}\right)$
$\displaystyle\quad+4(2R)^{3}\tan^{-1}\left(\frac{w^{\text{eq}}}{2R}\right)$
(33)

where

$w^{\text{eq}}=\sqrt{3(2R)^{2}-a^{2}}\quad\qquad\qquad\qquad\qquad\qquad$

and

$\displaystyle V_{3I}^{\text{rt}}(a)$
$\displaystyle=\frac{a^{2}w^{\text{rt}}}{6}-a\left(2(2R)^{2}-\frac{a^{2}}{6}\right)\tan^{-1}\left(\frac{w^{\text{rt}}}{a}\right)$
$\displaystyle\quad-\frac{a\sqrt{2}}{2}\left(2(2R)^{2}-\frac{a^{2}}{3}\right)\frac{\pi}{2}$
(34)
$\displaystyle\quad+\frac{4(2R)^{3}}{3}\left(2\tan^{-1}\left(\frac{w^{\text{rt}}}{2(2R)}\right)+\frac{\pi}{2}\right)$

where

$w^{\text{rt}}=\sqrt{4(2R)^{2}-2a^{2}}.\quad\qquad\qquad\qquad\qquad$

As explained in (30), the quadruple intersection volume may be computed simply
by

$V_{4I}(a)=2V_{3I}^{\text{rt}}(a)-V_{2I}(\sqrt{2}a).$ (35)

###  A.2 BCC lattice

In the BCC lattice the percolation threshold occurs at
$\eta^{\text{BCC}}_{\text{p}}\approx 0.19$, at which spheres are caged by 14
neighbors, located both at vertices of a cube as well as across cubic faces.
The void space is contained within 6 octahedra (regular square pyramids
mirrored at the base). At the leaky cell transition
$\eta^{\text{BCC}}_{\text{l}}\approx 0.19$, spheres can no longer escape
through the square octahedral midplanes, so that the cage is formed by the
eight neighbors on cube vertices. Unlike the other lattices, the leaky cell
model in the BCC lattice leads to lower free volumes since neighbors across
cubic faces occlude volume via their exclusion spheres leaking into the
Wigner-Seitz cell.

Upon making the identification $a(v)=\sqrt[3]{2v}$, $\mathcal{F}(v)$ may be
written as

$\mathcal{F}(v)=\begin{cases}0,&a<\frac{4}{\sqrt{3}R}\\\ \\\\[-14.45377pt]
a^{3}-V_{s}+3V_{2I}(a)+6V_{2I}(\sqrt{2}a)&\\\
\quad-12V_{3I}^{\text{eq}}(a)+3V_{4I}(a),&\frac{4}{\sqrt{3}R}<a<\frac{8}{3}R\\\
\\\\[-14.45377pt]
N_{o}V_{o}-N_{o}\left(4V_{s}\frac{2\Omega_{b}}{4\pi}+V_{s}\frac{\Omega_{a}}{4\pi}\right.&\\\
\quad\left.+4V_{2I}(a)\frac{2\alpha_{b}}{2\pi}+4V_{2I}(\frac{\sqrt{3}}{2}a)\frac{\alpha_{a}}{2\pi}\right.&\\\
\quad\left.+2V_{2I}(\sqrt{2}a)-2V_{3I}^{\text{iso}}(a,\frac{\sqrt{3}}{2}a)\right.&\\\
\quad\left.-4V_{3I}^{\text{rt}}(a)-2\widetilde{V}_{3I}^{\text{iso}}(\sqrt{2}a,\frac{\sqrt{3}}{2}a)\right.&\\\
\quad\left.+4\widetilde{V}^{\text{pyr}}_{4I}(a,\frac{\sqrt{3}}{2}a)\right)&\\\
\quad\left.+4V^{\text{sq}}_{4I}(a)-4\widetilde{V}^{\text{pyr}}_{5I}(a,\frac{\sqrt{3}}{2}a)\right),&\frac{8}{3}R<a<2\sqrt{2}R\\\
\\\\[-14.45377pt]
N_{o}V_{o}-N_{o}\left(4V_{s}\frac{2\Omega_{b}}{4\pi}+V_{s}\frac{\Omega_{a}}{4\pi}\right.&\\\
\quad\left.+4V_{2I}(a)\frac{2\alpha_{b}}{2\pi}+4V_{2I}(\frac{\sqrt{3}}{2}a)\frac{\alpha_{a}}{2\pi}\right.&\\\
\quad\left.-2V_{3I}^{\text{iso}}(a,\frac{\sqrt{3}}{2}a)\right),&2\sqrt{2}R<a<\frac{8\sqrt{2}}{3}R\\\
\\\\[-14.45377pt]
N_{o}V_{o}-N_{o}\left(4V_{s}\frac{2\Omega_{b}}{4\pi}+V_{s}\frac{\Omega_{a}}{4\pi}\right.&\\\
\quad\left.+4V_{2I}(a)\frac{2\alpha_{b}}{2\pi}+4V_{2I}(\frac{\sqrt{3}}{2}a)\frac{\alpha_{a}}{2\pi}\right),&\frac{8\sqrt{2}}{3}R<a<4R\\\
\\\\[-14.45377pt]
N_{o}V_{o}-N_{o}\left(4V_{s}\frac{2\Omega_{b}}{4\pi}+V_{s}\frac{\Omega_{a}}{4\pi}\right.&\\\
\quad\left.+4V_{2I}(\frac{\sqrt{3}}{2}a)\frac{\alpha_{a}}{2\pi}\right),&4R<a<\frac{8}{\sqrt{3}}R\\\
\\\\[-14.45377pt]
N_{o}V_{o}-N_{o}\left(4V_{s}\frac{2\Omega_{b}}{4\pi}+V_{s}\frac{\Omega_{a}}{4\pi}\right),&\frac{8}{\sqrt{3}}R<a.\end{cases}$
(36)

Here, $\Omega_{a}=2\pi/3$ and $\Omega_{b}=\pi/6$ are the solid angles
subtended by the apex and base, respectively, of the regular square pyramid,
and $\alpha_{a}=\cos^{-1}(-\frac{1}{2})$ and $\alpha_{b}=\pi/4$ are the
corresponding dihedral angles. Further, $V_{3I}^{\text{iso}}(x,y)$ is the
triple intersection between three spheres at vertices of an isosceles triangle
with edge lengths $x$, $y$, and $y$. To calculate $V_{3I}^{\text{iso}}$, we
use the general formula from (30):

$\displaystyle V_{3I}^{\text{iso}}(x,y)$
$\displaystyle=\frac{a^{2}w^{\text{rt}}}{6}-a\left(2(2R)^{2}-\frac{a^{2}}{6}\right)\tan^{-1}\left(\frac{w^{\text{rt}}}{a}\right)$
$\displaystyle\quad-\frac{a\sqrt{2}}{2}\left(2(2R)^{2}-\frac{a^{2}}{3}\right)\frac{\pi}{2}$
(37)
$\displaystyle\quad+\frac{4(2R)^{3}}{3}\left(2\tan^{-1}\left(\frac{w^{\text{rt}}}{2(2R)}\right)+\frac{\pi}{2}\right)$

where

$w^{\text{rt}}=\sqrt{4(2R)^{2}-2a^{2}}.\quad\qquad\qquad\qquad\qquad$

The remaining intersections appearing in (36) simplify upon inspection; one of
the intersections does not exclude any additional volume, and all the
intersections marked with tildes are equal to lower-order intersections. In
particular, the triple intersection
$\widetilde{V}_{3I}^{\text{iso}}(\sqrt{2}a,\frac{\sqrt{3}}{2}a)$ formed
between the spheres at apex of the regular right pyramid and vertices along a
diagonal at its base is in fact a double intersection, i.e.

$\widetilde{V}_{3I}^{\text{iso}}(\sqrt{2}a,\frac{\sqrt{3}}{2}a)=V_{2I}(\sqrt{2}a).$

Similarly, the quadruple intersection
$\widetilde{V}^{\text{pyr}}_{4I}(a,\frac{\sqrt{3}}{2}a)$ formed between the
sphere at the apex of the regular square pyramid and three spheres at its base
reduces to a triple intersection, so that

$\widetilde{V}^{\text{pyr}}_{4I}(a,\frac{\sqrt{3}}{2}a)=V_{3I}^{\text{rt}}(a).$

Finally, the quintuple intersection
$\widetilde{V}^{\text{pyr}}_{5I}(a,\frac{\sqrt{3}}{2}a)$ formed between the
sphere at the apex of the regular square pyramid and all four spheres at its
base is simply a quadruple intersection, i.e.

$\widetilde{V}^{\text{pyr}}_{5I}(a,\frac{\sqrt{3}}{2}a)=V^{\text{sq}}_{4I}(a).$

###  A.3 SC lattice

In addition to the percolation threshold at
$\eta^{\text{SC}}_{\text{p}}\approx 0.19$, at which diagonal neighbors on
square faces begin to intersect and spheres becomes caged by 26 neighbors,
there is a transition in the cubic lattice at
$\eta^{\text{SC}}_{\text{l}}\approx 0.29$ beyond which the sphere can no
longer escape through triangular faces, so that opposite corners of cubes are
lost and the number of neighbors decreases to 6.

Upon making the identification $a(v)=\sqrt[3]{v}$, $\mathcal{F}(v)$ may be
written as

$\mathcal{F}(v)=\begin{cases}0,\qquad&a<2R\\\
\frac{4}{3}a^{3}-6V_{s}\frac{\Omega_{o}}{4\pi}+12V_{2I}(\sqrt{2}a)\frac{\alpha_{o}}{2\pi}\\\
\quad-4V_{3I}^{\text{eq}}(\sqrt{2}a),\qquad&2R<a<\sqrt{6}R\\\
8a^{3}-7V_{s}+18V_{2I}(a)+36V_{2I}(\sqrt{2}a)\\\
\quad-60V_{3I}^{\text{rt}}(a)+12V_{4I}(a),\qquad&\sqrt{6}R<a<2\sqrt{2}R\\\
8a^{3}-7V_{s}+18V_{2I}(a),\qquad&2\sqrt{2}R<a<4R\\\
8a^{3}-7V_{s},\qquad&4R<a,\end{cases}$ (38)

where as before $\Omega_{o}=4\sin^{-1}(1/3)$ is the solid angle subtended by
the vertices of a regular octahedron, $\alpha_{o}=\cos^{-1}(-1/3)$ is the
corresponding dihedral angle, and the quadruple intersection may be computed
by $V_{4I}(a)=2V_{3I}^{\text{rt}}(a)-V_{2I}(\sqrt{2}a)$. Note that there is a
discontinuity in the second derivative of $\mathcal{F}(v)$ at the leaky cell
transition.

##  Appendix B  Asymptotic approximation of the free volume

Consider a general lattice, which is jammed at its densest packing
$\eta_{\text{cp}}$. We illustrate the situation in Fig. 9 using a square
lattice in 2D, where we show the centers of mass of the centre particle and
its neighbors (black dots), and the exclusion regions due to its nearest
neighbors (black circles). We approximate the void space by a polygon. That
is, we approximate the free volume of a particle at $0$ with nearest neighbors
at $a\mathbf{e}_{i}$, for $a$ a lattice spacing parameter and $\mathbf{e}_{i}$
unit vectors by replacing the exclusion regions
$\|a\mathbf{e}_{i}-\mathbf{x}\|<2r$ with half-spaces,
$\mathbf{e}_{i}\cdot\mathbf{x}>a-2r$. In the schematic Fig. 9, the boundaries
of these half spaces are given by the dashed lines, and the polygonal cavity
corresponds to the square contained within. First of all, we note that this
polygonal cavity will simply be a uniform scaling of the Voronoi cell of the
densest packing. We can see this because the bounding hyperplanes correspond
to the shared faces of the Voronoi cells of particles at $0$ and
$2(a-2r)\mathbf{e}_{i}$, which is simply a rescaled version of our lattice.
These “virtual point particles" bounding the Voronoi cells are given by grey
squares in Fig. 9. The characteristic width of this Voroni cell is thus
$2(a-2r)$, in comparison to $2r$ at the densest packing. Thus the scale factor
between the two Voronoi cells is given by $2\left(\frac{a}{2r}-1\right)$, and
thus the volume $\mathcal{F}$ of the polygonal cavity is given by

$\mathcal{F}(a)\approx
2^{d}\left(\frac{a}{2r}-1\right)^{d}v_{\text{cp}}=\frac{2^{d}\left(\frac{a}{2r}-1\right)^{d}}{\rho_{\text{cp}}},$
(39)

where $v_{\text{cp}}$ is the volume of the Voronoi cell in the densest
packing, and $\rho_{\text{cp}}$ is the number density at densest packing,
related by $\rho_{\text{cp}}^{-1}=v_{\text{cp}}$. The number density of a
lattice with spacing $a$ is readily given by

$\rho=\left(\frac{2r}{a}\right)^{d}\rho_{\text{cp}},$ (40)

by noting that $\frac{a}{2r}$ is the scale factor between the densest packing
and the lattice with parameter $a$. We may substitute this into our previous
expression to give

$\mathcal{F}(\rho)\approx\frac{2^{d}\left(\sqrt[d]{\frac{\rho_{\text{cp}}}{\rho}}-1\right)^{d}}{\rho_{\text{cp}}}.$
(41)

Remarkably, this only depends on the geometry of the lattice via the densest
packing parameter.

It is worth noting that as both $\mathcal{F}$ and this approximation,
$\mathcal{F}_{\text{approx}}$ converge to $0$ as
$\rho\nearrow\rho_{\text{cp}}$, this means that

$\lim\limits_{\rho\to\rho_{\text{cp}}^{-}}\frac{\mathcal{F}(\rho)}{\mathcal{F}_{\text{approx}}(\rho)}=\lim\limits_{\rho\to\rho_{\text{cp}}^{-}}\frac{\mathcal{F}^{\prime}(\rho)}{\mathcal{F}_{\text{approx}}^{\prime}(\rho)},$
(42)

by l’Hôpital’s rule, so that if $v_{ap}$ is a good approximation to 0th-order
(in the sense that the left-hand limit in the previous equation equals $1$),
it is automatically a good approximation to 1st-order (in the sense the right-
hand limit equals $1$). It is precisely this approximation which is used in
Fig. 2 of the Main Text.

Figure 9: Schematic of the polygonal approximation

##  Appendix C Exact free volumes on 2D lattices

In this appendix, we provide the analytic free volume formulas for lattice
packings in 2D using the same geometric argument applied in 3D.

###  C.1 Hexagonal lattice

For a hexagonal lattice of disks of radius $R$, the percolation threshold
$\eta^{\text{HE}}_{\text{p}}$ is reached when the lattice edge length $a$
satisfies $a=4R$, corresponding to $\eta^{\text{HE}}_{\text{p}}\approx 0.23$.
The close-packing limit $\eta^{\text{HE}}_{\text{cp}}$ is attained when
$a=2R$, corresponding to $\eta^{\text{HE}}_{\text{cp}}\approx 0.91$, at which
point the disk is in contact with its neighbor so that it has no free volume
to explore. Unlike any of the lattices investigated in the 3D case, there is
no leaky regime in the hexagonal 2D lattice; disks escape from their unit
cells at the percolation threshold, and their is no emergence of an expanded
cage that involves additional neighbors. The critical packing fractions of the
2D cell theory are summarized in Table 3.

Table 3: Percolation ($\eta_{\text{p}}$), leaky ($\eta_{\text{l}}$), and close-packed ($\eta_{\text{cp}}$) packing fractions in the hexagonal and square 2D lattices. 2D Lattice | $\eta_{\text{p}}$ | $\eta_{\text{l}}$ | $\eta_{\text{cp}}$
---|---|---|---
hexagonal (HE) | $\frac{\pi}{8\sqrt{3}}\approx 0.23$ | N/A | $\frac{\pi}{2\sqrt{3}}\approx 0.91$
square (SQ) | $\pi/16\approx 0.20$ | $\frac{\pi}{8}\approx 0.39$ | $\frac{\pi}{4}\approx 0.79$

Given a hexagonal lattice of $N$ disks in a 2D $L\times L$ domain, using the
identification $a(v)=\sqrt{\frac{2v}{\sqrt{3}}}$ that relates edge lengths and
local volumes in a hexagonal lattice, the free volume satisfies

$\mathcal{F}(v)=\begin{cases}0,\qquad&a<2R\\\
3v-2V_{s}+3V_{2I}(a),\qquad&2R<a<4R\\\ 3v-2V_{s},\qquad&4R<a,\end{cases}$ (43)

where $V_{s}=\pi(2R)^{2}$ is the area of a disk of radius $2R$ and where the
area of intersection of two such disks separated by a distance $a$ is given by

$V_{2I}(a):=2(2R)^{2}\cos^{-1}\left(\frac{a}{4R}\right)-\frac{a}{2}\sqrt{4(2R)^{2}-a^{2}}.$

###  C.2 Square lattice

Next, we provide the corresponding formula for the free area in a two-
dimensional square lattice. The percolation transition
$\eta^{\text{SQ}}_{\text{p}}$ at which $a=4R$ satisfies
$\eta^{\text{SQ}}_{\text{p}}\approx 0.20$. The close-packing fraction
$\eta^{\text{SQ}}_{\text{cp}}$ at which $a=2R$ occurs at
$\eta^{\text{SQ}}_{\text{cp}}\approx 0.79$. In addition, there is a leaky cell
transition $\eta^{\text{SQ}}_{\text{l}}\approx 0.39$ below which the number of
neighbors caging each disk expands from 8 to 4. Above the leaky cell fraction,
which corresponds to edge lengths of $a=2\sqrt{2}R$, disks becomes caged by
their orthogonal neighbors (i.e. the corners are lost).

Upon making the identification $a(v)=\sqrt{v}$ for a square lattice,
$\mathcal{F}(v)$ may be written as

$\mathcal{F}(v)=\begin{cases}0,\qquad&a<2R\\\
2v-V_{s}+2V_{2I}(\sqrt{2}a),\qquad&2R<a<2\sqrt{2}R\\\
4v-3V_{s}+4V_{2I}(a),\qquad&2\sqrt{2}R<a<4R\\\
4v-3V_{s},\qquad&a>4R.\end{cases}$ (44)

Note that, in contrast to the smooth free volume derived in the hexagonal
case, there is a discontinuity in the second derivative of $\mathcal{F}(v)$ at
the leaky cell transition.

Figure 10 shows the resulting compressibility factors computed through Eq.
(9). It is interesting to note that whereas the compressibility factor for the
hexagonal lattice is a smooth function of $\eta$, there is a kink in the
square lattice compressibility factor that arises from a discontinuity in the
second derivative of $\mathcal{F}(v)$ at the leaky cell transition as well as
a discontinuity at the percolation transition.

Figure 10: Comparison of 2D compressibility factors on square and hexagonal
lattices and the scaled particle theory (SPT) equation of state.

###  C.3 Comparison to 2D scaled particle theory

Finally, we compare the compressibility factors computed from the above
formulas to the 2D scaled particle theory equation of state:

$Z(\eta)=\frac{1}{(1-\eta)^{2}}.$ (45)

As shown in Fig. 10, we find that the hexagonal lattice compressibility factor
is in reasonable agreement with the scaled particle theory, whereas the square
lattice compressibility factor is significantly different and has distinctive
qualitative features, such as non-monotoniticity and a loss of smoothness at
$\eta\approx 0.4$ and $\eta\approx 0.2$.

##  Appendix D  Exact leaky model on 1D lattices: Tonks’ gas

Recall that for $N$ rods, each of length $\sigma$, confined within a segment
of length $\ell$, we call

$\lambda:=\frac{\ell}{N}$ (46)

the _specific length_ and

$\eta:=\frac{\sigma N}{\ell}=\frac{\sigma}{\lambda}$ (47)

is the packing fraction as before. We imagine the segment partitioned in $N$
equal permeable cells, each of length $\lambda$. We denote by $x_{j}$ the
coordinate of the center of mass of the $j$-th rod and, for a given parameter
$0\leqq\alpha\leqq 1$, we assume that

$\displaystyle(j-1)\lambda-\left(\alpha-\frac{1}{2}\right)\sigma\leqq x_{j}$
$\displaystyle\leqq j\lambda+\left(\alpha-\frac{1}{2}\right)\sigma,$
$\displaystyle\quad 2\leqq j\leqq N-1.$ (48)

We call $\alpha$ the _leak_ parameter, as $\alpha\sigma$ is the maximum length
each end-point of every rod can go outside the ideal restraining cell.

To compute the entropy $S_{N}$ of such a system, we start by computing the
partition function $Z_{N}$, which is defined by the following integral,

$\displaystyle Z_{N}:=$ $\displaystyle\int_{(N-1)\lambda-\beta\sigma}^{\lambda
N-\frac{\sigma}{2}}\int_{\lambda-\beta\sigma}^{\min\\{x_{3}-\sigma,2\lambda+\beta\sigma\\}}$
$\displaystyle\dots\int_{\frac{\sigma}{2}}^{\min\\{x_{2}-\sigma,\lambda+\beta\sigma\\}}\mathrm{d}x_{1}\mathrm{d}x_{2}\dots\mathrm{d}x_{N},$
(49)

where

$\beta:=\alpha-\frac{1}{2}.$ (50)

For a given configuration where $x_{2},\dots,x_{N}$ are frozen, we compute the
first of the nested integrals in (D) as

$\displaystyle f_{2}(x_{2})$
$\displaystyle:=\int_{\frac{\sigma}{2}}^{\min\\{x_{2}-\sigma,\lambda+\beta\sigma\\}}\mathrm{d}x$
$\displaystyle=\begin{cases}x_{2}-\frac{3}{2}\sigma&\lambda-\beta\sigma\leqq
x_{2}\leqq\lambda+(\beta+1)\sigma,\\\
\lambda+\left(\beta-\frac{1}{2}\right)\sigma&\lambda+(\beta+1)\sigma\leqq
x_{2}\leqq 2\lambda+\beta.\end{cases}$ (51)

It turns out that, provided that

$\lambda\geqq 2(\beta+1)\sigma,$ (52)

also $f_{n}$, the $(n-1)$-th of the nested integrals in (D), is a piece-wise
linear function of $x_{n}$,

$\displaystyle f_{n}(x)=$
$\displaystyle\begin{cases}A_{n}x+B_{n}&(n-1)\lambda-\beta\sigma\leqq
x\leqq(n-1)\lambda+(\beta+1)\sigma,\\\ C_{n}&(n-1)\lambda+(\beta+1)\sigma\leqq
x\leqq n\lambda+\beta\sigma,\end{cases}$ (53)

with

$C_{n}=A_{n}[(n-1)\lambda+(\beta+1)\sigma]+B_{n}.$ (54)

Moreover, the coefficients $A_{n}$ and $B_{n}$ satisfy the recurrence relation

$\begin{split}A_{n+1}&=C_{n}\\\ B_{n+1}&=(2\beta+1)\sigma
B_{n}+\frac{1}{2}(2\beta+1)\sigma[\sigma+2\lambda(n-1)]A_{n}\\\
&\quad-[(n-1)\lambda+(\beta+2)\sigma]C_{n}.\end{split}$ (55)

It follows from (D) that

$\begin{split}Z_{N}&=\int_{(N-1)\lambda-\beta\sigma}^{N\lambda-\frac{\sigma}{2}}f_{N}(x)\mathrm{d}x\\\
&=\frac{1}{2}(2\beta+1)\sigma[\sigma+2\lambda(N-1)]A_{N}+(2\beta+1)\sigma
B_{N}\\\
&\quad+\left[\lambda-\left(\beta+\frac{3}{2}\right)\sigma\right]C_{N}.\end{split}$
(56)

Combining (56), (55), and (54), we arrive at

$\begin{split}C_{N+1}&=Z_{N}+\left(\beta+\frac{1}{2}\right)\sigma C_{N},\\\
Z_{N+1}&=\left[\lambda+\left(\beta-\frac{1}{2}\right)\sigma\right]Z_{N}\\\
&\quad+\sigma\left[\left(\beta+\frac{1}{2}\right)\lambda-\left(2\beta+\beta^{2}+\frac{3}{4}\right)\sigma\right]C_{N},\end{split}$
(57)

which in matrix form reads as

$\begin{pmatrix}C_{N+1}\\\
Z_{N+1}\end{pmatrix}=\mathsf{A}\begin{pmatrix}C_{N}\\\ Z_{N}\end{pmatrix},$
(58)

where $\mathsf{A}$ is the matrix

$\mathsf{A}:=\begin{bmatrix}\left(\beta+\frac{1}{2}\right)\sigma&1\\\
\sigma\left[\left(\beta+\frac{1}{2}\right)\lambda-\left(2\beta+\beta^{2}+\frac{3}{4}\right)\sigma\right]&\lambda+\left(\beta-\frac{1}{2}\right)\sigma\end{bmatrix},$
(59)

which possesses two real positive eigenvalues,

$\displaystyle\mu_{1,2}=\frac{1}{2}\left[\lambda+(2\alpha-1)\sigma\pm\sqrt{[\lambda+(2\alpha-1)\sigma]^{2}-8\alpha^{2}\sigma^{2}}\right],$
(60)

where use has been made of (50). Letting

$\begin{pmatrix}v_{11}\\\
v_{12}\end{pmatrix}\quad\text{and}\quad\begin{pmatrix}v_{21}\\\
v_{22}\end{pmatrix}$ (61)

denote the corresponding eigenvectors, the solution of the recurrence relation
(58) with initial value

$\begin{pmatrix}C_{0}\\\
Z_{0}\end{pmatrix}=\frac{1}{2}\begin{pmatrix}\frac{1}{\sigma}\\\
1\end{pmatrix}$ (62)

is given by

$\begin{split}C_{N}&=c_{1}v_{11}\mu_{1}^{N}+c_{2}v_{21}\mu_{2}^{N},\\\
Z_{N}&=c_{1}v_{12}\mu_{1}^{N}+c_{2}v_{22}\mu_{2}^{N},\end{split}$ (63)

where $c_{1}$ and $c_{2}$ are solutions of the linear system

$\begin{split}c_{1}v_{11}+c_{2}v_{21}=\frac{1}{2\sigma},\\\
c_{1}v_{12}+c_{2}v_{22}=\frac{1}{2},\end{split}$ (64)

with the normalization $v_{11}=v_{21}=1$.

Since

$S_{N}=k_{B}T\ln Z_{N},$ (65)

in the limit as $N\gg 1$, $S_{N}$ is uniquely determined by the larger
eigenvalue $\mu_{2}$ of $\mathsf{A}$,  $\displaystyle S_{N}=k_{B}TN\ln\mu_{2}$
$\displaystyle=k_{B}TN\ln\frac{1}{2}\left[(\lambda+(2\alpha-1)\sigma)+\sqrt{[\lambda+(2\alpha-1)\sigma]^{2}-8\alpha^{2}\sigma^{2}}\right],$
(66) where use has also been made of (47). By (52), this formula is only valid
for

$0\leqq\eta\leqq\frac{1}{2\alpha+1}.$ (67)

For larger values of $\eta$, (D) is no longer valid, as $f_{n}$ is delivered
by a piece-wise quadratic function, which first becomes cubic, and then of
ever increasing order as the upper bound of $\eta$ approaches $1$. In such
cases the above computations become increasingly cumbersome.

## References

  * Torquato and Stillinger (2010) S. Torquato and F. H. Stillinger, “Jammed hard-particle packings: From Kepler to Bernal and beyond,” Reviews of Modern Physics 82, 2633 (2010).
  * Löwen (2000) H. Löwen, “Fun with hard spheres,” in  _Statistical Physics and Spatial Statistics_ (Springer, 2000) pp. 295–331.
  * Lee and Frenkel (2017) A. A. Lee and D. Frenkel, “Entropy, probability and packing,” in  _Mathematics and Materials_ , Vol. 23, edited by M. J. Bowick, D. Kinderlehrer, G. Menon, and C. Radin (American Mathematical Soc., 2017) pp. 103–165.
  * Santos, Yuste, and López de Haro (2020) A. Santos, S. B. Yuste, and M. López de Haro, “Structural and thermodynamic properties of hard-sphere fluids,” The Journal of Chemical Physics 153, 120901 (2020).
  * Santos (2016) A. Santos, “A concise course on the theory of classical liquids,” Lecture Notes in Physics 923 (2016).
  * Hoover, Hoover, and Hanson (1979) W. G. Hoover, N. E. Hoover, and K. Hanson, “Exact hard-disk free volumes,” The Journal of Chemical Physics 70, 1837–1844 (1979).
  * Hoover and Ree (1967) W. G. Hoover and F. H. Ree, “Use of computer experiments to locate the melting transition and calculate the entropy in the solid phase,” The Journal of Chemical Physics 47, 4873–4878 (1967).
  * Vieillard-Baron (1972) J. Vieillard-Baron, “Phase transitions of the classical hard-ellipse system,” The Journal of Chemical Physics 56, 4729–4744 (1972).
  * Chaikin  _et al._ (2006) P. Chaikin, A. Donev, W. Man, F. H. Stillinger, and S. Torquato, “Some observations on the random packing of hard ellipsoids,” Industrial & Engineering Chemistry Research 45, 6960–6965 (2006).
  * Ni  _et al._ (2012) R. Ni, A. P. Gantapara, J. de Graaf, R. van Roij, and M. Dijkstra, “Phase diagram of colloidal hard superballs: from cubes via spheres to octahedra,” Soft Matter 8, 8826–8834 (2012).
  * Carlsson  _et al._ (2012) G. Carlsson, J. Gorham, M. Kahle, and J. Mason, “Computational topology for configuration spaces of hard disks,” Physical Review E 85, 011303 (2012).
  * Sastry  _et al._ (1997) S. Sastry, D. S. Corti, P. G. Debenedetti, and F. H. Stillinger, “Statistical geometry of particle packings. I. algorithm for exact determination of connectivity, volume, and surface areas of void space in monodisperse and polydisperse sphere packings,” Physical Review E 56, 5524 (1997).
  * Sastry  _et al._ (1998) S. Sastry, T. M. Truskett, P. G. Debenedetti, S. Torquato, and F. H. Stillinger, “Free volume in the hard sphere liquid,” Molecular Physics 95, 289–297 (1998).
  * Kapfer  _et al._ (2012) S. C. Kapfer, W. Mickel, K. Mecke, and G. E. Schröder-Turk, “Jammed spheres: Minkowski tensors reveal onset of local crystallinity,” Physical Review E 85, 030301 (2012).
  * Maiti and Sastry (2014) M. Maiti and S. Sastry, “Free volume distribution of nearly jammed hard sphere packings,” The Journal of Chemical Physics 141, 044510 (2014).
  * Chen and Holmes-Cerfon (2017) E. R. Chen and M. Holmes-Cerfon, “Random sequential adsorption of discs on surfaces of constant curvature: Plane, sphere, hyperboloid, and projective plane,” Journal of Nonlinear Science 27, 1743–1787 (2017).
  * Noé  _et al._ (2019) F. Noé, S. Olsson, J. Köhler, and H. Wu, “Boltzmann generators: Sampling equilibrium states of many-body systems with deep learning,” Science 365, eaaw1147 (2019).
  * Eyring and Hirschfelder (1937) H. Eyring and J. Hirschfelder, “The theory of the liquid state,” Journal of Physical Chemistry 41, 249–257 (1937).
  * Hirschfelder, Stevenson, and Eyring (1937) J. Hirschfelder, D. Stevenson, and H. Eyring, “A theory of liquid structure,” The Journal of Chemical Physics 5, 896–912 (1937).
  * Lennard-Jones and Devonshire (1937) J. E. Lennard-Jones and A. Devonshire, “Critical phenomena in gases-I,” Proceedings of the Royal Society of London. Series A-Mathematical and Physical Sciences 163, 53–70 (1937).
  * Lennard-Jones and Devonshire (1938) J. E. Lennard-Jones and A. Devonshire, “Critical phenomena in gases. II. vapour pressures and boiling points,” Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 165, 1–11 (1938).
  * Buehler  _et al._ (1951) R. Buehler, R. Wentorf Jr, J. Hirschfelder, and C. Curtiss, “The free volume for rigid sphere molecules,” The Journal of Chemical Physics 19, 61–71 (1951).
  * Woodcock (1997) L. Woodcock, “Computation of the free energy for alternative crystal structures of hard spheres,” Faraday Discussions 106, 325–338 (1997).
  * Warshavsky, Ford, and Monson (2018) V. B. Warshavsky, D. M. Ford, and P. A. Monson, “On the mechanical stability of the body-centered cubic phase and the emergence of a metastable ci16 phase in classical hard sphere solids,” The Journal of Chemical Physics 148, 024502 (2018).
  * Rice (1944) O. Rice, “On the statistical mechanics of liquids, and the gas of hard elastic spheres,” The Journal of Chemical Physics 12, 1–18 (1944).
  * Kirkwood (1950) J. G. Kirkwood, “Critique of the free volume theory of the liquid state,” The Journal of Chemical Physics 18, 380–382 (1950).
  * Wood (1952) W. W. Wood, “Note on the free volume equation of state for hard spheres,” The Journal of Chemical Physics 20, 1334–1334 (1952).
  * Speedy and Reiss (1991) R. J. Speedy and H. Reiss, “Cavities in the hard sphere fluid and crystal and the equation of state,” Molecular Physics 72, 999–1014 (1991).
  * Momma and Izumi (2011) K. Momma and F. Izumi, “Vesta 3 for three-dimensional visualization of crystal, volumetric and morphology data,” Journal of Applied Crystallography 44, 1272–1276 (2011).
  * Gibson and Scheraga (1987a) K. D. Gibson and H. A. Scheraga, “Exact calculation of the volume and surface area of fused hard-sphere molecules with unequal atomic radii,” Molecular Physics 62, 1247–1265 (1987a).
  * Gibson and Scheraga (1987b) K. Gibson and H. A. Scheraga, “Volume of the intersection of three spheres of unequal size: a simplified formula,” Journal of Physical Chemistry 91, 4121–4122 (1987b).
  * Hoover and Ree (1968) W. G. Hoover and F. H. Ree, “Melting transition and communal entropy for hard spheres,” The Journal of Chemical Physics 49, 3609–3617 (1968).
  * Binder  _et al._ (2012) K. Binder, B. J. Block, P. Virnau, and A. Tröster, “Beyond the van der Waals loop: What can be learned from simulating Lennard-Jones fluids inside the region of phase coexistence,” American Journal of Physics 80, 1099–1109 (2012).
  * Ponce and Renon (1976) L. Ponce and H. Renon, “Analytical equation for the Helmholtz free energy of a pure fluid, using the perturbation theory and a square well potential,” The Journal of Chemical Physics 64, 638–640 (1976).
  * Mirzaeinia, Feyzi, and Hashemianzadeh (2017) A. Mirzaeinia, F. Feyzi, and S. M. Hashemianzadeh, “Equation of state and Helmholtz free energy for the atomic system of the repulsive Lennard-Jones particles,” The Journal of Chemical Physics 147, 214503 (2017).
  * Alder, Hoover, and Young (1968) B. Alder, W. Hoover, and D. Young, “Studies in molecular dynamics. V. high-density equation of state and entropy for hard disks and spheres,” The Journal of Chemical Physics 49, 3688–3696 (1968).
  * Colot, Baus, and Xu (1986) J.-L. Colot, M. Baus, and H. Xu, “The freezing of hard spheres: III. testing the approximations,” Molecular Physics 57, 809–823 (1986).
  * Rice (1938) O. Rice, “On communal entropy and the theory of fusion,” The Journal of Chemical Physics 6, 476–479 (1938).
  * Hoover and Alder (1966) W. G. Hoover and B. J. Alder, “Cell theories for hard particles,” The Journal of Chemical Physics 45, 2361–2367 (1966).
  * Tonks (1936) L. Tonks, “The complete equation of state of one, two and three-dimensional gases of hard elastic spheres,” Physical Review 50, 955 (1936).
  * David  _et al._ (1992) W. David, R. Ibberson, T. Dennis, J. Hare, and K. Prassides, “Structural phase transitions in the fullerene C60,” EPL (Europhysics Letters) 18, 219 (1992).
  * Wochner  _et al._ (1997) P. Wochner, X. Xiong, P. Chow, and S. Moss, “X-ray study of two-phase coexistence in a C60 single crystal through the cubic-cubic phase transition at 260 K,” Physical Review B 55, 5678 (1997).
  * Salsburg and Wood (1962) Z. Salsburg and W. Wood, “Equation of state of classical hard spheres at high density,” The Journal of Chemical Physics 37, 798–804 (1962).
  * Gujrati (2021) P. D. Gujrati, “The glass transition and the entropy crisis,” in  _Encyclopedia of Glass Science, Technology, History, and Culture_ , edited by P. Richet (Wiley Hoboken, NJ, 2021).
  * Barroso and Ferreira (2002) M. Barroso and A. Ferreira, “Solid–fluid coexistence of the Lennard-Jones system from absolute free energy calculations,” The Journal of Chemical Physics 116, 7145–7150 (2002).
  * Schilling and Schmid (2009) T. Schilling and F. Schmid, “Computing absolute free energies of disordered structures by molecular simulation,” The Journal of Chemical Physics 131, 231102 (2009).
  * Polson  _et al._ (2000) J. M. Polson, E. Trizac, S. Pronk, and D. Frenkel, “Finite-size corrections to the free energies of crystalline solids,” The Journal of Chemical Physics 112, 5339–5342 (2000).