dependent chemical potential $\Delta\mu(\sigma)$ (instead of fixing the parent
distribution $p(\sigma)$). This demonstration of qualitative agreement between
simulations and theoretical work (see Fig. 10 Wilding and Sollich (2010);
Sollich and Wilding (2011)), is supported by some evidence for fractionation
in experiments using light scattering Martin _et al._ (2003).

Intriguingly, recent swap-assisted EDMD simulations at fixed volume (see Sec.
V.1) report ordering of polydisperse hard spheres into complex crystals with a
large unit cell, such as Laves and Frank-Kasper phases, instead of
fractionating into fcc crystals Lindquist _et al._ (2018); Bommineni _et al._
(2019). These large unit cells, which incorporate small and large particles,
are reminiscent of what is observed in binary mixtures (see Sec. IX).
Similarly, polydisperse mixtures with a non-Gaussian size distribution have
also been observed to partially crystallize into an AlB2 structure Coslovich
_et al._ (2018). The spontaneous formation of these phases suggests that their
formation might be thermodynamically preferred over fractionation into crystal
(see Fig. 10), but this proposal has yet to be investigated using free-energy
calculations. Although this phenomenon has not yet been observed in
experimental colloidal hard spheres either, polydisperse charged silica
nanospheres do form comparably complex crystal phases Cabane _et al._ (2016).
The question of the true equilibrium phase diagram of strongly polydisperse
hard sphere mixtures therefore remains open.

###  VII.3 Equation of state

Equations of state (EoS) are some the most fundamental descriptions of
equilibrium systems. On the theoretical and numerical side, high-accuracy EoS
are now available for both liquid and crystal phases of hard spheres Pieprzyk
_et al._ (2019). Such a work caps decades of systematic improvements (see
e.g., Alder _et al._ (1968); Speedy (1997, 1998b); Almarza (2009); Santos _et
al._ (2020) and Secs. VI and II). Despite the quantitative success of these
efforts, physical insight largely emerges from approximation schemes, which
also suffice as reference and calibration for colloid experiments. The
expression obtained by Carnahan and Starling (1969) for the fluid phase (Sec.
VI, Eq. (24)), in particular, has proven particularly useful. For the crystal,
fitted data to computer simulation results, such as Hall’s fit Hall (1972) and
the Speedy equation of state Speedy (1998a), play a similar role.

Figure 11: Experimental determination of the equation of state, i.e.,
$Z(\phi_{\textrm{eff}})=\beta P(\phi_{\textrm{eff}})/\rho_{\textrm{eff}}$
versus $\phi_{\textrm{eff}}$, for $\sigma=0.720\mu\mathrm{m}$ PMMA particles
from the observed sedimentation-diffusion equilibrium: (a) supernatant,
$\phi_{\textrm{eff}}\approx 0$; (b) fluid phase, $0<\phi_{\textrm{eff}}\leq
0.492$; (c) sharp interface; (d) crystalline phase,
$0.543\leq\phi_{\textrm{eff}}\leq 0.74$. Lines denote the fluid and solid EoS
for true hard spheres. The inset shows a photograph of 1 cm of the sample.
Adapted from Rutgers _et al._ (1996).

Early experimental attempts at determining the fluid EoS were carried out by
Vrij _et al._ (1983), and marked methodological improvements were later
achieved by Piazza _et al._ (1993). Although both approaches obtained the EoS
by integrating the equilibrium density profile in a single sedimentation
experiment, the later work used a charge–stabilized system of poly-tetrafluoro
ethylene (a polymer colloid similar in some respects to polystyrene but with a
rather lower refractive index) with a Debye length much smaller than the
particle diameter ($\kappa^{-1}=2$ nm vs. $\sigma=146$ nm), thus yielding a
good hard-sphere approximation (see Sec. III.4). These later measurements
compared well with the Carnahan–Starling expression for the fluid and
reasonably well to Hall’s fit to simulation data for the crystal Hall (1972).
As shown in Fig. 11, subsequent work by Rutgers _et al._ using X-ray
densiometry Rutgers _et al._ (1996) also found reasonable agreement with
theory for the fluid branch as well as for the solid branch for $\phi$
approaching close packing.

What is _reasonable_ agreement in this context? As discussed in Sec. III.5,
although colloidal systems are excellent model systems to demonstrate a range
of phenomena of condensed matter, highly accurate standard measures, they are
not, at least compared to atomic and molecular systems. While the accuracy
limit has not been systematically characterized, it seems reasonable that an
error in volume fraction of about $\delta\phi\sim 0.01$ should result from
measurements that base their phase boundaries on the method of Paulin and
Ackerson (1990) discussed in Sec. VII.1.

At low $\phi$, however, the results of Piazza _et al._ (1993) exhibited a
systematic drift far in excess of the expected error. The measured volume
fraction dropped much slower than anticipated as a function of height. In
other words, the sedimentation profile was _extended_. The physical mechanism
that underlies such extension was later identified by van Roij (2003). In
these systems, colloids and small ions decouple under a gravitational
field—the latter being only negligibly affected on the experimental length
scale. The resulting macroscopic electric field then partially counters
gravity, thus extending the sedimentation profile. Subsequent experiments
confirmed that interpretation Rasa and Phillipse (2004); Royall _et al._
(2005).

Particle-resolved studies have enabled new developments in EoS measurements.
For example, Dullens _et al._ (2006a) performed numerical Widom particle
insertion on experimental particle configurations to determine the chemical
potential of colloidal hard spheres under the assumption of a hard-sphere pair
potential. The resulting thermodynamic properties, including the equation of
state are quite good, but become numerically inaccessible at higher volume
fraction ($\phi\gtrsim 0.43$). The relation between free volume measurements
and free energy has separately been explored in the supercooled fluid regime
($\phi\geq 0.54$) with reasonable quantitative agreement Zargar _et al._
(2013); Dang _et al._ (2022) despite the crudeness of the assumptions
underlying cell theory in this regime (Sec. VI.2). Furthermore, as elegant as
this approach can seem, the effects of polydispersity and tracking errors can
hamper such analyses based on particle coordinates. An analysis based on
cavity averages which compared with simulation data found errors in the
pressure to be 100% and up to 10kBT in the chemical potential Schindler and
Maggs (2015).

###  VII.4 Structure and dynamics of the hard-sphere fluid

_Structure._ – As discussed in Sec. VI, the radial distribution function
$g(r)$ and structure factor $S(k)$ of the hard-sphere fluid can be obtained
from the Percus-Yevick closure to the Ornstein-Zernike equation. The resulting
theoretical predictions have been extensively tested against simulation data
(e.g., Ree _et al._ (1966); Frenkel _et al._ (1986); Hansen and McDonald
(2013)). Experimental studies of colloidal fluids using static light
scattering (SLS, see Sec. IV.1), also provide a natural access to the pair
structure. Early SLS experiments Vrij _et al._ (1983); de Kruif _et al._
(1985b) found that $S(k)$ for hard-sphere–like colloidal suspensions was
consistent with the Percus-Yevick expression in the small wave vector limit
$k\to 0$. Because $S(k\to 0)$ is directly linked to the isothermal
compressibility of the fluid under consideration, as described by Eq. (11),
this correspondence also experimentally validated the hard-sphere EoS.
Subsequent small-angle neutron scattering De Kruif _et al._ (1988) and light
scattering experiments Moussaïd and Pusey (1999) extended this correspondence
over a range of $k$ that included the first peak of $S(k)$ (see Fig. 12).

Figure 12: Reciprocal-space pair structure of hard sphere fluids determined
using static light scattering Moussaïd and Pusey (1999). Experimental data is
shown as points, and compared with Percus-Yevick predictions (solid line) at
the specified volume fractions.

Access to the real-space structure of bulk colloidal fluids, made possible
with the advent of confocal microscopy, led to the first experimental
measurements of $g(r)$ of hard spheres van Blaaderen and Wiltzius (1995). The
example in Fig. 13 compares experimentally measured $g(r)$ with EDMD
simulation results. Note, however, that the measurements, especially the first
peak of $g(r)$, are sensitive both to polydispersity and to tracking errors
Ivlev _et al._ (2012); Royall _et al._ (2007b); Mohanty _et al._ (2014). Both
effects can be explicitly taken into account in simulations, in order to
improve the match with in silico predictions of the experimental $g(r)$. The
resulting pair correlation functions match to a high degree of accuracy. Very
recently, a range of solvents has been employed to obtain a behavior
indistinguishable from hard spheres when real–space $g(r)$ comparisons were
carried out. That is to say, the effects of softness (see Sec. III.3) were
found to be negligible Kale _et al._ (2023).

Figure 13: Radial distribution function of hard-sphere fluids in real space at
various (effective) volume fractions. Points are experimental data, lines are
simulations of hard sphere fluids with the experimental polydispersity.
Adapted from Royall _et al._ (2018b), with the permission of AIP Publishing.

Many-body correlation functions and structural features are also easier to
extract from real-space than from reciprocal-space information Royall and
Williams (2015). In particular, it has been possible to compare experiments
Taffs _et al._ (2013), simulations and theoretical predictions Robinson _et
al._ (2019, 2020) for the occurrence of clusters that are known to minimise
packing constraints Manoharan _et al._ (2003); Robinson _et al._ (2019). As
demonstrated in Fig. 14, theoretical predictions of the population of such
clusters show close agreement between the various approaches.

Figure 14: Higher–order structure in the bulk hard-sphere fluid. Populations
of higher-order structures $N_{c}/N$ containing $m=3$ to $m=13$ particles,
where $N_{c}$ is the total number of clusters of a certain topology. Lines are
from morphometric theory Robinson _et al._ (2019, 2020). Small data points are
from monodisperse Monte-Carlo simulations, except for the purple line for
$m=13$ (with $s=8$% polydispersity). Large data points are from confocal
microscopy experiment using particles with $\sigma=2.0\mu$m. Simulations and
experiments are analyzed with the topological cluster classification using a
simple “bond length” (or pair distance) of $1.2\sigma$ Malins _et al._ (2013).
Figure 15: Self-diffusion coefficient of colloidal hard spheres. (a) Short-
time diffusion coefficient obtained from various sources: experiments Segre
_et al._ (1995), Eq. (31) by Lionberger and Russel (1994), theoretical
prediction based on up to three-body hydrodynamic interactions Cichocki _et
al._ (1999), force multipole simulations by Ladd (1990), accelerated Stokesian
dynamics (ASD) simulations Banchio and Nägele (2008), and theoretical
predictions based on the theory of Beenakker and Mazur Beenakker and Mazur
(1983); Banchio and Nägele (2008). Figure taken from Banchio and Nägele
(2008). (b) Long-time diffusion coefficient from various sources: experiments
with hard-sphere like suspensions by van Megen _et al._ (1987)
($\triangledown$), Segrè _et al._ (1995a) ($\diamond$), Van Blaaderen _et al._
(1992) ($+,\triangle$), Imhof and Dhont (1995b) ($\times$), Kops-Werkhoven and
Fijnaut (1982) ($\bullet$) as well as the theory of Cohen _et al._ (1998)
(solid line) and of Medina-Noyola (1988) (dash-dotted line). Adapted from
Cohen _et al._ (1998).

_Dynamics._ – At very short times (prior to colliding with a significant
number of solvent molecules), colloidal motion is ballistic Franosch _et al._
(2011); Hammond and Corwin (2017). At slightly longer times, however, hard-
sphere colloids behave diffusively Tough _et al._ (1986); colloidal dynamics
then stems from Brownian motion. Due to the frequent collisions of colloids
with solvent molecules, inertia can typically be neglected, thus making their
motion overdamped, as in Eq. (14).

Before particles have had the opportunity to interact with their neighbors,
their diffusivity can be described by a short-time diffusion coefficient
$D_{S}$. In the dilute limit $\phi\to 0$, neglecting hydrodynamic interactions
between colloids, the short-time diffusion coefficient coincides with the free
diffusion coefficient, $D_{S}=D_{0}$, which for a single colloidal sphere in a
solvent with viscosity $\eta$ is accurately described by the Stokes-Einstein
relation

$D_{0}=\frac{k_{\text{B}}T}{3\pi\eta\sigma}.$ (30)

At higher volume fractions, the short-time diffusion coefficient is reduced by
hydrodynamic interactions. This slowdown has been demonstrated in a variety of
experimental hard-sphere realizations Ottewill and Williams (1987); van Megen
_et al._ (1987); Qiu _et al._ (1990); Van Megen and Underwood (1990); Zhu _et
al._ (1992); Segre _et al._ (1995), and reproduced in simulations that
specifically account for such interactions Ladd (1990); Segre _et al._ (1995);
Banchio and Nägele (2008). Various theoretical and empirical descriptions of
the relation between $D_{0}$ and $D_{S}$ have also been proposed Beenakker and
Mazur (1983, 1984); Pusey and Tough (1983); Lekerkerker and Dhont (1984);
Cohen _et al._ (1998); Cichocki _et al._ (1999). A common semi-empirical
expression was introduced by Lionberger and Russel (1994)

$\frac{D_{S}}{D_{0}}\simeq(1-1.56\phi)(1-0.27\phi).$ (31)

Figure 15(a) compares experimental results of $\frac{D_{S}}{D_{0}}$ with Eq.
(31), as well as with several other theoretical and simulation approaches.
Clearly, different simulation methods (force multipole simulations by Ladd
(1990) and accelerated Stokesian dynamics by Banchio and Nägele (2008))
accurately describe the experimental data, while the theoretical predictions
based on Beenakker and Mazur (1983) agree well up to reasonably high $\phi$.
Good agreement with the same theory up to $\phi\simeq 0.4$ was also obtained
by the experimental work of Orsi _et al._ (2012). The wave vector dependence
of $D_{S}$ has also been investigated and theoretical predictions verified van
Megen _et al._ (1985); Segrè _et al._ (1995a).

At still longer times, interactions with neighboring particles hinder colloid
mobility. This crowding leads to a subdiffusive regime (as shown by the slope
of the mean squared displacement at intermediate times) before motion
ultimately becomes diffusive once more. This long-time diffusive behavior is
now described by a different diffusion coefficient $D_{L}<D_{S}$ Tough _et
al._ (1986). As discussed in Sec. XII, $D_{L}$ is strongly suppressed at very
high $\phi$, in the regime of glassy dynamics. Here, we focus on the dynamics
in the equilibrium fluid phase alone.

The long-time diffusion coefficient and other transport properties of pure
hard spheres (i.e. in the absence of hydrodynamics) have been extensively
explored by means of theory and simulations. For instance, early simulations
considered the agreement between low-density fluid results with predictions
from Enskog’s kinetic theory of gasses Chapman and Cowling (1990); Alder and
Wainwright (1967); Alder _et al._ (1970); Easteal _et al._ (1983); Speedy
(1987). When appropriately rescaled, these results are essentially independent
of the choice of microscopic dynamics (i.e. Monte Carlo, Brownian, or
Newtonian dynamics) Scala _et al._ (2007); Sanz and Marenduzzo (2010). As a
result, for pure hard spheres the choice of simulation methods typically does
not significantly impact dynamical behavior.

In experimental colloidal hard-sphere fluids, the long-time dynamics is
influenced by collisions with neighbors as well as by hydrodynamics. While
hydrodynamic interactions typically only lead to an overall (density-
dependent) scaling of the dynamical time scales for systems close to
equilibrium, they can significantly influence non-equilibrium behavior (see
Sec. XI and XIII.4). Experimental measurements of the diffusion coefficient in
a variety of systems van Megen _et al._ (1987); Segrè _et al._ (1995a); Van
Blaaderen _et al._ (1992); Imhof and Dhont (1995b); Kops-Werkhoven and Fijnaut
(1982) have shown broadly consistent results for $D_{L}/D_{0}$ as a function
of $\phi$ (see Fig. 15b). Theoretical approaches based on the generalized
Langevin equation Medina-Noyola (1988), on the expected scaling of diffusion
near random close packing Brady (1994), or on a mean-field description Cohen
_et al._ (1998) provide a reasonably good description of this behavior. It
should be noted that although the experimental data in Fig. 15(b) represent a
variety of systems with different polydispersities and electrostatic charge,
these effects are largely hidden by the experimental noise. Since both charge
and polydispersity shift the effective hard sphere volume fraction, it is
natural to expect a corresponding shift in the dynamics. This effect can be
more clearly observed when considering viscosity measurements Segrè _et al._
(1995a); van der Werff and de Kruif (1989); de Kruif _et al._ (1985a); Jones
_et al._ (1991); Papir and Krieger (1970); Imhof _et al._ (1994), and is
closely matched by the theoretical treatment of Cohen _et al._ (1998). Charged
colloids also exhibit the same scaling and it has been suggested that the
collapse may be more coincidence than something of fundamental significance.

The reciprocal-space dynamics of hard sphere has also been extensively
studied. Intriguingly, experiments by Segré et al. on sterically stabilized
PMMA particles found an unexpected collapse of the dynamical behavior for
large wave vectors: both the wave vector-dependent diffusion coefficient
$D(k,t)$ and the intermediate scattering function $F(k,t)$ collapsed onto a
single curve for $k\sigma\gtrsim 5$ Segre and Pusey (1996, 1997). Later work
on charge-stabilized polystyrene spheres, however, could not reproduce that
collapse Lurio _et al._ (2000), thus suggesting that residual charges may have
played a role in that study even if their effect can be hardly distinguished
in the static structure factor. This interpretation is further supported by
the fact that both mode-coupling theory calculations on pure hard spheres
Fuchs and Mayr (1999) and additional experiments on sterically stabilized PMMA
spheres Orsi _et al._ (2012) reproduce the approximate collapse of $F(k,t)$
observed by Segré and coworkers.

While not strictly pertaining to _colloidal_ hard spheres, because the work
was carried out using MD simulations of hard spheres in a vacuum, we note _en
passant_ that some attention has been given to the tail of the velocity
autocorrelation function. Around the freezing volume fraction, this observable
turns from positive to negative a back to positive Williams _et al._ (2006);
Martinez _et al._ (2014), as a result of collisional back-scattering in high-
density fluids. Unlike some authors, however, we conclude based on the nature
of the phase transition and its equivalence in systems with other microscopic
dynamics that the emergence of this dynamical phenomenon is merely
coincidental with the thermodynamic freezing point and not reflective of a
deeper physical relationship.

To summarize, we have a clear and robust theoretical understanding and
experimental validation of the diffusion of colloidal hard-sphere fluids at
low and intermediate densities.

###  VII.5 Bulk hard spheres in two dimensions

_Equilibrium phase behavior of hard disks._ – When monodisperse hard spheres
are confined to a strictly two-dimensional set-up, they behave as hard disks,
a model system that has itself long been the focus of theoretical and
computational study. Interestingly, the phase behavior of hard disks is
profoundly different to that of hard spheres. As shown in Fig. 16, a
qualitatively new phase – the hexatic phase – then appears between the fluid
and the solid (which, strictly speaking, is not a crystal as it lacks periodic
positional order due to Goldstone modes fluctuations at small wavevector $k$).

Figure 16: Hard-sphere phase behavior in 2d as a function of area fraction.
Close packing is indicated at $\phi_{\mathrm{cp}}=0.907\ldots$. The first-
order fluid-hexatic transition is marked by a narrow coexistence region of
$\phi\approx 0.68--0.70$. The hexatic-solid transition is second order at
$\phi\approx 0.73$ as indicated by the shading. Phase boundaries are taken
from Bernard and Krauth (2011).

In the perfect crystal, each particle has six neighbors and at close packing
$\phi_{\text{cp}}=\pi/2\sqrt{3}=0.907\ldots$ of the area is covered. The
hexagonal nature is captured through the extent of spatial correlations of the
bond-orientation order (BOO) parameter $\psi_{6}$, which takes a value of
unity for perfect hexagonal order. The hexatic phase is characterized by an
algebraic decay of orientational correlations, in contrast to the exponential
decay and incomplete decay (i.e. long-range order) of these same correlations
in the fluid and the solid, respectively. There is an intimate link between
the XY model in 2d, for which the Kosterlitz-Thouless-Halperin-Nelson-Young
(KTHNY) theory predicts two continuous phase transitions Kosterlitz and
Thouless (1973); Halperin and Nelson (1978); Young (1979) as first dislocation
pairs appear and in a second step these dislocation pairs unbind and
proliferate, thus destroying orientational order in the fluid.

How closely hard disks follow this scenario was debated for decades.
Simulations finally settled the question about a dozen years ago Bernard and
Krauth (2011); Engel _et al._ (2013). The key difficulty is that finite-size
effects cut off the algebraic decay of correlations. Very large system sizes
are therefore required to assess their relevance. This requirement is
compounded by the slow decay of structural correlations in the regime of
interest. This dual hurdle was first overcome thanks to event-chain MC
simulations (Sec. V.1). The results clearly revealed that the solid-hexatic
transition is continuous, but that the hexatic-fluid transition is weakly
first order with a narrow coexistence region $0.700<\phi<0.716$, wherein
domains of long-range orientational correlations coexist with short-range
orientational correlations, see Table 2. In order to relate these predictions
to experimental systems, the effect of out-of-plane fluctuations was
investigated using event-driven molecular dynamics simulations Qi _et al._
(2014). These simulations showed that the two-stage melting scenario of hard
disks persists even for monolayers of hard spheres with out-of-plane buckling
as high as half a particle diameter. These simulations have further shown that
the solid-hexatic transition is of the Kosterlitz-Thouless type and occurs via
dissociation of bound dislocation pairs, whereas the hexatic-liquid transition
is driven by a spontaneous proliferation of grain boundaries.

Table 2:  Table summarizing prior reports on the phase behavior of hard disks in 2D using computer simulations. For each result, we list the area fractions of the coexisting fluid $\phi_{\mathrm{f}}$ and hexatic $\phi_{\mathrm{h}}$ phases, the area fraction at the hexatic-solid transition $\phi_{\mathrm{hs}}$, the coexistence pressure at the fluid-hexatic phase transition $\beta P^{\mathrm{fh}}_{\mathrm{coex}}\sigma^{2}$ as well as the system size considered. Note that the area fraction is defined as $\phi=\pi\sigma^{2}N/4A$ with $A$ the system area. Note also that the accuracy is to the last reported digit.  Source | $\phi_{\mathrm{f}}$ | $\phi_{\mathrm{h}}$ | $\phi_{\mathrm{hs}}$ | $\beta P^{\mathrm{fh}}_{\mathrm{coex}}\sigma^{2}$ | $N$
---|---|---|---|---|---
Bernard and Krauth (2011) | 0.700 | 0.716 | 0.720 | 9.185 | $1024^{2}$
Engel _et al._ (2013) | 0.700 | 0.716 | 0.720 | 9.185 | $1024^{2}$
Qi _et al._ (2014) | 0.699 | 0.717 | 0.725 | 9.183 | $1024^{2}$

Early experimental work aiming at elucidating the nature of the melting
transition of colloidal monolayers of PMMA spheres suggested a first-order
phase transition between the fluid and hexatic phases as well as a first-order
phase transition between the hexatic and solid phase Marcus and Rice (1997).
However, the authors argued that this observation could be linked to (non-
hard) attractive or repulsive interactions between particles, consistent with
an earlier simulation study Bladon and Frenkel (1995).

A definitive experimental elucidation of the nature of the hard-disk phase
diagram was later provided by experiments on a tilted monolayer of colloidal
hard spheres in sedimentation-diffusion equilibrium (Fig. 17) Thorneywork _et
al._ (2017b). By measuring the density profile at different tilts of the
plane, the full equation of state was extracted via the barometric law,
unveiling a discontinuity at $\phi=0.68$. By calculating the height-resolved
BOO correlation function in time and a modified Lindemann parameter, this
discontinuity was identified as a first-order fluid-hexatic phase transition,
and a study of the thermal fluctuations of the interface identified the width
of the fluid-hexatic coexistence region as $\Delta\phi\simeq 0.02$. In
addition, no finite discontinuities were found between the hexatic-solid
transition. In short, a reasonably good agreement with the simulation results
was observed (with essentially the same level of accuracy as that noted above
in Sec. VII.3). Later direct comparisons of correlations in translational and
bond-orientational orders in the various phases similarly demonstrated
excellent agreement between this experimental model system and simulations
Thorneywork _et al._ (2018a).

Figure 17: Experimental determination of the hard-sphere phase behavior in 2d.
(a) Experimental image of quasi-2d hard spheres in sedimentation-diffusion
equilibrium under a slight tilt. (b) The equation of state $P/\rho
k_{\text{B}}T$. The inset shows an expanded view of EoS around the
discontinuity. The solid red line is a prediction of scaled particle theory
for the fluid regime (Eq. 32). The solid blue line shows is a semiempirical
fit $P/\rho k_{\text{B}}T=a/(\phi_{\mathrm{cp}}-\phi)$. Reproduced from
Thorneywork _et al._ (2017b)

Interestingly, whether the hexatic phase is present or not depends strongly on
the details of the system. For instance, the fluid-hexatic phase transition
becomes metastable with respect to a first-order fluid-solid transition for
binary mixtures of large (L) and small (S) disks with diameter ratio
$q=\sigma_{L}/\sigma_{S}=1.4$, for molar fractions of small disks as small as
$1\%$ Russo and Wilding (2017); Russo _et al._ (2018). The two-stage melting
scenario of a continuous solid-hexatic and a first-order fluid-hexatic
transition also becomes metastable with respect to a first-order fluid-solid
transition for hard disks with less than $1\%$ of the particles pinned to a
triangular lattice Qi and Dijkstra (2015). These findings are corroborated
with a renormalization group analysis based on the KTHNY theory, which shows
that the renormalized Young’s modulus of the crystal increases in the presence
of pinned particles Qi and Dijkstra (2015). The melting transition was shown
to be qualitatively changed to the KTHNY scenario for polydisperse disks
Sampedro Ruiz _et al._ (2019), for which a significantly enlarged stability
range for the hexatic phase and re-entrant melting have been observed.

_Equation of state._ – The EoS of hard disks has been explored numerically
since the very introduction of the Metropolis algorithm in 1953 Metropolis _et
al._ (1953), and has been studied extensively ever since (see e.g. Li _et al._
(2022)). Unlike 3d hard spheres, however, the Percus-Yevick closure to the
Ornstein-Zernike equation cannot be solved analytically in 2d, thus
necessitating the use of numerical integration. As a result, a large number of
different approximations for the equation of state has been proposed (see e.g.
Helfand _et al._ (1961); Henderson (1975); Tejero and Cuesta (1993); Santos
_et al._ (1995); Mulero _et al._ (2009); Boublík (2011)). Scaled particle
theory offers a particularly simple expression Helfand _et al._ (1961),

$\frac{\beta P}{\rho}=\frac{1}{(1-\phi)^{2}},$ (32)

which agrees reasonably well with experimental results (see Fig. 17(b)). When
more accurate results are required, a common approach is to use the polynomial
fits of simulation data provided by Kolafa and Rottner (2006).

The experimental equation of state has been measured by Brunner _et al._
(2003) in the fluid regime and by Thorneywork _et al._ (2017b) for the whole
$\phi$ regime (see Fig. 17). Comparison using test-particle insertion agrees
remarkably well with simulations Stones _et al._ (2018). A very good agreement
between theory and simulations has also been reported for $g(r)$ at various
$\phi$ Marcus and Rice (1997); Santana-Solano and Arauz-Lara (2001); Brunner
_et al._ (2003); Thorneywork _et al._ (2014).

_Dynamics of the hard-disk fluid._ – The dynamics of 2d confined hard-sphere
colloids has been explored experimentally, with early work investigating the
decay of the short-time diffusion to its long-time limit Marcus _et al._
(1996). This study was followed by explorations of hydrodynamic effects on the
short–time dynamics and related these to the 3d system Carbajal-Tinoco _et
al._ (1997); Santana-Solano and Arauz-Lara (2001); Bonilla-Capilla _et al._
(2004).

Another way to form a 2d colloidal system is for particles to adhere to an
interface between two liquids. The dynamics in such a system provides a good
estimate of the drag coefficient of the particles Peng _et al._ (2009). Other
work has emphasised the importance of hydrodynamic interactions at short times
Thorneywork _et al._ (2015). Surprisingly, while higher $\phi$ lead to a lower
short-time diffusion coefficient $D_{S}$, as expected, the behavior of the
long-time diffusion coefficient $D_{L}$ agrees well with simulation results of
a pure hard-disk system. Hence, although hydrodynamics slows down the
individual motion of particles at short times, its effect is largely
compensated at longer time scales. Recently, the degree of confinement to a
plane has also been investigated using simulations incorporating hydrodynamic
interactions. The long-time diffusion was then found to be enhanced in weaker
confinement Tian _et al._ (2022).

##  VIII Interfaces in hard spheres

Hard-sphere systems spontaneously form interfaces as a result of fluid-crystal
phase separation as well as grain boundaries in the crystal phase. Interfacial
effects can also arise by bringing hard spheres in contact with a boundary.
This section will delve into these various phenomena, starting with the
conceptually important case of a single hard wall.

###  VIII.1 Fluid in contact with a wall

_Planar wall._ – A bulk fluid of hard spheres next to a single hard wall can
be considered as the simplest inhomogeneous system. Intuitively, the presence
of a flat wall induces some local structuring of the fluid, which in turn
promotes the formation of the crystal phase. Early simulations by Courtemanche
and van Swol (1992); Courtemanche _et al._ (1993) indeed demonstrated that
crystalline layers could spontaneously form near a planar hard wall, even at
pressures below bulk freezing. However, it is important to note that in
simulations, because of periodic boundary conditions, a single wall is
essentially a confined system with two walls at a plate separation determined
by the system size. Therefore, complete wetting of a hard-sphere crystal at a
single wall, also termed pre-freezing, is challenging to distinguish from
capillary freezing Kegel (2001); Groot _et al._ (1987); Dijkstra (2004).
Careful simulation studies have now demonstrated that (_i_) a (111) hard-
sphere crystal completely wets a hard wall-fluid interface at about $98.3\%$
of the bulk freezing density, via interfacial free-energy calculations Laird
and Davidchack (2007), and (_ii_) a crystalline film at a single wall,
independent of the plate separation, grows logarithmically, and is clearly
distinguishable from capillary freezing Dijkstra (2004).

The interface between a hard-sphere fluid and a hard wall also gives rise to
an interfacial free-energy $\gamma$, which depends on the fluid volume
fraction $\phi$. Although a virial-like expansion for this quantity is
possible, too few coefficients are known for it to be of much practical use
Charbonneau _et al._ (2010). Davidchack _et al._ (2016) numerically studied a
hard-sphere fluid close to a planar wall, and computed the density profiles,
$\gamma$ and the excess adsorption by means of computer simulations and
density functional theory (DFT) (see Sec. VI.3). They detected systematic
deviations between numerical and DFT results for $\phi>0.3$ (see Fig. 18).
Similar deviations were previously reported in a study of confined hard-sphere
fluids Deb _et al._ (2011), which suggested that nontrivial correlations
between fluid particles beyond nearest-neighbor shell might develop at higher
$\phi$. These correlations could not, however, be accurately captured by DFT.

Figure 18: Interfacial free energy of a hard-sphere fluid in contact with a
flat hard wall, as a function $\phi$. The lines are theoretical predictions
based on fundamental measure theory (FMT) and scaled particle theory (SPT).
The points are simulation data. Simulation and FMT data taken from Davidchack
_et al._ (2016); the SPT expression is based on Heni and Löwen (1999).

In experiments Hoogenboom _et al._ (2003b), the presence of a bottom wall was
found to initially induce layering in a sedimenting colloidal fluid; as
sediment thickness increased, crystallization occurred within these fluid
layers. Crystallization in the first layer appeared to proceed via a first-
order transition consistent with predictions from BD simulation and theory
Biben _et al._ (1994). Using grand-canonical Monte Carlo simulations including
a gravitational field, it was later shown that the fluid crystallized via a
first-order freezing transition at which several fluid layers close to the
bottom of the sample froze at the same chemical potential Marechal and
Dijkstra (2007). The number of such layers simultaneously freezing decreases
for higher gravitational field strengths.

The presence of a wall can also affect colloid dynamics. Evanescent dynamic
light scattering experiments have studied sterically stabilized PMMA particles
for various volume fractions over a range of scattering wave vectors. In the
dilute regime, Michailidou _et al._ (2009) observed that the near wall short-
time diffusion was slowed down due to particle-wall hydrodynamic interactions.
However, for a concentrated suspension, the wall effect progressively
diminished at all vectors $k$ and many-body hydrodynamic interactions became
less relevant Liu _et al._ (2015b); Michailidou _et al._ (2009).

_Patterned wall._ – The wetting behavior of the hard-sphere crystal has also
been studied for patterned walls Heni and Löwen (2000). In the case of a (111)
structured surface, complete wetting of the hard-sphere crystal sets in
already at $29\%$ of the bulk freezing density. Even crystal structures that
are unstable in bulk can be promoted by the wall surface pattern. As an
example, hcp (rather than fcc) has been epitaxially grown on a structured
template Hoogenboom _et al._ (2003a). By contrast, surface wetting by the
crystal below freezing can be completely suppressed by wall patterns that are
incommensurate with the equilibrium crystal structure Espinosa _et al._
(2019).

###  VIII.2 Fluid-solid interface

At coexistence, the fluid phase is separated from the ordered solid by a
(thermally fluctuating) fluid-solid interface. Understanding the properties of
this interface on a microscopic scale is fundamental to assess crystal
nucleation (see Sec. XIII). For a planar fluid-crystal interface at
coexistence, one can define the interfacial (or surface) free energy $\gamma$
as the (reversible) work needed to form a unit area of a flat interface. The
quantity $\gamma$ generally depends on the orientation of the interface normal
($\widehat{n}$) relative to the crystalline axes. When considering the
fluctuations of that same interface, the quantity of interest is instead the
interfacial stiffness $\tilde{\gamma}$. We here discuss efforts to quantify
both $\gamma$ and $\tilde{\gamma}$.

_Interfacial free energy._ – In a pioneering numerical simulation, Davidchack
and Laird (1998) characterized the face-centered cubic (100) and (111) fluid-
solid interfaces, thus demonstrating that the transition from crystal to fluid
occurred over few crystal planes made of domains of coexisting crystal and
fluid phases. They later extended the method developed by Broughton and Gilmer
Broughton and Gilmer (1986), based on a thermodynamic integration along a
reversible path defined by cleaving, to determine $\gamma$ for hard-sphere
fluids, whose value depended on the crystalline lattice structure Davidchack
and Laird (2000). Their results were then sufficiently precise to assess the
anisotropy of the fluid-solid interfacial free energy. This work was followed
by more simulation-based efforts to determine $\gamma$, using a wide variety
of approaches including the analysis of interfacial fluctuations using
capillary wave theory Mu _et al._ (2005); Davidchack _et al._ (2006); Härtel
_et al._ (2012), non-equilibrium work methods Davidchack (2010), tethered
Monte Carlo Fernández _et al._ (2012), classical nucleation theory Cacciuto
_et al._ (2003), mold integration Espinosa _et al._ (2014) and thermodynamic
integration Bültmann and Schilling (2020); Benjamin and Horbach (2015) (see
also Sec. V). Härtel _et al._ (2012) additionally used fundamental measure
theory (FMT, see Sec. VI), although this approach resulted in significantly
higher values than in simulations. An overview of results for $\gamma$ is
provided in Table 3, listing results for both individual interfaces and an
orientationally averaged $\bar{\gamma}$, which is a key quantity in crystal
nucleation (see Sec. XIII). Clearly, significant variations exist between
different methods, especially when compared to the relatively small error bars
reported. Schmitz _et al._ (2014); Schmitz and Virnau (2015) suggested that
systematic errors (related to finite-size effects) could be a possible
explanation for this disagreement. Benjamin and Horbach (2015) proposed a
theoretical tool to obtain reliable estimates for $\gamma$ in the
thermodynamic limit.

Table 3:  Interfacial free-energy $\gamma$ and stiffness $\tilde{\gamma}$, where ($hjk$) are the Miller indices of the interfacial crystal plane. For the stiffness we might also indicate the symmetry of the short (in plane) direction. $\bar{\gamma}$ is averaged over the (100), (110) and (111) (of an fcc crystal). Various approached have been employed for the determination: cleaving [CL], capillary wave [CW], non-equilibrium work [NEW], mold integration [MI], tethered Monte Carlo [TET], CNT analysis of the nucleation rate [CNT], thermodynamic integration [TI], ensemble switch [SW], fundamental measure DFT [FMT] (theory), and  experiments $[]_{\mathrm{exp}}$. Interfacial free energy ($\gamma$) |  Interfacial stiffness ($\tilde{\gamma}$) |
---|---|---
$\gamma$(111)$[k_{\text{B}}T/\sigma^{2}]$ |  $\tilde{\gamma}$(111)$[k_{\text{B}}T/\sigma^{2}]$ |
0.58 $\pm$ 0.01 Davidchack and Laird (2000) [CL] | 0.80 [$\bar{1}10$] Mu _et al._ (2005) [CW] |
0.61 $\pm$ 0.02 Mu _et al._ (2005) [CW] | 0.78 $\pm$ 0.04 [$\bar{1}10$] Davidchack (2010) [NEW] |
0.546 $\pm$ 0.016 Davidchack _et al._ (2006) [CW] |  1.025 $\pm$ 0.079  Härtel _et al._ (2012) [FMT] |
0.5416 $\pm$ 0.0031 Davidchack (2010) [NEW] | 0.810 $\pm$ 0.005 Härtel _et al._ (2012) [CW] |
0.636 $\pm$ 0.001 Härtel _et al._ (2012) [FMT] |  0.41 van Loenen _et al._ (2019) [CW]exp |
0.600$\pm$0.011 Härtel _et al._ (2012) [CW] |  |
0.556$\pm$ 0.003 Benjamin and Horbach (2015) [TI] |  |
0.544$\pm$ 0.008 Schmitz and Virnau (2015) [SW] |  |
0.554$\pm$0.006 Sanchez-Burgos _et al._ (2021) [MI] |  |
$\gamma$(100)$[k_{\text{B}}T/\sigma^{2}]$ |  $\tilde{\gamma}$(100)$[k_{\text{B}}T/\sigma^{2}]$ |
0.62 $\pm$ 0.01 Davidchack and Laird (2000) [CL] | 0.55 [$001$] Mu _et al._ (2005) [CW] |
0.64 $\pm$ 0.02 Mu _et al._ (2005) [CW] | 0.44 $\pm$ 0.03 [$001$] Davidchack _et al._ (2006) [CW] |
0.574$\pm$ 0.017 Davidchack _et al._ (2006) [CW] | 0.40 $\pm$ 0.02 [$001$] Davidchack (2010) [NEW] |
0.582 $\pm$ 0.002 Davidchack (2010) [NEW] | 0.49 $\pm$ 0.02 [$001$] Zykova-Timan _et al._ (2010) [CW] |
0.636 $\pm$ 0.011 Fernández _et al._ (2012) [TET] |  0.53$\pm$ 14 [$001$] Härtel _et al._ (2012) [FMT] |
0.687 $\pm$ 0.001 Härtel _et al._ (2012) [FMT] | 0.419$\pm$ 0.005 [$001$] Härtel _et al._ (2012) [CW] |
0.639 $\pm$ 0.0011 Härtel _et al._ (2012) [CW] |  |
0.586 $\pm$ 0.008 Espinosa _et al._ (2014) [MI] |  |
0.596$\pm$ 0.006 Benjamin and Horbach (2015) [TI] |  1.2 Hernandez-Guzman and Weeks (2009) [CW]exp |
0.581$\pm$ 0.003 Schmitz and Virnau (2015) [SW] |  1.3 $\pm$ 0.3 Ramsteiner _et al._ (2010) [CW]exp |
0.589 $\pm$ 0.001 Bültmann and Schilling (2020) [TI] |  0.47 van Loenen _et al._ (2019) [CW]exp |
0.586 $\pm$ 0.006 Sanchez-Burgos _et al._ (2021) [MI] |  |
$\gamma$(110)$[k_{\text{B}}T/\sigma^{2}]$ |  $\tilde{\gamma}$(110)$[k_{\text{B}}T/\sigma^{2}]$ |
0.64 $\pm$ 0.01 Davidchack and Laird (2000) [CL] | 0.71 [$001$] Mu _et al._ (2005) [CW] |
0.62 $\pm$ 0.02 Mu _et al._ (2005) [CW] | 0.49 [$\bar{1}10$] Mu _et al._ (2005) [CW] |
0.557$\pm$ 0.017 Davidchack _et al._ (2006) [CW] | 0.42 $\pm$ 0.03 [$\bar{1}10$] Davidchack _et al._ (2006) [CW] |
0.559 $\pm$ 0.002 Davidchack (2010) [NEW] | 0.70 $\pm$ 0.03 [$001$] Davidchack _et al._ (2006) [CW] |
0.616 $\pm$ 0.011  Härtel _et al._ (2012) [FMT] | 0.34 $\pm$ 0.05 [$\bar{1}10$] Davidchack (2010) [NEW] |
0.600$\pm$0.011 Härtel _et al._ (2012) [CW] | 0.74 $\pm$ 0.02 [$001$] Davidchack (2010) [NEW] |
0.577 $\pm$ 0.004 Benjamin and Horbach (2015) [TI] |  0.283 $\pm$ 0.035 [$\bar{1}00$] Härtel _et al._ (2012) [FMT] |
0.559 $\pm$ 0.001 Schmitz and Virnau (2015) [SW] | 0.401 $\pm$ 0.005 [$\bar{1}00$] Härtel _et al._ (2012) [CW] |
0.572 $\pm$ 0.007 Sanchez-Burgos _et al._ (2021) [MI] |  0.86 $\pm$ 0.14 [01$\bar{1}$] Härtel _et al._ (2012) [FMT] |
| 0.769 $\pm$ 0.005 [01$\bar{1}$] Härtel _et al._ (2012) [CW] |
|  1.0 $\pm$ 0.2 Ramsteiner _et al._ (2010) [CW]exp |
|  0.53 van Loenen _et al._ (2019) [CW]exp |
$\bar{\gamma}$ $[k_{\text{B}}T/\sigma^{2}]$ |  |
0.617 $\pm$ 0.006 Davidchack and Laird (2000) [CL] |  |
0.27 Gasser (2001) [CNT]exp |  |
0.616 $\pm$ 0.003 Cacciuto _et al._ (2003) [CNT] |  |
0.62 $\pm$ 0.02 Mu _et al._ (2005) [CW] |  |
0.56 $\pm$ 0.01 Davidchack (2010) [NEW] |  |
0.66 $\pm$ 0.01  Härtel _et al._ (2012) [FMT] |  |
0.5916 Sanchez-Burgos _et al._ (2021) [CNT] |  |

To establish whether the fcc or hcp crystal is favored due to a difference in
interfacial free energy, Sanchez-Burgos _et al._ (2021) computed $\gamma$ by
means of the mold integration method considering both ordered phases. The
authors concluded that the $(11\bar{2}0)$ plane of the hcp crystal had a
slightly higher $\gamma$ than any of the faces of the fcc polytype. They also
found a higher orientationally averaged $\bar{\gamma}$ for the hcp phase when
considering spherical nuclei using a seeding approach (see Fig. 19).

Figure 19: Interfacial free energy $\gamma$ for different polytypes at
coexistence. For planar interfaces (at the bulk melting pressure), an average
$\bar{\gamma}$ was obtained by averaging over different crystal orientations
for both polymorphs. The values at higher pressures correspond to
$\bar{\gamma}$ estimated for spherical nuclei (from the seeding technique),
using CNT Sanchez-Burgos _et al._ (2021).

Only more recently have experiments based on colloidal suspensions studied the
phenomenon via confocal microscopy Gasser (2001); Dullens _et al._ (2006b);
Hernandez-Guzman and Weeks (2009). Experimentally, $\gamma$ can be obtained
from indirect measurements of crystal nucleation rates using classical
nucleation theory (see Sec. XIII). Using this method on PMMA colloids, Gasser
(2001) obtained a surprisingly low $\bar{\gamma}$ (see Table 3), possibly due
to the colloids being softened by electrostatic interactions (see Sec. III.3).

_Interfacial stiffness._ – If the fluid-crystal interface is rough, capillary
wave theory (CWT) can be applied to define the interfacial stiffness
$\tilde{\gamma}$, that controls the resulting capillary waves. The tensorial
expression for $\tilde{\gamma}$ is then

$\tilde{\gamma}_{\alpha\beta}(\widehat{n})=\gamma(\widehat{n})+\frac{\partial^{2}\gamma(\widehat{n})}{\partial\widehat{n}_{\alpha}\partial\widehat{n}_{\alpha}},$
(33)

with $\widehat{n}_{\alpha}$ and $\widehat{n}_{\beta}$ two directions
orthogonal to $\widehat{n}$ Fisher and Weeks (1983). While the first
contribution ($\gamma(\widehat{n})$) describes the free-energy cost of
increasing the interfacial area, the second contribution accounts for the free
energy required to locally change the crystal orientation. Note that the
applicability of CWT requires the interface to be rough. This condition is
only valid above the (system-dependent) roughening transition temperature for
thermal systems Zykova-Timan _et al._ (2010); Schmid and Binder (1992), but is
always so for hard spheres.

If the anisotropy of the interfacial free energy is known, $\tilde{\gamma}$
can be calculated directly from Eq. (33), as in Härtel _et al._ (2012).
Alternatively, it can be obtained from interfacial fluctuations measured in
simulations using capillary wave theory Davidchack _et al._ (2006); Mu _et
al._ (2005); Härtel _et al._ (2012); Zykova-Timan _et al._ (2009, 2010).
Values from various different approaches are given in Table 3. Overall,
$\tilde{\gamma}$ shows significantly more anisotropy than $\gamma$, but there
are also significant discrepancies between different measurements.

Hernandez-Guzman and Weeks (2009) used confocal microscopy to study the
interfacial fluctuations of a solid-fluid interface and determined the
interfacial stiffness using capillary wave theory. They obtained a value about
twice as high as that obtained from computer simulations (see Table 3). This
discrepancy might be caused by slightly charged particles, but that remains to
be verified.

Figure 20: Rendering of particle coordinates determined from confocal
microscopy showing the fluid-crystal interface. Shown is a slice that is two
crystalline layers thick. Particles are colored according to the number of
ordered neighbors they have Hernandez-Guzman and Weeks (2009).

Similar values of the interfacial stiffness were obtained by Ramsteiner _et
al._ (2010). After sedimenting hard-sphere silica colloids onto (100) and
(110) oriented templates, they located the interface by confocal microscopy,
and again used capillary wave theory to determine $\tilde{\gamma}$. For all
three main crystal orientations, they noted that the Fourier amplitudes are
independent of the in-plane direction of the associated wave vectors in the
long wavelength limit. This result directly contradicts simulation results for
the (110) interface, from which a difference consistent with the stiffness
tensor of a cubic crystal is expected. Ramsteiner _et al._ (2010) suggested
that the relatively thin fluid layer and the small gravitational length could
explain this discrepancy.

Recently, van Loenen _et al._ (2019) studied the interfaces between hard-
sphere colloidal fluids and fcc crystals sedimented onto differently oriented
templates. The $\tilde{\gamma}$ they obtained from the capillary waves were
not only much closer to the simulation results than earlier estimates by
Hernandez-Guzman and Weeks (2009), but their order also was more consistent
with simulation predictions (see Table 3).

As these different results show, experimental measurements of the stiffness of
nominally the same system have yielded different results. A very recent study
by MacDowell (2023) suggests that the surface stiffness of colloidal hard
spheres, as measured by interface fluctuations, might include a gravity-
dependent contribution which could account for this discrepancy. The external
field dependence of the interfacial stiffness could then be explained by an
improved interface Hamiltonian that provides corrections to the capillarity
theory equation, thus reconciling experimental and theoretical results.

###  VIII.3 Grain boundaries

Figure 21: Grain boundaries in a polycrystalline 2d solid of colloidal hard
spheres. The left panel shows the particles (scale bar is 20$\mu$m). The right
panel shows the Voronoi cells of the particles, colored by the local
orientation of the BOO parameter, as encoded in the color bar, which allows
for easy distinction between the different domains. Figure adapted from
Lavergne _et al._ (2017).

Grain boundaries form spontaneously in crystalline materials at finite
temperature or, in the case of hard spheres, for volume fractions below close
packing, i.e., $\phi<\phi_{\mathrm{cp}}$. In a polycrystalline material,
crystalline grains with differing orientations are separated by an interface
made of amorphous grain boundaries. Most of the work on grain boundaries has
been performed in 2d Zhang _et al._ (2009); Gray _et al._ (2015); Lavergne _et
al._ (2015). Grain boundaries and other defects reduce long-range order of a
single crystal Zhang _et al._ (2009), deforming the underlying hexagonal order
(in 2d) due to a preference for five- and seven-fold coordinated particles to
be adjacent to impurities Gray _et al._ (2015).

Considering a polycrystalline monolayer of colloidal hard spheres (see Fig.
21), Lavergne _et al._ (2017) followed the grain growth process and detected
an anomalous slow growth of the BOO correlation length. The authors invoked
the curvature-driven coarsening of the large-angle grain boundaries at a rate
dependent on the grain boundary length to explain the effect. When dealing
with a polycrystalline monolayer of hard spheres with embedded impurities,
Lavergne _et al._ (2016) demonstrated that the size of the impurities
determined whether they behaved as interstitial or substitutional impurities
in the bulk crystal. Once formed, grain boundaries can also shrink.
Spontaneous shrinkage of circular grain boundaries has been studied in 2d
colloidal crystals by Lavergne _et al._ (2019), who demonstrated that the
shrinkage can be driven by three mechanisms: curvature-driven migration,
coupled grain boundary migration, and grain boundary sliding.

Grain boundary dynamics are driven by their local curvature and have been
thoroughly studied by means of point sampled surface analysis techniques
Lavergne _et al._ (2015). While the structure on length scales larger than the
grain boundaries distances strongly depends on the defects concentration,
local structural distortions close to a grain boundary occur only over short
distances compared to the grain boundary size, and are independent of the
defect concentration Gray _et al._ (2015). The kinetics of grain boundaries is
closely related to the topological constraints imposed on their dislocation
structure. As an example, Lavergne _et al._ (2018) showed that a local
rotational deformation of a 2d colloidal crystal with an optical vortex might
originate a grain boundary loop, thus underlining the relevance of defects in
the kinetics of grain boundaries.

In 3d, studies of hard-sphere grain boundaries are much more sparse. Recently,
Orr _et al._ (2021) introduced a new methodology to detect and characterize
grain boundaries in experimental data on colloidal grain boundaries, revealing
detailed misorientation distributions and grain boundary structures, but much
remains to be done.

##  IX Binary hard-sphere mixtures

Figure 22: Phase diagrams of binary hard-sphere mixtures displaying a spindle-
like phase behavior at size ratio $q=0.95$, an azeotropic phase diagram at
$q=0.92$, an azeotropic phase behavior with a solid-solid demixing transition
at $q=0.90$, and an eutectic phase diagram at $q=0.85$ as obtained from
simulations. Data taken from Kranendonk and Frenkel (1991).

The phase behavior of binary mixtures of large (L) and small (S) hard spheres
with diameters $\sigma_{\mathrm{L}}$ and $\sigma_{\mathrm{S}}$, respectively,
is remarkably rich. Even at close packing, complex binary crystals emerge as
the size ratio $q=\sigma_{\mathrm{S}}/\sigma_{\mathrm{L}}$ changes Dijkstra
(2014); Hopkins _et al._ (2012); Hudson and Harrowell (2011). In this section,
we review some of the salient thermodynamic features of these systems,
including the diversity of binary crystal structures they can form, the
structural crossover in binary fluids, fluid-fluid demixing and the quest for
its critical point as well as the behavior of sedimented monolayers of binary
hard-sphere mixtures.

###  IX.1 Crystal regime

In binary mixtures of hard spheres of two very similar sizes, the most stable
state at infinite pressure is phase separated into two separate fcc crystals:
one of large spheres and one of small spheres. Density functional theory,
computer simulations, and scaled particle theory have revealed further
enrichment of the crystal regime at finite pressure, away from close packing.
In the limit $q=\sigma_{\mathrm{S}}/\sigma_{\mathrm{L}}\rightarrow 1$, the
system reduces to a one-component hard-sphere system for which the phase
behavior is discussed in Sec. VII.1. As the two components become more
dissimilar in size, i.e., $q\lesssim 1$, the freezing transition changes from
a spindle-like via an azeotrope to an eutectic phase diagram Barrat _et al._
(1986, 1987); Zeng and Oxtoby (1990); Denton and Ashcroft (1990); Kranendonk
and Frenkel (1991); Kofke (1991); Cottin and Monson (1997), as shown in Fig.
22.

In the case of spindle-like phase behavior, the coexistence between a fluid
and a substitutionally disordered fcc crystal phase is narrow, in the sense
that only a small composition difference between the two phases develops. When
spheres become more dissimilar in size, the fluid-solid region broadens and an
azeotropic point appears at around $q=0.94$ Kranendonk and Frenkel (1991). At
higher packing fractions a coexistence region between two substitutionally
disordered fcc solids appears when the spheres become sufficiently dissimilar,
which shifts to lower densities as the size ratio is further increased. When
this miscibility gap in the solid phase interferes with the fluid-crystal
coexistence, the phase diagram becomes eutectic at $q\simeq 0.875$ Kranendonk
and Frenkel (1991). Note that only substitutionally disordered fcc phases
appear for the size ratio range $q\in(0.85,1)$, in contrast to what happens at
smaller $q$.

The phase behavior of more asymmetric mixtures was studied in experiments
using light scattering Bartlett _et al._ (1990, 1992); Pusey _et al._ (1994).
For size ratios $q=0.58$ and $0.62$, complex binary LS2 and LS13 superlattice
structures were found for sterically stabilized PMMA spheres Bartlett _et al._
(1990, 1992). The stability of these superlattice structures, which are
analogous to their atomic counterparts, AlB2 and NaZn13 as shown in Fig. 23,
was subsequently confirmed by computer simulations for size ratios
$0.42<q<0.625$ Eldridge _et al._ (1993b, c, a); Eldridge _et al._ (1995);
Dijkstra (2014) and by density functional approaches Xu and Baus (1992).

For size ratios $0.2<q<0.42$, a superlattice structure isostructural to the
rocksalt NaCl was predicted by computer simulations Trizac _et al._ (1997),
and verified experimentally Hunt _et al._ (2000); Vermolen _et al._ (2009)
(see Fig. 23). Surprisingly, the experimental observations of NaCl always
showed a crystal phase with many vacancies in the sublattice of the small
spheres in such a way that the L:S stoichiometry of the crystal is not 1:1 but
rather 1:$x$ with $x\leq 1$ Hunt _et al._ (2000); Vermolen _et al._ (2009).
Free-energy calculations from Monte Carlo simulations subsequently showed that
an interstitial solid solution constructed by filling the octahedral holes of
an fcc crystal of large spheres with small spheres is indeed the stable phase
Filion _et al._ (2011a). Upon increasing pressure, the fraction of octahedral
holes filled with a small sphere can be completely tuned from
zero—corresponding to a fcc phase of pure large spheres—to one—corresponding
to a NaCl phase. For less assymetric size ratios, non-equilibrium interstitial
solid solutions are found, which are long–lived on the experimental timescale
Rios de Anda _et al._ (2017).

For larger size ratios, another set of LS2 phases, known as Laves phases, were
found in computer simulations Hynninen _et al._ (2007, 2009); Bommineni _et
al._ (2020). Three Laves phase structures exist: hexagonal MgZn2, cubic MgCu2,
and hexagonal MgNi2. Each is characterized by the stacking of large-sphere
dimers in the crystal structure, as shown in Fig. 23. The MgCu2 phase is
particularly interesting because it consists of a diamond lattice of large
(Mg) spheres and a pyrochlore lattice of small (Cu) spheres. By selectively
removing one of the species, one can readily obtain either the diamond or
pyrochlore phase. Because both the diamond and pyrochlore phases have a
photonic band gap Hynninen _et al._ (2007), which makes them potential targets
for various applications in optics, Laves phases are strongly sought after.
However, these phases are notoriously difficult to self-assemble. At the high
packing fractions at which they become thermodynamically stable, the fluid
phase is dynamically sluggish Dasgupta _et al._ (2020) and may show
instabilities under sedimentation Milinković _et al._ (2011). Additionally,
the Laves phases are predicted to contain a significant concentration of
substitutional defects van der Meer _et al._ (2020), which can disrupt crystal
growth and hence hinder self-assembly Dasgupta and Dijkstra (2018).

Figure 23: Binary crystal phases of binary hard spheres. Schematic
representations of binary crystal structures known to be (nearly) stable in
binary hard-sphere mixtures.

”Laves phases have nevertheless been observed experimentally for binary
nanoparticle suspensions, which diffuse much faster than micron-scale colloids
Shevchenko _et al._ (2006); Evers _et al._ (2010); Cabane _et al._ (2016) (see
Sec. III.1), as well as sub-micron-sized spheres that interact with soft
repulsive potentials Nakagaki _et al._ (1983); Hasaka _et al._ (1984); Ma _et
al._ (1994); Gauthier _et al._ (2004). The observation of Laves phases in
softer spheres is consistent with recent simulation work showing that
interparticle softness pushes the glass transition in the binary fluid phase
to higher densities and hence enhances crystallization of Laves phases in
nearly hard spheres Dasgupta _et al._ (2020). However, Schaertl _et al._
(2018) have recently also demonstrated self-assembly of Laves phases in a
slightly off-stoichiometric mixture of (nearly) hard microgel particles, using
static light scattering. ”

Simulation and theoretical studies using a variety of methodologies have
explored other dense possible crystal structures of binary hard sphere
mixtures 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). These surveys have revealed a rich variety of dense packings
at different size ratios and compositions that are expected to be stable in
the limit of infinite pressures Hopkins _et al._ (2012). Interestingly, many
of them have yet to be observed in simulations or experiments at finite
pressure.

###  IX.2 Binary fluids: Structural crossover, demixing and dynamics

Figure 24: Structural crossover in binary hard-sphere fluids. Plots of
$\ln|rh_{ij}(r)|$ obtained by real-space experiment (top) and simulation
(bottom) for L-L, L-S, and S-S total correlation functions. The packing
fraction $\phi_{\mathrm{S}}$ (marked) is increasing from top (blue) to bottom
(red) for each case. Curves are shifted vertically for clarity. The black
lines in (a) and (b) are fits from which parameters controlling the crossover
are obtained Grodon _et al._ (2005). Simulations are for a binary hard-sphere
mixture with a size ratio $q=0.648$. Reproduced from Statt _et al._ (2016).

_Structure: crossover and demixing._ – DFT and computer simulations have
identified a _structural crossover_ line in the phase diagram of binary hard-
sphere fluids. This line marks a rather abrupt change in the wavelength that
dominates the asymptotic decay of the radial distribution function Grodon _et
al._ (2004, 2005). Experimental evidence for this effect was also found in
sedimented quasi-2d hard-sphere glasses Baumgartl _et al._ (2007) and
subsequently in 3d hard-sphere fluids where quantitative agreement with theory
and simulation was found (Fig. 24) Statt _et al._ (2016).

For smaller size ratios, it has long been debated whether a stable fluid-fluid
demixing transition exists or not. Such a spinodal instability in the fluid
mixture can be ascribed to the _depletion_ mechanism, which is known to drive
phase separation in colloid-polymer suspensions Asakura, S. and Oosawa, F.
(1954); Long _et al._ (1973); Vrij (1977); Lekkerkerker _et al._ (1992). The
depletion effect, first described by Asakura and Oosawa in 1954 Asakura, S.
and Oosawa, F. (1954); Oosawa (2021), induces an attractive interaction
between large colloids due to an unbalanced osmotic pressure of the small
(polymeric) spheres. Alternatively, this effect can be explained by the
increase in free volume available to the small spheres upon clustering of the
large spheres. The resulting entropy gain for the small spheres would then
drive phase separation of the colloids.

The effective depletion interaction between the larger spheres, due to the
smaller ones, has been calculated theoretically Attard (1989); Attard and
Patey (1990); Dickman _et al._ (1997); Mao _et al._ (1995); Götzelmann _et
al._ (1998, 1999); Roth _et al._ (2000) and extracted from computer simulation
Dickman _et al._ (1997); Biben _et al._ (1996); Götzelmann _et al._ (1999);
Ashton _et al._ (2011). The results exhibit a short-range attractive well
close to the surface of large hard spheres, followed by several oscillations
around zero at larger distances.

Figure 25: Radial distribution function of large colloids in binary hard
spheres at various concentrations of the smaller species for: (a) low
$\phi_{\mathrm{S}}\leq 0.25$; (b) high $\phi_{\mathrm{S}}>0.25$. Monte Carlo
simulation at small colloid volume fraction $\phi_{\mathrm{S}}$ are denoted by
solid lines, and compared to the experimental data (circles). Dashed lines
correspond to the relation $g(r)\approx\exp(-\beta u(r))$ where $u(r)$ is the
effective interaction between the large colloids consisting of the depletion
interaction induced by the smaller and a residual electrostatic contribution.
In addition to experimental resolution and polydispersity, MC simulations
account for residual colloid charge. Reproduced from Royall _et al._ (2007b).

Experimental investigations of the depletion interaction include tracking the
Brownian trajectory of a large sphere near a wall in a suspension of small
spheres with video microscopy Kaplan _et al._ (1994a) and using optical
tweezers to study the interaction between two large colloids Crocker _et al._
(1999). This latter work showed a discrepancy with subsequent theory Roth _et
al._ (2000). Later work which determined radial distribution functions from
which the interaction may be inferred showed better agreement with simulation
and theory (Fig. 25) Royall _et al._ (2007b); Roth _et al._ (2000).

Other experimental investigations probing fluid-fluid demixing have suggested
that it might be strongly coupled to freezing. For instance, Sanyal et al.,
using mixtures of polystyrene spheres with $q=0.2$, observed segregation in
regions rich in large spheres and rich in small spheres in the sediment at the
bottom of their samples Sanyal _et al._ (1992). However, when they suspended
the mixture in a density-matched solvent, neither sedimentation nor demixing
was seen. Van Duijneveldt _et al._ (1993) observed a phase instability in a
fairly narrow concentration range of small and large sterically stabilized
silica particles with $q=0.1667$, but sedimentation obscured whether the
transition corresponded to fluid-fluid demixing or to freezing. Experiments by
Kaplan et al. on mixtures of polystyrene particles with $0.069<q<0.294$
revealed the existence of a single homogeneous disordered phase, a coexistence
between two disordered phases, and a coexistence between one or two disordered
phases and a surface crystal on the sample wall Kaplan _et al._ (1994b).
However, bulk crystallization was not observed. A possible reason why no
surface crystallization was found in the experiments of van Duijneveldt et al.
is that silica spheres settle quickly compared to polystyrene particles. In
the experiments of Dinsmore et al. on mixtures of polystyrene particles with
$0.083<q<0.149$ Dinsmore _et al._ (1995), a phase separation into a disordered
fluid phase consisting primarily of small spheres and a crystalline solid of
large spheres permeated by a disordered fluid of small spheres was observed at
sufficiently high volume fractions. They additionally showed that the
crystallites on the surface of the sample cell had the same structure as the
bulk crystals, which they attributed to wetting of the bulk phase. Finally,
Imhof and Dhont observed a fluid-solid type of phase separation but no fluid-
fluid demixing in experiments on a binary mixture of silica spheres with
$q=0.1075$ Imhof and Dhont (1995a, b).

Conflicting results have also been reported from various theoretical
approaches. Integral equation theory approaches predict either that the
homogeneous fluid phase of a binary mixture of large and small hard spheres is
stable with respect to demixing Lebowitz and Rowlinson (1964); Mansoori _et
al._ (1971), that a spinodal instability occurs in the fluid phase Biben and
Hansen (1990, 1991a, 1991b); Rosenfeld (1994), or that the fluid-fluid
demixing transition is metastable with respect to a broad fluid-solid
transition by using a simple freezing criterion Caccamo and Pellicane (1997).
Free volume theory approaches have further predicted that the fluid-fluid
phase separation is metastable with respect to freezing Poon and Warren
(1994); Dinsmore _et al._ (1997). Furthermore, it has been shown that large
spheres could crystallize at the wall well below those corresponding to bulk
phase separation due to the presence of small spheres, in agreement with
experimental observations. A broad fluid-solid coexistence was also found in
DFT Xu and Barentin (1995), while theoretical approaches based on virial
coefficients predicted either a fluid-fluid phase separation in the limit of
highly asymmetric sizes Saija and Giaquinta (1996); Coussaert and Baus (1997);
Coussaert and Baus (1998b), a demixing transition metastable with respect to
either a broad or a narrow fluid-solid phase transition Coussaert and Baus
(1998a), a simple narrow fluid-solid phase transition in the limit
$q\rightarrow 0$ Vega (1998), or a complete absence of a fluid-fluid demixing
transition López de Haro _et al._ (2013). What this broad range of predictions
makes clear is that the phase behavior is too sensitive to the details of the
specific approximations used in the integral equation theories and virial
expansion approaches to reach a clear conclusion.

Dijkstra et al. followed a different approach by mapping the binary hard-
sphere mixture onto an effective one-component system by formally integrating
out the degrees of freedom of the small spheres in the partition function
Dijkstra _et al._ (1998); Dijkstra _et al._ (1999a). By using the two-body
(depletion potential) contribution to the effective Hamiltonian in
simulations, this effort revealed that the fluid-solid phase coexistence
region significantly broadens as $q$ becomes very small. In addition, it found
a stable isostructural solid-solid demixing transition for $q<0.05$ and fluid-
fluid demixing that remains metastable with respect to the fluid-solid
transition for $q<0.10$. These predictions were later validated by simulations
of the true binary hard-sphere mixture Dijkstra _et al._ (1999b, c). However,
only recently has the _critical point_ for the fluid-fluid demixing transition
been reported Kobayashi _et al._ (2021). Using a two-level simulation approach
based on a coarse-grained description with effective two- and three-body
interactions and the full (fine-grained) binary mixture, Kobayashi _et al._
(2021) managed that feat by matching the probability distribution for the
number of large particles to the 3d Ising universality class scaling form.
However, a non-trivial dependence on $q$ was also found.

In summary, it is now well-settled that a fluid-fluid demixing transition
exists in a binary mixture of hard spheres for sufficiently large size
asymmetries, but this transition is metastable with respect to a broad fluid-
solid transition.

_Dynamics in binary hard spheres._ – Dynamics in binary fluids, as probed with
multiple light scattering, has shown reasonable agreement with hydrodynamic
theory Kaplan _et al._ (1992). Later, Imhof and Dhont (1995b) measured the
diffusion in strongly asymmetric mixtures by photobleaching part of the sample
and examining the recovery of the fluoresence as the system rearranged. They
found good agreement with previous theory by Batchelor (1983) for the
diffusion of the small spheres in the dilute regime, but observed significant
deviations for the large spheres, which they attributed to the large size
disparity in their system. Significant experimental and simulation work has
also focused on the dense fluid regime Foffi _et al._ (2003); Marín-Aguilar
_et al._ (2020); Williams and Van Megen (2001); Götze and Voigtmann (2003),
thus revealing a sensitive dependence of the diffusivity behavior of both
species on the size ratio and composition of the chosen mixtures. This is of
particular importance when studying glassy systems (see Sec. XII), in which,
depending on $q$, one might find a single glass (where only the large
particles are dynamically arrested) or a double glass (where both species are
arrested) Voigtmann (2011); Lazaro-Lazaro _et al._ (2019); Laurati _et al._
(2018).

###  IX.3 Sedimented monolayers of binary hard spheres

Analogous to the monodisperse case, binary mixtures of hard spheres can be
confined to a quasi-2d setup by, e.g., allowing them to sediment into a
monolayer on a substrate. In this case, the system can be mapped onto a binary
hard-disk model, which typically has to be non-additive in order to account
for the centers of the small spheres being positioned below those of the large
spheres Thorneywork _et al._ (2014). Thorneywork et al. performed a detailed
characterization of the radial distribution function Thorneywork _et al._
(2014), structure factor Thorneywork _et al._ (2018b), and self-diffusion
coefficient Thorneywork _et al._ (2017a) of these systems, and found their
experimental results to be in good agreement with simulations and theory of
model mixtures. An intriguing consequence of the depletion interaction in
sedimented binary hard spheres is that the large particles are in fact
_repelled_ from the edge of a raised surface Dinsmore _et al._ (1996).

As in the bulk 3d case, the introduction of a second size of particles vastly
increases the phase diagram complexity. In 2d hard disks, fluid-fluid demixing
is generally expected only for positively non-additive hard disks Sillrén and
Hansen (2010), which would be difficult to achieve in colloidal experiments
with hard spheres. In the crystal regime, however, binary hard disks can
stabilize extremely diverse behavior. Considering only structures in the limit
of infinite pressure, Likos and Henley identified a wide range of distinct
periodic crystals as well as the possibility of a quasicrystal with 12-fold
symmetry Likos and Henley (1993). For extreme size ratios, the densest
possible crystal phases are expected to consist of a hexagonal packing of
large disks with an increasing number of small particles inserted into the
triangular voids Uche _et al._ (2004). For the specific case of hard spheres
sedimented onto a substrate, simulations interestingly predict the spontaneous
self-assembly of both 12-fold and 8-fold quasicrystals Fayen _et al._ (2023).

It should be noted that while monolayers of hard-sphere mixtures closely
approximate true 2d systems, simulations incorporating hydrodynamics on binary
mixtures confined between two plates suggest that even small amounts of
vertical freedom can significantly speed up dynamics of quasi-2d binary
mixtures Tian _et al._ (2022).

##  X Confinement

As detailed in Sec. VII, the entropy-driven equilibrium phase diagram of
_bulk_ hard-sphere systems is fairly well understood. The insertion of one or
more confining walls in the fluid, however, decreases the number of possible
configurations and generally alters the system structure and phase behavior
Bechinger (2002); Araújo _et al._ (2023). In this section, we specifically
consider the interfacial behavior of hard-sphere fluids between two hard
walls, under quasi-1d and quasi-2d conditions as well as within cylindrical,
spherical, and flexible immurement. Note that the case of a single planar
wall, which results in an interface albeit not in confinement proper is
presented with in Sec. VIII.1.

###  X.1 Quasi-2d confinement

The effect of confinement between two planar hard walls was first investigated
in Pieranski _et al._ (1983)using colloidal polystyrene spheres. They observed
a sequence of layered solid structures with triangular ($\triangle$) and
square ($\square$) symmetry, $1\triangle\rightarrow 2\square\rightarrow
2\triangle\rightarrow 3\square\rightarrow 3\triangle\ldots\rightarrow
n\square$ upon increasing the plate separation, where $n$ denotes the number
of layers. Using MC simulations and cell theory (Sec. VI.2), Schmidt and Löwen
later mapped out the phase diagram for plate separations ranging from one to
two particle diameters, and found additional buckled and rhombic phases
Schmidt and Löwen (1996, 1997). The sequence of high-density structures was
determined more accurately in subsequent experiments Neser _et al._ (1997);
Fontecha _et al._ (2005), further identifying prism phases with both square
and triangular symmetry.

Figure 26: a) Equilibrium phase diagram of hard spheres confined between two
parallel hard walls with plate separation $H$ versus packing fraction
representation. The white, shaded (yellow) and dotted regions indicate the
stable one-phase region, the two-phase coexistence region, and the forbidden
region, respectively. (b,c) Optical microscopy image of charged spheres in an
aqueous system confined in a wedge geometry. Shown are both a 2$\Delta$ phase
(b) and a $2\square$ phase (c). Reproduced with permission from (a) Fortini
and Dijkstra (2006) and (b,c) Fontecha _et al._ (2005). Copyright IOP
Publishing. All rights reserved.

Using extensive free-energy calculations in MC simulations, the full phase
diagram for plate separations from one to five hard-sphere diameters was
mapped out as a function of packing fraction Fortini and Dijkstra (2006), see
Fig. 26. These results identify a first-order fluid–solid transition,
corresponding to either capillary freezing or capillary melting depending on
the plate separation, with the coexisting solid phases consisting of
crystalline layers with either triangular or square symmetry. At high
densities, prism, buckled, and rhombic phases are found to be
thermodynamically stable in agreement with experiments Neser _et al._ (1997);
Fontecha _et al._ (2005); Oğuz _et al._ (2012).

Curk _et al._ (2012) used simulations to investigate hard spheres in soft
quasi-2d confinement with a parabolic potential along one dimension. Like in
hard confinement, a sequence of confined hexagonal and square-symmetric
packings was found, but none of the intervening ordered phases were observed,
the system undergoing phase separation instead.

###  X.2 Cylindrical confinement

Hard spheres perfectly confined to a 1d line form a _Tonks gas_. This model,
which has a rich theoretical history Lieb and Mattis (1966), captures
reasonably well the behavior of various condensed-matter systems, such as
linear chains of mercury in Hg3-δAsF6 Spal _et al._ (1980); Chaikin, P. M. and
Lubensky, T. C. (1995)101010Note that the analysis of Refs. Spal _et al._ ,
1980; Chaikin, P. M. and Lubensky, T. C., 1995 did not take advantage of the
exact 1d solution. This analysis is left as an exercise to the reader., and
chains of cobalt atoms confined by platinum steps Gambardella _et al._ (2002).
(Its quantum mechanical relative, the Tonks-Girardeau model, has also been
experimentally realized Bloch _et al._ (2008).)

The Tonks gas can be solved by a transfer matrix (TM) scheme. A rich set of
structural and thermodynamic observables can thus be obtained at minimal
computational cost. With appropriate discretization schemes, the TM scheme can
also be extended to quasi-1d systems Barker (1962). Both hard disks between
hard lines Kofke and Post (1993); Gurin and Varga (2013); Godfrey and Moore
(2014); Robinson _et al._ (2016); Hu _et al._ (2018); Hicks _et al._ (2018);
Zhang _et al._ (2020) and hard spheres within a hard cylinder Kamenetskiy _et
al._ (2004); Hu _et al._ (2018) (as well as more exotic hard shapes Kantor and
Kardar (2009); Gurin and Varga (2015); Gurin _et al._ (2016); Gurin _et al._
(2017a, b)) have thus been considered using this approach.

Figure 27: Close packing fraction of hard spheres confined in cylinders of
diameter $D$ along with some of the structures. These packings (a) are
phyllotactic-like for $D/\sigma\lesssim 2.8$, and (b) exhibit a core–shell
structure beyond that point. Reproduced from Fu _et al._ (2016) with
permission from the Royal Society of Chemistry.

The TM size, however, jumps markedly as the number of interacting particles
increases, which limits the maximal cylinder width that can be computationally
resolved. At present, it is nevertheless possible to study systems with up to
next-nearest-neighbor interactions. Spheres in a cylinder with a diameter of
$2\sigma$ (and for disks a line spacing of up to $5\sigma/2$) can thus also be
solved this way Hu _et al._ (2018). Despite the limitations of this size
regime, it suffices for zig-zag and helical order to emerge at high density.

The analysis of these systems provides key insight into quasi-1d ordering. For
instance, although 1d systems with finite-range pairwise interactions have
long been known to be unable to exhibit a phase transition van Hove (1950)
(see also (Ruelle, 1999, Thm. 5.6.7)), various numerical simulations of
cylindrically confined hard spheres have identified somewhat abrupt structural
changes as density increases Gordillo _et al._ (2006); Varga and Ball[́o ;
Koga and Tanaka (2006); Duran-Olivencia and Gordillo (2009). Similarly, it has
been proposed that confined disks could exhibit a Kosterlitz-Thouless–type
phase transition Huerta _et al._ (2020); Hu and Charbonneau (2020). The TM
scheme neatly resolves these apparent paradoxes. First, their consideration of
infinite and fully equilibrated systems removes all hints of a thermodynamic
singularity (as given by the largest TM eigenvalue). Second, the marked
structural changes find a microscopic explanation. Because the correlation
length that describes the spatial decay of structural order is given by the
ratio of the two largest (magnitude-wise) TM eigenvalues,
$\xi\sim\ln(\lambda_{1}/\lambda_{2})$, crossing of sub-dominant TM eigenvalues
can give rise to marked structural changes without any thermodynamic
singularity Hu _et al._ (2018).

In order to observe richer order types, larger cylinder diameters need to be
considered. The diameter dependence indeed does not result in continuous
structural changes, as is most obvious at close packing. (The effect is also
notable at finite pressure Mon and Percus (2000).) The associated
morphological richness was first studied by Pickett _et al._ (2000), who noted
that chiral order spontaneously develops for certain diameter ratios. Mughal
et al. Mughal _et al._ (2011, 2012) later noted that for cylinders of diameter
up to $2\geq D/\sigma=1+1/\sin(\pi/5)\approx 2.70$ close-packing is described
using a phyllotactic construction because all spheres then coat the cylinder
wall. Beyond this diameter regime, not all spheres touch the cylinder wall,
which eventually results in a separation between core and shell particles Fu
_et al._ (2016). Exotic arrangements, complex helices and limit periodic
structures then follow. More or less systematic numerical exploration of these
structures ends around $D/\sigma=4$, but it is conceivable that larger
diameters might accommodate even more exotic structures.

Considerable interest has also been paid to the transport dynamics of these
systems. In particular, for the hard-sphere case, Mon and Percus identified a
crossover from single-file to Fickian diffusion for sufficiently wide
cylinders, i.e., $D\geq 2\sigma$ Mon and Percus (2002, 2003). A transition
state theory description of the hopping mechanism that enables sphere passing
was later proposed Wanasundara _et al._ (2012), and the impact of microscopic
dynamics Sané _et al._ (2010); Flomenbom (2010) as well as size dispersity
Wanasundara _et al._ (2012) on these has been considered.

From the experimental standpoint, hard spheres in quasi-1d confinement have
notably been used to rationalize the packing behavior of C60 in nanotubes
Mickelson _et al._ (2003); Khlobystov _et al._ (2004, 2005); Troche _et al._
(2005), of confined nanoparticles Tymczenko _et al._ (2008); Sanwaria _et al._
(2014), and of vacuolated cells in an embryonic structure Norman _et al._
(2018). A few colloid-based realizations have also been achieved Lohr _et al._
(2010); Jiang _et al._ (2013); Liu _et al._ (2015a); Fu _et al._ (2017), as
have granular-scale ones Bogomolov _et al._ (1990). In all cases, a good
correspondence between theory and experiments is obtained, although the
assembly pathway may need to be taken into account to rationalize the observed
packings Mughal _et al._ (2012); Fu _et al._ (2017). Nevertheless, relatively
few of the above numerical and theoretical predictions have been tested
directly in experiments.

###  X.3 Quasi-1d confinement

A lower-dimensional generalization of cylindrical confinement considers hard
disks between parallel lines a distance $D$ apart. These models have largely
been studied for the insight they offer into the physics of disordered
systems. Given that the geometry of these systems is a lot simpler than that
of spheres in cylindrical confinement, only fairly unremarkable structures
form the densest close packing as $D$ increases. Subtleties in these packings
nevertheless lead to the emergence of nontrivial local features Ashwin and
Bowles (2009); Zhang _et al._ (2020). Such irregularities have allowed these
models to find particular use in clarifying the physics of low-dimensional
disordered systems.

Given the robustness of the jamming phenomenology down to 2d systems, such as
the algebraic scaling of weak forces and small gaps (see Sec. XII.8), an
interesting target has been to consider the structural criticality of quasi-1d
systems. For $D<1+\sqrt{3}\sigma/2$, jammed states are then isostatic and have
a nonzero complexity Ashwin and Bowles (2009); Bowles and Ashwin (2011);
Ashwin _et al._ (2013); Godfrey and Moore (2018), but do not exhibit a
critical structure. For $1+\sqrt{3}/2<D<2$, however, it is possible to
consider states that are critical. Then, depending on the system details,
different critical exponents have been reported Ikeda (2020); Zhang _et al._
(2020). This sensitivity of the structural criticality likely follows from
these systems being below the lower critical dimension for jamming.

These models have also been used to study dense fluids Godfrey and Moore
(2014); Robinson _et al._ (2016); Hicks _et al._ (2018), so as to assess how
some of the mean-field theory predictions then fare. Unsurprisingly for such
low-dimensional models, a qualitatively different physics is observed. The
findings should nevertheless caution against possible confounding physical
factors in the study of higher-dimensional fluids, such as the importance of
local structure and the possibility of crossovers at fairly large system
sizes.

Given the rich physics of such simple models, a seemingly open area of
research is their experimental study. Microfluidic channels with colloids Mark
_et al._ (2010), for instance, might be a promising approach to consider.

###  X.4 Spherical confinement

Arguably the most natural way of confining colloidal spheres in three
dimensions is inside a larger sphere. As a natural toy model system for small
many-body systems, even extremely small systems of only a few hard spheres in
spherical confinement have attracted significant theoretical attention
exploring, e.g., the effect of the thermodynamic ensemble on the observed
structure González _et al._ (1998, 1997), and the thermodynamic properties of
the confined fluid Urrutia (2011); Urrutia and Pastorino (2014); Urrutia and
Castelletti (2012).

As one would expect, hard-sphere fluids confined to larger spherical cavities
are known to show structuring Macpherson _et al._ (1987); Chui (1991); Zhou
and Stell (1989); Calleja _et al._ (1991) near the cavity walls. The dynamics
inside the cavity then depends strongly on whether the cavity walls are rough
or smooth Németh and Löwen (1999), with rough walls more strongly inducing
dynamical arrest. In particular, when a glassy fluid of hard spheres is
confined in a spherical cavity with rough walls consisting of pinned
particles, the dynamics approaches that of a bulk hard-sphere glass for large
cavities, providing a route to probe dynamical correlations in these systems
Németh and Löwen (1999); Zhang and Cheng (2016).

Experimentally, self-assembly of hard-sphere colloids in spherical confinement
can be achieved by confining the colloids inside emulsion droplets Yi _et al._
(2002); Manoharan _et al._ (2003); Manoharan (2006). When the droplets shrink
(due to either evaporation or Ostwald ripening Schmitt _et al._ (2016)), the
slowly increasing density of the colloids can eventually lead to self-
assembly. The resulting structure depends sensitively on the number of
colloidal particles within the cluster. Clusters containing only a handful of
colloids form small polyhedral clusters whose geometry can often be understood
as the densest configuration possible in the circumstances. However, in these
small clusters, capillary forces during the final stages of solvent
evaporation also play a significant role, and in many cases drive the cluster
to minimize the second moment of its mass distribution Manoharan _et al._
(2003); Cho _et al._ (2005a); Yi _et al._ (2004); Lauga and Brenner (2004).
For hard spheres, this set of minimal clusters has been studied extensively by
Sloane _et al._ (1995).

For larger clusters of monodisperse spheres in spherical confinement, the
natural tendency of hard spheres to crystallize into an fcc crystal competes
with the curvature of the surface. Although sufficiently large clusters
($N\gtrsim 10^{5}$) form a simple fcc structure, intermediate cluster sizes
($100\lesssim N\lesssim 10^{5}$) spontaneously form clusters with icosahedral
symmetry de Nijs _et al._ (2015). The core of these icosahedral clusters
consists of tetrahedral-shaped domains of distorted fcc crystal, with one
hexagonal plane from each domain forming the faces of an icosahedron, often
called Mackay clusters Mackay (1962). For sufficiently large clusters, these
domains can be capped by additional surface layers, resulting in a family of
“anti-Mackay” sphere packings with varying surface reconstructions de Nijs _et
al._ (2015); Wang _et al._ (2018b) (see Fig. 28). Free-energy calculations
based on computer simulations have shown that these clusters are indeed
thermodynamically stable configurations for pure hard spheres in spherical
confinement de Nijs _et al._ (2015), with clusters containing certain “magic
numbers” of spheres corresponding to defect-free clusters being particularly
stable Wang _et al._ (2018b); Wang _et al._ (2019).

Figure 28: Icosahedral clusters of spheres self-assembled in spherical
confinement. The top row shows both an electron microscopy image of a cluster
of nanoparticles and a reconstruction of its surface layer. The bottom row
shows (left) a cluster obtained in experiments on micron-sized colloids and
(right) one obtained in simulations of hard spheres. Figure adapted with
permission from Ref. de Nijs _et al._ (2015).

Naturally, this already complex behavior can be further tuned by considering,
e.g., binary mixtures of spheres in spherical confinement. For small clusters,
complex anisotropic colloidal supraparticles can result Cho _et al._ (2005b).
For larger clusters, Wang et al. demonstrated the formation of clusters
consisting of the thermodynamically stable MgZn2 Laves phase when there is an
excess of small spheres Wang _et al._ (2021b), but the clusters turn into an
icosahedral cluster consisting of tetrahedral domains of the less stable MgCu2
Laves phase when the fluid composition is iso-stoichiometric with Laves phase
Wang _et al._ (2021a).

As the self-assembly behavior of colloidal spheres under spherical confinement
is robust on both the nanometer and micrometer scale de Nijs _et al._ (2015),
and for a variety of materials Yi _et al._ (2004), it provides a versatile
route for tuning material structure across various length scales.
Supraparticles created via this technique have been proposed as building
blocks for additional self-assembly steps, e.g., for creating materials with
structural order on two different length scales Bai _et al._ (2007).
Additionally, the icosahedral nature of intermediate-sized clusters of
colloidal spheres may lead to intriguing optical effects, including
(iridescent) structural color patterns Wang _et al._ (2020a) that can be used
to track the self-assembly dynamics inside the clusters in real time.

###  X.5 Flexible confinement

Thus far in this section, we have seen that a large variety of packing
geometries can be obtained by confining spheres in a container of fixed shape.
A different question concerns the most efficient packing of spheres in a
natural bin, which is the smallest convex hull (volume-wise) that can enclose
a certain number of spheres. Counterintuitively, the result is not always a
compact cluster of spheres. For a packing of spheres with up to $n=55$, along
with $n=57,58,63,64$, the linear conformation in which the spheres lie on a
straight line, also called a sausage, is denser than a cluster-like or plate-
like configuration Gandini and Wills (1992). The optimal packing for $n=56$ is
not fcc like, and the exact configuration remains unknown. (By contrast, in
4d, this sudden transition from a sausage packing to a cluster shape is
conjectured to happen at $n=375769$ spheres, and is therefore referred as the
“sausage catastrophe” Henk and Wills (2021).) Hard spheres in a flexible
container can be used to model colloids in a fluctuating vesicle, which was
studied theoretically and in simulations by Maibaum _et al._ (2001). Recently,
such a system of non-close-packed colloids in a flaccid lipid vesicle has been
experimentally realized Marín-Aguilar _et al._ (2023). Through a combined
experimental and simulation study, the authors obtained a state diagram that
includes linear, planar, and cluster conformations of spheres, as well as
bistable states, which alternate between cluster-plate and plate-linear
conformations due to membrane fluctuations. Additionally, Marín-Aguilar _et
al._ (2023) have identified truncated polyhedral packings of $56\leq N\leq 70$
spheres (excluding $N=57$ and 63) that pack more efficiently than linear
arrangements.

In experiment, flexible confinement may be induced by the use of optical
tweezers. To the best of our knowledge, in the context of hard spheres, using
lasers to induce structuring may be traced to the work of colloid pioneers
Bruce Ackerson and Noel Clark and coworkers in their use of lasers to induce
freezing in quasi–2d systems Chowdhury _et al._ (1985). A more recent approach
was to confine a system inside a ring of particles held in optical traps. Such
traps have well–defined potentials, which means that the osmotic pressure and
therefore the EoS could be measured. At low density, the results were found to
agree with the bulk EoS. At higher density, by contrast, the fluid takes on a
layered structure, thus reflecting the confinement. Interestingly, at still
higher density, a bi-stable state between the layered fluid and a hexagonal
structure was observed, reminiscent of the bulk hexatic phase Williams _et
al._ (2013).

##  XI Out-of-equilibrium phenomena

In this section, we discuss how hard spheres have shed some light on material
phenomena away from thermal equilibrium. When studying out-of-equilibrium
suspensions it is useful to distinguish linear response around thermal
equilibrium, which can be treated perturbatively within the framework of
statistical mechanics, from phenomena that occur due to strong driving beyond
the linear regime. Non-equilibrium effects often appear as fluxes of
macroscopic observables and transport phenomena. Because particle motion
remains governed by thermal fluctuations, the relative magnitude of this
transport compared to diffusion, as captured by the various dimensionless
Péclet numbers (Table 4), e.g. for gravitation, shear flow Russel _et al._
(1989), or other fields such as optical tweezers Williams _et al._ (2016), is
a crucial measure how strongly a system is driven.

| general force | gravitation | drift speed | shear flow | glassy flow
---|---|---|---|---|---
drive | $F$ | $g$ | $v$ | $\dot{\gamma}$ | $\dot{\gamma}$
$\tau_{x}$ | $\frac{\sigma}{2\mu F}$ | $\frac{\sigma}{2\mu m_{b}g}$ | $\frac{\sigma}{2v}$ | $\dot{\gamma}^{-1}$ | $\dot{\gamma}^{-1}$
$\operatorname{P\kern-0.80002pte}_{x}$ | $\frac{F\sigma}{2k_{\text{B}}T}$ | $\frac{\pi\sigma^{4}\delta\rho g}{12k_{\text{B}}T}$ | $\frac{2\sigma v}{D}$ | $\dot{\gamma}\tau_{\text{B}}$ | $\dot{\gamma}\tau_{\alpha}$
Table 4: Péclet numbers
$\operatorname{P\kern-0.80002pte}_{x}=\tau_{\text{B}}/\tau_{x}$ relate the
diffusive time $\tau_{\text{B}}$, as defined in Eq. (1) to the time scale
$\tau_{x}$ of a directed transport process using the particle radius
$\sigma/2$ as reference length. $\mu=D_{0}/k_{\text{B}}T$ is the particle
mobility. Sedimentation (second column) is due to the gravitational force
$m^{\ast}g$ with the gravitational $g$ and buoyant mass
$m_{b}=(\pi\sigma^{3}/6)\delta\rho$, where $\delta\rho$ is the density
difference between particle and solvent. In this case the Péclet number can
also be written $\operatorname{P\kern-0.80002pte}_{g}=\sigma/(2\xi)$ with
gravitational length $\xi$. In the case of shear flow, the timescale of
directed transport is $\dot{\gamma}^{-1}$. For normal fluids, the Brownian
time $\tau_{\text{B}}$ is appropriate, but for glassy systems the Weissenberg
number $\dot{\gamma}\tau_{\alpha}$ is often used (see Sec. XII.10).

Thus far in this review, the solvent has largely been considered as a
structureless quiescent medium. Out of equilibrium, however, this
simplification is often no longer appropriate. Solvent dynamics has to be
taken into account. In this section, we therefore shift our focus away from
hard spheres as a model system and consider instead the colloidal material as
a complex fluid. In that context, the full dynamics of a particle-laden
Newtonian solvent is described by the Navier-Stokes equation together with a
suitable (typically no-slip) boundary condition on the particle surfaces. By
eliminating the solvent, forces between suspended particles effectively couple
through the mobility tensor which can, in principle, be derived from the Oseen
tensor Dhont (1996). In bulk systems, hydrodynamic coupling is long-ranged due
to momentum conservation of the solvent. Taking these hydrodynamic
interactions into account in simulations is feasible but computationally
costly (see Sec. V.4 for available methods). A recurring theme is the
importance of colloidal forces over hydrodynamic coupling. If the latter can
be neglected, one can resort to computationally cheaper BD simulations (Sec.
V.3).

###  XI.1 Non-equilibrium sedimentation in hard-sphere colloids

Figure 29: Non-equilibrium sedimentation of hard spheres on the particle
scale. (a) Illustration of a system of sterically stabilized PMMA particles
under the influence of gravity $g$ and vertically confined between two walls
separated a distance $L$. (b-d) Time series of confocal micrographs taken in
the (vertical) $xz$ plane at times (b) $t=3$, (c) $26$, and (d)
$200\tau_{\text{B}}$. The scale bars denote 20 $\mu$m; the horizontal lines
indicate the position of the walls. For (b-d) the sedimentation Péclet number
is $\operatorname{P\kern-0.80002pte}_{g}=0.625$. (e) Time evolution of the
sedimentation profile $\phi(z,t)$ for a system with
$\operatorname{P\kern-0.80002pte}_{g}=1.11$. Solid lines are experimental data
from particle-resolved studies, dashed lines are from dynamical DFT. No fit
parameters are used. Reprinted with permission from Royall _et al._ , Phys.
Rev. Lett. 98, 188304 (2007) Royall _et al._ (2007a). Copyright (2007) by the
American Physical Society. Figure 30: Rayleigh-Taylor–like instability in
hard-sphere colloids. (a-c) Time series of images taken with a confocal
microscope, for $1.43$, $5.5$ and $11.2$ $\tau_{\text{B}}$. $\lambda$ denotes
the characteristic growing wavelength of the Rayleigh-Taylor instability. (d)
Image in $xy$ plane at the height of the dashed line in (c). Reproduced from
Wysocki _et al._ (2009) with permission from the Royal Society of Chemistry.

Let us first consider what is arguably the simplest non-equilibrium situation,
namely the gravitational settling of a suspension. As discussed in Sec. VII.3,
waiting for the suspension to relax yields a density profile from which the
equation of state may be inferred. While such equilibrium profiles are
straightforward to calculate, the approach to equilibrium, from an out-of-
equilibrium starting point is a challenging problem to address.

As alluded to above, a useful way to distinguish different regimes is through
the _gravitational_ (or sedimentation) Péclet number

$\operatorname{P\kern-0.80002pte}_{g}=\frac{\tau_{\text{B}}}{\tau_{g}}=\frac{\sigma/2}{\xi_{g}}=\frac{\pi\sigma^{4}\delta\rho
g}{12k_{\mathrm{B}}T}$ (34)

defined as the ratio of the Brownian diffusive time $\tau_{\text{B}}$ [Eq.
(1)] to the time $\tau_{g}$ for a sphere to sediment its own radius, where
$\delta\rho$ is the mass density difference between colloidal particles and
the solvent. Alternatively it can be expressed through the gravitational
length $\xi_{g}=k_{\text{B}}T/m^{\ast}g$ with $m^{\ast}$ the buoyant mass.
Note that without loss of generality, sedimentation is assumed to occur in the
$z$ direction and has no dependence on the $x,y$ coordinates. Equation 34
describes the sedimentation of a single colloid, and so holds in the dilute
limit.

For more concentrated systems, when the density varies on lengthscales
sufficient for local packing effects to be neglected, and for which
gravitational settling is slow, i.e., $\operatorname{P\kern-0.80002pte}_{g}\ll
1$, the time evolution of the sedimentation profile $\phi(z,t)$ largely
follows a _batch settling_ process Russel _et al._ (1989). Under these
conditions, the volume fraction at any height and at any time may be captured
by a set of relatively simple coupled equations. (For more strongly driven
systems, in the granular regime, strong swirls due to hydrodynamic coupling
can be observed Segrè _et al._ (2001, 1997).)

For smaller systems, for which packing effects can be more noticeable (see
Sec. X), classical DFT (see Sec. VI.3) is needed. Using dynamical DFT, it is
indeed possible to propagate these density profiles forward in time, as shown
in Fig. 29(d,e). The time-dependent density profiles predicted by dynamical
DFT agree remarkably well with particle-resolved results from BD simulations.
With the inclusion of a simple treatment for the volume fraction-dependent
slowdown in dynamics due to hydrodynamic interactions (treated with the
Hayakawa-Ichiki method), the time evolution of the density profiles of
sedimenting hard spheres in experiment could be accurately described (see Fig.
29(a-e)) Royall _et al._ (2007a); Schmidt _et al._ (2008).

This description of sedimentation presumes translational invariance in the
$(x,y)$ plane, which holds (empirically) for a starting configuration that is
approximately homogeneous Royall _et al._ (2007a); Schmidt _et al._ (2008). An
inhomogeneous configuration such as that shown in Fig. 30(a) presents quite a
different–and at first sight surprising–proposition. Density fluctuations with
a characteristic length scale (Fig. 30(b)) then develop, leading to a behavior
reminiscent of the catastrophic Rayleigh-Taylor instability when two
immiscible liquids are prepared with the denser liquid above. What is
particularly surprising here is that the hard-sphere fluid, which of course is
a single phase, then behaves like a phase-separated system. These complex
time-dependent patterns can further be accurately captured by multi-particle
collision dynamics simulations Wysocki _et al._ (2009, 2010). Therefore, even
phenomena quite far-from-equilibrium can be accurately captured by theory and
simulation.

Adding a second colloidal species, whose gravitational Péclet number can be
tuned independently, provides a further means of controlling the structural
properties of the swirls. Interestingly, these swirls depend mainly on the
relative magnitudes of the Péclet numbers of the two species, and much less on
the composition of the mixture which would be the case closer to equilibrium
(see Sec. IX) Milinković _et al._ (2011).

###  XI.2 Rheology, flow, and shear-induced order

The most prevalent way to probe mechanical properties is to subject a sample
to external forces and measure its deformation (rate) Wagner (2022); Larson
(1998); Chen _et al._ (2010b). In addition to the packing fraction $\phi$, the
strain $\gamma$ (strain rate $\dot{\gamma}$, not to be confused with the
interfacial free energy in Sec. VIII.2) with conjugate stress $\tau$ are then
needed to characterize a hard-sphere system. The absence of cohesive forces
between hard spheres implies that their solids deform easily. The elastic
moduli, which encode how a solid linearly deforms in response to an applied
force, nevertheless diverge as $\propto
k_{\text{B}}T/[\sigma^{3}(\phi_{\text{cp}}/\phi-1)^{2}]$ upon approaching
close packing Stillinger and Salsburg (1967); Farago and Kantor (2000).

Similarly, a hard-sphere fluid starts flowing in response to shear forces,
thus entering a dissipative non-equilibrium state. In certain geometries
(e.g., Couette or cone-plate shear cells, cf. Fig. 31(a)), a uniform strain
rate $\dot{\gamma}$ can be achieved, thus defining viscosity
$\eta=\tau/\dot{\gamma}$. How strongly a suspension is sheared is described
through the dimensionless shear Péclet number
$\operatorname{P\kern-0.80002pte}_{\dot{\gamma}}=\dot{\gamma}\tau_{\text{B}}$,
again using the Brownian time $\tau_{\text{B}}$ [Eq. (1)] (Table 4),
characterizing the importance of advection over diffusion. Figure 31(b) shows
typical flow curves $\tau(\dot{\gamma};\phi)$ for dense hard spheres. For the
fluid ($\phi=0.52$) one observes shear thinning with a stress that is smaller
than the initial linear increase. The denser ($\phi\geqslant 0.59$) no-slip
samples behave as Herschel–Bulkley fluids following the empirical flow curve
$\tau-\tau_{0}\propto\dot{\gamma}^{n}$ with some exponent $n$. They approach a
finite yield stress $\tau_{0}$ in the limit $\dot{\gamma}\to 0$, thus
indicating that the sample behaves as a (disordered) solid. Rheology therefore
provides a mechanical route to probe the glass transition, to which we return
in Sec. XII.10.

Figure 31:  (a) Sketch of a cone-plate rheometer together with a confocal
microscope. (b) Measured shear stress $\tau$ as a function of applied strain
rate $\dot{\gamma}_{\text{a}}$ for several packing fractions $\phi$ below and
above the glass transition. Open symbols indicate uncoated smooth plates while
full symbols show data for coated plates with no-slip boundary conditions. (c)
The excess stress $\tau-\tau_{0}$ is linear with slope $\eta_{\text{eff}}$
defining an effective viscosity. Reprinted with permission from Ballesta _et
al._ , Phys. Rev. Lett. 101, 258301 (2008) Ballesta _et al._ (2008). Copyright
(2008) by the American Physical Society.

Traditional rheological studies probe macroscopic volumes through small-
amplitude oscillatory shear, so as to only weakly perturb the material from
thermal equilibrium. In this linear response regime with
$\operatorname{P\kern-0.80002pte}_{\dot{\gamma}}\ll 1$, the material settles
in a periodic state with stress

$\tau(t)=\int_{-\infty}^{t}ds\;G(t-s)\dot{\gamma}(s)$ (35)

leading to the storage modulus $G^{\prime}(\omega)$ and loss modulus
$G^{\prime\prime}(\omega)$, the real and imaginary part of the complex shear
modulus $G(\omega)$, as a function of external frequency $\omega$.

For a dilute hard-sphere suspension, Einstein showed that (to linear order in
$\phi$) $\eta/\eta_{\text{sol}}=1+\tfrac{5}{2}\phi$ Einstein (1906); Mewis and
Wagner (2011). In other words, viscosity increases with respect to the solvent
viscosity $\eta_{\text{sol}}$ due to volume excluded by the suspended
particles. The quadratic order correction in $\phi$ arises from the
hydrodynamic coupling between spheres. The exact value of the coefficient,
however, depends on the analytical approach (Batchelor (1977), for instance,
found $6.2\phi^{2}$).

At higher strain rates beyond linear response (for
$\operatorname{P\kern-0.80002pte}_{\dot{\gamma}}\approx 1$), shear thinning
sets in. Viscosity then decreases upon approaching a plateau $\eta_{\infty}$
van der Werff and de Kruif (1989), as has been observed in BD simulations
without hydrodynamics Strating (1999). The phenomenon can therefore be linked
to a change in the local arrangement of particles Xu _et al._ (2013) through
the excess shear stress

$\Delta\tau=-\frac{1}{2}k_{\text{B}}T\rho^{2}\int
d^{3}\bm{r}\;\frac{xy}{r}g(\bm{r})\delta(r-\sigma).$ (36)

Importantly, the pair distribution function $g(\bm{r})$ is then no longer
isotropic Lin _et al._ (2013); Lin _et al._ (2014a), because an external force
or flow defines a preferred direction. In particular, its contact value,
$g^{+}(\sigma,\theta)$, using polar coordinates within the $xy$-plane through
the particle center, varies with orientation and derived quantities like
pressure become anisotropic with the off-diagonal component determining the
shear stress, cf. Eq. (36). In linear response, Brady (1993) determined the
Brownian stress as the equilibrium $g(\sigma^{+})$ at contact divided by the
short‐time self‐diffusivity.

The deformation at intermediate Péclet numbers is shown in Fig. 32(a-c). The
relationship between macroscopic material properties and microscopic structure
is at the heart of theoretical approaches. Much effort has been devoted to
predict the deformation of the pair distribution in response to external
forces and flows Squires and Brady (2005), allowing to calculate flow curves
$\tau(\dot{\gamma})$ from first principles. Coupling of confocal microscopy
(of sufficient frame rate) and rheology provides a powerful means to access
the local structure in hard-sphere fluids far from equilibrium Besseling _et
al._ (2009); Lin _et al._ (2014b).

Figure 32: The microstructure $g(\bm{r})$ deforms in shear flow at (a)
$\operatorname{P\kern-0.80002pte}_{\dot{\gamma}}=0.1$, (b)
$\operatorname{P\kern-0.80002pte}_{\dot{\gamma}}=1$, and (c)
$\operatorname{P\kern-0.80002pte}_{\dot{\gamma}}=10$. Reprinted with
permission from Foss and Brady, J. Rheol. 44, 629 (2000) Foss and Brady
(2000a). Copyright 2000, The Society of Rheology. (d) Sketch of the different
flow regimes of hard sphere suspensions. Lines indicate crossovers, with shear
thinning being predominately due to colloidal forces and shear thickening
having a hydrodynamic origin. Shear thickening has also been observed in
dilute suspensions at high strain rates Bergenholtz _et al._ (2002). At high
packing fractions and high strain rates shear-induced ordering is possible.

At even higher strain rates (or shear stresses) and moderate to high packing
fractions ($\phi>0.3$), some experiments report an abrupt (sometimes
discontinuous) viscosity increase as the stress is increased, i.e., shear
thickening d’Haene _et al._ (1993); Bender and Wagner (1996); Meeker _et al._
(1997); Cheng _et al._ (2002). A comprehensive understanding of this somewhat
counter-intuitive behavior (for instance, it is absent in atomistic and
molecular liquids) has been an open challenge Brown and Jaeger (2014). For
Brownian hard spheres it is now broadly accepted to be a hydrodynamic
phenomenon caused by the (reversible) formation of “hydroclusters” held
together through lubrication forces Wagner and Brady (2009). Confocal
microscopy evidence Cheng _et al._ (2011) has also been used to investigate
the hydrodynamic and contact force contributions to shear thickening Lin _et
al._ (2015). Shear thickening extends to dilute suspensions ($\phi<0.1$),
although only at very high strain rates Bergenholtz _et al._ (2002). Figure
32(d) summarizes these different uniform flow regimes. Shear thickening has
also been proposed as the mechanism that links the rheology of Brownian hard
spheres ($\sigma\lesssim 1\;\mu$m) to the physics of non-Brownian (granular)
hard spheres ($\sigma\gtrsim 50\;\mu$m) Guy _et al._ (2015).

While the shear thickening behavior is associated with an abrupt increase in
the hydrodynamic stresses between the colloidal spheres, which locks them
together giving rise to the large hydroclusters resistant to flow,
hydrodynamic interactions between smooth surfaces predict a continuous shear-
thickening transition Ball and Melrose (1995); Foss and Brady (2000b); Bender
and Wagner (1996); Melrose and Ball (2004). The discontinuous appearance of
shear thickening instead is found to depend on the small scale surface
asperities of the colloidal particles. These asperities can break the
lubrication layer, and can give rise to a frictional contribution. The
relative motion of the particles is then akin to a stick-slip scenario
triggered by the breakage of the lubrication layer Mari _et al._ (2015);
Hermes _et al._ (2016); Morris (2018); Kawasaki and Berthier (2018); Jamali
and Brady (2019). Simulations have explicitly shown a transition from
continuous to discontinuous shear thickening by progressively increasing the
surface roughness of the particles Wang _et al._ (2020b). The importance of
surface characteristics has also been confirmed in a large variety of
experimental works Fernandez _et al._ (2013); Hsu _et al._ (2018); Schroyen
_et al._ (2019); Hsiao _et al._ (2017). In the absence of inertia, two types
of discontinuous shear thickening are predicted to occur depending on whether
the frictional particles are below or above their jamming point, where the
thickened phase either flows smoothly or is completely jammed Wyart and Cates
(2014). Curiously, in experiments an additional high-frequency and low-
amplitude shear orthogonal to the primary shearing flow has been shown to
suppress shear thickening Lin _et al._ (2016b).

The very high $\phi$ and high strain rate regime leads to even more surprising
observations. Experiments first revealed Ackerson and Pusey (1988); Ackerson
(1990), and BD simulations later confirmed Strating (1999), that particles
then form layers perpendicular to the shear gradient that slide over each
other. The system therefore keeps flowing with a viscosity that suddenly drops
at $\operatorname{P\kern-0.80002pte}_{\dot{\gamma}}\approx 10$. Hard spheres
under shear can also exhibit non-uniform flow profiles, so-called shear
banding. One mechanism is the formation of an arrested band due to small
variations of the local packing fraction that trigger the arrest of a much
bigger region of the flow Besseling _et al._ (2010). Other exotic phenomena
include the formation of twinned fcc crystals and sliding layers Haw _et al._
(1998), strings of particles in hard-sphere fluids Cheng _et al._ (2012) and
novel configurations which optimized packing Cohen _et al._ (2004). The non-
equilibrium phase behavior of a fluid of colloidal hard spheres under
oscillatory shear was investigated in real space with experiments on PMMA
colloidal suspensions and BD simulations as a function of the frequency of the
oscillations and $0\leq\operatorname{P\kern-0.80002pte}_{\dot{\gamma}}\leq
15$, displaying a shear-induced oscillating twinned fcc phase, a sliding layer
phase, a string phase and a tilted layer phase Besseling _et al._ (2012).

Shearing hard-sphere crystals opens yet another range of phenomena. These
include melting Wu _et al._ (2009), shear banding Cohen _et al._ (2004, 2006);
Dhont and Briels (2008), and in the case of confined crystals buckling
phenomena Schall _et al._ (2004). It is even possible to infer information
about stresses between defects in hard-sphere crystals through careful
analysis of particle trajectories Lin _et al._ (2016a).

###  XI.3 Microrheology

Trapping a colloidal probe with optical tweezers and imaging its stochastic
motion provides insights into the mechanical properties of the host material
Wilson _et al._ (2009); Puertas and Voigtmann (2014). This microrheological
technique is particularly important for (biological) materials that are
difficult to prepare in amounts sufficiently large for conventional
rheological studies Wilhelm (2008). The approach can also be used to resolve
local mechanical properties in inhomogeneous (soft) materials.

Figure 33: The microstructure $g(\bm{r})$ deforms around a driven probe with
(a) $\operatorname{P\kern-0.80002pte}_{v}=0.15$, (b)
$\operatorname{P\kern-0.80002pte}_{v}=0.5$, and (c)
$\operatorname{P\kern-0.80002pte}_{v}=1.5$. Reprinted from Squires and Brady,
Phys. Fluids 17, 073101 (2005) Squires and Brady (2005), with the permission
of AIP Publishing.

Depending on whether the probe is forced or not, one distinguishes active from
passive microrheology. The latter exploits the fluctuation-dissipation
theorem, which greatly simplifies the analysis, but analyzing the data from
driven probes requires a pre-existing model of the host material. In this
context, hard spheres have emerged as a particularly useful reference. Squires
and Brady (2005) provide a comprehensive analysis of the deformation of the
microstructure (Fig. 33) around a hard probe with diameter
$\sigma_{\text{probe}}$ forced through a bath of hard spheres with drift
Péclet number
$\operatorname{P\kern-0.80002pte}_{v}=v(\sigma+\sigma_{\text{probe}})/(2D_{0})$.
Here, $v$ is the probe speed. Analytical expressions for the viscosity are
derived in the limits of small and large $\operatorname{P\kern-0.80002pte}$
yielding an accurate extrapolation to intermediate Péclet numbers. Other than
fluid media, forcing a probe through a crystal Vossen (2004); Dullens and
Bechinger (2011) and a colloidal glass Habdas _et al._ (2004); Gazuz _et al._
(2009); Gruber _et al._ (2016) have been studied.

###  XI.4 Other out-of-equilibrium phenomena

Time-resolved confocal microscopy has been used to study other out-of-
equilibrium phenomena. One such example is the dynamics of colloidal particles
in externally created energy landscapes, such as optical potentials generated
by interfering laser beams111111Due to the external forcing in these systems,
precise control over interparticle interactions to be hard is often
unnecessary, but interactions nevertheless need to be sufficiently short-
ranged.. Other significant examples include ordering of polystyrene spheres in
quasiperiodic patterned potentials Mikhael _et al._ (2008), driving colloidal
monolayers through time-dependent fields Bohlein _et al._ (2012); Brazda _et
al._ (2018), diffusion in random landscapes Evers _et al._ (2013), and the
transmission of forces through dense colloidal aggregates Williams _et al._
(2016).

Following the trajectories of single colloidal particles has also been
instrumental in experimentally verifying fluctuation theorems Wang _et al._
(2002); Carberry _et al._ (2004); Blickle _et al._ (2006), and in verifying a
class of exact relations from nonequilibrium statistical physics starting with
Jarzynski’s seminal work relation Jarzynski (1997). (See Seifert (2012) for a
comprehensive theoretical review.) In a nutshell, thermodynamic (current-like)
quantities like work and heat can be extended to single stochastic
trajectories, which are therefore described by probability distributions
constrained by the (near-universal) fluctuation theorems. These relations have
also been tested in colloidal suspensions Gomez-Solano _et al._ (2015), and
recent experimental results, including with colloidal particles, have been
reviewed by Ciliberto (2017).

More recently, “active” colloidal particles have moved into focus. A wealth of
interfacial phoretic mechanisms between molecular and colloidal solutes can be
exploited to generate self-sustained gradients that move with the colloidal
particles and lead to directed motion (revealed as a correlation in
displacements absent in passively diffusing particles). One phenomenon where
simulations of active hard disks (and spheres) have been, and still are,
instrumental is the coexistence of dense and dilute regions following
_motility-induced phase separation_ Cates and Tailleur (2015). Given the
absence of cohesive forces, these processes are genuine far-from-equilibrium
phase transitions. In light of the many reviews that have been devoted to this
rich and fast moving field Elgeti _et al._ (2015); Bechinger _et al._ (2016);
Marchetti _et al._ (2013); Janssen (2019), we will not, however, delve further
here.

##  XII Hard-sphere glasses and their formation

At its core, the glass problem consists of understanding how the equilibrium
dynamics of a liquid that exhibits no obvious structural change grows to be so
sluggish as to freeze particles in place Berthier and Biroli (2011). In a
materials context, while glasses are traditionally associated with cohesive
supercooled liquids, dense hard spheres have also long played a key conceptual
role. The geometrical simplicity of the constituent particles, in particular,
enables a variety of material and theoretical approaches to explore the
underlying physics. Therefore, since the 1976 numerical simulations of Gordon
_et al._ (1976) and Woodcock (1976), hard spheres have been fueling the debate
about the nature and location of the glass transition, initially with computer
simulations Frenkel and McTague (1980); Woodcock (1981); Woodcock and Angell
(1981); Speedy (1987), and later with experimental advances, such as dynamic
light scattering (Sec. IV.1) Pusey and van Megen (1987); van Megen _et al._
(1991b); van Megen and Pusey (1991); van Megen _et al._ (1991a), real-space
analysis (Sec. IV.2) van Blaaderen and Wiltzius (1995); Weeks _et al._ (2000);
Kegel and Van Blaaderen (2000); Leocmach and Tanaka (2012); Hallett _et al._
(2018), and rheology Mason and Weitz (1995); Bonn _et al._ (2017).

In this section, we review the contribution to our understanding of the glass
transition made by work on hard spheres. We begin by providing some context
(Sec. XII.1), before discussing the early work which mainly used light
scattering (Sec. XII.2), developments made possible by real-space imaging
(Sec. XII.3) and more recent work which has approached rather longer time
scales than some of the earlier studies (Sec. XII.4). We also consider the
effects of confinement (Sec. XII.5) as well as its 2d limit (Sec. XII.6). We
then move on to the related phenomenon of jamming (Sec. XII.8, and to specific
properties of glasses, such as their vibrational behavior (Sec. XII.7), aging
(Sec. XII.9) and rheology (Sec. XII.10).

Inevitably, this section—like many others in this review—is limited in scope.
We here exclusively discuss contributions to our understanding of glasses
(which has certain universal qualities Berthier and Biroli (2011)) that result
from hard-sphere studies. We therefore do not include important work that used
other systems. For a more complete picture, we refer the reader to reviews on
the glass transition Berthier and Biroli (2011) and on specific aspects such
as mean field/high dimensional work Charbonneau _et al._ (2017), dynamical
heterogeneity Berthier _et al._ (2011), dynamical and structural length scales
Karmakar _et al._ (2014), local structure Royall and Williams (2015), jamming
Torquato and Stillinger (2010); Liu and Nagel (2010); van Hecke (2010);
Charbonneau _et al._ (2017); Arceri _et al._ (2022), and aging Arceri and
Corwin (2020) and rheology Bonn _et al._ (2017).

In this section, we also make a departure from the nomenclature used in the
rest of this review. Hard spheres, of course, have no liquid phase (Fig. 3,
Sec. VII.1). However, because the glass transition in molecular systems is
typically driven by cooling a liquid below its melting point, such that it is
_supercooled_ , the analogous behavior in hard spheres is to compress the
system beyond its freezing point $P_{f}$ (or $Z_{f}$, for the reduced pressure
$Z=\beta P/\rho$) as the control parameter Berthier and Witten (2009).
Therefore, in keeping with much of the hard-sphere glass literature, and to
emphasize the analogy with molecular systems, we refer to fluid state points
thus compressed as supercooled liquids (instead of supercompressed fluids).
Similarly, we shall refer to state points far beyond $Z_{f}$ as being “deeply
supercooled”. While the (reduced) pressure $Z$ emphasizes the link with
molecular systems Berthier and Witten (2009), some of the literature expresses
state points in terms of the volume fraction. Here we prefer $Z$, but when a
particular reference uses $\phi$, we will also often use it.

###  XII.1 Historical theoretical developments and persistent challenges

Figure 34: Roadmap to the glass transition in hard spheres. Configurational
entropy $S_{\mathrm{conf}}$ as a function of the inverse pressure (osmotic
pressure in the case of colloidal systems). Typically, $S_{\mathrm{conf}}$ of
molecular liquids decreases faster than that of crystals as a function of
temperature. The configurational entropy of the crystal is non-zero due to
defects. At some pressure $P_{k}$, the (supercooled) liquid configurational
entropy would fall below that of the crystal. $P_{g}$ is the operational glass
transition, which is mapped from that of molecules where the structural
relaxation time $\tau_{\alpha}$ typically exceeds that of the liquid by
$10^{14}$ (corresponding to 100s for molecules). Further compression (fast
relative to $\tau_{\alpha}$) leads to non-equilibrium states (colored lines).
$P_{\mathrm{mct}}$ is the MCT crossover; $P_{\text{f}}$ is the freezing point.
We plot this figure in terms of pressure to emphasise the connection with
molecular systems Cavagna (2009); Kauzmann (1948).

Understanding the glass transition has been an active area of research for
over a century Berthier and Biroli (2011), but as it pertains to hard spheres,
the 1980s saw the independent emergence of two major microscopic theories of
the glass transition: a density-functional description of amorphous solids
Singh _et al._ (1985); Wolynes (1985); Baus and Colot (1986); Lowen (1990),
and kinetic-theory based Dorfman _et al._ (2021) mode-coupling theory (MCT) of
glasses Bengtzelius _et al._ (1984); Leutheusser (1984); Kirkpatrick (1985);
Barrat _et al._ (1989); Fuchs _et al._ (1992); Goetze (2009). Using spin glass
models as inspiration, the random first-order theory (RFOT) of glasses was
then proposed to unify the two descriptions Kirkpatrick and Wolynes (1987);
Lubchenko and Wolynes (2007). In short, this framework posits the existence of
a dynamical (or mode-coupling) transition at which ergodicity is lost–at
reduced pressure $Z_{d}$ (or $Z_{\mathrm{mct}}$)–due to phase space then
subdividing into a number of metastable states so large that it is accompanied
by a thermodynamic entropy contribution, i.e., a configurational entropy (or
complexity).

While this description is now understood to hold in high dimensions
Charbonneau _et al._ (2017); Parisi _et al._ (2020), in 2d and 3d, the
situation is more subtle. This largely reflects the fact that this ideal
scenario neglects _activated_ processes that can restore ergodicity beyond the
MCT transition (thus turning it into a crossover) perhaps even up to the point
at which the configurational entropy $S_{\mathrm{conf}}=S-S_{\mathrm{vib}}$ of
the supercooled liquid is extrapolated to become equal to that of the crystal
Kauzmann (1948), where $S_{\mathrm{vib}}$ is the vibrational contribution to
the entropy Berthier _et al._ (2019c). In the case of molecular liquids, this
would happen at a finite (Kauzmann) temperature. The analogous scenario for
hard spheres is sketched in Fig. 34 with the configurational entropy of the
supercooled liquid becoming equal to that of the crystal at a certain pressure
$P_{k}$. (Here we use pressure to connect with the crystal in the spirit of
work with molecular systems Cavagna (2009); Kauzmann (1948).) In Sec. XII.4,
we see that measurements of the configurational entropy using advanced
computer simulation methods confirm this picture Berthier _et al._ (2017).

Now the drop in configurational entropy is understood to be accompanied by an
increase in a structural length scale Lubchenko and Wolynes (2007), which
would diverge at $P_{k}$ (or $Z_{k}$) (although the transition may be avoided
Stillinger _et al._ (2001); Royall _et al._ (2018a)). Such a length scale may
be accessed by _point-to-set_ correlations, which measure the impact of a
frozen, disordered boundary (often spherical) on the fluid behavior away from
it (often the center of the sphere) Biroli _et al._ (2008). Cooperative
relaxation at deep supercooling (for $P>P_{\mathrm{mct}}$ or
$Z>Z_{\mathrm{mct}}$) leads to _dynamical heterogeneity_ with some regions
relaxing slower or faster than one another in supercooled liquids Berthier _et
al._ (2011). The size of these dynamically heterogeneous regions is then
characterized by a dynamical length scale.

Despite its very considerable success and broad epistemic reach, RFOT is far
from being the only theory of the glass transition to have been formulated and
used to interpret results from colloidal hard spheres. A well-studied
structural mechanism for dynamical arrest has its roots in the work of Frank
(1952), who hypothesized that local five-fold symmetric arrangements of
particles would inhibit crystallization. This line of thought was later
extended into a dynamical theory built on geometric frustration which imagines