# Colloidal Hard Spheres: Triumphs, Challenges and Mysteries

C. Patrick Royall<EMAIL_ADDRESS>Gulliver UMR CNRS 7083, ESPCI
Paris, Université PSL, 75005 Paris, France.  School of Physics, HH Wills
Physics Laboratory, University of Bristol, Tyndall Avenue, Bristol, BS8 1TL,
UK.  School of Chemistry, University of Bristol, Cantock’s Close, Bristol, BS8
1TS, UK.  Patrick Charbonneau  Department of Chemistry, Duke University,
Durham, North Carolina 27708, USA.  Department of Physics, Duke University,
Durham, North Carolina 27708, USA.  Marjolein Dijkstra  Soft Condensed Matter,
Debye Institute for Nanomaterials Science, Utrecht University, Princetonplein
1, 3584 CC, Utrecht, The Netherlands.  John Russo  Department of Physics,
Sapienza University of Rome, Piazzale Aldo Moro 5, 00185 Rome, Italy.  Frank
Smallenburg  Université Paris–Saclay, CNRS, Laboratoire de Physique des
Solides, 91405 Orsay, France.  Thomas Speck  Institut für Theoretische Physik
IV, Universität Stuttgart, Heisenbergstr. 3, 70569 Stuttgart, Germany.
Chantal Valeriani  Departamento de Estructura de la Materia, Fisica Termica y
Electronica, Facultad de Ciencias Fisicas, Universidad Complutense de Madrid,
28040 Madrid, Spain.

###### Abstract

The simplicity of hard spheres as a model system is deceptive. Although the
particles interact solely through volume exclusion, that nevertheless suffices
for a wealth of static and dynamical phenomena to emerge, making the model an
important target for achieving a comprehensive understanding of matter. In
addition, while real colloidal suspensions are typically governed by complex
interactions, Pusey and Van Megen [_Nature_ 320 340–342 (1986)] demonstrated
that suitably tuned suspensions result in hard-sphere like behavior, thus
bringing a valuable experimental complement to the celebrated theoretical
model. Colloidal hard spheres are thus both a material in their own right and
a platform upon which phenomena exhibited by simple materials can be explored
in great detail. These various purposes enable a particular synergy between
experiment, theory and computer simulation. Here we review the extensive body
of work on hard spheres, ranging from their equilibrium properties such as
phase behavior, interfaces and confinement to some of the non–equilibrium
phenomena they exhibit such as sedimentation, glass formation and nucleation.

###### Contents

  1. I Introduction
  2. II Historical background
    1. II.1 History of hard-sphere models
    2. II.2 The foretelling: How PMMA colloids came to be
  3. III Realizing hard-sphere systems
    1. III.1 Experimental methods to prepare hard-sphere-like systems
    2. III.2 Synthesizing hard-sphere-like colloidal particles
    3. III.3 Interactions in real hard-sphere-like systems
    4. III.4 Mapping soft spheres to hard spheres
    5. III.5 Accuracy in hard-sphere experiments
  4. IV Measuring hard spheres _in vitro_
    1. IV.1 Light scattering
    2. IV.2 Optical microscopy and particle-resolved studies
    3. IV.3 Size _does_ matter: colloid dynamics in different experiments
  5. V Hard-sphere systems _in silico_
    1. V.1 Monte Carlo methods
    2. V.2 Event-driven simulations
    3. V.3 Brownian dynamics simulations
    4. V.4 Simulations including hydrodynamics
    5. V.5 Advanced sampling and free-energy methods
    6. V.6 Simulation models
  6. VI Hard-sphere systems in theory
    1. VI.1 Integral equation theory
    2. VI.2 Cell theory
    3. VI.3 Classical density functional theory
  7. VII Bulk equilibrium hard spheres
    1. VII.1 Equilibrium phase behavior of monodisperse hard spheres
    2. VII.2 Impact of polydispersity on the equilibrium phase diagram
    3. VII.3 Equation of state
    4. VII.4 Structure and dynamics of the hard-sphere fluid
    5. VII.5 Bulk hard spheres in two dimensions
  8. VIII Interfaces in hard spheres
    1. VIII.1 Fluid in contact with a wall
    2. VIII.2 Fluid-solid interface
    3. VIII.3 Grain boundaries
  9. IX Binary hard-sphere mixtures
    1. IX.1 Crystal regime
    2. IX.2 Binary fluids: Structural crossover, demixing and dynamics
    3. IX.3 Sedimented monolayers of binary hard spheres
  10. X Confinement
    1. X.1 Quasi-2d confinement
    2. X.2 Cylindrical confinement
    3. X.3 Quasi-1d confinement
    4. X.4 Spherical confinement
    5. X.5 Flexible confinement
  11. XI Out-of-equilibrium phenomena
    1. XI.1 Non-equilibrium sedimentation in hard-sphere colloids
    2. XI.2 Rheology, flow, and shear-induced order
    3. XI.3 Microrheology
    4. XI.4 Other out-of-equilibrium phenomena
  12. XII Hard-sphere glasses and their formation
    1. XII.1 Historical theoretical developments and persistent challenges
    2. XII.2 Reciprocal-space picture: early studies of the hard-sphere glass transition
    3. XII.3 Real-space picture: Local structure and dynamic heterogeneity
    4. XII.4 Deeper supercooling: Beyond the mode-coupling (dynamical) crossover
    5. XII.5 The hard-sphere glass transition under confinement
    6. XII.6 The glass transition in 2d hard spheres
    7. XII.7 Vibrational properties of hard-sphere glasses
    8. XII.8 The jamming transition and its influence
    9. XII.9 Aging in hard-sphere glasses
    10. XII.10 Hard-sphere glasses far from equilibrium
  13. XIII Nucleation and growth
    1. XIII.1 A primer on classical nucleation theory
    2. XIII.2 Early nucleation experiments
    3. XIII.3 Light scattering versus real-space experiments
    4. XIII.4 The nucleation rate discrepancy
    5. XIII.5 Homogeneous nucleation in binary hard-sphere mixtures
    6. XIII.6 Heterogeneous and seeded nucleation
    7. XIII.7 Crystal growth
    8. XIII.8 Challenges to the fundamental assumptions of CNT
    9. XIII.9 Beyond the CNT regime: Crystallization at high volume fractions
    10. XIII.10 Nucleation in external fields
    11. XIII.11 Melting
  14. XIV Summary and Outlook
    1. XIV.1 What have we learned?
    2. XIV.2 Open challenges

##  I Introduction

Polymers, liquid crystals, surfactants, and colloidal dispersions are key
material pillars of soft matter. While physicists have long focused on the
study of the first three, colloids did not receive comparable attention until
relatively recently. Only thanks to significant advancements in the synthesis
of colloidal particles in the second half of the 20th century – motivated by
both industrial interest Vincent (2012); Tadros (2018); Vanderhoff _et al._
(1956) and scientific curiosity Fijnaut _et al._ (1978) – did the field come
to light. The potential of hard–sphere–like colloidal dispersions as physical
models of “simple liquids” exhibiting a fluid-crystal phase transition Poon
_et al._ (1996); Frenkel (2002); Evans _et al._ (2019) was subsequently
developed by Hachisu and co–workers Hachisu and Kobayashi (1974); Kose and
Hachisu (1974) and by Vrij _et al._ (1978). Hard-sphere particles, which are
impenetrable at distances less than their diameter $\sigma$ and otherwise
exhibit no interaction (see Fig. 1), had long been the object of theoretical
enquiry, but had until then largely remained but an experimental fantasy (see
Sec. II).

During the 1960s, ’70s and ’80s, advancements in colloid synthesis in industry
Barrett (1973); Walbridge and Waters (1966) and by pioneering colloid
scientist Ron Ottewill Cairns _et al._ (1976); Antl _et al._ (1986), combined
with the application of light scattering by soft matter physicists Peter Pusey
and Bill van Megen, led to the development of a well-controlled colloidal
model system (Fig. 2) van Megen _et al._ (1985); Pusey and Van Megen (1986);
van Megen _et al._ (1987). Pusey and Van Megen then not only convincingly
demonstrated hard-sphere freezing (Fig. 3), thus experimentally validating
theoretical Kirkwood and Monroe (1940) and computer simulation Alder and
Wainwright (1957); Wood and Jacobson (1957) predictions from a generation
prior, but also—and perhaps more importantly—elegantly demonstrated the
epistemic potential of colloidal suspensions.

Figure 1: Hard-sphere pair interaction potential $u_{\mathrm{hs}}(r)$ as a
function of the center-to-center distance $r$ between particles of diameter
$\sigma$.  Figure 2: The colloidal hard spheres observed by Pusey and Van
Megen are depicted (a) immediately, (b) one, and (c) four days after mixing.
After four days, the system is presumed to have completed its phase
separation. Volume fraction increases from left to right, with effective
values (determined by reference to the phase behavior)
$\phi=0.491,0.517,0.525,0.542,0.568,0.593,0.611,0.637$ and $0.654$ for samples
2–10 respectively. The samples range from a fluid phase to a fluid coexisting
with an iridescent crystal; at slightly larger volume fractions the whole
sample is crystalline while at yet higher volume fractions glassy amorphous
states are encountered. These states initially “coexist” with the crystal
until the entire sample ultimately becomes glassy. Reproduced from Pusey _et
al._ (2009) with permission. Figure 3: Equation of state and phase diagram of
hard spheres. The pressure $\beta P\sigma^{3}$ as a function of volume
fraction $\phi$ (solid blue) is (approximately) given by the Carnahan-Starling
expression Carnahan and Starling (1969) for the fluid and by that of Hall Hall
(1972) for the crystal. The metastable extension of the fluid branch (dashed
blue) is where slow dynamics and the hard–sphere supercooled liquids and glass
transition can be found, see Sec XII. Indicated on the $x$–axis are the fluid
and crystal phases, along with freezing volume fraction of about
$\phi_{f}=0.492$ and melting volume fraction $\phi_{m}=0.543$ and phase
coexistence at pressure $\beta P\sigma^{3}=11.6$ (see Sec. VII.1). Also
indicated are random close packing $\phi_{\mathrm{rcp}}\approx 0.64$ and
crystal close packing $\phi_{\mathrm{cp}}=\pi/(3\sqrt{2})=0.740\ldots$.
Approximate state points of the confocal microscopy images are denoted by grey
arrows.

Underlying the newfound physical interest is the fact that colloidal
particles—like atoms and molecules—exhibit thermal motion which allows them to
explore configurational space, and self-assemble into different phases, such
as _colloidal_ crystals, liquids, and gases (Fig. 3). Of the plethora of
possible colloidal systems, hard spheres have naturally emerged as the
benchmark. Despite the challenge of synthesizing perfectly hard colloids,
(micron-sized) hard spheres have four key strengths: (_i_) they have a single
control parameter, the volume fraction $\phi$; (_ii_) they are _classical_ in
nature, thus enabling accurate comparison with a broad array of theoretical
predictions and large–scale computer simulation; (_iii_) they diffuse their
own radius on the order of seconds, thus making their dynamics readily
accessible in experiment; (_iv_) their size makes them amenable to optical
techniques such as light scattering and confocal microscopy, thus enabling
accurate measurement of spatial and dynamical correlations and even particle
positions, without the need for large facilities such a synchotrons. Colloidal
hard spheres are therefore prized as model systems. Quantitative tests and
validations with theory nevertheless remain challenging. In this review, we
specifically take a critical look at what has been achieved and what
challenges remain to be faced in using colloids to verify theory. We also
emphasize the experimental observations that have yet to be given an accurate
theoretical description.

In order to manage the scope of this review, we largely restrain our
consideration to experiments with colloids. Only where exceptionally relevant
do we mention work with nanoparticles (see Boles _et al._ (2016) for a
dedicated review) or with granular matter (for which reviews are available in
Liu and Nagel (2010); van Hecke (2010); Charbonneau _et al._ (2017); Torquato
and Stillinger (2010); Arceri _et al._ (2022); Forterre and Pouliquen (2008);
Bi _et al._ (2015)). Where possible, we have referenced relevant review papers
from these other fields, but we humbly beg for grace from the reader regarding
important material that we may have missed. We nevertheless hope that our
review conveys that similar behavior can be observed in very different
systems, and that the transfer of ideas between these fields can be a fruitful
source of future research objectives.

Even after taking these thematic restrictions into account, the table of
contents hints at the wide reach of the remaining scope. Because at first
sight going through this whole review may seem daunting, we note that Secs.
II-VI deal with the historical perspective and methodology, with the results
from hard sphere colloids being from Sec. VII.1 onwards. We also emphasize a
few highlights:

  * •

phase behavior—colloidal experiments confirm the entropy-driven fluid-crystal
transition in hard-sphere systems (see Sec. VII);

  * •

fluid structure—experimental measures validate pair structure estimates and
track the development of even higher–order structures (see Sec. VII);

  * •

interfaces—the interfacial free energy underlies the barrier to nucleation,
and grain boundaries in hard-sphere systems form a fundamental model for a key
failure mechanism in crystalline materials (see Sec. VIII);

  * •

binary mixtures—colloidal suspensions have fueled the exploration of the rich
phase behavior of multi–component systems (see Sec. IX);

  * •

confinement—the stabilization of a wealth of different structural arrangements
of hard spheres has been achieved in both experiments and simulations (see
Sec. X);

  * •

far-from-equilibrium behavior—hard spheres have been instrumental in
disentangling hydrodynamic interactions from distortions of the local
structure to, _inter alia_ , shed light onto the phenomenon of shear
thickening (see Sec. XI);

  * •

glass transition—colloidal hard spheres have markedly advanced the
understanding of glasses, thanks to their simplicity and the possibility to
image them in real space, thus providing access to a host of properties that
are hard to access in molecular systems (see Sec. XII);

  * •

nucleation—the ability to directly image hard spheres enabled the first direct
observation of a critical nucleus, which makes it possible to test more
stringently the approximations that underlie classical nucleation theory (see
Sec. XIII).

We trust that this selection will motivate at least a few of the more hesitant
readers to continue through.

Before embarking on the bulk of this review, a few essential quantities must
be defined. We have already encountered the first one: the hard-sphere
diameter $\sigma$ (see Fig. 1). Although in theory it is a well defined
quantity, in practice it might refer to the effective diameter, the
hydrodynamic diameter, or the mean diameter in a suspension of slightly size
polydisperse colloidal particles, see Sec. III.4. The second has also been
alluded to above: the number density of $N$ colloids within a volume $V$,
i.e., $\rho=N/V$. More often, we will employ the dimensionless volume (or
packing) fraction $\phi=Nv_{1}/V=\pi\sigma^{3}\rho/6$ with
$v_{1}=\pi\sigma^{3}/6$ the volume of a single three-dimensional (3d) sphere.
The conjugate variable to density is the pressure $P$, our third quantity. It
is often beneficial to consider the dimensionless _reduced pressure_ $Z=\beta
P/\rho$ (also known as the compressibility factor), which combines pressure
with density and the inverse temperature $\beta=(k_{\text{B}}T)^{-1}$ where
$k_{\text{B}}$ is Boltzmann’s constant and $T$ is temperature. Note that i the
case of colloidal hard spheres, the relevant quantity is the _osmotic_
pressure. The final quantity is the (Brownian) time taken by a single sphere
to diffuse by its own radius

$\tau_{\text{B}}=\frac{(\sigma/2)^{2}}{D_{0}}=\frac{3\pi\eta\sigma^{3}}{4k_{\text{B}}T}$
(1)

where $D_{0}$ is the (bare) diffusion coefficient and $\eta$ is the solvent
viscosity.

As Fig. 1 makes clear, the potential energy in a hard sphere system is zero.
This means that only entropy contributes to the phase behavior. Furthermore,
provided the temperature is non-zero, ie the system is thermal, temperature
plays no role in the behavior of the system beyond scaling the timescale as
shown in Eq. 1. It is this independence of temperature that leads to density
or volume fraction or pressure being the sole control parameter for hard
sphere.

##  II Historical background

Given the extended and intricate history of hard spheres before they became an
object of study experimentally, a _longue-durée_ overview of the topic is in
order. This section summarizes the distinct histories of hard spheres as a
model of matter and as colloidal particles that mimick that model.

###  II.1 History of hard-sphere models

Describing atoms as hard elastic spheres – akin to billiard balls – finds its
origin in the early days of the kinetic theory of gases Brush (2003). The laws
of impact derived in the 17th century underlie the first such theory, which
was formulated by Bernoulli in the early 18th century. That approach, however,
was deemed disputable at the time and did not have much immediate scientific
impact. Further formalization by Maxwell in the mid-1800s, which built
explicitly on “systems of particles acting on each other only by impact”, was
more kindly received, and gave rise to a large body of work on kinetic
theories Brush (2003). In the words of science historian Stephen Brush, this
_exercise in mechanics_ “not only helped establish the theory, but laid the
foundations for modern statistical mechanics”.

It is therefore unsurprising that Van der Waals’ attempt at a microscopic
model of condensation, a couple of decades later, used a model effectively
consisting of _hard spheres_ with weak, fairly short-range attractive forces
Nairn and Kilpatrick (1972) (Brush, 1976, p. 407)111We now know the Van der
Waals equation of state to be exact only for a one-dimensional system with
infinitely long range and weak attraction Kac _et al._ (1963); Hansen and
McDonald (2013); Niss (2018).. Interestingly, Van der Waals’ description gave
the correct second virial coefficient, $B_{2}$, for the reference model Brush
(1976). Virial coefficients, which characterize how the equation of state for
the pressure $P$ as a function of density $\rho$ deviates from the ideal gas
reference, i.e.,

$Z=\frac{\beta P}{\rho}=1+\sum_{i=1}B_{i+1}\rho^{i},$ (2)

quickly became an important tool to characterize models of intermolecular
interactions. (An example is shown in Fig. 3. See Secs. VI.1 and VII.3 for
further discussion.) For hard spheres, $B_{3}$ and $B_{4}$, were hence
determined by Boltzmann before the end of the 19th century Nairn and
Kilpatrick (1972). Beyond that point, however, the algebra is intractable. It
therefore took more than fifty years and the advent of electronic computers
for estimates of higher-order coefficients to be obtained. The numerical
determination of $B_{5}$ by Rosenbluth and Rosenbluth (1954) renewed the quest
for accurate estimates, which currently extend up to $B_{12}$ Clisby and McCoy
(2006); Wheatley (2013).

Hard spheres also formed the basis for studying liquid _structure_. Early
experiments with lead shot in the 1920s Smith _et al._ (1929) and with hard
gelatin balls in the 1930s Morrell and Hildebrand (1936) were followed after
the war by Bernal’s extensive work with ball bearings Bernal and Mason (1960);
Finney (2013) and by more controlled theoretical descriptions, as further
discussed in Sec. XII.8. Hard spheres also contributed to establishing the
scaled particle theory of liquids Reiss _et al._ (1959), and were considered
by Wertheim as a solution of the Percus-Yevick approximation to the
Ornstein–Zernike equation that formally describes the liquid structure
Wertheim (1963). Hard spheres were subsequently a natural object for other
integral equation approaches Hansen and McDonald (2013); Rowlinson (2005) (see
also Sec. VI.1).

A parallel research effort into hard spheres stemmed from Kirkwood’s
speculation Kirkwood (1939), which he quickly substantiated with his postdoc
Elizabeth Monroe Kirkwood and Monroe (1940, 1941); Kirkwood (1951), that a
hard-sphere fluid becomes unstable to crystallization at high density
Charbonneau (2022); Hoddeson _et al._ (1992) (see Fig. 3). The somewhat
surprising proposal that the ordered phase could be entropically more
favorable than the dense (and disordered) liquid was independently formulated
by Fisher (1955) a decade later, and supported by early simulations using
molecular dynamics by Alder and Wainwright (1957) and Monte Carlo sampling by
Wood and Jacobson (1957)222It is important to note that the contributions of
coders for many of these works were not acknowledged through authorship
Battimelli _et al._ (2020). The author list of these works therefore does not
fully reflect the full scope of intellectual contributions that went into
their realization.. The ensuing debate about the physicality (and eventual
validation) of such a phenomenon led to various computational advances, as was
recently reviewed by Battimelli _et al._ (2020).

The study of the crystal state of hard spheres also has a long history. (That
of disordered jammed packings of hard spheres is somewhat more recent; see
Sec. XII.8.) The Kepler conjecture for their densest packing – their infinite-
pressure state with a packing fraction $\phi_{\text{cp}}=\pi/(3\sqrt{2})$ –
dates back to the 17th century and was proven by Hales only relatively
recently Aste and Weaire (2008); Hales _et al._ (2017). Finite-pressure
descriptions of hard-sphere crystals are less ancient, but have reasonably
rich antecedents as well. The familiar description in terms of free volume,
known as cell theory (Sec. VI.2) takes root in various attempts to
characterize the _liquid state_ in the 1930s Rowlinson (2005). The free-volume
description only became a useful reference for hard-sphere crystals in the
1950s Wood (1952); Salsburg and Wood (1962) and persists to this day in the
pedagogical literature (see, e.g., Barrat and Hansen (2003); Kamien (2007)).
However, because the associated corrections are less well controlled and are
computationally harder to evaluate than for the virial series (see Charbonneau
_et al._ (2021b) and references therein), cell theory has not been as
quantitatively influential.

###  II.2 The foretelling: How PMMA colloids came to be

Despite the extended theoretical importance of hard spheres, the system long
remained but a distant abstraction as far as experiments were concerned. It is
only through a gradual increase in the control over colloidal experiments that
the state of affairs has changed. In this section, we focus specifically on
the developments that have led to the sterically stablized poly-methyl
methacrylate (PMMA) system. Although other systems have been utilized to model
hard spheres Russel _et al._ (1989); Hunter (2001); Royall _et al._ (2013)
(see Sec. III.2), PMMA-based colloids have predominated since their
introduction. During the first 10-15 years of its use, 333Specifically before
the density matching of larger PMMA particles for confocal microscopy , see
Sec. III.3 around the millenium. and as the field burgeoned,any particle
softness is attributed to the steric stabilization layer, which is much
shorter than the particle diameter. Because PMMA particles have a bulk modulus
that is orders of magnitude larger, the tiny stablization layer only slightly
perturbs the system away from true hard sphere behavior (see Sec. III.3).

The origins of modern colloid science may be traced back to the latter part of
the nineteenth century with the emergence of physical chemistry as a distinct
discipline Vincent (2012, 2018), and early discoveries, such as the first
synthetic suspension of (optically active) gold nanoparticles by Faraday
(1859) and the use of colloidal sedimentation profiles to determine
Boltzmann’s constant by Perrin Perrin (1913). A major challenge to the use and
application of colloidal suspensions is that they typically aggregate
irreversibly when concentrated. Understanding of this phenomenon took a great
leap forward with the independent development—during the Second World War—by
Derjaguin and Landau Derjaguin and Landau (1941) in the Soviet Union and by
Verwey and Overbeek Verwey and Overbeek (1948) in the Netherlands of the DLVO
theory for the interactions between charged colloids. DLVO theory notably
provided a first-principles explanation for the long-noted _salting out_
phenomenon for aqueous suspensions of both colloids and proteins, in which
salt addition leads to aggregation and precipitation.

Colloid synthesis further progressed in the postwar years Vanderhoff _et al._
(1956); Ihler (1979) not least as certain industrial products, especially
paints and coatings, called for novel stabilization schemes Barrett (1973). It
had long been known that addition of (colloidal) spheres to a solution
increases viscosity Batchelor (1967), hence when in the 1960s control was
sought over the rheology of paints and coatings, adding colloids seemed like a
natural solution. However, these largely oil–based products —this being long
before legislation promoted the use of water–based coatings which release
fewer volatile organic compounds as they dry Derksen _et al._ (1995)—called
for stabilizing colloids in non–aqueous solvents. In these systems,
electrostatic charging was assumed to be negligible and therefore incapable of
stabilizing colloidal dispersions. (It was not until the turn of the third
millennium that electrostatic charging in non–aqueous solvents would be
investigated thoroughly (see Sec. III.3) Yethiraj and Van Blaaderen (2003);
Roberts _et al._ (2007); Royall _et al._ (2013, 2021).) Another method for
stabilization was therefore sought out.

Inspiration came from the burgeoning developments in polymer physics, notably
from the work of Flory de Gennes (1979). The loss of configurational entropy
of polymers grafted to surfaces was proposed as a colloidal stablization
mechanism. Early work included the stablization of carbon black through
physisorption of aromatic hydrocarbons Van Der Waarden (1950); Overbeek
(1966). In the 1960s, steric stabilization was first combined with a
(relatively) monodisperse nonaqueous colloidal system at Imperial Chemical
Industries (ICI), a major paint and coatings manufacturer of the time, by
coating polymethyl methacrylate (PMMA) particles with polyhydroxy steric acid
(PHSA) Walbridge and Waters (1966); Barrett (1973). In the UK, at that time
leading academic research in colloids took place in Cambridge University,
University College London and the University of Bristol Vincent (2012). In
Bristol, in particular, Ottewill had strong links with ICI from where the
sterially stablized PMMA was initially obtained. Once synthesized in Bristol
Cairns _et al._ (1976), these particles played a central role in the hard-
sphere story. They not only closely approximated the model of interest, but
they were amenable to solvents which match their refractive index and mass
density. As we will see, the former is essential for light scattering (Sec.
IV.1), and both are necessary for confocal microscopy (Sec. IV.2).

In 1974, Ottewill met Pusey who was then based at Malvern, close to Bristol,
and the two became interested in exploring light scattering as a means to
study strongly interacting colloidal systems Pusey _et al._ (1974); Brown _et
al._ (1975) as a means to study the PMMA system Brown _et al._ (1975). With
the arrival in Bristol of Van Megen from the Royal Melbourne Institute of
Technology and through a series of extended back-and-forth visits, Ottewill’s
ability to synthesize the hard-sphere-like PMMA and Pusey’s interest in these
system naturally led to the series of seminal papers that gave rise to Fig. 2,
see e.g., Pusey and Van Megen (1986); Pusey and van Megen (1987); van Megen
_et al._ (1985, 1987).

##  III Realizing hard-sphere systems

In this section, we briefly discuss the main experimental methods that have
been brought to bear in the study of colloidal hard spheres more specifically.
For a more general treatment of colloidal synthesis, we refer the reader to
appropriate reviews Russel _et al._ (1989); Hunter (2001); Ihler (1979). It is
worth emphasizing that, in experiments, truly hard spheres do not exist Poon
_et al._ (2012); Royall _et al._ (2013). In order to understand how far real
experimental systems deviate from hard spheres, different mappings to hard
spheres are discussed in Sec. III.4.

###  III.1 Experimental methods to prepare hard-sphere-like systems

Experimental techniques for colloidal hard spheres into three main categories:
synthesis, characterization and observation. The first has a chemical flavor
and is typically carried out by chemists, or by specially-trained physicists.
The second and third, by contrast, are more typically practised by physicists
(and are discussed in Sec. IV). These tasks can be conducted independently, as
it is common to buy colloidal particles (especially silica and polystyrene)
from commercial suppliers. However, it is worth noting that there is no
commercial supplier for the (canonical) sterically stabilized PMMA suspensions
discussed in Sec. II.2 444Worse, the PHSA sythesis has proven to be
challenging to challenging to reproduce, leading to very polydisperse, or
colloidally unstable, PMMA particles. At present, the global stocks of PHSA
are less than 10kg.. (A short review of PMMA synthesis is also given in
Dullens (2006).) Here, we mainly focus on the characterization and observation
of colloidal hard-sphere suspensions.

Aside from identifying the “hardness” of the particles (see Sec. III.3), a key
property of these suspensions is their polydispersity $s$. This is typically
expressed as the standard deviation of the particle diameter normalized by the
mean diameter expressed as a percentage. For PMMA synthesis, a “good” level of
polydispersity is considered to be $s\lesssim 5\%$ such that the system
crystallizes as readily as in the monodisperse case. (In Sec. VII.2 we will
report a detailed discussion of the impact of polydispersity on the phase
behavior.) Colloids produced from other materials, such as polystyrene, can
exhibit significantly lower polydispersity, down to $s\approx 1\%$ Russel _et
al._ (1989); Vanderhoff _et al._ (1956).

It is worth emphasizing that polydispersity distributions in experimental
systems cannot be prescribed. While in numerical simulations one may easily
change from Gaussian, to top hat or lognormal distributions, in experiment the
distribution is at the mercy of the particular synthesis. Even bimodal
distributions can sometimes arise (due to secondary nucleation events, in
which a second wave of colloids form during the reaction) Kawaguchi and Ito
(2005).

###  III.2 Synthesizing hard-sphere-like colloidal particles

The method adopted for PMMA synthesis was _dispersion polymerisation_
Kawaguchi and Ito (2005); Barrett (1973). In this approach, the methyl-
methacrylate monomers and methacrylic-acid are dispersed in a suitable solvent
(typically an alkane). Under the action of an _initiator_ , polymerisation
begins at a somewhat elevated temperature (around 80∘C) and in principle
proceeds until the monomer is exhausted. Relatively early in the
polymerization process, the initially transparent solution acquires a blueish
tint, due to the growing polymers falling out of solution and condensing into
colloids that scatter predominately blue light. Polyhydroxy-stearic acid
(PHSA) is added to the synthesis and is typically physisorbed to provide
steric stabilization. In this way, the strerically stabilized polymethyl
methacrylate (PMMA) system was born Walbridge and Waters (1966); Barrett
(1973); Cairns _et al._ (1976); Antl _et al._ (1986). The size of the
particles can be controlled with the monomer concentration, enabling for
example binary systems with a given size ratio to be synthesized, or a system
whose size optimized density matching (easier with small particles) against
higher-quality imaging (better with large particles) Bosma _et al._ (2002);
Poon _et al._ (2012).

Developments in synthesis include (_i_) partial control over the particle size
by the monomer concentration Poon _et al._ (2012); Bosma _et al._ (2002),
(_ii_) locking the stabilizer, such that it is covalently bonded to the PMMA
polymer backbone, and (_iii_) crosslinking the PMMA chain. The latter prevents
particles swelling and dissolving in a good solvent. Using such crosslinked
particles as a core is important for confocal microscopy because it enables a
PMMA _shell_ to be grown around the particle Dullens _et al._ (2003); Dullens
(2006)555The particles used by Kose and Hachisu (1974) had little to do with
modern PMMA particles, except their chemical nature. In particular, they were
not sterically stabilized.: if the initial step (the _core_) is labeled with a
fluorescent dye and the shell is unlabeled, then tracking particle locations
with confocal microscopy is much easier, thanks to the clean separation
between different fluorescent cores (see Sec. IV.2) van Blaaderen _et al._
(1992); Ivlev _et al._ (2012).

While many key early hard-sphere experiments were carried out with sterically
stabilized PMMA, and indeed the majority of the work we shall consider used
PMMA colloids, it is far from being the only material of interest. Silica may
also be stabilized in a non–aqueous solvent van Helden and Vrij (1980).
Aqueous systems of polystyrene and silica colloids form reasonable
approximations to hard spheres, as do microgel particles Royall _et al._
(2013) and emulsion droplets Dong _et al._ (2022). Because the diameter of
microgel particles can be tuned _in-situ_ , the effective volume fraction to
be changed at will, and therefore a single sample can access multiple state
points Yunker _et al._ (2014).

Two promising systems have been developed more recently. One involves the use
of an aqueous copolymer of fluorinated methacrylate and methacrylate with a
tunable density and refractive index, as demonstrated by Kodger _et al._
(2015). The other involves the use of 3-(trimethoxysilyl)propyl methacrylate
(TPM) particles Liu _et al._ (2016, 2019). Impressive results also have been
obtained with nanoparticles, whose small size enables their faster diffusion
than for the larger colloids and their self–assembly into larger structures
Boles _et al._ (2016) (see Sec. IV.3). Early work in the same size range used
microemulsions, which are thermodynamically stable droplets of _nanometer_
dimensions (despite their name).

###  III.3 Interactions in real hard-sphere-like systems

Central to identifying the state point of a hard-sphere-like system is
determining the (effective) volume fraction $\phi_{\mathrm{eff}}$, as reviewed
in Poon _et al._ (2012); Royall _et al._ (2013). The notion of effective hard-
sphere diameter will be properly discussed in Sec. III.4, but we here first
discuss interactions commonly used in experiments to mimic hard spheres. We
consider two sizes, $\sigma=200$ nm and $\sigma=2000$ nm, which roughly
reflects the range of colloid sizes used in the work discussed here. We note
the impact of the colloid size upon the system dynamics (Sec. IV.3) and their
applicability for particular experimental techniques (see Sec. IV).

While we have discussed the origins of some hard-sphere experimental systems,
we should assess how close to hard spheres these truly are Royall _et al._
(2013); Poon _et al._ (2012). Here, we provide a summary for the purposes of
this review. Practical hard-sphere–like systems fall into three broad
categories, as illustrated in Fig. 4. Colloidal stabilisation can be achieved
either via (a) steric stabilization or (b) charge stabilization. Because (c)
microgel particles consist of densely crosslinked polymers in a swelling
solvent, they do not explicitly need stabilization Schneider _et al._ (2017);
Lyon and Fernandez-Nieves (2012). The hardness of microgel colloids depends on
the quality of the solvent and the degree of crosslinking and can be density
dependent Royall _et al._ (2013). Emulsion droplets are an intermediate case.
They use a molecularly thin layer of surfactant for stabilization. Although
these particles are liquid, even a _microscopic_ ($\lesssim 1$nm2) change in
surface area leads to an interfacial free energy cost of order
$k_{\text{B}}T$. The softness of mesoscopic-scale droplets is, therefore, very
limited Dong _et al._ (2022); Morse and Witten (1993); Lacasse _et al._
(1996).

Figure 4: Schematic representation of various models for hard-sphere colloids.
(a) A sterically-stabilized particle has surface “hairs” (not to scale) of
average thickness $\bar{\delta}$, resulting in core-stabilizer diameter
$\sigma_{\mathrm{cs}}=\sigma_{c}+2\bar{\delta}$. (b) A charged colloid has an
electrical double layer (shaded) that gives rise to an effective diameter
$\sigma_{\mathrm{eff}}$. (c) A microgel particle is a heavily cross-linked
polymer. Figure reproduced from Royall _et al._ (2013).

_Softness due to stabilization._ – A sterically stabilized particle is shown
schematically in Fig. 4(a). One of us Royall _et al._ (2013) has previously
reviewed ways to obtain an effective hard-sphere diameter for these particles
via mapping various hard-sphere properties (see Sec. III.4). A mapping is also
possible via direct measurement of colloidal interactions. This procedure has
been applied to poly-hydroxy-steric acid-stabilized PMMA particles, whose
interactions were measured with the surface force apparatus and were found to
be quite well described by an inverse power law potential with energy scale
$\epsilon_{\mathrm{ipl}}$ Bryant _et al._ (2002)666While an inverse power law
of course is long-ranged, unlike the steric stablization, here we follow the
work of Bryant _et al._ (2002).:

$u_{s}(r)\approx\epsilon_{\mathrm{ipl}}\left(\frac{\sigma}{r}\right)^{n},$ (3)

where $\sigma$ denotes the particle diameter and $r$ the interparticle
distance. The relative range of $u_{s}$ depends on the particle size. For
example, $n=170$ was determined for particles with a diameter of $\sigma=200$
nm. Likewise, the strength of the interactions also depends on the particle
size, with $u_{s}(\sigma)=146$ $k_{\text{B}}T$ was reported for $\sigma=200$
nm, as shown by the pink line in Fig. 5(a). For larger particles
($\sigma=2000$ nm), $n=1800$ and $u_{s}(\sigma)\approx 1800$ $k_{\text{B}}T$,
were determined Bryant _et al._ (2002) [Fig. 5(b), pink line]. The results of
Bryant _et al._ Bryant _et al._ (2002) quantify what is intuitively obvious,
namely, that for a fixed length of stabilizing “hairs”, larger particles are
relatively harder. Similar conclusions were obtained in the rheological study
of Mewis _et al._ (1989), which varied the thickness of the stabilizing layer.

Figure 5: Estimation of effective colloid-colloid interactions for various
hard-sphere-like scenarios for particle diameters (a) $\sigma=200$ nm and (b)
$\sigma=2000$ nm. Shown is steric stabilization (pink line), strong
electrostatic interaction (solid blue line) and weak electrostatic interaction
(dashed blue line). See text for additional details. Figure reproduced from
Royall _et al._ (2013).

_Electrostatics._ – It is now generally accepted that immersion of a colloid
in a liquid medium systematically gives rise to some degree of charging
Yethiraj and Van Blaaderen (2003); Leunissen (2006); Roberts _et al._ (2007);
Royall _et al._ (2013). This charge is a potential source of softness that
should always be considered in experiments. In determining the degree of
softness due to electrostatic interactions, the linearized Poisson-Boltzmann
(PB) theory is incorporated into the Derjaguin-Landau-Verwey-Overbeek (DLVO)
theory Verwey and Overbeek (1948) to describe the interactions between charged
colloids. The original DLVO potential consists of Van der Waals (VdW) and
electrostatic components. That contribution often becomes very small due to
refractive index matching between colloids and solvent as in the case of light
scattering or (3d) microscopy. In the case of quasi-2d systems, or techniques
which do not require index matching such as rheology, electrostatic
stabilization ensures that particles do not come sufficiently close together
for VdW interactions to become important. Therefore, we shall henceforth
neglect VdW interactions except where explicitly stated.

We are often interested in systems in which sterically-stabilized particles
become charged. In this case, we can consider an inter-particle potential
consisting of a steric repulsion, $u_{s}(r)$, and an electrostatic
interaction, which in linearized Poisson–Boltzmann theory has a Yukawa form,
$u_{\mathrm{yuk}}(r)$:

$\displaystyle u(r)$ $\displaystyle=$ $\displaystyle
u_{s}(r)+u_{\mathrm{vdW}}(r)+u_{\mathrm{yuk}}(r),$ (4) $\displaystyle
u_{\mathrm{yuk}}(r)$ $\displaystyle=$
$\displaystyle\epsilon_{\mathrm{yuk}}\frac{\exp\left[-\kappa(r-\sigma)\right]}{r/\sigma},$
(5)

where $\kappa$=$\sqrt{4\pi\lambda_{B}\rho_{\mathrm{ion}}}$ is the inverse
Debye screening length with $\rho_{\mathrm{ion}}$ the number density of
monovalent ions, and $\epsilon_{\mathrm{yuk}}$ denotes the contact value of
the Yukawa interaction given by

$\beta\epsilon_{\mathrm{yuk}}=\frac{Z_{e}^{2}}{(1+\kappa\sigma/2)^{2}}\frac{\lambda_{B}}{\sigma},$
(6)

with $Z_{e}$ the number of electronic charges on the colloid. The Bjerrum
length

$\lambda_{B}=\frac{e^{2}}{4\pi\epsilon_{0}\epsilon_{r}k_{\text{B}}T},$ (7)

is the distance at which the interaction energy between two elementary charges
equals $k_{\text{B}}T$, $e$ is the electronic elementary charge,
$\epsilon_{0}$ the permittivity of the vacuum, and $\epsilon_{r}$ the relative
dielectric constant.

We first consider the case for aqueous systems in which the charging is rather
strong. In 3d, most work with aqueous systems has used rather small colloids.
We provide an estimation of the electrostatic interaction potential by
plotting the solid blue line in Fig. 5(a), where we evaluate the Yukawa
interaction Eq. (5) as follows. We assume that salt is added such that the
Debye length is 4 nm. It has been suggested that an upper bound for the
effective charge in many systems Alexander _et al._ (1984) can be described by
the rule of thumb $Z_{e}\lambda_{B}/\sigma\approx 6$. For higher values of the
colloid charge, the electric field at distances beyond the Stern layer reduces
to the electric field corresponding to the effective charge due to ion
condensation Alexander _et al._ (1984); Russel _et al._ (1989). With water as
the solvent, this leads to an (effective) colloid charge of $Z_{e}=1500$ for
the colloid diameter of interest, $\sigma=200$ nm. The resulting electrostatic
interaction (represented by the solid blue line in Fig. 5(a)) demonstrates a
noticeable degree of softness.

Sterically stabilized PMMA is used in non-aqueous solvents. For example, the
work of Pusey and Van Megen used a mixture of cis-decalin and carbon
disulphide with dielectric constant $\epsilon_{r}=2.64$ Pusey and Van Megen
(1986), which leads to a Bjerrum length $\lambda_{B}\approx 21$ nm. With ionic
sizes in the range of $\lesssim 1$ nm, one expects strong coupling between
oppositely charged ions and consequently very little dissociation of surface
groups. It was therefore a long-held assumption that electrostatics could be
safely neglected. However, more recent work has demonstrated that some
electrostatic charging is _always_ present Yethiraj and Van Blaaderen (2003);
Roberts _et al._ (2007); Klix _et al._ (2013); Leunissen (2006).

To evaluate electrostatic interactions in solvents with such low dielectric
constants, we again use the $Z_{e}\lambda_{B}/\sigma\approx 6$ criterion with
a Debye length of $\kappa^{-1}=5\mu$m, which is consistent with measurements
of the very small ionic strength found in such solvents Klix _et al._ (2013).
We can then evaluate Eq. (5) for particles with diameters of $\sigma=200$ nm
and $\sigma=2000$ nm, which are represented by blue dashed lines in Fig. 5(a)
and Fig. 5(b), respectively. We see that the electrostatic repulsion changes
very little over the chosen range. It is worth noting that for the smaller
particles ($\sigma=200$ nm), these parameters correspond to an electrostatic
charge of just $Z_{e}=2$. Although small, comparable values have been measured
in experiments Klix _et al._ (2013). For this case, it is highly probable that
the spherically symmetric DLVO approach would not be accurate. However, the
effects seem small enough that one can reasonably neglect electrostatic
charging for this combination of solvent dielectric constant and particle
size. Larger colloids, however, require density matching to suppress
sedimentation. Since two of the denser solvents used in experiment, namely
tetrachloroethylene (TCE) and carbon tetrachloride, can be readily absorbed by
PMMA, both the density and the refractive index of the particles change,
resulting in a turbid system characterized by substantial Van der Waals
interparticle attractions. The hard-sphere-like behavior of the particles is
then lost Ohtsuka _et al._ (2008); Royall _et al._ (2013). Very recent work
has nevertheless found excellent hard-sphere behavior using a solvent mixture
with TCE Kale _et al._ (2023).

Finally, we consider stronger electrostatic charging in the PMMA system often
applied in studies using confocal microscopy for larger colloids with
$\sigma=2000$ nm. The solvents used for the smaller PMMA particles used for
light scattering have a relative dielectric constant $\epsilon_{r}\approx 2$.
Relative to these, the density matching mixture of cyclo-hexyl bromide (CHB)
and _cis_ -decalin has a dielectric constant of $\epsilon_{r}=5.4$, which is
due to the CHB component with $\epsilon_{r}=7.9$. While non-aqueous, the
dielectric constant is nevertheless much higher than that mentioned in the
preceding paragraph and tetrabutyl ammonium bromide (or a similar salt) is
typically used to screen the electrostatic interactions Yethiraj and Van
Blaaderen (2003). CHB advantageously appears to be less aggressive towards the
PMMA particles than does TCE for example, with rather less absorption and
swelling. If we again use the same criterion $Z_{e}\lambda_{B}/\sigma\approx
6$ we arrive at $Z_{e}\approx 500$. The Debye length corresponding to a
saturated solution of this salt is around $\kappa^{-1}\approx 100$ nm
Leunissen (2006); Royall _et al._ (2006, 2003), Evaluating Eq. (5) then gives
the solid blue line in Fig. 5(b), which exhibits a very considerable degree of
softness. Therefore for the larger PMMA particles, which require density
matching, one can either accept some softness and add salt, or risk
attractions due to solvent absorption (in the case of TCE for example).

Imaging quasi-2d or strongly confined systems – obtained, for example, by
sedimenting particles onto a substrate Ivlev _et al._ (2012) – is somewhat
relieved from such drawbacks. For example, polystyrene colloids in an aqueous
solvent with a suitable amount of salt added are rather hard as the Debye
length of a few nm is very much less than the micron scale particle diameter
Royall _et al._ (2013). Particles of size 3-5 $\mu$m are readily imaged and
their gravitational length 777The gravitational length is a measure of
vertical movement due to thermal energy, see Sec. XI.1. can be $\approx
0.01\sigma$ so that out–of–plane motion is rather small Williams _et al._
(2013, 2015).

###  III.4 Mapping soft spheres to hard spheres

As discussed above, real colloids inevitably display some degree of softness.
It is therefore important to be able to map their behavior to that of perfect
hard spheres for the purpose of comparison with theory and simulation results.
Mapping entails finding an effective hard-sphere diameter $\sigma_{\rm eff}$,
such that one may translate the experimentally-controllable particle number
density $\rho$ to an effective hard-sphere volume fraction $\phi_{\rm eff}$
using $\phi_{\rm eff}=\pi\sigma_{\rm eff}^{3}\rho/6$. Two distinct approaches
can be used to determine $\sigma_{\rm eff}$: directly from the interaction
potential and indirectly from the observed behavior of the system. The latter
method rests on the assumption that the weight fraction (which may be
accurately determined) at freezing may be taken to correspond to the freezing
volume fraction. This provides a calibration of the volume fraction through an
accurately known quantity. Developments of this method to address the effects
of graviational settling are discussed in Sec. VII.

If the interaction potential is known, several analytical routes to estimate
an effective hard-sphere diameter have been proposed. Arguably the simplest
one is to set an effective hard-sphere diameter $\sigma_{\mathrm{kT}}$ such
that the inter-particle repulsive energy at this centre-to-centre separation
between two particles is equal to the thermal energy $k_{\text{B}}T$. A more
sophisticated approach entails taking into account the functional behavior of
the pair interaction as proposed by Barker and Henderson (BH) Barker and
Henderson (1967),

$\sigma_{\mathrm{BH}}=\intop_{0}^{\infty}dr\left[1-\exp\left(-\beta
u(r)\right)\right].$ (8)

An alternative mapping due to Andersen, Weeks, and Chandler Andersen _et al._
(1971) takes into account the structure of a hard-sphere reference fluid at
the same density, making the effective diameter density-dependent. This
approximation has shown to be effective for mapping the structural relaxation
time of soft spheres near the glass transition to the hard-sphere model
Schmiedeberg _et al._ (2011). For the commonly used Weeks-Chandler-Andersen
potential Weeks _et al._ (1971), both the Barker-Henderson and Andersen-Weeks-
Chandler mapping have been shown to provide excellent phase behavior
predictions at sufficiently low temperatures Attia _et al._ (2022); Dasgupta
_et al._ (2020). However, for higher temperatures – or equivalently, particles
that are less hard-sphere-like – significant deviations emerge. As a result,
care is needed when interpreting the results of this type of mapping in
colloidal systems Royall _et al._ (2013). The situation is illustrated for a
model system of charged colloids in Fig. 6, where the fluid-crystal phase
boundaries from simulations Hynninen and Dijkstra (2003) are mapped to hard
spheres using $\sigma_{\mathrm{kT}}$ and $\sigma_{\mathrm{BH}}$. Both
approaches show significant deviations from the true hard spheres.
Importantly, the fluid-crystal coexistence gap is also narrowed. Therefore, as
corroborated by simulation work showing that nucleation barriers continue to
be sensitive to the softness even at screening lengths on the order of a few
percent of the particle diameter de Jager and Filion (2022), simply obtaining
$\sigma_{\mathrm{eff}}$ may not suffice to confidently reproduce hard-sphere
behavior.

Figure 6: Mapping the phase behavior of weakly charged hard-core colloids to
pure hard spheres. The hypothetical particles have an electrostatic charge
$Z=500$, and the solvent has $\sigma=2\mu$m with $\kappa^{-1}=100$nm, which is
typical for non–aqueous solvents Royall _et al._ (2005). The fluid-solid
coexistence gap of the hypothetical charged-colloid system as a function of
packing fraction $\phi$ is taken from simulations Hynninen and Dijkstra (2003)
(red, “sim”), and is mapped to pure hard spheres assuming that the effective
hard-sphere diameter corresponds to the interaction potential $\beta
u(\sigma_{\mathrm{kT}})=1$ (lilac, $\sigma_{\mathrm{kT}}$), and using the
mapping of Eq. (8) (blue, $\sigma_{\mathrm{BH}}$). The phase diagram of pure
hard spheres (grey, HS) is shown for comparison.

The second option is to measure quantities that depend on volume fraction.
Structural quantities such as the radial distribution function may be compared
with simulation or theory (see Sec. VII.4) Royall _et al._ (2007b, 2018b);
Thorneywork _et al._ (2014); Kale _et al._ (2023). However, for dense hard
spheres, these pair correlations vary rather slowly with $\phi$.
Alternatively, one can choose a quantity which strongly depends on $\phi$,
such as the relaxation time at high volume fraction (see Fig. 36, Sec. XII) or
higher–order structural observables (see Fig. 14, Sec. VII.4) Royall _et al._
(2018b); Pinchaipat _et al._ (2017). However, it is important to note that
there is absolutely no guarantee that different methods give the same mapping.
For example, higher-order structure measurements are highly sensitive to
missing – or “ghost” particles (see Sec. IV.2), while measuring dynamics in
real space is more sensitive to errors in particle positions.

In short, nearly-hard spheres, whether realized in experiments or in
simulations, can closely approach the behavior of true hard spheres, but
subtle differences remain. Choosing an effective diameter is an essential step
in assigning an effective packing fraction to an experimental system, and
hence in comparing its behavior to hard-sphere simulations or theory. However,
no unambiguous definition of such an effective packing fraction exists, and
different approaches to obtain one result in different estimates.

###  III.5 Accuracy in hard-sphere experiments

The above discussion makes clear that, with care, an effective hard sphere
volume fraction may be arrived at for an experimental system. However, a
degree of uncertainty will persist, caused by (_i_) softness, due to
electrostatics, stabilization or other sources; (_ii_) attraction, due to Van
der Waals interactions or other sources; and (_iii_) polydispersity, which is
typically represented via the second moment of the distribution. (An unbiased
way to characterize the distribution remains to be found. Even the laborious
method of sizing every particle by electron microscopy likely samples from a
biased set of particles, due to inhomogeneities during evaporation in
preparing the sample.) These qualities pertain to the particles, but,
especially for non-aqueous systems, the parameters which might effect the
state point include the solvent choice, its age and purity as well as the
ambient temperature and humidity (which are often not controlled in colloid
laboratories).

Even if these uncertainties could somehow be mitigated, there remains a
fundamental aspect of colloidal systems that sets them apart from most
materials in that _model colloids are a synthetic material_. As noted above,
PMMA, for instance, cannot be purchased commercially; it needs to be
synthesized, typically in small batches of 10-1000g. Therefore, even if the
quantities above were somehow characterized for a given synthesis (which would
entail a very large amount of work), once that had been exhausted, the
operation would need to be carried out again for the next synthesis.

As a consequence, _colloidal systems are not as well–characterized as atomic
and molecular systems_. A container of argon of sufficient purity has
identical characteristics to another. However, like home-made cakes,
_colloidal hard spheres from two different syntheses will never have identical
characteristics_. That is to say, the concept of an accurate measurement of a
physical property fundamentally differs. For example, the triple point
temperature of argon is 83.8058 K, i.e. known to six significant figures,
drawing on many successive and increasingly accurate measurements of the same
material. Such an accuracy is essentially impossible to achieve in colloids.
This brings us to the question of just how much accuracy we can expect in
colloids. Poon _et al._ (2012) concluded that a _relative error_ in volume
fraction around the freezing transition of 6% might be reasonable. Our opinion
is that when careful comparisons are made, as described in Sec. III.4, an
absolute error of 1% in volume fraction seems achievable, albeit not always
realized. In the remainder of this review, we shall revisit the impact of
colloidal–level accuracy where appropriate.

A more far-reaching interpretation of accuracy in the context of this review
is to consider what is meant by _hard_ spheres. This matter is not as
straightforward as it might first seem, as whether a system is sufficiently
similar to hard spheres or not depends on the physical property considered.
For example, it has long been known that even quite soft repulsive particles
have a fluid structure similar to that of hard spheres Weeks _et al._ (1971).
By contrast, phase boundaries are highly sensitive to the interaction details.
For the purposes of this review, we can clearly state that systems that are
_not_ hard-sphere–like are (_i_) so soft that they form a bcc crystal, i.e.,
their phase behavior differs qualitatively from that of hard spheres; (_ii_)
attract each other significantly, by which we mean a well depth of 0.5
$k_{\text{B}}T$ or above; (_iii_) use temperature rather than volume fraction
as control parameter. Beyond these three criteria, the degree of hardness
required for our purposes is therefore context–dependent. The interactions of
the particles in the work that we discuss in this review then is sufficiently
hard for appropriate hard sphere like behavior to be observed in the context
in question. We discuss certain instances where the degree of hardness is
particularly relevant.

##  IV Measuring hard spheres _in vitro_

We now turn to techniques for analyzing colloidal hard spheres in experiments.
The configurations adopted by the particles and their dynamics are the main
quantities of interest. The experimental work that we review primarily used
direct observation through microscopy or scattering methods. (Rheological
approaches have been carefully reviewed elsewhere Wagner (2022); Larson
(1998); Jacob _et al._ (2019).) While some scattering studies of hard spheres
have used neutron Cebula _et al._ (1981); De Kruif _et al._ (1988) and small-
angle X-ray scattering Petukhov _et al._ (2002) and more exotic X-ray
scattering methods Wochner _et al._ (2009), most have used light scattering.
This last method is therefore the focus of our discussion.

###  IV.1 Light scattering

Characterization of materials via the scattering of electromagnetic radiation
dates back almost to the discovery of X-rays, whose wavelength is comparable
to length scales relevant to atomic and molecular systems. The larger length
scales of colloids correspond to the wavelength of visible light, which is
scattered by spatial and dynamical fluctuations in the refractive index of the
material. A typical light scattering set–up is illustrated in Fig. 7, by a
global view and by a close up of the scattering volume with the scattering
wavevector $\mathbf{k}$ defined in terms of the wavevector of the incident
$\mathbf{k}_{i}$ and scattered $\mathbf{k}_{s}$ light, together with the
scattering angle $\theta$.

Figure 7: Schematic of a light-scattering set-up. (a) Global view of the set-
up. In the expanded view, scattering though an angle $\theta$ is indicated in
a small volume $dV$. (b) Close up of the scattering volume with the scattering
wavevector $\mathbf{k}$ defined in terms of the wavevector of the incident
$\mathbf{k}_{i}$ and scattered $\mathbf{k}_{s}$ light, and also the scattering
angle $\theta$. Inspired by Pusey (2002).

The key observables in light scattering are the intensity of the scattered
light $I(\mathbf{k},t)$ and its time correlation. Here $\mathbf{k}$ is the
wavevector through which the light is scattered, $t$ is the time and
$\mathbf{r}$ is the particle position (Fig. 7). The scattered intensity is
then Pusey (2002)

$\displaystyle I(\mathbf{k},t)$
$\displaystyle=\frac{E_{0}^{2}}{r^{2}}\sum_{j=1}^{N}\sum_{k=1}^{N}b_{j}(\mathbf{k},t)b_{k}^{*}(\mathbf{k},t)\cdot$
(9)
$\displaystyle\exp\left[-i\mathbf{k}\cdot\left(\mathbf{r}_{j}(t)-\mathbf{r}_{k}(t)\right)\right],$

where $E_{0}$ is the amplitude of the incident light field, $b_{i}$ is the
scattering length of the $i$th particle and is closely related to the form
factor.

_Static light scattering._ – In the case of static light scattering (SLS)
temporal quantities are averaged over time. Then from Eq. (9) for a given
$b(\mathbf{k},t)$ (which can be readily calculated for spheres made of a
material with a fixed refractive index) we can obtain the structure factor
Pusey (2002)

$S(\mathbf{k})=\frac{1}{N}\sum_{j=1}^{N}\sum_{k=1}^{N}\langle\exp\left[-i\mathbf{k}\cdot\left(\mathbf{r}_{j}-\mathbf{r}_{k}\right)\right]\rangle,$
(10)

and therefore its Fourier transform, the radial distribution function $g(r)$.
This quantity also provides a link to the thermodynamics of the system, since
in the long-wavelength limit it is related to the isothermal (osmotic)
compressibility $\chi_{T}$

$\lim_{k\to 0}S(k)=\rho
k_{\text{B}}T\chi_{T}=k_{\text{B}}T\left(\frac{\partial\rho}{\partial
P}\right)_{T}.$ (11)

Light scattering has been applied to soft-matter systems since the 1930s Doty
and Edsall (1951). For example, Doty and Steiner (1952) studied solutions of
Bovine Serum Albumin (BSA) protein, and soon were followed by famed
crystallographer Rosalind Franklin and coworkers Klug _et al._ (1959) who
studied colloidal-type crystals of charged viruses. However, the approach only
really took off after 1962, when Spectra Physics _lasers_ became commercially
available. The intense light produced by these sources revolutionized the
technique by dramatically increasing the signal-to-noise ratio, making it much
more practical and widely available. For typical scattering angles $\theta$
and a HeNe laser wavelength of 633 nm, sub-micron sized colloids fall within
the ideal size range for light scattering. By the mid-1970s, the measurement
technology was ready for the emergence of experimental hard spheres from the
point of view of two-point spatial correlation functions.

Pusey was an early champion of light scattering, moving from the study of
critical phenomena in molecular liquids to the study of viruses as colloidal
particles Pusey _et al._ (1972). Because synthetic colloids held the promise
of being more robust than viruses, he arranged to receive polystyrene samples
from the Bristol colloid group via Ottewill and Vincent (see Sec. II.2). This
enabled his pioneering work on determining static structure factors from the
mid 1970s onward, starting with charged polystyrene colloids Brown _et al._
(1975), which complemented earlier work on colloidal crystals Hiltner and
Krieger (1969); Williams and Crandall (1974). When using light scattering for
denser suspensions (by which we mean denser than the dilute limit,
$\phi\gtrsim 0.01$), precise refractive index matching between colloids and
solvent is essential as noted in Secs. II.2 and III.1. For polystyrene,
however, this is challenging due to its high refractive index (relative to the
solvent). Early work, which explicitly matched the refractive index of the
colloids and solvent, was carried out by Vrij _et al._ (1983).

Analysis of light scattering data in dense suspensions can be complicated by
the occurrence of multiple scattering, whereby incident photons are scattered
more than once before exiting the sample. For both static and dynamic light
scattering, however, this issue can be circumvented by using scattering set-
ups that isolate the signal from singly scattered light by cross-correlating
the intensity of two different wave lengths of light Segre _et al._ (1995);
Schatzel (1991); Moussaïd and Pusey (1999).

_Dynamic light scattering._ – While the discussion above outlines early
efforts to extract structural information such as the static structure factor
from SLS, as Eq. (9) makes clear, determining time correlations in the
fluctuations in refractive index – in the scattered form of the so–called
_speckle pattern_ – can reveal dynamical behavior, i.e., _dynamic light
scattering_ (DLS). In the 1960s, the early pioneers in the field, Benedek and
coworkers Clark _et al._ (1970) and Cummins and Swinney (1970), used a
spectral technique to analyze the speckle pattern, which was experimentally
tedious. The now dominant technique, time correlation of detected photons
(frequently called “photon correlation spectroscopy”), was pioneered by Pike
and coworkers Foord _et al._ (1970).

The key quantity measured in DLS is the normalized time correlation function

$g^{(2)}(k,\tau)\equiv\frac{\langle
I(\mathbf{k},\tau)I(\mathbf{k},0)\rangle}{\langle I(\mathbf{k})\rangle^{2}}.$
(12)

Importantly, $g^{(2)}(k,\tau)$ is intimately related to a key dynamical
property, the intermediate scattering function $F(k,\tau)$, via the Siegert
relation

$F(k,\tau)=\sqrt{g^{(2)}(k,\tau)-1},$ (13)

which constituted the mainstay of dynamical measurements of hard-sphere
systems until dynamical analysis using confocal microscopy became available
around the turn of the millennium.

A significant advantage of DLS over confocal microscopy is the possibility to
access a broader range of time scales. Microscopy is limited to a relatively
small time window. Fast time scales ($<10^{-2}$ s) and very long ones
($>10^{4}$ s) are inaccessible, due to limitations in sample stability. By
contrast, DLS can probe much shorter time scales. (See also the discussion
concerning the consequences of particle size in Sec. IV.3.)

_More exotic light scattering: small angles, multiple length scales and
multiple color._ – While the essential methods of light scattering were
established by the mid–1970s, the technique has since enjoyed considerable
development. As Fig. 7 makes clear, going to very small angles (small-angle
light scattering) probes smaller wave vectors and larger real–space length
scales Cipelletti and Weitz (1999). Mounting multiple detectors, such as one
corresponding to a “standard” wave vector around $2\pi/\sigma$ and one at a
smaller angle $\ll 2\pi/\sigma$, then enables multiple length scales to be
probed simultaneously. This approach has since been implemented to great
effect Franke _et al._ (2014); Tamborini and Cipeletti (2012). Other
developments include two-color techniques which enable rather more turbid
samples ot be studied Segrè _et al._ (1995b). More recent advancements such as
multi-speckle correlation spectroscopy enable simultaneous measurement of
dynamical and spatial heterogeneities Golde _et al._ (2013, 2016).

###  IV.2 Optical microscopy and particle-resolved studies

Figure 8: The principle of confocal microscopy in the conventional _Epi–_
illumination mode. In _Epi–_ illumination, the objective and condenser lenses
are the same, and the dichroic mirror allows the illuminating beam to pass
through, while reflecting the returning beam to the detector. Light is focused
into a point in the sample by the condenser lens. The confocal pinhole rejects
all light except that from the point in focus (dotted lines show reflection of
light from out-of-focus regions). In this way, only one point in a 3d sample
is in focus and the incident beam is scanned relative to the sample to
generate a 3d image. _Epi–_ illumination records fluorescent light emitted by
the sample, so the distribution of fluorescent dye is imaged.

Unlike atomistic systems, colloids are amenable to _in-situ_ observation via
light microscopy techniques. Analysis of such real-space images of colloidal
hard-sphere systems dates back at least to the early 1970s, with the work of
Hachisu and coworkers Kose _et al._ (1973); Kose and Hachisu (1976); Okamoto
and Hachisu (1977), but antecedents can be found in the Bragg scattering study
of “colloidal” crystals of viruses of Franklin and coworkers Klug _et al._
(1959). Microscopy work may be broadly classified into two categories–one in
which individual particles are resolved (so–called particle–resolved studies
Ivlev _et al._ (2012)) and the other in which particles are too small to be
resolved.

In the case of particle-resolved studies of 3d systems, some kind of 3d
imaging is clearly important. Here, the development of _confocal_ microscopy
would prove to be transformative. A schematic of a confocal microscope is
shown in Fig. 8. Confocal microscopy was proposed in the mid-20th century with
Minsky developing an operational device in 1955 Minsky (1957). Like light
scattering, confocal microscopy was transformed by the advent of lasers. By
the mid-1980s, stage controlled systems, which scan the stage relative to the
beam, were commercially available from Oxford Optoelectronics, later acquired
by Biorad, the company which became synonymous with early confocal
microscopes.

Since the signal recorded in (_epi–_ illumination) confocal microscopy is the
distribution of fluorescent dye, some form of labeling of colloidal particles
(or solvent) is necessary. Novel fluorescent particles suitable for confocal
microscopy had been developed by Van Blaaderen _et al._ van Blaaderen _et al._
(1992). The additional synthetic step of growing an undyed shell around the
fluorescently labeled core contributed greatly to the ability to distinguish
particles (see Sec. III.2). To our knowledge, the first 3d imaging of
hard–sphere-like systems with confocal microscopy was that of Van Blaaderen
_et al._ van Blaaderen _et al._ (1992), who observed colloidal crystals and
glasses. The field, which has since grown markedly, has also been extensively
reviewed Ivlev _et al._ (2012); Hunter and Weeks (2012). Confocal microscopy
has since been developed to incorporate a shear stage as discussed in Sec. XI
and XII.10 Besseling _et al._ (2009); Lin _et al._ (2014b) and optical
tweezers (see below).

_Particle tracking._ – In a seminal paper, Crocker and Grier (1996) developed
an algorithm for tracking the coordinates of colloidal particles from
microscopy data. This approach interpreted bright pixels as potential particle
centres. In the Crocker and Grier method, the trial centers are subsequently
refined through a series of criteria, such as being the centre of an object
weighted by pixel intensity and being far enough apart from other objects, so
that other particles do not overlap. It is worth noting that the majority of
work on 2d hard-sphere systems has been carried out using microscopy, as such
systems are less amenable to scattering techniques.

While the algorithms introduced by Crocker and Grier for quasi-2d systems
Crocker and Grier (1996) have undoubtedly formed the backbone of most work
since, it is important to note that more sophisticated approaches are
available. Among the pioneers of the field, van Blaaderen and Wiltzius (1995)
sought to improve the tracking accuracy in the axial direction by fitting a
Gaussian to the integrated intensities in each plane that constituted a given
colloid. The early work examined hard-sphere glasses and crystals, i.e.,
colloidal solids in which diffusion could be neglected van Blaaderen and
Wiltzius (1995). Even though this method made imaging easier, the slow scan
rates needed to acquire images with low noise left dynamical information
beyond reach. For a quasi–2d system, time–resolved tracking was carried out by
the mid-1990s Marcus _et al._ (1996). In 3d systems, two breakthrough papers,
published simultaneously in the year 2000, obtained dynamical information by
tracking particle coordinates over time Weeks _et al._ (2000); Kegel and Van
Blaaderen (2000). This local information led to a breakthrough experimental
measurement of dynamical heterogeneity, a particularly important
characteristic of glass physics (see Sec. XII).

More recently, Jenkins and Egelhaaf (2008) pushed the limits of the technique
as they sought to identify contacts between colloids in a dense sediment
through ultra-high precision coordinate location. Such precision was achieved
by first determining an empirical image of a colloid, which could then be
compared to the original image. Further systematic improvements were made by
Gao and Kilfoil (2009). Other, more specialized, developments include that of
Lu _et al._ Lu _et al._ (2006), whose algorithmic advances for image
processing made it possible to follow experiments in real time, rather than
via post facto analysis. This distinction is particularly significant. The
spatial region of interest can often drift out of the field of view, which is
usually quite small given the large file sizes generated by time-resolved 3d
imaging, but real-time tracking allowed feedback to be applied to the stage,
such that the region of interest could be followed irrespective of drift.

Ever-present polydispersity is typically neglected, except in the simple case
of binary systems wherein two-color imaging enables the two species to be
differentiated Royall _et al._ (2007b). More recently, algorithms which
identify the size of _each_ particle have been developed Kurita and Weeks
(2012); Leocmach and Tanaka (2013). Ideas to extract further information than
possible using conventional techniques have also been implemented, such as
comparing simulated trial coordinates and experimental data Statt _et al._
(2016) as well as maximal precision methods Bierbaum _et al._ (2017). One
might imagine that machine-learning techniques could be applied to the vexing
problem of obtaining particle coordinates from confocal microscopy images, as
proposed by Bailey _et al._ (2022) for a dilute suspension of out-of-
equilibrium (active) colloids. However, this is a challenging problem in the
case of concentrated colloidal dispersions in 3d, not least because of the
paucity of the “ground truth” (other than synthetic data) against which a
neural network can be trained Kawafi (2023).

When comparing experimental particle-resolved data for hard spheres with
numerical simulations, the effects of polydispersity and particle tracking
errors should never be omitted. _High-quality_ data nevertheless implies an
error of $0.05\sigma$ in particle position and a few percent missing or
additional “ghost” particles Ivlev _et al._ (2012); Royall _et al._ (2007b).
This problem may be mitigated to some extent by mimicking such features into
computer simulation data Royall _et al._ (2007b); Pinchaipat _et al._ (2017),
and through other techniques, as extensively reviewed Ivlev _et al._ (2012);
Lu and Weitz (2013); Hunter and Weeks (2012).

_Optical tweezers._ – The development of lasers had yet another important
consequence for colloid experiments. A tightly focused laser beam enables
trapping of colloidal particles and their manipulation Ashkin (1970). These
optical tweezers have great potential to enable _in-situ_ observation of
phenomena which depend upon exquisite control of constituent particles Grier
(1997, 2003). Like microscopy, 3d implementation of optical tweezers carries
challenges, mainly related to the requirement of refractive index matching
colloids and solvent for imaging while simultaneously _mismatching_ refractive
indexes for tweezing. The first fully 3d system for colloids was introduced by
Vossen _et al._ (2004) and single-lens 3d tweezer methods have since been
developed Curran _et al._ (2014). It is important to note, however, that the
maximum number of particles which can be tweezed is about 40 and that tweezing
leads to a relatively weak confinement, of the order of a few $k_{\text{B}}T$
at most Williams _et al._ (2013).

###  IV.3 Size _does_ matter: colloid dynamics in different experiments

Light scattering and microscopy techniques have typically been used for rather
different sized colloidal particles. While of course the fundamental behavior
of the system remains unaltered, there are significant consequences concerning
the ability to access phenomena of interest. In particular, light scattering
can operate with particles in the 100s of nanometer size range, while
microscopy typically uses particles larger than a micron.

As Eq. (1) shows, the Brownian time scales as the volume of the colloidal
particle, $\tau_{B}\sim\sigma^{3}$. In other words, a factor of ten drop in
particle size means that the Brownian time falls by a factor of 1000, with
major experimental implications. For instance, a phenomenon observed on the
time scale of a day in a light scattering set-up with 200 nm diameter colloids
might take a year to observe with confocal microscopy with 2 $\mu$m diameter
particles. In practice, that would make it unobservable. This distinction is
particularly important for rare events, such as crystal nucleation for which
the “discrepancy regime” is inaccessible to confocal microscopy (see Sec.
XIII), and for phenomena with long time scales, such as the glass transition
(see Sec. XII).

##  V Hard-sphere systems _in silico_

From a numerical simulation perspective, hard spheres might seem to be as
simple as it gets. Their interaction is pairwise and short-ranged, and is also
trivial to calculate. These features should result in efficient simulations,
as is indeed the case for Monte Carlo (MC) Frenkel and Smit (2002) and event-
driven molecular dynamics (EDMD) Rapaport (2004) schemes. Both of these
methods perfectly sample the structure and phase behavior of a pure hard-
sphere model, and EDMD simulations faithfully reproduce the dynamics of
perfect hard spheres in a vacuum. In practice, however, experiments using
colloidal spheres inevitably also include a solvent, which has profound
effects on dynamics especially in out-of-equilibrium systems (see Secs. VII.4
and XI). These hydrodynamic effects can be included at different levels of
detail, with Brownian dynamics (BD) simulations being the most straightforward
implementation that reproduces the diffusive short-time motion of colloidal
particles. Simulating hard-sphere colloidal systems away from thermal
equilibrium, however, requires more complex methods that explicitly model a
coarse-grained version of the background solvent. It is nevertheless important
to note that for dense fluids close to equilibrium, the difference between
different simulation methods in terms of the dynamics is rather small Berthier
and Kob (2007); Montero de Hijes _et al._ (2017). Detailed comparison of EDMD
with experimental data, for instance, found no meaningful difference Royall
_et al._ (2018b).

This section provides an overview of the simulation techniques most commonly
used for hard-sphere models. Because not all of these simulation methods are
easily applicable to perfectly hard spheres, we additionally briefly discuss
commonly used models of _nearly_ hard spheres.

###  V.1 Monte Carlo methods

Since MC simulations were first used to extract the radial distribution
functions and pressure of hard-sphere fluids Rosenbluth and Rosenbluth (1954),
the scheme has remained instrumental for evaluating much of the physical
behavior of the model. A key strength of the approach is its adaptability.
Different thermodynamic ensembles, boundary conditions, or biased sampling can
be considered Frenkel and Smit (2002). Although standard implementations of
these methods are often straightforward, DL_MONTE Brukhno _et al._ (2021), a
general purpose code with advanced sampling methods, and the HOOMD-blue
package Anderson _et al._ (2020), which includes (CPU- or GPU-based)
parallelization for fast simulation of large systems Anderson _et al._ (2013),
should nevertheless be noted.

The simplest local Metropolis algorithm for hard spheres consists of single-
particle attempted displacements, which are accepted only if no overlap is
created. (Interestingly, in the limit of small step sizes, this algorithm
approximates BD Cichocki and Hinsen (1990); Sanz and Marenduzzo (2010).) In
many scenarios, hard-sphere MC simulations can be made more efficient by
incorporating clever collective displacements, such as cluster moves Dress and
Krauth (1995); Almarza (2009); Ashton _et al._ (2011); Dijkstra and van Roij
(1997); Buhot and Krauth (1998); Liu and Luijten (2004), or by employing non-
local sampling schemes. One particularly effective example is to exchange the
diameter of particles in dense mixtures of hard spheres of different sizes
Berthier _et al._ (2016); Ninarello _et al._ (2017). In dense mixtures of not-
too-dissimilar hard spheres, this _swap_ MC approach leads to a major
efficiency gain in configurational sampling. In glass forming hard-sphere
mixtures, in particular, swap MC can probe the equilibrium fluid structure at
packing fractions far beyond the point at which conventional sampling would be
essentially arrested (see Sec. XII.4) Berthier _et al._ (2016); Ninarello _et
al._ (2017).

Another effective strategy for speeding up configurational sampling of hard
spheres is the rejection-free event-chain MC moves introduced in Bernard _et
al._ (2009). In event-chain MC, particle displacements involve a chain of
neighboring particles. A particle is displaced until it collides with a
neighbor, which then becomes the next particle to be displaced. These
rejection-free moves can accelerate configurational sampling by over an order
of magnitude compared to single-particle MC in dense systems of hard spheres
in 2d or 3d Bernard _et al._ (2009); Klement and Engel (2019); Li _et al._
(2022). A variant of this approach called Newtonian event-chain MC, which
incorporates aspects of Newtonian dynamics into the event-chain moves, was
shown to lead to a speed-up of yet another factor of about five for dense
hard-sphere systems Klement and Engel (2019).

###  V.2 Event-driven simulations

Historically, MC simulations were soon followed by EDMD simulations, which
evolve the positions of elastic hard spheres in vacuum Alder and Wainwright
(1957). In other words, each particle interacts purely via instantaneous
elastic collisions, and moves at constant velocity between such collisions
Rapaport (2004). The efficiency of these simulations hinges strongly on the
efficiency of the algorithms used to predict and schedule future collisions
Lubachevsky (1991); Marín _et al._ (1993); Paul (2007); Donev _et al._ (2005),
and can rival event-chain MC Klement and Engel (2019); Smallenburg (2022) for
equilibrating dense systems. As a result, they are vastly more efficient than
MC methods based on single-particle moves. From a technical standpoint, EDMD
simulations are more challenging to implement and adapt, but publicly
available codes have greatly helped spread the approach Donev _et al._ (2005);
Bannerman _et al._ (2011); Smallenburg (2022). EDMD simulations have also been
adapted to, e.g., reproduce Brownian motion Scala _et al._ (2007), to include
swap moves Bommineni _et al._ (2019), and to incorporate isotropic compression
Donev _et al._ (2005).

###  V.3 Brownian dynamics simulations

Suspensions of colloidal hard spheres embedded in a solvent effectively
undergo Brownian motion. The stream of collisions with the much smaller
solvent molecules results in the spheres experiencing random forces. The
ensuing equation of motion is

$\frac{d\mathbf{r}}{dt}=\frac{D_{S}}{k_{\text{B}}T}\mathbf{f}(t)+\sqrt{2D_{S}}\mathbf{R}(t),$
(14)

where $\mathbf{r}(t)$ is the colloid position at time $t$, $D_{S}$ is the
short-time self-diffusion coefficient, and $\mathbf{R}(t)$ is a Gaussian
random process with zero mean and unit variance. The total force
$\mathbf{f}(t)$ on the particle includes its interactions with neighbors as
well as any external forces present in the system.

Simulating Brownian dynamics entails numerically integrating Eq. (14). In
practice, this integration is often achieved using a fixed small time step.
For hard spheres, however, this approach does not work, because any fixed time
step inevitably leads to overlaps. Since the event-driven approach mentioned
above is not particularly efficient Scala _et al._ (2007) – frequent
stochastic updates require re-predicting collisions for all particles after
each velocity change – BD simulations are most commonly performed on model
systems that interact via a continuous approximation of the hard-sphere
potential. (Variants are presented in Sec. V.6.) The implementation approach
for BD simulations closely follows standard molecular dynamics (MD) schemes.
Each time step consists of calculating the forces on all particles, followed
by updating particle positions. As a result, MD codes can be readily adapted
to perform BD simulations, and MD packages such as LAMMPS Thompson _et al._
(2022), HOOMD-Blue Anderson _et al._ (2020), and Espresso Weik _et al._ (2019)
provide efficient implementations of BD, including options for parallelization
or GPU computation.

###  V.4 Simulations including hydrodynamics

All simulation methods described thus far neglect hydrodynamic coupling, which
is (necessarily) present for colloids moving in a solvent. While hydrodynamic
interactions do not modify static equilibrium properties, they can play a
significant role in the dynamics of colloidal hard spheres (see Sec. VII.4),
especially when out-of-equilibrium phenomena are considered (an issue to which
we return in Sec. XI). A number of different numerical methods have therefore
been developed to account for hydrodynamic effects in hard-sphere simulations.

Stokesian dynamics Brady and Bossis (1988) assumes an incompressible fluid and
a low Reynolds number, both of which are typically valid for colloidal hard
spheres. Hydrodynamic interactions between colloids are then taken into
account via a many-particle mobility matrix that models the solvent flow as a
superposition of Stokeslets centered on each colloid. Due to the slow
algebraic decay of this effect, the matrix effectively couples the motion of
all particles in the system, making the method computationally expensive,
particularly when Brownian forces are included. Significantly better
performance can be obtained using the so-called accelerated Stokesian dynamics
(ASD) Sierou and Brady (2001); Banchio and Brady (2003). Publicly available
implementations can be found in the simulation package Espresso Weik _et al._
(2019), as a plugin Fiore and Swan (2019) for HOOMD-Blue Anderson _et al._
(2020), as well as in specialized codes Yan _et al._ (2020).

Lattice Boltzmann (LB) methods Ladd and Verberg (2001); Cates _et al._ (2004);
Dünweg and Ladd (2009) divide the solvent into lattice sites, each carrying a
set of local solvent densities associated with a discrete set of permitted
solvent velocities. These methods reproduce the dynamics of a Newtonian liquid
with a given shear viscosity, and recover the relevant hydrodynamic variables
as moments of the one-particle velocity distribution functions Succi (2001).
The effect of thermal fluctuations can also be included directly into the
discrete Boltzmann equation Adhikari _et al._ (2005); Gross _et al._ (2011).
The fluid dynamics then emerges from the evolution of the one-particle
distribution function. During each simulation time step, these densities are
updated by considering both the streaming of the solvent to neighboring
lattice sites and a collision step, which relaxes the distribution of
velocities at each lattice site towards the equilibrium distribution. The
lattice-based solvent is then coupled to the colloidal particles by imposing
boundary conditions at their surface, thus allowing momentum exchange between
the solvent and the particles. While explicit treatment of the solvent
inherently increases the computational cost, LB methods scale well with
particle size and can be readily parallelized for large systems. Simulation
codes are also publicly available Desplat _et al._ (2001); Weik _et al._
(2019); Bonaccorso _et al._ (2020).

Multi-particle collision dynamics (MPCD) simulations Malevanets and Kapral
(1999); Gompper _et al._ (2009); Howard _et al._ (2019); Kapral (2008) are
similar to LB methods in that they explicitly incorporate the solvent as a
simplified fluid evolved via a streaming and a collision update. In MPCD
approaches, however, the solvent is modeled as discrete _effective particles_
moving in continuous space. At each time step, all particles first move
forward according to their instantaneous velocities (streaming), and then
exchange momentum with other nearby solvent particles during a collision step.
In the most commonly used collision algorithm, the system is divided into
cells, and the velocities of all particles within a cell are mixed via a
collective rotation that conserves their total momentum. Due to this rotation,
MPCD methods using this update are also known as stochastic rotational
dynamics (SRD) Malevanets and Kapral (1999); Ihle and Kroll (2001). Coupling
to the colloids can be done either by including the colloidal particles as
collision partners in the collision step, or via direct interactions between
solvent particles and colloids during the streaming step. It should be noted
that MPCD methods model compressible solvents. In terms of implementation,
MPCD methods are closely connected to standard MD, and public (parallellized)
implementations are available in, e.g., LAMMPS Thompson _et al._ (2022) and
HOOMD-Blue Anderson _et al._ (2020); Howard _et al._ (2018).

Dissipative particle dynamics Hoogerbrugge and Koelman (1992); Espanol and
Warren (1995); Groot and Warren (1997); Espanol and Warren (2017) (DPD) is
similar to MPCD in spirit but introduces soft repulsive interactions that
include both a dissipative and a stochastic term to replace the collision
step. These interactions are idealized in such a way that the dynamics of the
solvent particles can be resolved with time steps much larger than those used
for the (steeper) colloid interactions, thus significantly lowering the
computational cost compared to full-solvent simulations. Coupling to the
colloids is achieved by associating each colloid particle with one or more
solvent particles that overlap with it, and interact with the remaining
solvent. Similar to the MPCD approach, DPD can be seen as a variation on
classical MD, and efficient implementations can be found in several simulation
packages, including LAMMPS Thompson _et al._ (2022), Espresso Weik _et al._
(2019), and HOOMD-Blue Anderson _et al._ (2020).

Fluid particle dynamics Tanaka and Araki (2000) (FPD) is based on a direct
numerical simulation of the Navier-Stokes equations, resolving the problem of
coupling the fluid and the colloidal particles by approximating each colloidal
particle as highly viscous particle with a smooth viscosity profile at its
interface with the fluid. The biggest advantage of the method is that it
avoids the difficulties associated with moving solid-fluid boundary conditions
and that it allows the fluid to couple to additional fields. The fluid itself
is also incompressible. Thermal fluctuations can additionally be considered by
integrating the fluctuating hydrodynamics equations Furukawa _et al._ (2018).
A method akin to FPD is the _smooth profile method_ (SPM) Nakayama and
Yamamoto (2005), which avoids introducing a large viscosity in the particle
domain by separating the calculation of the hydrodynamic forces from that of
the boundary conditions. The approach has found application to a large variety
of passive and active soft-matter systems Yamamoto _et al._ (2021).

###  V.5 Advanced sampling and free-energy methods

For the purpose of studying more detailed aspects of hard-sphere behavior,
such as phase boundaries, defects, interfacial properties and nucleation, a
wide variety of advanced simulation techniques have been combined with the
general approaches outlined above.

Arguably the most immediate way of obtaining phase boundaries is to simulate
the direct coexistence between the two phases of interest, e.g., a hard-sphere
fluid and a hard-sphere crystal. In practice, however, this approach suffers
from large finite-size effects, in addition to being sensitive to statistical
error as well as to details in the set-up of the initial configuration
Espinosa _et al._ (2013). Moreover, this approach is usually unsuitable for
determining crystal-crystal coexistence, due to its slow equilibration
dynamics. As a result, studies of equilibrium phase behavior in hard spheres
often make use of free-energy calculations. Although free energies can
typically not be directly measured in simulations, thermodynamic integration
can be used to determine free-energy differences with respect to some known
reference state Frenkel and Smit (2002). The series of simulations that define
the integration path between the system of interest and a reference system can
be chosen in a variety of ways. Common choices are the ideal gas for the fluid
phase Frenkel and Smit (2002) and an Einstein crystal for an ordered phase
Frenkel and Ladd (1984); Vega and Noya (2007), but a variety of other
reference states have been employed for hard-sphere systems, including cell
models Hoover and Ree (1968); Nayhouse _et al._ (2011) and pinning schemes
Schilling and Schmid (2009); Moir _et al._ (2021). Note that the hard-sphere
system itself has also often been used as a reference system for more complex
systems (e.g. Dijkstra _et al._ (1999a); Bolhuis and Frenkel (1997); Patey and
Valleau (1973); Smallenburg and Sciortino (2013)). After obtaining the free
energy of all relevant phases for a range of state points, equilibrium phase
boundaries can be determined by identifying which phase (or coexistence of
phases) has the lowest free energy at any given state point. As an alternative
approach, accurate phase coexistence conditions can be extracted by mapping
out the free-energy landscape connecting the fluid and crystal phases Martin-
Mayor _et al._ (2011); Fernández _et al._ (2012). Note that all of these
methods require initial knowledge of the potential phases of interest, which
may require an extensive search for possible crystal structures Filion _et
al._ (2009); Filion and Dijkstra (2009); Kummerfeld _et al._ (2008); O’Toole
and Hudson (2011); Hudson and Harrowell (2011); Hopkins _et al._ (2012),
especially in multicomponent mixtures.

Free-energy integration methods also provide routes to calculate the effect of
interfaces and defects on the free energy of hard-sphere systems. In
particular, the equilibrium concentration of point defects in hard-sphere
crystals has been determined by setting up a thermodynamic integration path
between systems with and without a defect Pronk and Frenkel (2001); van der
Meer _et al._ (2017). Similar approaches can also be used to determine the
fluid-solid interfacial free energy of hard spheres (Sec. VIII)—by
constructing an integration path between a direct fluid-solid coexistence and
a state of two separate bulk phases Davidchack and Laird (2000); Davidchack
(2010); Schmitz and Virnau (2015)—or between a pure fluid and a fluid in
coexistence with a slab of crystal Espinosa _et al._ (2014); Sanchez-Burgos
_et al._ (2021). Alternatively, interfacial free energies can be determined
from interfacial fluctuations Davidchack _et al._ (2006), or via thermodynamic
integration Bültmann and Schilling (2020). The free-energy difference between
face-centered-cubic (fcc) and hexagonal-close-packed (hcp) crystals and the
hcp-fcc interfacial free energy can be computed using the lattice-switch
method Bruce _et al._ (1997); Pronk and Frenkel (1999); Wilding and Bruce
(2000).

For studying rare events, such as nucleation, it is common to use sampling
schemes that purposefully bias the simulation towards sampling rare
configurations. In particular, umbrella sampling Torrie and Valleau (1977);
Kästner (2011); Filion _et al._ (2010) has been used extensively to determine
the barrier to nucleating hard-sphere crystals by calculating the free-energy
cost of creating crystalline nuclei of different sizes. This approach,
however, assumes that the cluster is in local equilibrium at each cluster
size, and therefore ignores the possibility of any non-equilibrium dynamics.
Hence, to explore nucleation trajectories that take into account the dynamics
of the system, more specialized methods such as forward-flux sampling Allen
_et al._ (2005); Allen _et al._ (2006b, a); Allen _et al._ (2009); Hussain and
Haji-Akbari (2020), transition path sampling Bolhuis _et al._ (2002), or
transition interface sampling Moroni _et al._ (2004) are required. Biasing
methods have also found uses in other areas of hard-sphere simulations, such
as with potentials that suppress crystallization to study hard-sphere glasses
Valeriani _et al._ (2011); Taffs and Royall (2016).

###  V.6 Simulation models

The general strategy to approximate the hard-sphere interaction is to consider
a harshly-repulsive continuous interaction potential. Perhaps the simplest of
such forms is the inverse power-law potential of exponent $n$, as given in Eq.
(3). In the limit $n\rightarrow\infty$, this potential converges to the true
hard-sphere potential with diameter $\sigma$. Similar to hard spheres, the
thermodynamic behavior of the inverse power-law model with a given exponent
$n$ is effectively controlled by a single parameter, since a change in the
interaction strength $\epsilon$ and the particle size $\sigma$ both have the
same trivial effect of scaling the total energy of the system. Experimental
measurements of the interactions between sterically stablized colloids have
further been shown to agree well with power-law potentials Bryant _et al._
(2002) (see Sec. III.1). However, the inverse power law does not naturally
capture the short-range nature of the soft interaction that results from
steric stabilization.

Another commonly used continuous approximation is the Weeks-Chandler-Andersen
(WCA) potential. Originally introduced to consider the separate roles of
attraction and repulsion on the structure of simple liquids Weeks _et al._
(1971), the WCA potential corresponds to a purely repulsive variant of the
Lennard-Jones potential

$u(r)=\begin{cases}4\epsilon\left(\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}+\frac{1}{4}\right)&r\leq
2^{1/6}\sigma\\\ 0&r>2^{1/6}\sigma\end{cases},$ (15)

with energy parameter $\epsilon$. This interaction form smoothly approaches
zero at the cutoff distance $r_{c}=2^{1/6}\sigma$. At sufficiently low
temperatures (a common choice is $\epsilon/k_{\text{B}}T=40$), the WCA
potential has been shown to map well to hard spheres in terms of, e.g., the
equation of state and nucleation rates Dasgupta _et al._ (2020); Richard and
Speck (2018a); Filion _et al._ (2011b). Variations of Eq. (15), with different
(typically higher) exponents have been designed to match the hard-sphere
behavior even more closely. A notable example of this is the so-called pseudo-
hard-sphere potential Jover _et al._ (2012), for which the freezing transition
lies very close to that of true hard spheres Espinosa _et al._ (2013). This
class of potentials has a long history as a substitute for hard spheres, and
has the advantage of being purely short-ranged and being continuous, thus
allowing for an easy BD implementation. However, the thermodynamic behavior of
this class of models depends on both the chosen temperature and density.

Note that regardless of the exact functional form chosen, a trade off must
generally be made between simulation efficiency and the accuracy of the
approximation of the true hard-sphere potential. Generally, although a sharper
interaction potential provides a better approximation of hard spheres, it also
results in a more rapid variation of the interaction forces as two particles
approach. Since MD and BD simulations numerically integrate the force
experienced by each particle over time, these faster variations then
necessitate a smaller integration time step to maintain accurate results, thus
slowing simulations down.

##  VI Hard-sphere systems in theory

As noted in Secs. I and II, hard spheres came into existence as a minimal
model system for exploring condensed matter, and have been extensively
reviewed Hansen and McDonald (2013); Santos _et al._ (2020); Roth (2010);
Tarazona and Evans (1984). Here we only cover a couple of key theories which
are mentioned in later sections, and which we believe had a significant impact
in the success of the hard-sphere system.

###  VI.1 Integral equation theory

Integral equation theory lies at the heart of descriptions of the liquid state
Hansen and McDonald (2013); Barker and Henderson (1976). In the context of
hard-sphere colloids as model atoms and molecules, it provides an elegant and
simple theoretical framework for obtaining predictions that may be directly
probed in experiments. For our purposes, integral equation theory boils down
to finding approximate solutions (or closures) for the Ornstein–Zernike (OZ)
relation

$h(r)=c(r)+\rho\int c(r)h(r)d\mathbf{r},$ (16)

which relates the total correlation function $h(r)=g(r)-1$ (with $g(r)$ the
radial distribution function) to the direct correlation function $c(r)$.

In the framework of classical density functional theory (see Sec. VI.3), the
direct correlation function is related to the second functional derivative of
the excess part of the free-energy functional with respect to the one-particle
density of the system Hansen and McDonald (2013). A homogeneous and isotropic
system at density $\rho$ and temperature $T$ interacting with a pairwise
additive potential $u(r)$ is therefore _uniquely_ defined by its $h(r)$ and
thus by its $g(r)$, which is measurable in experiments (see Sec. III).
Measuring $g(r)$ – or its Fourier transform $S(k)$ – is therefore an important
means of characterizing the structure and thermodynamics of an experimental
system.

The OZ relation can also be solved using a second (closure) relation. A
relatively simple relation, first proposed by Percus and Yevick (PY) Percus
and Yevick (1958), is given by

$c(r)=(\exp[-\beta u(r)]-1)(h(r)-c(r)+1).$ (17)

This relationship is highly accurate for strongly repulsive and short-range
interactions, such as for hard spheres. In 1963, Wertheim and Thiele
independently showed that the PY closure to the OZ equation for a fluid of
hard spheres with diameter $\sigma$ and volume fraction
$\phi=\pi\sigma^{3}\rho/6$ yields for the direct correlation function Wertheim
(1963); Thiele (1963)

$c(r)=\\\
\begin{cases}\frac{-(1+2\phi)^{2}+6\phi(1+\frac{1}{2}\phi)^{2}\frac{r}{\sigma}-\frac{1}{2}\phi(1+2\phi)^{2}(\frac{r}{\sigma})^{3}}{(1-\phi)^{4}}&r\leq\sigma\\\
0&r>\sigma.\end{cases}$

$c(r)$ can be analytically Fourier transformed, from which the structure
factor $S(k)$ follows using the Fourier transform of the Ornstein-Zernike
equation. The pair correlation function $g(r)$ can subsequently be obtained by
a numerical Fourier transform of $S(k)$. The theoretical predictions agree
well with computer simulation results for volume fractions $0\leq\phi\lesssim
0.5$, i.e., over the entire stable fluid regime. Examples of experimental
comparisons against the PY approximation to the OZ equation are shown in Fig.
12 in Sec. VII.

Given $g(r)$, three independent routes can then be followed to extract
thermodynamic quantities:

$\displaystyle P$ $\displaystyle=$ $\displaystyle\rho
k_{\text{B}}T-\frac{\rho^{2}}{6}\int d{\bf r}ru^{\prime}(r)g(r),$ (18)
$\displaystyle\frac{E}{V}$ $\displaystyle=$ $\displaystyle\frac{3}{2}\rho
k_{\text{B}}T+\frac{\rho^{2}}{2}\int d{\bf r}u(r)g(r),$ (19) $\displaystyle
k_{\text{B}}T\left(\frac{\partial\rho}{\partial P}\right)_{T}$
$\displaystyle=$ $\displaystyle 1+\rho\int d{\bf r}\big{(}g(r)-1\big{)},$ (20)

which are the virial, caloric, and compressibility routes, respectively. An
important exact sum rule for hard spheres can be obtained straightforwardly
from Eq. (18)

$\frac{\beta P}{\rho}=1+\frac{2\pi\rho\sigma^{3}g(\sigma^{+})}{3},$ (21)

which relates the pressure of the hard-sphere fluid to the contact value of
$g(r)$. Using $g(r)$ from the PY approximation to the OZ equation in Eq. (20)
yields the compressibility pressure $P_{c}$

$\frac{\beta P_{c}}{\rho}=\frac{1+\phi+\phi^{2}}{(1-\phi)^{3}},$ (22)

while the virial pressure $P_{v}$ from Eq. (18) gives

$\frac{\beta P_{v}}{\rho}=\frac{1+2\phi+3\phi^{2}}{(1-\phi)^{2}}.$ (23)

The compressibility equation of state overestimates the pressure obtained from
simulations, whereas the virial equation of state underestimates the
simulation results Hansen and McDonald (2013). It turns out that the Carnahan
and Starling (CS) linear combination $P_{\mathrm{CS}}=(2P_{c}+P_{v})/3$
Carnahan and Starling (1969),

$\frac{\beta
P_{\mathrm{CS}}}{\rho}=\frac{1+\phi+\phi^{2}-\phi^{3}}{(1-\phi)^{3}}.$ (24)

is essentially indistinguishable from simulations up to $\phi\simeq 0.5$
Hansen and McDonald (2013). It should be noted, however, that Carnahan and
Starling did not use simulation results but rather the pressure from the
virial expansion for the equation of state Eq. (2) to motivate their
expression Carnahan and Starling (1969). Using the analytical expressions for
the virial coefficients $B_{2},B_{3}$ and $B_{4}$, and the numerical
expressions for $B_{5}$ and higher, the equation of state for hard spheres
becomes

$\displaystyle\frac{\beta P}{\rho}$ $\displaystyle=$ $\displaystyle
1+4\phi+10\phi^{2}+18.365\phi^{3}+28.225\phi^{4}+39.74\phi^{5}$ (25)
$\displaystyle+53.5\phi^{6}+70.8\phi^{7}+\cdots.$

Approximating coefficients as $B_{i+1}=(\pi\sigma^{3}/6)^{i}(i^{2}+3i)$ then
gives Eq. (24) Carnahan and Starling (1969).

While integral equation theory is successful in predicting the structure and
thermodynamics of equilibrium hard-sphere fluids, the closures are somewhat
uncontrolled and not systematically improvable. A clear manifestation of this
effect is that predictions for binary mixtures of hard spheres turn out to be
extremely sensitive to the details of the approximation scheme. A more robust
approach in this sense is density functional theory, and especially
fundamental measure theory, described in Refs. Tarazona _et al._ (2008); Roth
(2010) and Sec. VI.3.

###  VI.2 Cell theory

The crystal branch of the equation of state of hard spheres is less amenable
to integral equation theory, due to the broken translational and rotational
symmetry in the system. One early mean-field theory that can be used as an
approximate analytical prediction was originally introduced by Lennard-Jones
and Devonshire (1937) to estimate fluid free energies, and is now commonly
known as cell theory or free volume theory when applied to solids and
especially crystals Wood (1952); Salsburg and Wood (1962). In this approach,
it is assumed that the Helmholtz free energy $F(N,V,T)$ of the crystal phase
can be subdivided into contributions from each particle, and that each can be
computed by assuming that all other particles are fixed exactly at their
lattice site. In that approximation, the free energy $f_{1}$ of a single
particle in an FCC lattice is given by

$\beta f_{1}(\rho)=-\log\frac{V_{\mathrm{free}}(\rho)}{\Lambda^{3}},$ (26)

where $V_{\mathrm{free}}$ is the free volume available to the particle
confined in the cage of its neighbors, and $\Lambda$ is the thermal De Broglie
wavelength. This free volume is then typically approximated as a sphere with
diameter $a-\sigma$, with $a$ the nearest-neighbor distance, which is given by
$a=(\rho_{\mathrm{max}}/\rho)^{1/3}$. Here, the close-packed density of
crystal lattices of monodisperse hard spheres is
$\rho_{\mathrm{max}}\sigma^{3}=\sqrt{2}$. Hence, cell theory gives the
approximate total free energy

$\frac{\beta F}{N}\simeq\beta
f_{1}\simeq-\log\left[\frac{4\pi\sigma^{3}}{3\Lambda^{3}}\left(\left(\frac{\rho_{\mathrm{max}}}{\rho}\right)^{1/3}-1\right)^{3}\right].$
(27)

Taking the derivative with respect to the volume gives the pressure

$\frac{\beta
P}{\rho}=\left(1-\left(\rho/\rho_{\mathrm{max}}\right)^{1/3}\right)^{-1}.$
(28)

Although this expression significantly underestimates the pressure (on the
order of $\sim 1k_{\text{B}}T/\sigma^{3}$), cell theory provides a quick and
physically intuitive estimate for the crystal free energy and equation of
state. As a result, it is commonly used for pedagogical purposes Barrat and
Hansen (2003); Kamien (2007), as well as for estimating phase transitions.
Cell theory has been extended in a variety of ways (see, e.g., Rudd _et al._
(1968); Koch _et al._ (2005); Charbonneau _et al._ (2021b)), including for
binary mixtures of hard spheres Cottin and Monson (1993, 1995); van der Meer
_et al._ (2020) and glassy systems. It should be noted that when cell theory
is applied to disordered systems, such as glasses, it cannot be expected to
take into account the configurational entropy of the system (see Sec. XII).
Furthermore, in the fluid regime, particles are no longer trapped in cells in
any meaningful way, and hence the description should be expected to break down
there as well.

###  VI.3 Classical density functional theory

Classical density functional theory (DFT) provides an exact theoretical
framework for describing the thermodynamic and structural properties of
interacting many-body systems, starting from a microscopic description of the
interparticle interactions. The approach is based on the observation that the
grand potential of a specified inhomogeneous fluid is a functional
$\Omega[\rho(\mathbf{r})]$ of the variational one-body density profile
$\rho(\mathbf{r})$, with the properties that (_i_) the equilibrium density
profile $\rho_{0}(\mathbf{r})$ minimizes the functional
$\Omega[\rho(\mathbf{r})]$, and (_ii_) this minimum equals the equilibrium
grand potential $\Omega_{0}[\rho_{0}(\mathbf{r})]$. From $\Omega_{0}$, all
thermodynamic properties then follow. For instance, one can obtain the
homogeneous system pressure

$P=-\frac{\Omega_{0}}{V}$ (29)

Functional derivatives of the Helmholtz free-energy functional additionally
provide correlation functions. In particular, the OZ two-body direct
correlation function is related to the second functional derivative of the
functional from which the two-body structure of the system follows. However,
because the exact Helmholtz free-energy functional for a given interaction
potential is unknown, DFT depends on approximate free-energy density
functionals. The efficacy of a given DFT approximation is therefore commonly
tested by comparing DFT predictions with computer simulations.

Exploiting the wealth of available information on the thermodynamics and
structure of hard-sphere systems, high-accuracy approximations of that
functional have been constructed for hard-sphere mixtures starting with
Rosenfeld’s Fundamental Measure Theory (FMT) in 1989 Rosenfeld (1989). In FMT,
the Helmholtz free-energy functional is based on weighted densities that are
convolutions of the density profiles with weight functions depending on the
geometrical properties of the spheres. Various subsequent extensions and
modifications have been proposed to improve the description of various
inhomogeneous systems of hard spheres. Further details on this approach can be
found in the comprehensive and excellent reviews by Tarazona Tarazona _et al._
(2008) and Roth Roth (2010).

##  VII Bulk equilibrium hard spheres

Given the wide array of experimental, numerical and theoretical methods
described in the previous sections, we can now consider the equilibrium
behavior of hard-sphere systems. This behavior not only serves as a reference
point in nearly all subsequent sections of this review, but is also an
important scientific matter in its own right.

Naturally, the first hard-sphere systems studied by computer simulation were
monodisperse Alder and Wainwright (1957); Wood and Jacobson (1957), but only
_quasi_ –monodisperse in colloidal experiments. However, because–as we will
see later in this section–a mean size polydispersity $s\lesssim 5\%$ barely
changes the equilibrium phase diagram Wilding and Sollich (2010), both
experimental and numerical work nevertheless focused on equivalent systems. In
this context, a natural starting topic is their bulk phase behavior. Other
topics considered include the equation of state as well as fluid structure and
dynamics for both 3d and (quasi–)2d hard spheres.

###  VII.1 Equilibrium phase behavior of monodisperse hard spheres

Since hard spheres do not interact beyond contact, their (Helmholtz) free
energy only comprises an entropic contribution. Temperature plays but a
trivial role, and hence the hard-sphere equilibrium phase diagram depends
solely on the volume fraction $\phi$ (density) of the system. At low $\phi$
the system is naturally in a fluid phase, whereas at high $\phi$ the system is
in a crystal phase. (In 3d, the two are separated by a first-order phase
transition.) Put differently, at sufficiently high $\phi$, the crystal entropy
is higher than that of the fluid; the regular arrangement of spheres on the
crystal lattice provides each particle with more local free volume to move
around than it would have in a fluid at the same density.

As mentioned in Sec. II, the first numerical validations of this theoretical
prediction were obtained in the 1950s by studying the melting of a face-
centered cubic (fcc) crystal Alder and Wainwright (1957); Wood and Jacobson
(1957). (The crystal phase with the fcc structure is narrowly
thermodynamically favored over competing close packed structures, as discussed
later in this section.) Since then, considerable efforts have been devoted to
pinning down the details of this transition, including characterizing finite-
size effects Polson _et al._ (2000). Table 1 provides an overview of numerical
results for the coexistence packing fractions, pressures, and chemical
potentials. Apart from a few outliers, often associated with small system
sizes, results for these quantities are in good agreement, and have
established that the freezing and melting volume fractions of monodisperse
hard spheres are $\phi_{f}\simeq 0.492$ and $\phi_{m}\simeq 0.543$,
respectively. Similar predictions have also been obtained from various
theoretical treatments, notably fundamental measure theory Roth (2010).
Experimental measurements on hard-sphere colloids, such as those shown in Fig.
2 Pusey and Van Megen (1986), agree reasonably well with this result, provided
a suitable $\sigma_{\mathrm{eff}}$ is chosen to account for any residual
softness (see Sec. III.3).

Note that while these studies do not typically take into account the presence
of vacancies in the equilibrium crystal phase, its effect on the coexistence
pressure is expected to be smaller than the typical error bars considered here
($\beta\Delta P\sigma^{3}\simeq-0.0026$) Pronk and Frenkel (2001). Simulation
results have revealed that the equilibrium fcc crystal is indeed nearly
flawless. The equilibrium fraction of defective sites at melting has been
computed to be approximately $10^{-4}$ for vacancies and $10^{-8}$ for
interstitials Bennett and Alder (1971); Pronk and Frenkel (2001, 2004).
However, comparing these predictions to experiments is inherently difficult
beyond the usual concerns about coexistence determination, given (_i_) the low
concentrations involved, (_ii_) their sensitivity to the crystal packing
fraction, and (_iii_) the possibility that defects become kinetically trapped
during crystal formation (Sec. III.3) Poon _et al._ (2012); Royall _et al._
(2013).

Interestingly, various properties of the bulk crystal phase, including its
elastic constants Runge and Chester (1987); Pronk and Frenkel (2003); Sushko
and van der Schoot (2005) and defect diffusivity van der Meer _et al._ (2017),
have been studied in detail in simulation but remain to be systematically
investigated in experiments.

Table 1:  Table summarizing prior reports of hard sphere fluid-fcc phase coexistence using computer simulations (black) and DFT (green). For each result, we list the main method used to determine the free energy or stability of the solid phase (see also Sec. V), the freezing and melting packing fractions, the coexistence pressure and chemical potential, and the system size considered. A system size of $\infty$ indicates that results were extrapolated to the thermodynamic limit $N\rightarrow\infty$, but the system sizes considered and the extrapolation scheme vary between studies. The final consensus averages results from five studies with high reported accuracy (marked by †), namely Frenkel and Smit (2002); Fortini and Dijkstra (2006); Fernández _et al._ (2012); Pieprzyk _et al._ (2019); Moir _et al._ (2021) Source | Method | $\phi_{f}$ | $\phi_{m}$ | $\beta P_{\mathrm{coex}}\sigma^{3}$ | $\beta\mu_{\mathrm{coex}}$ | $N$
---|---|---|---|---|---|---
Hoover and Ree (1968) | Single-occupancy cell | 0.494 | 0.545 | 11.70 |  | $\infty$
Speedy (1997) | Unspecified | 0.491(1) | 0.543(1) | 11.55(11) |  |
Davidchack and Laird (1998) | Direct coexistence | 0.491 | 0.543 | 11.55(5) |  | 10752
Wilding and Bruce (2000) | Phase switch |  |  | 11.49(9) |  | $\infty$
†Frenkel and Smit (2002) | Ladd-Frenkel | 0.4917 | 0.5433 | 11.567 | 16.071 | $\infty$
Roth _et al._ (2002) | FMT (White Bear I) | 0.489 | 0.536 |  |  |
Errington (2004) | Phase switch |  |  | 11.43(2) |  | $\infty$
†Fortini and Dijkstra (2006) | Ladd-Frenkel | 0.4915(5) | 0.5428(5) | 11.57(10) | 16.08(10) |
Vega and Noya (2007) | Einstein molecule | 0.492 | 0.543 | 11.54(4) | 16.04 | $\infty$
Noya _et al._ (2008) | Direct coexistence | 0.491(1) | 0.543(2) | 11.54(4) |  | 5184
Odriozola (2009) | Parallel tempering | 0.492(4) | 0.545(4) | 11.43(17) |  | $\infty$
Zykova-Timan _et al._ (2010) | Direct coexistence | 0.492 | 0.545 | 11.576(6) |  | $160000$
Nayhouse _et al._ (2011) | Single-occupancy cell | 0.4912(6) | 0.5425(5) | 11.502(19) |  | $\infty$
Oettel _et al._ (2010) | FMT (White Bear II) | 0.495 | 0.545 | 11.89 | 16.40 |
de Miguel (2008) | Ladd-Frenkel | 0.489(1) | 0.541(2) | 11.43(11) | 15.92(11) | 864
†Fernández _et al._ (2012) | Free-energy landscape | 0.49188(2) | 0.54312(5) | 11.5727(10) |  | $\infty$
Ustinov (2017) | Kinetic MC | 0.491 | 0.543 | 11.534 | 16.054 | 3000
†Pieprzyk _et al._ (2019) | Literature values | 0.4917(5) | 0.5433(5) | 11.5712(10) | 16.0758(20) | $\infty$
†Moir _et al._ (2021) | Tethered particle | 0.49161(3) | 0.54305(4) | 11.550(4) | 16.053(4) | $\infty$
Consensus |  | 0.4917(2) | 0.5431(2) | 11.56(2) | 16.06(2) | $\infty$

_Phase behavior in experiment._ – As noted in Sec. II, interest in hard
spheres received a tremendous boost from the experiments of Pusey and Van
Megen (see Fig. 2) Pusey and Van Megen (1986), but earlier work by Hachisu and
coworkers Kose and Hachisu (1974); Hachisu and Kobayashi (1974) had already
demonstrated that hard sphere-like colloids crystallize. Additionally, results
for the strongly screened charged stabilized aqueous suspension of the
Japanese group were consistent with theoretical and simulation predictions
Hachisu and Kobayashi (1974). When it came to non-aqueous systems (for which
electrostatic effects were expected to be minimal), Kose and Hachisu (1974)
used what would now be called microgel colloids of crosslinked PMMA dispersed
in benzene 888It is worth noting that benzene is now typically avoided, given
its toxicity.. Given the difficulty of determining effective volume fractions
for such particles (see Sec. III.3) Poon _et al._ (2012); Royall _et al._
(2013), however, a quantitative comparison was challenging for this system
Kose and Hachisu (1974). Later work indeed found that the softness of microgel
particles (which can be tuned via the degree of crosslinking used Lyon and
Fernandez-Nieves (2012)) indeed can have a profound effect on the phase
coexistence Paulin and Ackerson (1996), as may be inferred from the discussion
of Fig. 6. Experimental efforts to determine the osmotic pressure at freezing
were pioneered by Takano and Hachisu Takano and Hachisu (1977); Hachisu and
Takano (1982), who, after drying and measuring the total mass of particles
from the fluid phase in a sedimented system, found reasonable agreement with
simulation results.

Figure 9: Determining the equilibrium coexistence volume fraction from the
time–evolution of a sedimenting suspension. The height of a system in the
fluid-crystal coexistence region undergoing sedimentation is recorded as a
function of time. After a short waiting period of two days, sharp interfaces
form between the clear supernatant (A), the fluid phase (B), the crystal phase
(C), and a dense polycrystalline sediment (D). Extrapolating the fluid-crystal
(B-C) interface (red line) back to time $t=0$ provides an estimate of the
phase boundary prior sedimentation settling in, i.e., in the absence of
gravity. Reproduced with permission from Paulin and Ackerson (1990).

The subsequent development of sterically stabilized PMMA (see Secs. II.2 and
III.2) brought better controlled hard-sphere-like experimental systems to the
field. A key hurdle to achieving quantitative accuracy, however, has been to
properly determine the (effective) $\phi$ in experiment (see Sec. III.3) Poon
_et al._ (2012); Royall _et al._ (2013). One might think that even if
particles have some degree of softness (as is always the case), phase
boundaries should nevertheless lead to volume fractions that can be mapped
onto the predicted values. In fact, as discussed in Sec. III.3 and Fig.6, it
is not that simple, as softness effects the relative coexistence gap between
fluid and crystal. Furthermore, as noted by Pusey and Van Megen, sedimentation
complicates matters Pusey and Van Megen (1986).

One clever method to address the effects of sedimentation was proposed by
Paulin and Ackerson (1990). As shown in Fig. 9, the fluid-crystal phase
boundary (in vials such as those shown in Fig. 2) tends to move upwards as the
fluid phase slowly crystallizes. While the system starts as a metastable
fluid, with an effective volume fraction that falls within the coexistence
regime (see Fig. 3), the final sedimentation-diffusion equilibrium state is a
crystal with only a thin fluid phase above. (The thickness of the fluid phase
is of the order of the gravitational length defined in Sec. XI, which in this
case is very much less than the container size.) In Fig. 9, the top of the
sediment is characterized by the supernatant (A)–colloidal fluid (B)
interface. Sedimentation proceeds over the course of the experiment–some two
months–but after approximately 27 days, the colloidal fluid (B) is no longer
visible to the naked eye. The fluid “disappearance” reflects the approach to
sedimentation–diffusion equilibrium. The colloidal fluid–crystal (C) phase
separation, however, took only about one day. (An additional, distinct layer
in the form of a polycrystalline sediment (D) also formed, akin to that
observed in Fig. 2(c), but plays no role in the subsequent analysis.)

Phase diagram determination from these observations requires a key assumption:
that the boundary in the sedimenting system extrapolated back to $t=0$
reflects the phase-separated system prior to any effects from sedimentation.
In other words, $t=0$ should correspond to a phase-separated system in the
_absence_ of gravity. Using different samples of varying initial $\phi$,
application of the lever rule then gives the equilibrium fluid-crystal
coexistence region. This method, therefore requires a separation of timescales
between that of sedimentation (1 month) and that of fluid-crystal phase
separation (1 day). The resulting phase behavior was later found to be in
reasonable agreement with the experimental equation of state Rutgers _et al._
(1996); Phan _et al._ (1996). (Section VII.3 discusses what _reasonable_
agreement means in the context of colloidal experiments.)

_Polymorphism._ – At sufficiently high $\phi$, one expects the stable crystal
to be one of the close-packed structures, known as Barlow packings, considered
in the Kepler conjecture Hales _et al._ (2017), which consist of stacked
layers of spheres arranged at the vertices of a hexagonal (or triangular)
lattice. The non-uniqueness of the stacking follows from the fact that each
additional layer has two different positions in which it can be placed with
respect to the previous one. Computer simulations combined with free-energy
calculations have demonstrated that the fcc structure (with an ABCABC stacking
order) is slightly entropically favored over the hexagonal close-packed (hcp)
structure (with an ABAB stacking order) Frenkel and Ladd (1984); Bolhuis _et
al._ (1997); Mau and Huse (1999); Polson _et al._ (2000). Given that the free
energy difference is only about $10^{-3}k_{\text{B}}T$ per particle, however,
spontaneously formed hard-sphere crystals often contain defects in the fcc
stacking sequence, resulting in a random-hexagonal-close-packing (rhcp)
structure. This polytype999Polytypes are polymorphs whose symmetry only
differs in one direction. is unlikely to anneal out to the equilibrium fcc
crystal phase over any realistic time scale in simulation Pronk and Frenkel
(1999). However, annealing of rhcp into fcc crystal structures has been
observed in X-ray crystallographic experiments of suspensions of quite small
particles (600nm) Kegel and Dhont (2000); Martelozzo _et al._ (2002). Small
external biases can also play a significant role. For instance,
crystallization resulting from sedimentation under gravity seems to more
strongly favor the formation of fcc Hoogenboom _et al._ (2002); Marechal _et
al._ (2011), whereas experiments in microgravity result in rhcp crystals Zhu
_et al._ (1997); Cheng _et al._ (2001b). A possible explanation is that
mechanical instability occurring during the packing process of the former
selects the fcc structure Heitkam _et al._ (2012).

###  VII.2 Impact of polydispersity on the equilibrium phase diagram

As mentioned in Sec. III, preparing a sample of colloidal spheres necessarily
leads to a size distribution $p(\sigma)$ of diameters instead of the idealized
identical-sphere model. The standard deviation of $p(\sigma)$ relative to the
mean diameter defines the polydispersity $s$, which typically quantifies this
effect.

The phase behavior of systems with small polydispersity, $s\leqslant 5$%, is
essentially indistinguishable ($\lesssim 2\%$ difference in coexistence
$\phi$) from that of monodisperse hard spheres. Going to larger
polydispersities, however, strongly suppresses crystallization. The existence
of a “terminal polydispersity” ($s_{\text{t}}\approx 5-12\%$ depending on the
higher moments of the distribution), beyond which the crystal cannot be
stabilized was therefore initially suggested Barrat and Hansen (1986); Pusey
(1987); Bolhuis and Kofke (1996); Phan _et al._ (1998); Bartlett and Warren
(1999); Ackerson and Schätzel (1995).

Figure 10: Phase diagram of polydisperse hard spheres. Data are from Sollich
and Wilding’s free-energy calculations for a top hat distribution Fasolo and
Sollich (2004). Beyond 6%, the single solid phase fractionates into coexisting
solids of different size distributions with smaller polydispersity. The larger
the polydispersity and packing fraction, the larger the number of coexisting
phases needed to satisfy the polydispersity constraint on individual crystal
phasesSollich and Wilding (2010). Reproduced with permission from Fasolo and
Sollich (2004).

Follow-up simulations, however, revealed that polydisperse suspensions
circumvent this barrier through _fractionation_ (see Fig. 10). The solid phase
then splits into coexisting solid phases with smaller polydispersity (each
being limited to $s_{\ast}\approx 6\%$) Kofke and Bolhuis (1999). In other
words, instead of a single crystal with widely different particle sizes
(implying large local strains), the system favors the formation of separate
crystalline domains of small and large particles. This proposal has since been
put onto solid theoretical ground through (approximate) free-energy
calculations Fasolo and Sollich (2003, 2004).

Support for fractionation can be obtained from computer simulations, albeit
with particular methodological care. In standard fixed-volume simulations,
accessing equilibrated fractionated solids is obfuscated by the extremely slow
nucleation dynamics (see Sec. XIII) and the slow dynamics of domain
boundaries. Sollich and Wilding (2010) have side-stepped these difficulties by
employing constant pressure Monte-Carlo simulations that prescribe a size-