# Cavity volume and free energy in many-body systems

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

###### Abstract

Within this work we derive and analyse an expression for the free energy of a
single-species system in the thermodynamic limit in terms of a generalised
cavity volume, that is exact in general, and in principle applicable to
systems across their entire range of density, as well as to particles within a
general coordinate space. This provides a universal equation of state, and can
thus relate the cavity volume to classical results, such as Mayer’s cluster
expansions. Through this we are able to provide some insight into the
connections between cavity volume and free energy density, as well as their
consequences. We use examples which permit explicit computations to further
probe these results, reclaiming the exact results for a classical Tonks gas
and providing a novel derivation of Onsager’s free energy for a single
species, isotropic system. Given the complexity of the problem we also provide
a local lattice ansatz, exact in one dimension, with which we may approximate
the cavity volume for hard sphere systems to provide an accurate equation of
state in the cases of hard disks and spheres in both dilute regimes as well as
beyond the freezing transition.

##  1 Introduction

The hard particle system is an attractive toy model of complex systems, as the
analysis of the complex energy landscape reduces to purely geometric
considerations. Despite the apparent simplicity of such systems, they are
capable of demonstrating a rich variety of thermodynamic behaviour such as
phase transitions and phase separation [5, 19, 29, 36]. Results that are both
exact and explicit are however generally unavailable, with scarce exceptions
such as the exactly solvable one-dimensional Tonks gas [42]. To this end, the
predominant tools for analysing hard particle systems are simplified
theoretical models and either Monte Carlo or molecular dynamics simulations
[2, 3, 4]. One particularly successful approach is the virial expansion about
the vacuum state, which presumes that the equation of state may be written as
a power series in the number density, and thus studying the thermodynamic
behaviour reduces to finding, either exactly or numerically, the values of the
Taylor coefficients [8, 12, 43]. Owing to the fact that the virial expansion
is a Taylor series about the density zero state, while it may perform
incredibly well at low density, at higher densities, or beyond a phase
transition, its ability to accurately describe the equation of state
deteriorates, and other techniques must be used to describe the system.

One such method to treat the dense regime is the free volume theory, which
considers the amount of space accessible via continuous movements to an
existing particle in the system, [10, 21, 39]. The properties of free volume
can be related to thermodynamic variables, and have been investigated
numerically, typically via the analysis of equilibrium Monte Carlo or
molecular dynamics simulations [7, 22, 37]. Calculating averages of free
volumes is a many-body problem, making explicit closed-form solutions
generally unobtainable. One way of vastly simplifying these calculations is
the cell theory, where the system is presumed to be well approximated by a
lattice, whose lattice parameters can be derived as a function of the number
density, with the consequence that complex multi-particle interactions may be
reduced to investigating the local environment of a single particle within the
lattice [13, 14, 20, 27, 28]. The cell theory is sufficiently simple to
provide closed-form results and generally performs well in the dense regime,
as dense hard particle systems are lattice-like in two and three dimensions.
As this structural assumption is less accurate in dilute systems, however, at
lower densities the predictive capability of the model falters.

Within this work we will focus on the cavity volume, which is defined to be
the amount of weighted phase space accessible to a new particle within the
entire system. In the case of hard particle systems, this simply reduces to
the accessible volume, whose complement is the excluded volume, and all
quantities of interest become purely geometric in nature. While we shall
predominantly focus on hard particle systems, the general framework can
equally be applied to soft, long-ranged interactions with sufficient decay at
large particle separation. In this case, the weighting of phase space used to
calculate the cavity volume corresponds to a Boltzmann factor-type weighting
according to the interaction energy. Our goal is to understand how the cavity
volume relates to the free energy and equation of state. Geometric properties
of the cavity volume for hard spheres in $d$ dimensions have previously been
related to thermodynamic variables by Speedy and Reiss, who provided the
equation of state

$\frac{P}{\rho kT}=1+\frac{\sigma}{2d}\frac{\langle s\rangle}{\langle
v\rangle},$ (1)

where $\sigma$ is the particle diameter, and $\langle s\rangle,\langle
v\rangle$ are the average surface area and volume of a cavity, respectively
[40]. Furthermore, geometric properties of the cavity volume have been
experimentally measured in [6, 11, 41].

Within this work we consider two aspects of the problem. In Section 2 we
present a general study on the cavity volume. By performing thermodynamic
integration with respect to the number density, we may relate the cavity
volume to the free energy density, equation of state, and virial expansion. We
have in mind the case of purely steric interactions between spheres, which
will also form the bulk of our later analysis of concrete systems; however,
the results of this section are applicable to more general systems, provided
the integral defined in (7) exists. This includes cases in which the pair
potential depends on internal molecular degrees of freedom such as orientation
or conformation, or “soft” interactions.

In Section 3 we apply these relationships to several model systems, where for
definiteness we mostly limit ourselves to familiar cases of hard spheres in
one, two and three dimensions. Within this framework we reclaim exact known
results for the 1D Tonks gas (Section 3.0.1), a Flory-Huggins-type entropy of
mixing which can alternatively be viewed as a linear constitutive assumption
on the cavity volume (Section 3.0.2), and a new derivation of the Onsager free
energy from an un-correlated, rather than dilute, assumption (Section 3.0.3).
However, in general we will not be so fortunate as to obtain exact results,
and thus we are obliged to invoke simplifying assumptions. In light of this,
we propose a fluctuating lattice model, based on an approximation by a local
lattice structure that may vary across the domain, which we believe
qualitatively captures the significant features of the cavity volume and thus
the equation of state, and is exact in one dimension (Section 3.1). Using the
fluctuating lattice model, we are able to approximate the cavity volume via a
one-dimensional integral that can be computed numerically. This leads to an
equation of state which is accurate both in dilute regimes and beyond the
freezing transition in the exemplary systems of hard spheres and disks. In
particular, this offers advantages over the limitations of the virial
expansion and cell theory, which are only valid in dilute and dense regimes,
respectively. Comparisons of the obtained equations of state with other models
can be found in Figure 3 and Figure 6 for the two-dimensional and three-
dimensional equation of state, respectively. This is also compared to Monte
Carlo estimation of the cavity volume, as outlined in Appendix A. However, in
both our theoretical and numerical treatment of the problem we observe a rapid
decay of the cavity volume at even moderate densities. This is problematic at
the level of evaluating the equation of state, as the logarithm of the cavity
volume is the key quantity of interest, meaning that small (absolute) errors
in the Monte Carlo scheme can produce catastrophically large errors in the
equation of state when the cavity volume is small. Consequently we only
consider moderate densities for Monte Carlo experiments.

##  2 The cavity volume

###  2.1 The free energy density and equation of state

Let $F_{N,V}$ denote the configurational Helmholtz free energy of a system of
$N$ particles with a state space $\Gamma$ of volume $V$, and interacting via a
pair potential $U:\Gamma^{2}\to[0,\infty]$. $\Gamma$ may be simply the
positional coordinates of particles, or it may contain internal degrees of
freedom such as orientation. Explicitly, we may write the configurational
partition function of $N$ particles in the volume $V$, denoted
$\mathcal{Z}_{N,V}$, as

$\mathcal{Z}_{N,V}=\frac{1}{N!}\int_{\Gamma^{N}}e^{-\frac{1}{kT}\sum\limits_{1\leq
i<j\leq N}U(q_{i},q_{j})}\,d^{N}\bar{q}_{N},$ (2)

where $\bar{q}_{N}=(q_{1},q_{2},...,q_{N})\in\Gamma^{N}$. We then have the
configurational Helmholtz free energy given as

$F_{N,V}=-kT\ln\mathcal{Z}_{N,V}.$ (3)

with $k$ the Boltzmann constant and $T>0$ the absolute temperature. It is
known, and a cornerstone of the Widom insertion method [44] (see also [24]),
that

$\begin{split}&\mathcal{Z}_{N,V}\\\
&=\frac{1}{N}\mathcal{Z}_{N-1,V}\left\langle\int_{\Gamma}e^{-\frac{1}{kT}\sum\limits_{i=1}^{N-1}U(q,q_{i})}\,dq\right\rangle_{P_{N-1,V}},\end{split}$
(4)

where the average is taken over particles $\bar{q}_{N-1}$ according to the
Gibbs distribution $P_{N-1,V}$. In this work, we express Widom’s relationship
using a different but equivalent set of variables. For a system of $N-1$
particles in a volume $V$ corresponding to the interaction $U$, $P_{N-1,V}$ is
given as

$\begin{split}&P_{N-1,V}(\bar{q}_{N-1})\\\
=&\frac{1}{\mathcal{Z}_{N-1,V}}e^{-\frac{1}{kT}\sum\limits_{1\leq i<j\leq
N-1}U(q_{i},q_{j})}.\end{split}$ (5)

Following (4), we have the recurrence relation

$\begin{split}&F_{N,V}-F_{N-1,V}\\\ &=kT\ln
N-kT\ln\left\langle\int_{\Gamma}e^{-\frac{1}{kT}\sum\limits_{i=1}^{N-1}U(q,q_{i})}\,dq\right\rangle_{P_{N-1,V}}.\end{split}$
(6)

For brevity, we will denote

$V_{f}(N,V)=\left\langle\frac{1}{V}\int_{\Gamma}e^{-\frac{1}{kT}\sum\limits_{i=1}^{N-1}U(q,q_{i})}\,dq\right\rangle_{P_{N-1,V}},$
(7)

which we interpret to be the generalised cavity volume fraction, with a
possible weighting corresponding to soft interactions. In the case of purely
steric interactions, where $U(q_{i},q_{j})=+\infty$ if particles with
coordinates $q_{i},q_{j}$ intersect, and $0$ otherwise, this is precisely the
volume fraction available to a new particle, randomly sampled according to the
integration measure on $\Gamma$.

This allows us to calculate the total free energy density
$\mathcal{F}_{N,V}=\frac{F_{N,V}}{V}$ by an incremental scheme, as

$\begin{split}&\mathcal{F}_{N,V}\\\
=&\frac{1}{V}F_{0,V}+\sum\limits_{n=1}^{N}\frac{1}{V}(F_{n,V}-F_{n-1,V})\\\
=&\mathcal{F}_{0,V}+kT\sum\limits_{n=1}^{N}\frac{1}{V}\left(\ln n-\ln
VV_{f}(n-1,V)\right).\end{split}$ (8)

We presume that $V_{f}$ is well defined in the thermodynamic limit, that is,
there exists some function $\mathcal{V}_{f}$ so that for $N,V$ large,
$V_{f}(N,V)\approx\mathcal{V}_{f}\left(\frac{N}{V}\right)$. If
$\ln\mathcal{V}_{f}$ is Riemann integrable, we may approximate the sum (8) as
a Riemann sum with step size $\Delta x=\frac{1}{V}$, $x_{n}=n\Delta x$.
Furthermore taking the vacuum energy $F_{0,V}=0$, without loss of generality,
and replacing $V_{f}$ with $\mathcal{V}_{f}$, this reduces to

$\begin{split}\frac{\mathcal{F}_{N,V}}{kT}\approx&\sum\limits_{n=1}^{N}\frac{1}{V}\left(\ln\frac{n}{V}-\ln\mathcal{V}_{f}\left(\frac{n}{V}\right)\right)\\\
\approx&\int_{0}^{\frac{N}{V}}\ln(x)-\ln\mathcal{V}_{f}(x)\,dx\\\
=&\rho\ln\rho-\rho-\int_{0}^{\rho}\ln\mathcal{V}_{f}(x)\,dx,\end{split}$ (9)

where $\frac{N}{V}=\rho$ denotes the number density. Following this, we define
the free energy density $\mathcal{F}$ at number density $\rho$ as

$\mathcal{F}(\rho):=kT\left(\rho\ln\rho-\rho-\int_{0}^{\rho}\ln\mathcal{V}_{f}(x)\,dx\right).$
(10)

We expect the function $\mathcal{V}_{f}$ to satisfy $\mathcal{V}_{f}(0)=1$ and
$\mathcal{V}_{f}(\rho^{*})=0$ where $\rho^{*}$ is the highest possible number
densit. We do not rule out the possibility of non-monotonicty of
$\mathcal{V}_{f}$ between these regimes. In the case of purely steric
interactions, where $U$ attains only the values zero (no penetration of
particles) or $\infty$ (penetration of particles), $\mathcal{V}_{f}$ only
depends on the number density and particle shape. However in the case of soft
interactions, in which $U$ attains values other than $0$ and $+\infty$,
$\mathcal{V}_{f}$ will also be dependent on the temperature, which is
immediate from (7).

Equation (10) shows the significance of the cavity volume fraction in
calculating the energy, and we can provide a universal equation of state for
such a system as

$\begin{split}\frac{P}{kT}=&\frac{1}{kT}\left(-\mathcal{F}(\rho)+\rho\frac{\partial\mathcal{F}}{\partial\rho}(\rho)\right)\\\
=&\rho-\rho\ln\mathcal{V}_{f}(\rho)+\int_{0}^{\rho}\ln\mathcal{V}_{f}(x)\,dx.\end{split}$
(11)

By integration by parts if $\mathcal{V}_{f}$ is differentiable, we may also
write that

$\frac{P}{\rho
kT}=1-\frac{1}{\rho}\int_{0}^{\rho}\frac{x\mathcal{V}_{f}^{\prime}(x)}{\mathcal{V}_{f}(x)}\,dx.$
(12)

This formula is particularly interesting as it has a superficial similarity to
the Speedy and Reiss formula (1) for hard spheres. In dilute systems when the
excluded volume is made entirely of disjoint exclusion spheres,
$\mathcal{V}_{f}(\rho)=1-\rho\sigma^{d}\omega_{d}$, where $\omega_{d}$ is the
volume of a ball of radius $1$ in $\mathbb{R}^{d}$. Then the surface area of
this excluded volume is $d\rho\sigma^{d-1}\omega_{d}$, so that
$\frac{\sigma}{2d}\langle
s\rangle=\frac{1}{2}\rho\sigma^{d}\omega_{d}=-\frac{\rho}{2}\mathcal{V}_{f}^{\prime}(\rho)$.
Roughly speaking, this suggests that (12) may be a kind of integral
formulation of the Speedy formula.

It is possible to invert the equation (12) to provide the cavity volume as a
function of pressure, if the latter is known. It is immediate that by
rearranging and taking a derivative, we have the relationship

$1-\frac{\partial}{\partial\rho}\frac{P(\rho)}{kT}=\frac{\rho\mathcal{V}_{f}^{\prime}(\rho)}{\mathcal{V}_{f}(\rho)}=\rho\frac{\partial}{\partial\rho}\ln\mathcal{V}_{f}(\rho).$
(13)

By further rearrangement and integration, we have that

$\mathcal{V}_{f}(\rho)=\exp\left(\int_{0}^{\rho}\frac{1}{x}\left(1-\frac{1}{kT}\frac{\partial
P}{\partial x}\right)\,dx\right).$ (14)

The cavity volume fraction may also be expressed in terms of the
compressibility factor,

$Z:=\frac{P}{\rho kT},$ (15)

which is typically defined in terms of the packing fraction $\eta$, defined by

$\eta:=\rho v_{0},$ (16)

with $v_{0}$ the single particle volume, giving

$\mathcal{V}_{f}(\eta)=\exp\left(\int_{0}^{\eta}\frac{1}{x}\left(1-\frac{\partial}{\partial
x}(xZ(x))\right)\,dx\right).$ (17)

For example, by inserting into (17) the equation of state for Carnahan-
Starling (see, e.g., [35, Page 76]), we obtain

$\mathcal{V}_{f}(\eta)=\exp\left(\frac{\eta(3\eta^{2}-9\eta+8)}{(1-\eta)^{3}}\right)$
(18)

Furthermore, we can infer from (10) that if $\mathcal{V}_{f}$ is $k$-times
differentiable and non-zero, $\mathcal{F}$ is consequently $k+1$-times
differentiable. More so, if $\mathcal{V}_{f}$ is continuous and monotonically
decreasing, then
$\frac{\partial\mathcal{F}}{\partial\rho}=kT\ln\rho-1-kT\ln\mathcal{V}_{f}(\rho)$
is necessarily continuous and monotonically increasing, implying that
$\mathcal{F}$ is convex, and thus there can be no phase separation.

This approach is essentially a form of thermodynamic integration [17, Section
7.1], where the free energy difference between two states is calculated by
integrating a tractable derivative of the energy between the states. In our
case, we are calculating the energy difference between the vacuum and a given
state by integrating the chemical potential as a function of number density in
the thermodynamic limit. A standard application of this technique is to
consider two systems with different interactions, parametrise a path between
the two types of interaction, and integrate the derivative of the free energy
across this path to evaluate the energy difference of the two systems.
Typically thermodynamic integration is performed numerically, with the
derivative of the free energy estimated via simulation methods. The advantage
of such a scheme is that the quadrature scheme implicitly averages the errors
of constituent simulations when evaluating the energy difference, which
consequently makes the final results typically more accurate than the
individual simulations used to obtain them.

###  2.2 Virial expansion and cavity volume

The virial equation of state for a system is

$\frac{1}{kT}P=\sum\limits_{n=1}^{\infty}\textup{B}_{n}\rho^{n},$ (19)

where the constants $\textup{B}_{n}$ are known as the virial coefficients.
There is a rich history of attempts to compute and approximate the virial
coefficents [12, 25]. In this discussion we will relate the virial
coefficients to the cavity volume in an exact and general way, expressing the
virial coefficients in terms of an analogous series expansion for the cavity
volume fraction.

The coefficients $\textup{B}_{n}$ have the same units of $\rho^{1-n}$ and
$v_{0}^{n-1}$, thus we introduce the reduced virial coefficients
$\text{b}_{n}$, defined as $\text{b}_{n}=v_{0}^{1-n}\textup{B}_{n}$. The
reduced virial coefficients $\text{b}_{n}$ are thus dimensionless constants
that satisfy

$\frac{P}{\rho kT}=Z=\sum\limits_{n=0}^{\infty}\text{b}_{n+1}\eta^{n},$ (20)

with $\eta=\rho v_{0}$ the packing fraction as before.

For the following it is more convenient to work with the excluded volume than
the cavity volume, with the excluded volume fraction $\bar{\mathcal{V}}_{f}$
defined as

$1-\bar{\mathcal{V}}_{f}(\rho)=\mathcal{V}_{f}(\rho)$ (21)

Assume that the excluded volume is analytic at the dilute regime, so that it
may be written as

$\bar{\mathcal{V}}_{f}(\rho)=\sum\limits_{n=1}^{\infty}\frac{v_{n}\rho^{n}}{n!}$
(22)

for $\rho$ sufficiently close to $0$ and coefficients $v_{n}$. We note that as
$\mathcal{V}_{f}(0)=1$, $\bar{\mathcal{V}}_{f}(0)=0$, hence the sum is over
$n\geq 1$.

In the following we will use $B_{n,k}$ to denote the partial exponential Bell
polynomials. For an introduction and discussion on Bell polynomials the reader
is directed to [9, Section 3.3]. These may be defined in numerous ways, one of
which is via a generating function type equality, so that they satisfy the
relationship

$\frac{1}{k!}\left(\sum\limits_{j=1}^{\infty}x_{j}\frac{t^{j}}{j!}\right)^{k}=\sum\limits_{n=k}^{\infty}B_{n,k}(x_{1},x_{2},...,x_{n-k+1})\frac{t^{n}}{n!}$
(23)

for real $t$, integer $k$ and sequences $(x_{j})_{j=1}^{\infty}$ so that the
sums are absolutely convergent. Immediately following the definition, we see
that if $y_{j}=\frac{x_{j}}{\alpha^{j}}$, then $y_{j}(\alpha
t)^{j}=x_{j}t^{j}$. This implies a homogeneity condition

$\alpha^{n}B_{n,k}(y_{1},y_{2},...,y_{n-k+1})=B_{n,k}(x_{1},x_{2},...x_{n-k+1}).$
(24)

For notational brevity, we consider the excess free energy
$\bar{\mathcal{F}}$, defined as

$\bar{\mathcal{F}}(\rho)=\mathcal{F}(\rho)-kT(\rho\ln\rho-\rho).$ (25)

This is simply the non-ideal part of the free energy. $\bar{\mathcal{F}}$ may
then be written as

$\begin{split}&\frac{1}{kT}\bar{\mathcal{F}}=\\\
&-\int_{0}^{\rho}\ln\big{(}1-\bar{\mathcal{V}}_{f}(t)\big{)}\,dt\\\
=&\int_{0}^{\rho}\sum\limits_{k=1}^{\infty}\frac{\bar{\mathcal{V}}_{f}(t)^{k}}{k}\,dt\\\
=&\int_{0}^{\rho}\sum\limits_{k=1}^{\infty}(k-1)!\frac{1}{k!}\left(\sum\limits_{n=1}^{\infty}\frac{v_{n}t^{n}}{n!}\right)^{k}\,dt\\\
=&\int_{0}^{\rho}\sum\limits_{k=1}^{\infty}(k-1)!\sum\limits_{n=k}^{\infty}B_{n,k}(v_{1},...,v_{n-k+1})\frac{t^{n}}{n!}\,dt\\\
=&\int_{0}^{\rho}\sum\limits_{n=1}^{\infty}\left(\sum\limits_{k=1}^{n}(k-1)!B_{n,k}(v_{1},...,v_{n-k+1})\right)\frac{t^{n}}{n!}\,dt\\\
=&\sum\limits_{n=1}^{\infty}\left(\sum\limits_{k=1}^{n}(k-1)!B_{n,k}(v_{1},...,v_{n-k+1})\right)\frac{\rho^{n+1}}{(n+1)!}.\end{split}$
(26)

This gives the virial equation of state as

$\begin{split}&\frac{P}{\rho kT}-1=\\\ &-\frac{1}{\rho
kT}\bar{\mathcal{F}}+\frac{1}{kT}\frac{\partial\bar{\mathcal{F}}}{\partial\rho}\\\
=&\sum\limits_{n=1}^{\infty}\left(\sum\limits_{k=1}^{n}(k\shortminus
1)!B_{n,k}(v_{1},...,v_{n\shortminus
k+1})\right)\frac{n\rho^{n}}{(n+1)!},\end{split}$ (27)

which in turn gives the virial coefficients as

$\textup{B}_{n}=\frac{(n\shortminus
1)}{n!}\left(\sum\limits_{k=1}^{n-1}(k\shortminus 1)!B_{n\shortminus
1,k}(v_{1},...,v_{n\shortminus k})\right).$ (28)

We note that this expression is valid for $n\geq 2$, as the ideal gas term
corresponding to $n=1$ is not considered.

Dividing through by $v_{0}^{n-1}$ to give the reduced virial coefficients,
using the homogeneity condition (24) gives

$\text{ b}_{n}=\frac{n\shortminus 1}{n!}\left(\sum\limits_{k=1}^{n\shortminus
1}(k\shortminus 1)!B_{n\shortminus 1,k}(u_{1},...,u_{n\shortminus k})\right).$
(29)

where $u_{j}=v_{0}^{-j}v_{j}$, so that for $\eta=\rho v_{0}$,

$\bar{\mathcal{V}}_{f}(\rho)=\sum\limits_{i=1}^{\infty}\frac{v_{i}\rho^{i}}{i!}=\sum\limits_{i=1}^{\infty}\frac{u_{i}\eta^{i}}{i!}.$
(30)

In principle, the coefficients $v_{k}$ can be computed numerically. This is
because for $V$ sufficiently large, we know that
$\mathcal{V}_{f}\left(\frac{N}{V}\right)\approx V_{f}\left(N,V\right)$, and
$V_{f}$ is a quantity defined for finite systems. Using this, we can
approximate derivatives at $0$ by finite differences with step sizes
$h=\frac{1}{V}$, giving

$\begin{split}v_{k}=&\left.\frac{d^{k}}{d\rho^{k}}\bar{\mathcal{V}}_{f}(\rho)\right|_{\rho=0}\\\
\approx&V^{k}\sum\limits_{i=1}^{k}\binom{k}{i}\bar{\mathcal{V}}_{f}\left(\frac{i}{V}\right)(-1)^{k-i}\\\
\approx&V^{k}\sum\limits_{i=1}^{k}\binom{k}{i}\left(1-V_{f}\left(i,V\right)\right)(-1)^{k-i}\\\
\end{split}$ (31)

The first of the approximations should be exact as $V\to\infty$ if
$\mathcal{V}_{f}$ is sufficiently smooth, and the second if $V_{f}(i,V)$ is
sufficiently well approximated by $\mathcal{V}_{f}\left(\frac{i}{V}\right)$.
This implies that to find the virial coefficient $\text{B}_{n}$, in principle
it suffices to have a very good estimate of $V_{f}(i,V)$ for $i=1,...,n-1$,
which amounts to estimating how much space is excluded in a system of $n-1$
particles.

In terms of the dimensionless coefficients $u_{k}$, the first few reduced
virial coefficients are thus given by

$\begin{split}{\text{b}}_{2}=&\frac{u_{1}}{2},\\\
{\text{b}}_{3}=&\frac{1}{3}\left(u_{1}^{2}+u_{2}\right),\\\
{\text{b}}_{4}=&\frac{1}{8}\left(2u_{1}^{3}+3u_{1}u_{2}+u_{3}\right),\\\
{\text{b}}_{5}=&\frac{1}{30}\left(6u_{1}^{4}+12u_{1}^{2}u_{2}+3u_{2}^{2}+4u_{1}u_{3}+u_{4}\right).\end{split}$
(32)

Using $V_{i}=V_{f}(i,V)$ for brevity, we can thus give the virial coefficients
in terms of $i$-particle excluded volumes as

$\begin{split}{\text{b}}_{2}=&\frac{V}{2v_{0}}(1-V_{1}),\\\
{\text{b}}_{3}=&\frac{V^{2}}{3v_{0}^{2}}\left(V_{1}^{2}-2V_{2}\right),\\\
{\text{b}}_{4}=&\frac{V^{3}}{8v_{0}^{3}}\left(-2V_{1}^{3}+3V_{1}V_{2}-V_{3}\right),\\\
{\text{b}}_{5}=&\frac{V^{4}}{30v_{0}^{4}}\left(6V_{1}^{4}-12V_{1}^{2}V_{2}+3V_{2}^{2}+4V_{1}V_{3}-V_{4}\right).\end{split}$
(33)

There will be three sources of error in using such representations to
numerically obtain virial coefficients, which are finite size effects in a
sampled system, sampling and integration error in evaluating $V(i,V)$, and
error in the approximations of the derivatives.

###  2.3 Comparison with the free volume theory

The theory presented within this work concerns cavity volume, which is
related, although distinct from, the well studied theory of free volume [7,
13, 14, 27, 28]. Given a configuration of hard particles, the free volume of a
particular particle is the volume that it may access by continuous motion
while holding all other particles in place, without particle interpenetration.
We emphasise that the free volume, in contrast to the cavity volume, is not a
globally defined property of the system; rather, it is a region of space
defined with respect to a probe particle within the system. Regions of free
volume with respect to different probe particles may be disjoint, or have a
non-empty intersection. In a dilute regime, where typically particles are
typically far apart, the particle can access almost any space outside of the
exclusion spheres of other particles, which are disjoint. This implies that
the free volume is an extensive quantity roughly equal to $V-Nv_{ex}^{0}$,
where $v_{ex}^{0}$ is the pairwise excluded volume of an average particle. In
dense systems however, particles are expected to be caged by their neighbours,
in which case the free volume is an intensive quantity, and the size of the
free volume should be comparable to the size of the particle itself. The
qualitative change of the free volume from an extensive to intensive quantity
is associated with the percolation transition. The advantage of the free
volume theory is that it is relatively straightforward to analyse in dense
systems, and the geometric properties of the free volume can be related to
various thermodynamic quantities. As noted by Sastry et al. [37], however, the
free volume is linked with the cavity volume in an intimate way. We replicate
the key points from their discussion here for comparison.

We define a connected cavity to be a connected component of the cavity volume,
and we will denote the volume of a connected cavity by $u$. Furthermore, given
a particular particle we denote its free volume as $v$. Then in a system of
$N+1$ particles, the free volume of the $N+1$-th particle corresponds to a
connected cavity in a system of $N$ particles, and the $N+1$-th particle will
be added to any particular connected cavity with a probability proportional to
its size, $u$. This argument can be used to show that the average free volume
and average cavity volume can be related by

$\langle v^{-1}\rangle^{-1}=\langle u\rangle.$ (34)

That is, the arithmetic mean of the connected cavity size is equal to the
harmonic mean of the free volume. If we have $N_{c}$ connected cavities in a
typical system, we may then relate the total cavity volume with the harmonic
mean of the free volume by the relationship

$V_{f}(N,V)=1-\frac{N_{c}}{V}\langle u\rangle=1-\frac{N_{c}}{V\langle
v^{-1}\rangle}.$ (35)

As $V$ is extensive and $V_{f}$ is intensive, this implies that
$\frac{N_{c}}{\langle v^{-1}\rangle}$ is also an extensive quantity, and thus
at the percolation transition, when $v$ undergoes a qualitative change from an
extensive to intensive quantity, we must have that $N_{c}$ does the reverse,
going from an intensive quantity (equal to $1$ in the dilute limit for hard
sphere systems in dimension greater than 1) to an extensive quantity. This
illustrates that while there is a link between the two theories, there is a
further quantity, $N_{c}$, which obscures the relationship between the two and
remains elusive. Furthermore, as both $N_{c}$ and $v$ undergo an
extensive/intensive exchange as the concentration increases and decreases
(respectively), any attempt to define them in the thermodynamic limit for all
concentrations would prove difficult. This particular problem is not present
in the cavity volume approach, as the cavity volume fraction is always an
intensive variable and can thus be defined without complication in the
thermodynamic limit. In particular, this allows us to use the cavity volume
approach across the entire range of number densities without the need to
impose different ad-hoc definitions in different regimes.

###  2.4 $\lambda$-function

In [31] a constitutive equation for the excluded volume fraction was proposed,
of the form

$1-\mathcal{V}_{f}(\rho)=\lambda(\rho)v_{ex}^{0}\rho,$ (36)

where $v_{ex}^{0}$ is the pairwise excluded volume given by

$v_{ex}^{0}=\int_{\mathbb{R}^{d}}1-\exp\left(-\frac{1}{kT}U(q_{1},q)\right)\,dq,$
(37)

which is implicitly assumed to be independent of $q_{1}$ on symmetry grounds.
The introduction of the function $\lambda$ was mostly ad-hoc, and aimed to act
as an effective correction for the over-counting of pairwise excluded volume.
Within this section we aim to revisit $\lambda$ and its interpretation.

Consider the case of a pair potential $U$, taken to be non-negative without
loss of generality. Recall we may write the (weighted) cavity volume fraction
in a single species system as

$V_{f}(N,V)=\frac{1}{V}\left\langle\int_{\Gamma}e^{-\frac{1}{kT}\sum\limits_{i=1}^{N-1}U_{i,N}}\,dq\right\rangle_{P_{N-1}}.$
(38)

We consider arbitrary particles with generalised coordinates $q\in\Gamma$.
Following the traditions of cluster expansions [30] define the new variable
$f_{i}$, analogous to a Mayer function, as

$f_{i}=1-\exp\left(-\frac{1}{kT}U_{i,N}\right),$ (39)

where for brevity we denote $U(q_{i},q_{N})=U_{i,N}$. Then we may re-write the
cavity volume fraction as

$\begin{split}&1-V_{f}(N,V)\\\
=&1-\frac{1}{V}\int_{\Gamma}\left\langle\prod\limits_{i=1}^{N-1}(1-f_{i})\right\rangle\,dq\\\
=&\frac{1}{V}\int_{\Gamma}\left\langle{\textstyle\sum\limits_{i_{1}}}f_{i_{1}}-\sum\limits_{i_{1},i_{2}}f_{i_{1}}f_{i_{2}}\right.\\\
&\hskip
28.45274pt\left.+{\textstyle\sum\limits_{i_{1},i_{2},i_{3}}}f_{i_{1}}f_{i_{2}}f_{i_{3}}+...\right\rangle\,dq\\\
=&{\textstyle\sum\limits_{i_{1}}}\frac{1}{V}\int_{\Gamma}\left\langle
f_{i_{1}}\left(1-{\textstyle\sum\limits_{i_{2}}}f_{i_{2}}+{\textstyle\sum\limits_{i_{2},i_{3}}}f_{i_{2}}f_{i_{3}}+...\right)\right\rangle\,dq\end{split}$
(40)

The sums are taken over $i_{j}=1,...,N-1$, and exclude cases where
$i_{j}=i_{k}$ for $j\neq k$, and averages $\langle\cdot\rangle$ are with
respect to $P_{N-1}$. As $f_{i_{1}}P_{N-1}$ is non-negative, if the series in
the brackets is bounded above and below we may find a dimensionless value
$\lambda_{i_{1}}>0$ so that

$\begin{split}&\int_{\Gamma}\left\langle
f_{i_{1}}\left(1-{\textstyle\sum\limits_{i_{2}=1}^{N-1}}f_{i_{2}}+{\textstyle\sum\limits_{i_{2},i_{3}=1}^{N-1}}f_{i_{2}}f_{i_{3}}-...\right)\right\rangle\,dq\\\
=&\lambda_{i_{1}}\int_{\Gamma}\left\langle
f_{i_{1}}\right\rangle_{P_{N-1}}\,dq\end{split}$ (41)

If the energy $U$ is continuous, we may apply the mean value theorem for
products of integrals to conclude that there exists some
$(q_{1},q_{2},...,q_{N-1})\in\Gamma^{N-1}$, $q\in\Gamma$ such that the scalar
$\lambda_{i}$ precisely satisfies

$\begin{split}\lambda_{i_{1}}&=\left(1-{\textstyle\sum\limits_{i_{2}=1}^{N-1}}f_{i_{2}}+{\textstyle\sum\limits_{i_{2},i_{3}=1}^{N-1}}f_{i_{2}}f_{i_{3}}-...\right),\\\
f_{j}=&1-\exp\left(-\frac{1}{kT}U(q_{j},q)\right).\end{split}$ (42)

Following this we may interpret $\lambda_{i_{1}}$ as a “typical” re-counting
factor for overlaps, and in the discontinuous case we may still view it as a
kind of effective recount. By symmetry, $\lambda_{i_{1}}$ may only be a
function of number density. That is, $\lambda_{i_{1}}=\lambda(\rho)$. Hence we
may write that

$\begin{split}&V_{f}(N,V)\\\
=&1-\lambda(\rho)\sum\limits_{i_{1}=1}^{N-1}\frac{1}{V}\int_{\Gamma}\langle
f_{i_{1}}\rangle_{P_{N-1}}\,dq\\\
=&1-\lambda(\rho)\frac{N-1}{V}v_{ex}^{0}.\end{split}$ (43)

This then gives the cavity volume fraction in the thermodynamic limit as

$\mathcal{V}_{f}(\rho)=1-\lambda(\rho)\rho v_{ex}^{0}.$ (44)

Of course $\lambda$ can be readily claimed from a known expression for
$\mathcal{V}_{f}(\rho)$ provided $v_{ex}^{0},\rho$ are known, but we claim its
interpretation as an effective correction for overlapping exclusion zones
makes it a worthwhile quantity to consider in itself.

In the case of hard spheres in $n$-dimensions, $v_{ex}^{0}=2^{n}v_{0}$, so the
fact that $\mathcal{V}_{f}\in[0,1]$ implies that $\lambda(\rho)\in[2^{-n},1]$,
where we expect $\lambda(0)=1$.

##  3 Applications and results

We now turn towards the application of our previous general results to
concrete systems. Our main novel contribution will be our fluctuating lattice
model, which is defined and analysed in Section 3.1. Before addressing this
however, we aim to use the following subsection to illustrate uses of the
framework developed throughout Section 2 within more familiar contexts to
demonstrate its efficacy in obtaining several classical results. Explicitly,
we may calculate exactly the cavity volume to rederive the equation of state
for the Tonks gas in one dimension (Section 3.0.1), use a linear constitutive
equation on the cavity volume to obtain the Flory-Huggins entropy of mixing
(Section 3.0.2), and take an ansantz of an uncorrelated system to obtain the
virial expansion to second order (Section 3.0.3).

####  3.0.1 Tonks Gas

A one-dimensional system of hard rods on a line is an exactly solvable system
[42]. Here we verify that classical results may be re-obtained via calculation
of the cavity volume and (10). More so, given that many exact computations are
possible for this one-dimensional system, we can use this as a test case to
understand in more depth the nature of cavity volume. Speedy and Reiss [40]
performed similar calculations for evaluating a distinct quantity, the free
volume. As the calculations we perform here are near-equivalent, we direct the
reader to their work for the details and provide only the key steps. We first
note a system of $N$ hard rods of length $\sigma$ on a periodic line (or
equivalently, a circle) of length $V$ is equivalent to $N$ non-interacting
points $(x_{i})_{i=1}^{N}$ on a periodic line of length $V-N\sigma$. Given
some probe point $x$, the periodic distances from $x$ are then distributed
uniformly on $\frac{V-N\sigma}{2}$. Thus the probability that all particles
are at a distance at least $r>0$ from $x$ is given by

$\left(\frac{V-N\sigma-2r}{V-N\sigma}\right)^{N}.$ (45)

In particular, taking $r=\frac{\sigma}{2}$, we see that the probability of
being able to insert a new particle of length $\sigma$ on the line of length
$V-N\sigma$ is

$\left(\frac{V-N\sigma-\sigma}{V-N\sigma}\right)^{N}=\left(1-\frac{\sigma}{V-N\sigma}\right)^{N}.$
(46)

This then gives the cavity volume fraction on a line of length $V$ as

$V_{f}(N,V)=\left(1-\frac{\sigma}{V-N\sigma}\right)^{N}\left(\frac{V-N\sigma}{V}\right),$
(47)

where the latter bracket corresponds to the probability of a probe point not
lying within an existing particle. Introducing the number density
$\rho=\frac{N}{V}$, we see this has a readily computable limit as
$N,V\to\infty$ with $\frac{N}{V}=\rho$ given by

$\begin{split}V_{f}(N,V)=&\left(1-\frac{\sigma}{V-N\sigma}\right)^{N}\left(\frac{V-N\sigma}{V}\right)\\\
=&\left(1-\frac{\rho\sigma}{N(1-\rho\sigma)}\right)^{N}\left(1-\rho\sigma\right)\\\
\to&\exp\left(-\frac{\rho\sigma}{1-\rho\sigma}\right)\left(1-\rho\sigma\right).\end{split}$
(48)

This expression decays to zero faster than any polynomial as
$\rho\to\frac{1}{\sigma}$, in fact it is a classical example of a function
where every (left-)derivative vanishes at $\rho=\frac{1}{\sigma}$, but the
function is non-zero for $\rho<\frac{1}{\sigma}$. In particular, a truncated
series expansion of the cavity volume at $\rho=\frac{1}{\sigma}$ would only
produce $\mathcal{V}_{f}(\rho)=0$, giving no information.

We may then calculate

$\begin{split}&\int_{0}^{\rho}\ln\mathcal{V}_{f}(x)\,dx\\\
=&\int_{0}^{\rho}\ln\left(\exp\left(-\frac{x\sigma}{1-x\sigma}\right)(1-x\sigma)\right)\,dx\\\
=&\int_{0}^{\rho}\frac{-x\sigma}{1-x\sigma}+\ln(1-x\sigma)\,dx\\\
=&\rho\ln(1-\rho\sigma).\end{split}$ (49)

By (10), we obtain a logarithmic singularity in the free energy density,

$\frac{1}{kT}\mathcal{F}=\rho\ln\rho-\rho-\rho\ln(1-\rho\sigma),$ (50)

and the classical equation of state for a Tonks gas,

$P=-\mathcal{F}+\rho\frac{\partial\mathcal{F}}{\partial\rho}=\frac{kT\rho}{1-\rho\sigma}.$
(51)

By means of the free energy we can observe that there are no phase transitions
as the energy is analytic in the number density for $0<\rho<\frac{1}{\sigma}$,
and we can rule out phase separation as the energy is strictly convex, as
verified from its second derivative

$\frac{1}{kT}\frac{\partial^{2}}{\partial\rho^{2}}\mathcal{F}=\frac{1}{\rho(1-\rho\sigma)^{2}}>0.$
(52)

We can reclaim $\lambda$, as per Section 2.4, from the cavity volume as

$\begin{split}\lambda(\rho)=&\frac{1-\exp\left(-\frac{\rho\sigma}{1-\rho\sigma}\right)(1-\rho\sigma)}{2\rho\sigma},\end{split}$
(53)

and for which we have included a plot in Figure 1. We observe a monotonic,
smooth decay from $\lambda(0)=1$ to
$\lambda\left(\frac{1}{\sigma}\right)=\frac{1}{2}$.

Figure 1: $\lambda(\rho)$ versus $\rho$ for $\sigma=1$.

We can consider the total number of cavities, $N_{c}$, as defined in Section
2.3. It should be noted that $N_{c}$ is of a very different nature in one
dimension, as a single particle disconnects space. Because of this, $N_{c}$
has limiting behaviour $N_{c}\to N$ in the dilute regime, and is always
extensive. This is in contrast to the behaviour of spheres in $d>1$
dimensions, which would have $N_{c}\to 1$. Using the same argument as used to
obtain (46), we see that a cavity exists between a particle and its neighbour
if and only if they are separated by a distance of at least $2r=\sigma$, which
has probability

$\left(\frac{V-(N+1)\sigma}{V-N\sigma}\right)^{N}.$ (54)

The number of cavities per volume is then simply this quantity multiplied by
$\rho=\frac{N}{V}$, giving

$\begin{split}&\frac{1}{V}N_{c}\\\
=&\frac{N}{V}\left(\frac{V-(N+1)\sigma}{V-N\sigma}\right)^{N}\\\
\to&\rho\exp\left(-\frac{\rho\sigma}{1-\rho\sigma}\right).\end{split}$ (55)

Furthermore, we can then infer the average cavity size and harmonic mean of
the free volume as

$\langle u\rangle=\langle
v^{-1}\rangle^{-1}=\frac{V\mathcal{V}_{f}}{N_{c}}=\frac{1-\rho\sigma}{\rho}$
(56)

In the dense limit, as $\rho\nearrow\rho^{*}=\frac{1}{\sigma}$, we thus see
that the leading contribution to the decay of the cavity volume is not from
the size of individual cavities, but their rarity, as $\langle u\rangle$
converges linearly in $\rho^{*}-\rho$ to zero, while $N_{c}$ decays faster
than any polynomial.

####  3.0.2 Linear Depletion and hole/particle equivalence

Outside of some toy systems, we generally cannot expect to evaluate the cavity
volume exactly, so we will need to invoke an ansatz or constitutive equation
on the cavity volume. The simplest such ansatz we may take is that the cavity
volume is an affine function of number density. To agree with exact results in
the dilute limit, we take the constitutive equation

$\mathcal{V}_{f}(\rho)=1-\rho v_{ex}^{0},$ (57)

where $v_{ex}^{0}$ is the two-particle generalised excluded volume, as in
(37). This constitutive equation can be viewed in one of two manners; the
first being that it is a linear approximation to the excluded volume of a
general system, the second is that this is the exact behaviour of a system of
finitely many discrete “bins” that can each hold at most one particle. In
terms of the previously defined $\lambda$-function, this corresponds to the
constitutive equation $\lambda(\rho)=1$. In the case of the constitutive
equation (57), the free energy can be evaluated explicitly, as

$\begin{split}\int_{0}^{\rho}\ln(\mathcal{V}_{f}(x))\,dx=&-\rho-\frac{1-\rho
v_{ex}^{0}}{v_{ex}^{0}}\ln(1-\rho v_{ex}^{0}).\end{split}$ (58)

Introducing a rescaled number density $\xi=\rho v_{ex}^{0}$, we then have the
free energy via (10) as

$\begin{split}&\frac{1}{kT}\mathcal{F}\\\ =&\rho\ln\rho+\frac{1-\rho
v_{ex}^{0}}{v_{ex}^{0}}\ln(1-\rho v_{ex}^{0})\\\
=&\frac{1}{v_{ex}^{0}}\left(\xi\ln\xi+(1-\xi)\ln(1-\xi)\right)-\xi\frac{\ln
v_{ex}^{0}}{v_{ex}^{0}}.\end{split}$ (59)

The system is saturated as $\rho\to\frac{1}{v_{ex}^{0}}^{-}$, and the energy
is bounded in this limit, with

$\lim\limits_{\rho\to\frac{1}{v_{ex}^{0}}^{-}}\mathcal{F}=-kT\frac{\ln
v_{ex}^{0}}{v_{ex}^{0}}.$

This free energy gives the equation of state

$P=-\frac{kT}{v_{ex}^{0}}\ln(1-\rho v_{ex}^{0}),$ (60)

exhibiting logarithmic divergence of the pressure

In the case of discrete bins, $v_{ex}^{0}=1,\rho=\xi$, in which case the
energy (59) reduces to a simple Flory-Huggins-type entropy of mixing [16, 23].
In this case as the energy is symmetric under inversion $\rho\mapsto 1-\rho$,
the theory of “particles” and the theory of “holes” are equivalent.

####  3.0.3 Uncorrelated systems and the Onsager approximation

In this section the state space $\Gamma$ will be general, but we have in mind
the case of a positional coordinate, potentially with internal degrees of
freedom such as orientation or conformation. The term “volume” will refer to
the measure on $\Gamma$. We presume that particles are uncorrelated, which we
express mathematically as the approximation

${P_{N}(\bar{q})\approx\frac{1}{|\Gamma|^{N}}=\text{Un}_{N}(\bar{q})},$ (61)

where $\Gamma$ is the generalised state space and $\text{Un}_{N}$ denotes the
uniform distribution on $N$ particles in the state space $\Gamma$. Thus we may
approximate

$\begin{split}&V_{f}(N,V)\\\
&\approx\left\langle\frac{1}{|\Gamma|}\int_{\Gamma}\exp\left(-\frac{1}{kT}\sum\limits_{i=1}^{N-1}U(q,q_{i})\right)\,dq\right\rangle_{\text{Un}_{N-1}}.\end{split}$
(62)

As we only consider steric interactions, we may employ a simple probabilistic
argument. If the average excluded volume of a single particle in $\Gamma$ is
then $v^{0}_{ex}$, then the probability of not intersecting a particular
particle is $1-\frac{v^{0}_{ex}}{|\Gamma|}$, where we assume $v_{ex}^{0}$ to
be independent of the sample size. We may thus write the probability of a
probe particle not intersecting any particle as

$\left(1-\frac{v_{ex}^{0}}{|\Gamma|}\right)^{N}=\left(1-\frac{v^{0}_{ex}\rho}{N}\right)^{N},$
(63)

as all the particles are independently distributed.111This argument may be
readily generalised in a straightforward manner to soft interactions with
sufficient decay at a large distance. This gives the cavity volume fraction in
the thermodynamic limit as $N,|\Gamma|\to\infty$ as

$\mathcal{V}_{f}(\rho)=\exp(-v^{0}_{ex}\rho).$ (64)

We may substitute this into the free energy density equation (10) to obtain

$\begin{split}\frac{1}{kT}\mathcal{F}(\rho)=&\rho\ln\rho-\rho-\int_{0}^{\rho}\ln\left(\exp(-v^{0}_{ex}x)\right)\,dx\\\
=&\rho\ln\rho-\rho+v^{0}_{ex}\int_{0}^{\rho}x\,dx\\\
=&\rho\ln\rho-\rho+\frac{v^{0}_{ex}}{2}\rho^{2},\end{split}$ (65)

which readily obtains the equation of state

$\frac{1}{kT}P=\rho+\frac{v_{ex}^{0}}{2}\rho^{2}.$ (66)

This produces Onsager’s acclaimed free energy for a single species isotropic
liquid [32], and correctly reclaims the second virial coefficient. Perhaps
curiously, the factor of $\frac{1}{2}$ before the excluded volume does not
arise as a traditional ad-hoc tool to avoiding double counting, but in a much
more indirect manner via the integration step. Furthermore, this approach does
not use an explicit mean-field, as is usually tradition in justifications of
Onsager (see for example [33], and Onsager’s original work), but instead uses
a purely geometric ansatz on the systems considered. While we do not attempt
to estimate the error, it is immediate that the error of the approximation is
precisely controlled by the validity of the approximation
$P_{N}\approx\text{Un}_{N}$, and in this sense we may view Onsager as an
uncorrelated limit of the system, rather than a dilute limit, although the two
interpretations are not independent.

###  3.1 A fluctuating lattice approximation

As exact solutions are generally unavailable, and detailed numerical
simulations are expensive and noisy, here we propose an alternative method of
approximating the cavity volume, and thus the free energy and equation of
state, that is efficient to compute and based on capturing the salient aspects
of the real physical system. The key idea is to presume that the cavity volume
behaves as if the local environment of each particle were a lattice, although
the lattice parameters may be subject to spatial fluctuations according to a
probability distribution. As such we will refer to the model as the
fluctuating lattice model. In particular, this contrasts to cell theories
which approximate the entire system by a single lattice, defined globally.

Consider a system of identical hard particles where the densest packing of $N$
particles forms a known lattice of volume $Nv_{c}$, where $v_{c}$ is the
volume of the Voronoi cell in this lattice. If a particular configuration has
total volume $V$, we define the excess volume to be $V-Nv_{c}$. We then
propose a method for approximating the cavity volume fraction according to two
assumptions:

  * •

We can attribute to each particle an individual excess volume $v_{e}$, such
that $N\langle v_{e}\rangle=V-Nv_{c}$, and $v_{e}$ corresponds to the amount
of excess volume living within the particle’s Voronoi cell. $v_{e}$ is a
quantity that is exchanged freely and reversibly between particles as the
system evolves.

  * •

Given a particular particle with excess volume $v_{e}$, its local environment
is well approximated by a uniform dilation of the densest packed lattice that
gives the same excess volume per particle.

The first assumption implies that $v_{e}$ should be distributed according to a
Boltzmann distribution (see [34, Section 1]), so
$\mathbb{P}(v_{e})=\nu\exp(-\nu v_{e})$, with $\langle
v_{e}\rangle_{\mathbb{P}}=\nu^{-1}$. As the total excess volume is $V-Nv_{c}$,
$\nu=\left(\frac{V}{N}-v_{c}\right)^{-1}=\frac{\rho^{*}\rho}{\rho^{*}-\rho}$
where $\rho=\frac{N}{V}$ is the number density and $\rho^{*}=\frac{1}{v_{c}}$
is the number density at the densest packing.

The true distribution of Voronoi cell sizes in hard particle systems has been
investigated previously by Senthil Kumar and Kumaran [38] and empirically
shown to follow a 2- or 3-parameter $\Gamma$ distribution. Within their work
they also review the case for Poisson distributed points (corresponding to the
dilute limit), showing there are still conflicts within the literature even in
this simpler case. By using the more naive Boltzmann distribution, we give an
explicit and analytic form for the probability distribution function
$\mathbb{P}$, and capture salient aspects of the problem.

The second assumption means that (locally) we have a dimensionless lattice
parameter $a>1$, so that the local environment of a particle resembles a
dilation of the densest lattice by factor $a$. We can relate $a$ to
$v_{c},v_{e}$ as $v_{e}+v_{c}=a^{d}v_{c}$ where $d$ is the dimension of space,
which can be rearranged to give $a=\sqrt[d]{\frac{v_{e}}{v_{c}}+1}$.

Given a lattice described by spacing $a$, we denote $v_{f}(a)$ to be the
amount of cavity volume per Voronoi cell. We may use $v_{f}$ to approximate
the cavity volume fraction of the entire system as

$\begin{split}&V_{f}(N,V)\\\
\approx&\frac{N}{V}\int_{0}^{\infty}v_{f}\left(a(v_{e})\right)\mathbb{P}(v_{e})\,dv_{e}\\\
=&\frac{\rho^{*}\rho^{2}}{\rho^{*}-\rho}\int_{0}^{\infty}v_{f}\left(\sqrt[d]{\frac{v_{e}}{v_{c}}+1}\right)e^{-\frac{\rho^{*}\rho}{\rho^{*}-\rho}v_{e}}\,dv_{e}\end{split}$
(67)

This equation is in principle applicable to general systems, but from here we
limit ourselves to hard spheres in $\mathbb{R}^{d}$. Before doing any
numerical computations, we make a few qualitative statements about the
constitutive equation.

It is immediate that for $0<\rho<\rho^{*}$ the integral in
(LABEL:eqFluctuatingLatticeApprox) is an analytic function of $\rho$, and
consequently the free energy will be smooth over its entire domain. Any
discontinuity in the free energy or its derivatives would necessarily have to
be accompanied by corresponding discontinuities in the map
$\rho\mapsto\mathbb{P}$. This statement is in accordance with the work of
Senthil Kumar and Kumaran [38], which shows a sharp transition in the
distribution of Voronoi cell sizes as the freezing transition is reached.
Phase separation is more complicated to describe, and would have to arise from
non-monotonicity of $\mathcal{V}_{f}$, leading to non-convexity in
$\mathcal{F}$. The map
$\rho\to\frac{\rho^{*}\rho^{2}}{\rho^{*}-\rho}\exp\left(-\frac{\rho^{*}\rho}{\rho^{*}-\rho}v_{e}\right)$
is bell-shaped for any $v_{e}>0$, which makes it unclear if $\mathcal{V}_{f}$
is monotonic in general, and consequently unclear if we may rule out phase
separation.

Different constitutive equations on $v_{f}$ and $\mathbb{P}$ may be employed
in (LABEL:eqFluctuatingLatticeApprox) to provide alternative models. The
simplest such change would be to permit different lattice structures.
Employing a broader family of probability distributions, perhaps based on
theoretical approximations or fitting to numerical experiments, would open
other avenues for investigation as well. For the sake of this work however we
will limit ourselves to the simplest possible assumptions that can provide
tangible results. As we demonstrate in the sequel, these incredibly coarse
assumptions are sufficient to describe the equation of state both in the high-
density and dilute regimes, within a single framework. This puts it at an
advantage over theories such as the virial expansion, only valid in dilute
regimes, and cell theories, which are only valid in dense regimes.

In particular, we can estimate the leading order behaviour in these extreme
regimes, with the calculations deferred to Appendix B. We find that

$\begin{split}\mathcal{V}_{f}\approx&1-2^{d}\rho v_{0},\end{split}$ (68)

as $\rho\to 0$, in accordance with the virial expansion. Futhermore, if there
exists some $v_{1}>0$ such that for $v_{e}<v_{1}$, the cavity volume is zero,
then for every $\epsilon>0$ there exist constants
$c_{\epsilon},C_{\epsilon}>0$ so that

$c_{\epsilon}\nu\exp(-\nu(1+\epsilon)v_{1})<\mathcal{V}_{f}(\rho)<C_{\epsilon}\nu\exp(-\nu(1-\epsilon)v_{1})$
(69)

as $\rho\to\rho^{*}$, where $\nu=\frac{\rho^{*}\rho}{\rho^{*}-\rho}$ as
before. In particular, we may estimate the free energy and pressure as
$\rho\to\rho^{*}$ by

$-c_{1}\ln(\rho^{*}-\rho)\leq\mathcal{F}(\rho)\leq-c_{2}\ln(\rho^{*}-\rho),$
(70) $\frac{c_{3}}{\rho^{*}-\rho}<\frac{P}{kT}<\frac{c_{4}}{\rho^{*}-\rho},$
(71)

for appropriate positive constants $c_{i}$.

####  3.1.1 1D Tonks gas

We will consider the Tonks gas as it is a simple enough system that the
fluctuating lattice approximation may be calculated exactly and explicitly, as
well as compared to the exact free energy density and equation of state, as in
Section 3.0.1. We have that $v_{c}=2r$, and $v_{f}(a)=\max(0,a-2r)$. Thus the
cavity volume fraction can readily be computed with a change of variables
$x=v_{e}-2r$ as

$\begin{split}\mathcal{V}_{f}=&\rho\int_{0}^{\infty}\max(0,v_{e}-2r)\nu\exp(-\nu
v_{e})\,dv_{e}\\\ =&\rho\int_{2r}^{\infty}\max(0,v_{e}-2r)\nu\exp(-\nu
v_{e})\,dv_{e}\\\ =&\rho\int_{0}^{\infty}x\nu\exp(-\nu(x+2r)\,dx\\\
=&\rho\exp(-2r\nu)\int_{0}^{\infty}\nu x\exp(-\nu x)\,dx\\\
=&\frac{\rho\exp(-2r\nu)}{\nu}=\frac{\rho^{*}-\rho}{\rho^{*}}\exp\left(-2r\frac{\rho\rho^{*}}{\rho^{*}-\rho}\right).\end{split}$
(72)

Introducing $\sigma=2r=\frac{1}{\rho^{*}}$ to make explicit the volume of the
particles, this can be written as

$\mathcal{V}_{f}=(1-\rho\sigma)\exp\left(-\frac{\rho\sigma}{1-\rho\sigma}\right),$
(73)

which is the exact result in one dimension as seen in Section 3.0.1. Therefore
any and all derived quantities (pressure, free energy) will also be correct
using this approach. This is however a truly exceptional situation for the
one-dimensional system. In a Tonks gas, our assumption of simple pairwise
exchange of excess volume is more accurate, as the excess volume can be
exchanged between next-nearest neighbours via translating the particle between
them. In higher dimensional systems where the geometry of a cavity is no
longer trivial, translations of a single particle would be expected to
redistribute the excess volume in a highly complicated way between both its
various neighbours and itself.

####  3.1.2 2D Hexagonal lattice

(a) The cavity volume that remains in a relatively dilute hexagonal lattice
($a=4.2$)

(b) The cavity volume that remains in a denser hexagonal lattice ($a=3.85$)

Figure 2:  The figures above show the cavity volume within a single Voronoi
cell. The black circles correspond to the particles in the system, the white
circle surrounding the central particle is its exclusion region, and the blue
region surrounding it is the cavity volume within the Voronoi cell of the
central particle.

In the case of an hexagonal lattice, we can explicitly find the cavity volume
function, which is done in Section C.1. We evaluate the resulting equation of
state by numerically integrating (LABEL:eqFluctuatingLatticeApprox) to obtain
$\mathcal{V}_{f}$, and then performing a second numerical integration in the
equation of state (11). This is then compared with the Monte Carlo simulation
data of Kolafa and Rottner [26] in the fluid regime, and Alder, Hoover and
Young [1] in the solid regime in Figure 3. We constrast our results with the
two-dimensional analogue of the Carnahan-Starling equation of state, which
gives the compressibility factor in terms of the packing fraction $\eta$ as
$Z=(1-\eta)^{2}$, which is expected to be accurate before the freezing
transition at around $\eta\approx 0.706$. Furthermore, we make comparison with
the leaky cell theory for an hexagonal lattice [14], which is an extension of
the classical cell theory [7], to dilute regimes. The cell theory approaches
are expected to be better approximations in dense regimes. We see good
agreement with the fluid regime up to moderate densities, with accuracy
becoming worse as we approach the freezing transition. However once we enter
the solid regime, we see very good agreement between the fluctuating lattice
model and the Monte Carlo data. We include the relative error between the
measurements of Alder, Hoover and Young and the fluctuating lattice model in
Figure 4(a), and we see that the relative accuracy is improving towards the
dense limit, and furthermore is within 5% beyond $\eta=0.8$. As the quality of
approximation of the equation of state is directly related to the quality of
approximation of the cavity volume, we show in Figure 4(b) a comparison of the
cavity volume against an ad-hoc Monte Carlo experiment up to modest densities
($\eta=0.6$, details in Appendix A), in which we see a consistent over-
estimate. Due to the expectation of vanishingly small cavity volume in higher-
density regimes, where incredibly large systems and long runs would be
required to obtain stable results, we do not probe the high-density regime in
our Monte Carlo methods.

Given our asumptotic result (69) in the dense regime, we expect that
$\mathcal{V}_{f}$ should behave as
$c_{1}\exp\left(-\frac{c_{2}}{\rho^{*}-\rho}\right)$ as $\rho\to\rho^{*}$ for
some constants $c_{1},c_{2}$. We perform a linear least-squares regression of
$\frac{-1}{\ln\mathcal{V}_{f}}$ versus $\eta^{*}-\eta$. We see a good linear
fit, as shown in Figure 4(c), so that $\frac{-1}{\ln\mathcal{V}_{f}}\approx
0.498(\eta^{*}-\eta)$, which would give the leading order approximation as
$\mathcal{V}_{f}\approx ce^{-\frac{2.08}{\eta^{*}-\eta}}$.

Figure 3: Comparison of the data of Kolafa and Rotter (KR), and Alder, Hoover
and Young (AHY) against the obtained equation of state for the fluctuating
lattice model (FL) for the two-dimensional hard disk system. The two-
dimensional Carnahan-Starling (CS) and Leaky Cell Theory (LCT) equations of
state are included for comparison.

(a) Quotients $Q=P_{HAY}/P$ with the pressure $P$ obtained from the Leaky Cell
Theory (LCT) and Fluctuating Lattice (FL) against the numerical data $P_{HAY}$
of Alder, Hoover and Young in the dense regime.

(b) Comparison of the cavity volume fraction from the fluctuating lattice
approximation with ad-hoc Monte Carlo runs, and the inferred cavity volume of
the 2-dimensional Carnahan-Starling, obtained by (14), using a logarithmic
scale.

(c) Linear fit of $-(\ln\mathcal{V}_{f})^{-1}$ against $\alpha(\eta^{*}-\eta)$
in the dense regime with optimal constant found to be $\alpha\approx 0.498$

Figure 4: Further comparisons of the fluctuating lattice model for the two-
dimensional hard disk system

####  3.1.3 3D FCC lattice

We now consider the case of a 3-dimensional FCC lattice. We defer calculations
of the cavity volume function $v_{f}$ to Section C.2. We compare the packing
fraction $\eta=\rho\frac{4\pi}{3}r^{3}$ against the compressibility
$Z=\frac{P}{\rho kT}$ in Figure 6. Furthermore, as in the two-dimensional
case, we compare the obtained equation of state against the classical
Carnahan-Starling equation of state given in terms of the packing density as
$Z=\frac{1+\eta+\eta^{2}-\eta^{3}}{(1-\eta)^{3}}$, the Leaky Cell Theory for
an FCC lattice, and Monte Carlo data, in this case from Wu and Sadus [45].

(a) The cavity volume that remains in a relatively dilute FCC lattice
($a=4.5$)

(b) The cavity volume that remains in a denser FCC lattice ($a=3.85$)

Figure 5:  The figures above show the cavity volume within a single Voronoi
cell. Only the part of the cavity volume with $x_{3}>0$ is shown for
visibility. The black spheres correspond to the particles in the system. The
region around the central particle is the cavity volume within the Voronoi
cell of the central particle.

Figure 6: Comparison of the fluctuating lattice (FL) equation of state for
three dimensional hard spheres, with the commonly used Carnahan-Starling (CS),
the Leaky Cell Theory (LCT), and comparison to Monte Carlo simulation data of
Wu and Sadus (WS), with the fluid branch and solid branch of their data
separated for comparison. Note that for $\eta$ greater than $\approx 0.502$,
the fluid branch is only metastable.

As expected from our analysis, there is generally good agreement with the
Monte-Carlo data and the Carnahan-Starling equations of state in the dilute
regime. However in intermediate, fluid, regimes we start to see stronger
disagreement. We see that for the entirety of the solid branch we have a
modest agreement with the fluctuating lattice equation of state. In Figure
7(a), we see a monotonic increase in relative accuracy of the fluctuating
lattice equation of state in the solid phase. Furthermore, we see that the
fluctating lattice equation of state always provides an underestimate of the
compressibility factor. Similarly to the two-dimensional case, we perform a
linear fit of $-(\ln\mathcal{V}_{f})^{-1}$ against $\alpha(\eta^{*}-\eta)$,
shown in Figure 7(c), which shows good agreement with $\alpha=0.413$. This
implies an asymptotic approximation $\mathcal{V}_{f}\approx
c\exp\left(-\frac{2.32}{\eta^{*}-\eta}\right)$ in the dense limit. We can
directly compare the calculated cavity volume against Monte Carlo data in
intermediate regimes, which we show in Figure 7(b).

(a) The quotients $Q=P_{WS}/P$ of the Leaky Cell Theory (LCT) and Fluctuating
Lattice (FL) pressures $P$ with the data of Wu and Sadus $P_{WS}$ on the solid
branch.

(b) Comparison of the fluctating lattice prediction of cavity volume fraction
against Monte Carlo simulation, and the inferred cavity volume of the
3-dimensional Carnahan-Starling, obtained by (14), using a logarithmic scale.

(c) Linear fit of $-\ln(\mathcal{V}_{f})^{-1}$ against $\alpha(\eta^{*}-\eta)$
in the dense regime with optimal constant found to be $\alpha\approx 0.413$

Figure 7:  Further comparisons of the fluctuating lattice model for the three-
dimensional hard sphere system

##  4 Concluding remarks

Within our work we have used the notion of cavity volume to provide
expressions for the free energy and equation of state in single-species
systems by performing a type of thermodynamic integration with respect to
number density from the vacuum state. The framework we have used, owing to its
generality, permits us to reclaim known results such as the Onsager
approximation and exact solution of the Tonks gas in one dimension. Our
strategy offers new expressions for the virial coefficients given in terms of
the cavity volume in the dilute regime, showing that if the cavity volumes of
systems of up to $n-1$ particles are known, in principle, we may evaluate the
reduced virial coefficient $\text{b}_{n}$.

This has strengthened our confidence in the theoretical setting envisioned
here. A number of novel conclusions come from our proposed approach. We
systematised a heuristic expression for the cavity volume, which had already
proven predictive in [31]. We profited from the generality we can afford by
identifying a constitutive function $\lambda$ that we have newly interpreted
as an effective correction for overlapping exclusion zones between interacting
particles.

In addition, we proposed a fluctuating lattice model, which is obtained by
making simplifying assumptions in order to calculate the cavity volume
fraction in hard particle systems. The model has shown that with a small
number of assumptions, we can capture the salient aspects of the equation of
state for hard rods, disks and spheres, in different space dimensions, in both
the dilute and dense regimes. This offers an advantage over the virial
expansion and Carnahan-Starling equation of state, which are based on series
expansions in the dilute regime and therefore cannot be expected to provide
accurate results beyond a phase transition; as well as advantages over the
cell theories, which perform well near dense packings but however fail to
capture the correct behaviour in dilute regimes.

The main limitation of the fluctuating lattice model as presented in our work
is that it does not capture phase transitions, presumably as a consequence of
the independence of the cell cavity volume function $v_{f}$ and the density,
which ignores qualitative changes in the geometry of the Voronoi cell as the
density varies. Similarly, the assumptions that lead to a Boltzmann
distribution of the excess cavity volume, while being convenient for the
arithmetic and analysis of the free energy, contradict known results even in
the dilute regime (see [15, 38]). Improving the structural assumptions on the
distribution of excess cavity volume and the function $v_{f}$ remains an
avenue for future investigation, whereby permitting competition between
different lattice structures could lead to a model capable of describing phase
transitions. Furthermore, while only hard interactions have been considered in
the fluctuating lattice model within this work, the same methodology would be
extendable to long-range, soft interactions, which may or may not exist in
tandem with hard-core repulsion at short range. As well as broadening the
nature of particle interactions, we believe this framework will be extendable
to systems with more complex structures, namely anisotropic and multi-species
systems. We have in mind liquid crystalline systems as a future avenue, where
suitable adjustments to the model could provide a classical density functional
theory for spatially inhomogeneous liquid crystalline systems, similar to the
approach as seen in [31].

##  5 Acknowledgments

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”.

##  Appendix A MC simulations

Approximations and constitutive equations of the excluded volume can tell us
certain information relating to the system, but ultimately in 2 or more
dimensions we are constrained to numerical methods if we wish to obtain
quantitatively accurate data for the excluded volume of even modest density
systems. For this we perform a simple Monte Carlo study to provide results for
comparison. The methodology is as follows:

  1. 1.

Equilibrate a configuration by moving each particle at least 50 successful
Monte Carlo steps.

  2. 2.

Estimate the cavity volume by testing if 150 randomly chosen points lie in any
exclusion regions.

  3. 3.

Generate a new realisation by further Monte Carlo methods on the previous
realisation, until all particles have been successfully moved at least 10
times.

  4. 4.

Go back to (2), unless the error of the logarithm of the cavity volume is
sufficiently small.

We consider the error of the logarithm as this is the quantity appearing in
the energy, and this means that in denser systems a significantly larger
number of runs is required as $\mathcal{V}_{f}$ is close to zero, meaning
small fluctuations lead to large fluctuations in $-\ln\mathcal{V}_{f}$. We
estimate the error of the logarithm as follows. Given $n$ test points and $s$
points in the cavity space, we estimate a 95% confidence interval using the
Wilson score interval,

$p=\frac{2s+z^{2}}{2n+2z^{2}}\pm\frac{z}{n+z^{2}}\sqrt{\frac{s(n-s)}{n}+\frac{z^{2}}{4}}.$
(74)

Here $z=1.96$ corresponds to the $95\%$ confidence interval, $s$ is the number
of successful trials and $n$ is the total number of trials. Then we estimate
the error of the logarithm as the width of the logarithm of this interval.

##  Appendix B Limiting behaviour in the fluctuating lattice model

###  B.1 Dilute systems

Let the local lattice spacing be $a$. We note that if $a=1$ we are at the
densest packing and no particles overlap, and are thus spaced at least by
twice their radius $r$. When $a>2$, all particles must be at distance at least
$4r$, meaning no exclusion regions overlap. In particular, this means that
$v_{f}(a)=a^{d}v_{c}-2^{d}v_{0}$, where $v_{0}$ is the volume of one particle
and $2^{d}v_{0}$ is the pairwise excluded volume. This can be simplified as
$v_{f}(a(v_{e}))=v_{e}+v_{c}-2^{d}v_{0}$. This means we can write the cavity
volume fraction, writing $\nu=\frac{\rho^{*}\rho}{\rho^{*}-\rho}$ and $v_{1}$
to be the solution of $a(v_{1})=2$ as

$\begin{split}&\mathcal{V}_{f}\\\ =&\rho\left(\int_{0}^{v_{1}}\nu\exp(-\nu
v_{e})\left(v_{f}(a(v_{e}))-a(v_{e})^{d}v_{c}+2^{d}v_{0}\right)\,dv_{e}\right.\\\
&\hskip
14.22636pt\left.+\int_{0}^{\infty}(a(v_{e})^{d}v_{c}-2^{d}v_{0})\nu\exp(-\nu
v_{e})\,dv_{e}\right)\\\
\approx&\rho\left(c\nu+\int_{0}^{\infty}(v_{c}-2^{d}v_{0}+v_{e})\nu\exp(-\nu
v_{e})\,dv_{e}\right)\\\
=&\rho\left(c\nu+v_{c}-2^{d}v_{0}+\frac{\rho^{*}-\rho}{\rho^{*}\rho}\right),\end{split}$
(75)

where $c$ is an order 1 quantity as $\rho\to 0$, or equivalently $\nu\to 0$.
Thus in the dilute regime $\rho\approx 0$, we have that the leading order
contribution to $\mathcal{V}_{f}$ is

$\begin{split}\mathcal{V}_{f}\approx&\rho v_{c}-\rho
2^{d}v_{0}+1-\frac{\rho}{\rho^{*}}\\\ =&1-2^{d}\rho v_{0},\end{split}$ (76)

recalling $v_{c}=\frac{1}{\rho^{*}}$. We thus reclaim the expected first order
behaviour of the cavity volume. It should be noted that this argument depends
only on the fact that $\mathbb{P}(v_{e})\to 0$ for $0<v_{e}<v_{1}$ as $\rho\to
0$, and would extend to more general systems without issue if the constitutive
probability distributions $\mathbb{P}$ satisfied this.

###  B.2 Dense systems

Let $\rho\approx\rho^{*}$. Assume that $v_{f}(a(v_{e}))=0$ for $v_{e}<v_{2}$,
and $v_{2}>0$. That is, in a nearly-dense packing there is no cavity volume,
which is certainly reasonable for spheres. Then for any fixed $1>\delta>0$,
$\nu$ sufficiently large, and $v_{e}>v_{2}$, we have that $\nu
v_{e}>(\delta\nu-1)v_{e}>\delta\nu v_{2}-v_{e}$. Then this means that

$\begin{split}&\mathcal{V}_{f}\\\ =&\int_{v_{2}}^{\infty}\rho
v_{f}\left(\sqrt[d]{\frac{v_{e}}{v_{c}}+1}\right)\nu\exp(-\nu
v_{e})\,dv_{e}\\\ \leq&\nu\exp(-\delta\nu v_{2})\int_{v_{2}}^{\infty}\rho
v_{f}\left(\sqrt[d]{\frac{v_{e}}{v_{c}}+1}\right)\exp(-v_{e})\,dv_{e}.\end{split}$
(77)

We note that the integral is finite and independent of $\nu$. Therefore we
have a bound

$\mathcal{V}_{f}\leq c_{\delta}\nu\exp(-\delta
v_{2}\nu)=\frac{c_{\delta}\rho\rho^{*}}{\rho^{*}-\rho}\exp\left(-\frac{\delta\rho^{*}\rho}{\rho^{*}-\rho}v_{2}\right).$
(78)

for some $c_{\delta}>0$ and $0<\delta<1$, for sufficiently large $\nu$
(equivalently, $\rho$ sufficiently close to $\rho^{*}$).

By a similar argument, if $\gamma>1$, using that $\nu
v_{e}<(\gamma\nu+1)v_{e}\leq\gamma\nu v_{e}+\gamma^{2}v_{2}$ for
$v_{2}<v_{e}<\gamma v_{2}$,

$\begin{split}&\mathcal{V}_{f}\\\ =&\int_{v_{2}}^{\infty}\rho
v_{f}\left(\sqrt[d]{\frac{v_{e}}{v_{c}}+1}\right)\nu\exp(-\nu
v_{e})\,dv_{e}\\\ \geq&\nu\exp(-\nu\gamma^{2}v_{2})\int_{v_{2}}^{\gamma
v_{2}}\rho
v_{f}\left(\sqrt[d]{\frac{v_{e}}{v_{c}}+1}\right)\exp(-v_{e})\,dv_{e}.\end{split}$
(79)

This gives a similar scaling, that for any $\gamma>1$, for $\rho$ sufficiently
close to $\rho^{*}$,

$\begin{split}&\mathcal{V}_{f}\\\
\geq&c_{\gamma}\nu\exp(-\nu\gamma^{2}v_{2})\\\
=&\frac{c_{\gamma}\rho\rho^{*}}{\rho-\rho^{*}}\exp\left(-\frac{\gamma^{2}\rho\rho^{*}}{\rho^{*}-\rho}v_{2}\right).\end{split}$
(80)

This means that the cavity volume tends to zero faster than any polynomial,
and the upper bound is a near-classical example of a function with vanishing
derivatives of all orders as $\rho\to\rho^{*}$, despite being a non-zero
function. Written in another way, we can say that for sufficiently small
$\epsilon>0$ there exists constants $c_{\epsilon},C_{\epsilon}>0$ so that

$\begin{split}&c_{\epsilon}\nu\exp(-\nu(1+\epsilon)v_{1})\\\
&<\mathcal{V}_{f}(\rho)\\\
&<C_{\epsilon}\nu\exp(-\nu(1-\epsilon)v_{1}).\end{split}$ (81)

The argument we have presented exemplifies a key point that we have seen in
the Tonks gas, that the leading cause for the loss of cavity volume is not
that particular cavities are exceptionally small, but that they are
exceptionally rare, in dense systems. This is quantified in the preceding
argument by noting that moderately dense lattices still have zero cavity
volume, and thus we may estimate the total cavity volume by the “tail” of the
distribution $\mathbb{P}$, which is rapidly decaying.

Substituting these inequalities into (10) provides the estimates

$-c_{1}\ln(\rho^{*}-\rho)\leq\mathcal{F}(\rho)\leq-c_{2}\ln(\rho^{*}-\rho),$
(82) $\frac{c_{3}}{\rho^{*}-\rho}>\frac{P}{kT}>\frac{c_{4}}{\rho^{*}-\rho}$
(83)

as $\rho\to\rho^{*}$.

##  Appendix C Cavity volume calculations

###  C.1 Cavity volume of a 2D Hexagonal lattice

As the configuration is highly symmetric it suffices to consider half a unit
cell, which at spacing $a>0$ has three relevant particles sitting on vertices
of an equilateral triangle of side length $2ar$, where $r$ is the particle
radius. We see there are three regimes. When $a<\sqrt{3}$ the centre point is
contained in all exclusion regions and there is zero cavity volume. When $a>2$
all particles are at a distance of at least $4r$ and no exclusion zones
overlap. thus it suffices to evaluate the overlaps in the intermediate regime.
For $\sqrt{3}<a<2$, we must account for the overlap of the exclusion regions.
We may thus compute the cavity volume within the triangle as

$\frac{\sqrt{3}}{4}(2ar)^{2}-3\frac{1}{6}4\pi r^{2}+\frac{3}{2}V_{2}(a).$ (84)

The first term corresponds to the area of the triangle, the second to the
removal of the exclusion regions of the three particles, of which only a sixth
is inside the triangle. The final term accounts for double overlaps, where
$V_{2}(a)$ is the intersection area of two disks of radius $2r$ at a distance
$2ar$, given explicitly by
$V_{2}(a)=8r^{2}\cos^{-1}\left(\frac{a}{2}\right)-2ar^{2}\sqrt{4-a^{2}}$. The
densest packing occurs at $\rho^{*}=\frac{1}{2\sqrt{3}r^{2}}$.

###  C.2 Cavity volume function of a 3D FCC lattice

Here we calculate the cavity volume function for an FCC lattice of spacing
$a>1$. We perform this analysis by considering a unit cell, which contains the
equivalent of four particles (three halves, and eight quarters). The unit cell
can be decomposed into $N_{t}=8$ tetrahedra and $N_{o}=4$ octahedra, where
each vertex is the centre of mass of a particle. We note that the following
calculations also apply to an HCP structure, as its unit cell can be
deconstructed into equivalent blocks. Thus it suffices to find the cavity
volume in a particular tetrahedron and octohedron. We introduce the following
notation. $V_{2}(a)$ will denote the intersection volume of two spheres of
radius $2$ at distance $2a$. $V_{3}^{eq}(a)$ denotes the intersection volume
of three spheres on the vertices of an equilateral triangle of side length
$2a$, and $V_{3}^{rt}(a)$ denotes the intersection volume of three spheres on
the vertices of a right-angled triangle of hypotenuse length $2\sqrt{2}a$. We
needn’t consider intersections of higher order, as we will have hard particles
of radius $1$, and in an FCC lattice the intersection of four exclusion
spheres only appears when the cavity volume is zero. We can evaluate the
corresponding functions using formulae from [18] as

$\begin{split}V_{2}(a)=&\left\\{\begin{array}[]{c
l}\frac{\pi}{12}(4-2a)^{2}(8+2a)&a<2\\\ 0&a\geq 2\end{array}\right.\\\
V_{3}^{eq}(a)=&\left\\{\begin{array}[]{c
l}\frac{2a^{2}}{3}\sqrt{12-4a^{2}}+32\tan^{-1}\left(\frac{\sqrt{12-4a^{2}}}{2}\right)-3a\left(8-\frac{2a^{2}}{3}\right)\tan^{-1}\left(\frac{\sqrt{12-4a^{2}}}{a}\right)&a<\sqrt{3}\\\
0&a\geq\sqrt{3}\end{array}\right.\\\ V_{3}^{rt}(a)=&\left\\{\begin{array}[]{c
l}\frac{2a^{2}}{3}\sqrt{16-8a^{2}}-\frac{a\pi(24-a^{2})}{3\sqrt{2}}+\frac{32}{3}\left(\frac{\pi}{2}+2\tan^{-1}\left(\frac{\sqrt{16-8a^{2}}}{4}\right)\right)\\\
-2a\left(8-\frac{2a^{2}}{3}\right)\tan^{-1}\left(\frac{\sqrt{16-8a^{2}}}{2a}\right)&a<\frac{2\sqrt{6}}{3}\\\
0&a\geq\frac{2\sqrt{6}}{3}\end{array}\right.\end{split}$ (85)

Furthermore, we introduce
$\Omega_{t}=\cos^{-1}\frac{23}{27},\Omega_{o}=4\sin^{-1}\frac{1}{3}$ as the
solid angles of the vertices of the tetrahedron and octahedron, respectively,
and $\alpha_{t}=\cos^{-1}\frac{1}{3},\alpha_{0}=\cos^{-1}-\frac{1}{3}$, as the
corresponding dihedral angles. The volumes of the corresponding tetrahedron
and octahedron respectively are
$V_{t}=\frac{\sqrt{2}(2a)^{3}}{6},V_{o}=\frac{\sqrt{2}(2a)^{3}}{3}$. Then we
can evaluate the cavity volumes of a single tetrahedron and octahedron,
$v_{f}^{t},v_{f}^{o}$ respectively as

$\begin{split}v_{f}^{t}(a)=&V_{t}-4\frac{\Omega_{t}}{4\pi}\frac{32\pi}{3}+6\frac{\alpha_{t}}{2\pi}V_{2}(a)-4\frac{1}{2}V_{3}^{eq}(a)\\\
v_{f}^{o}(a)=&V_{t}-6\frac{\Omega_{o}}{4\pi}\frac{32\pi}{3}+12\frac{\alpha_{o}}{2\pi}V_{2}(a)-8\frac{1}{2}V_{3}^{eq}(a)-4V_{3}^{rt}(a)\\\
\end{split}$ (86)

when $a>\frac{2\sqrt{3}}{6}$ and zero otherwise. These formulae arise simply
from adding even-numbered overlaps and subtracting odd-numbered overlaps, and
factors $\frac{\Omega_{X}}{4\pi}$, $\frac{\alpha_{X}}{2\pi}$ account for the
ratio of the sphere contained within the polyhedra, with the integer
coefficients count how many relevant (combinations of) spheres need to be
considered. This gives the cavity volume per particle as

$v_{f}(a)=\frac{1}{4}\left(8v_{f}^{t}(a)+4v_{f}^{o}(a)\right),$ (87)

where the prefactor $\frac{1}{4}$ accounts for the fact there are four
particles per unit cell, and the coefficients $8,4$ correspond to the number
of tetrahedra and octahedra in each cell.

## References

  * [1] Alder, B., Hoover, W., and Young, D. Studies in molecular dynamics. V. High-density equation of state and entropy for hard disks and spheres.  The Journal of Chemical Physics 49, 8 (1968), 3688–3696.
  * [2] Alder, B. J., and Wainwright, T. E. Phase transition for a hard sphere system.  The Journal of Chemical Physics 27, 5 (1957), 1208–1209.
  * [3] Alder, B. J., and Wainwright, T. E. Studies in molecular dynamics. I. general method.  The Journal of Chemical Physics 31, 2 (1959), 459–466.
  * [4] Allen, M. P., Evans, G. T., Frenkel, D., and Mulder, B. Hard convex body fluids.  Advances in chemical physics 86 (1993), 1–166.
  * [5] Bowick, M. J., Kinderlehrer, D., Menon, G., and Radin, C. Mathematics and Materials, vol. 23.  American Mathematical Soc., 2017.
  * [6] Bowles, R. K., and Speedy, R. J. Cavities in the hard sphere crystal and fluid.  Molecular Physics 83, 1 (1994), 113–125.
  * [7] Buehler, R., Wentorf Jr, R., Hirschfelder, J., and Curtiss, C. The free volume for rigid sphere molecules.  The Journal of Chemical Physics 19, 1 (1951), 61–71.
  * [8] Clisby, N., and McCoy, B. M. Ninth and tenth order virial coefficients for hard spheres in D dimensions.  Journal of Statistical Physics 122, 1 (2006), 15–57.
  * [9] Comtet, L. Advanced Combinatorics: The art of finite and infinite expansions.  Springer Science & Business Media, 2012.
  * [10] Corti, D. S., and Bowles, R. K. Statistical geometry of hard sphere systems: exact relations for additive and non-additive mixtures.  Molecular Physics 96, 11 (1999), 1623–1635.
  * [11] Debenedetti, P. G., and Truskett, T. M. The statistical geometry of voids in liquids.  Fluid phase equilibria 158 (1999), 549–556.
  * [12] Dymond, J., Marsh, K., Wilhoit, R., and Wong, K. Virial coefficients of pure gases and mixtures, vol. 21.  Springer Landord-Bornstein, 2003.
  * [13] Eyring, H., and Hirschfelder, J. The theory of the liquid state.  Journal of physical chemistry 41, 2 (1937), 249–257.
  * [14] Fai, T. G., Taylor, J. M., Virga, E. G., Zheng, X., and Palffy-Muhoray, P. Leaky cell model of hard spheres.  The Journal of Chemical Physics 154, 10 (2021), 104505.
  * [15] Ferenc, J.-S., and Néda, Z. On the size distribution of poisson voronoi cells.  Physica A: Statistical Mechanics and its Applications 385, 2 (2007), 518–526.
  * [16] Flory, P. J. Thermodynamics of high polymer solutions.  The Journal of Chemical Physics 10, 1 (1942), 51–61.
  * [17] Frenkel, D., and Smit, B. Understanding molecular simulation: from algorithms to applications, vol. 1.  Elsevier, 2001.
  * [18] Gibson, K., and Scheraga, H. A. Volume of the intersection of three spheres of unequal size: a simplified formula.  Journal of Physical Chemistry 91, 15 (1987), 4121–4122.
  * [19] Hansen, J.-P., and McDonald, I. R. Theory of simple liquids.  Elsevier, 1990.
  * [20] Hirschfelder, J., Stevenson, D., and Eyring, H. A theory of liquid structure.  The Journal of Chemical Physics 5, 11 (1937), 896–912.
  * [21] Hoover, W. G., and Grover, R. Exact dynamical basis for a fluctuating cell model.  In Bulleting of the American Physical Society (1972), vol. 17, p. 80.
  * [22] Hoover, W. G., Hoover, N. E., and Hanson, K. Exact hard-disk free volumes.  The Journal of Chemical Physics 70, 4 (1979), 1837–1844.
  * [23] Huggins, M. L. Solutions of long chain compounds.  The Journal of Chemical Physics 9, 5 (1941), 440–440.
  * [24] Jackson, J., and Klein, L. Potential distribution method in equilibrium statistical mechanics.  The Physics of Fluids 7, 2 (1964), 228–231.
  * [25] Kamerlingh Onnes, H. Expression of the equation of state of gases and liquids by means of series.  Koninklijke Nederlandse Akademie van Wetenschappen Proceedings Series B Physical Sciences 4 (1901), 125–147.
  * [26] Kolafa, J., and Rottner, M. Simulation-based equation of state of the hard disk fluid and prediction of higher-order virial coefficients.  Molecular Physics 104, 22-24 (2006), 3435–3441.
  * [27] Lennard-Jones, J. E., and Devonshire, A. Critical phenomena in gases-I.  Proceedings of the Royal Society of London. Series A-Mathematical and Physical Sciences 163, 912 (1937), 53–70.
  * [28] Lennard-Jones, J. E., and Devonshire, A. Critical phenomena in gases. II. vapour pressures and boiling points.  Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences (1938), 1–11.
  * [29] Löwen, H. Fun with hard spheres.  In Statistical physics and spatial statistics. Springer, 2000, pp. 295–331.
  * [30] Mayer, J. E. The statistical mechanics of condensing systems. I.  The Journal of Chemical Physics 5, 1 (1937), 67–73.
  * [31] Nascimento, E. S., Palffy-Muhoray, P., Taylor, J. M., Virga, E. G., and Zheng, X. Density functional theory for dense nematic liquid crystals with steric interactions.  Physical Review E 96, 2 (2017), 022704.
  * [32] Onsager, L. The effects of shape on the interaction of colloidal particles.  Annals of the New York Academy of Sciences 51, 4 (1949), 627–659.
  * [33] Palffy-Muhoray, P., Virga, E. G., and Zheng, X. Onsager’s missing steps retraced.  Journal of Physics: Condensed Matter 29, 47 (2017), 475102.
  * [34] Pitaevskii, L., and Lifshitz, E. Physical Kinetics: Volume 10, vol. 10.  Butterworth-Heinemann, 2012.
  * [35] Santos, A. A concise course on the theory of classical liquids.  Lecture Notes in Physics 923 (2016).
  * [36] Santos, A., Yuste, S. B., and López de Haro, M. Structural and thermodynamic properties of hard-sphere fluids.  The Journal of Chemical Physics 153, 12 (2020), 120901.
  * [37] Sastry, S., Truskett, T. M., Debenedetti, P. G., Torquato, S., and Stillinger, F. H. Free volume in the hard sphere liquid.  Molecular Physics 95, 2 (1998), 289–297.
  * [38] Senthil Kumar, V., and Kumaran, V. Voronoi cell volume distribution and configurational entropy of hard-spheres.  The Journal of Chemical Physics 123, 11 (2005), 114501.
  * [39] Speedy, R. J. Statistical geometry of hard-sphere systems.  Journal of the Chemical Society, Faraday Transactions 2: Molecular and Chemical Physics 76 (1980), 693–703.
  * [40] Speedy, R. J., and Reiss, H. Cavities in the hard sphere fluid and crystal and the equation of state.  Molecular physics 72, 5 (1991), 999–1014.
  * [41] Sturgeon, K. S., and Stillinger, F. H. Cavities in the hard-disk crystal: A Monte Carlo simulation study.  The Journal of Chemical Physics 96, 6 (1992), 4651–4654.
  * [42] Tonks, L. The complete equation of state of one, two and three-dimensional gases of hard elastic spheres.  Physical Review 50, 10 (1936), 955.
  * [43] Wheatley, R. J. Calculation of high-order virial coefficients with applications to hard and soft spheres.  Physical Review Letters 110, 20 (2013), 200601.
  * [44] Widom, B. Some topics in the theory of fluids.  The Journal of Chemical Physics 39, 11 (1963), 2808–2812.
  * [45] Wu, G.-W., and Sadus, R. J. Hard sphere compressibility factors for equation of state development.  AIChE journal 51, 1 (2005), 309–313.