# Diffusion, phase behavior and gelation in a two-dimensional layer of
colloids in osmotic equilibrium with a polymer reservoir

Sam E. Griffiths  School of Physics and Astronomy, The University of
Edinburgh, Edinburgh EH9 3FD, United Kingdom  Nick Koumakis  School of Physics
and Astronomy, The University of Edinburgh, Edinburgh EH9 3FD, United Kingdom
Aidan T. Brown  School of Physics and Astronomy, The University of Edinburgh,
Edinburgh EH9 3FD, United Kingdom  Teun Vissers  School of Physics and
Astronomy, The University of Edinburgh, Edinburgh EH9 3FD, United Kingdom
Patrick B. Warren  Hartree Centre, Science and Technology Facilities Council
(STFC), Sci-Tech Daresbury, Warrington, WA4 4AD, United Kingdom  Wilson C. K.
Poon  School of Physics and Astronomy, The University of Edinburgh, Edinburgh
EH9 3FD, United Kingdom

###### Abstract

The addition of enough non-adsorbing polymer to an otherwise stable colloidal
suspension gives rise to a variety of phase behavior and kinetic arrest due to
the depletion attraction induced between the colloids by the polymers. We
report a study of these phenomena in a two-dimensional layer of colloids. The
three-dimensional phenomenology of crystal-fluid coexistence is reproduced,
but gelation takes a novel form, in which the strands in the gel structure are
locally crystalline. We compare our findings with a previous simulation and
theory, and find substantial agreement.

††preprint: AIP/123-QED

##  I INTRODUCTION

Adding non-adsorbing polymers to a suspension of hard-sphere colloids (radius
$R$, volume fraction $\phi_{\mathrm{c}}$) induces a depletion attraction
between the particles.Asakura and Oosawa (1954); Vrij (1976); Lekkerkerker and
Tuinier (2011) Exclusion of polymers from the space between two nearby
particles leaves an unbalanced osmotic pressure pushing them together. The
depth and range of the depletion attraction $U_{\rm dep}(r)$ between two
particles with centre-to-centre distance of $r$ is proportional to the polymer
activity, $a_{p}$, and the polymer’s radius of gyration, $r_{\mathrm{g}}$,
respectively. When the size ratio $\xi=r_{\mathrm{g}}/R\lesssim 0.30\pm
0.05$,Lekkerkerker  _et al._ (1992); Ilett  _et al._ (1995) the equilibrium
$(\phi_{\mathrm{c}},a_{p})$ phase diagram displays an expanded region of
fluid-crystal (F-X) coexistence at $a_{p}\gtrsim O(10^{-1})$,Poon (2002) which
occurs for $0.494<\phi_{\mathrm{c}}<0.545$ at $a_{p}=0$.

‘Buried’ within the equilibrium F-X coexistence region is a metastable gas-
liquid (G-L) coexistence binodal that terminates at a critical point. A
homogeneous sample inside this binodal should in principal first phase
separate into metastable coexisting G-L phase, before further separating into
equilibrium F-X coexistence.Evans, Poon, and Cates (1997) This scenario is,
however, seldom observed, because kinetic arrest intervenes.

For $0.1\lesssim\phi_{\mathrm{c}}\lesssim 0.3$, samples inside the metastable
G-L coexistence region will phase separate by spinodal decomposition into a
bicontinuous structure. With time, this structure coarsens in length scale,
and the concentration difference between the two phases increases, until the
concentration of the liquid phase crosses the ‘attractive glass transition’.
linePham  _et al._ (2002) The texture arrests and the system becomes a
gel.Zaccarelli (2007),111Note that if the colloids are too polydisperse to
crystallise, this scenario still holds, except that, now, G-L coexistence is
no longer metastable, but is the sole equilibrium thermodynamic phase
transition in the system.

Such ‘depletion gels’ have been intensively studied for some time, but
mysteries remain, perhaps especially how they age with time. In some cases, a
depletion gel can undergo sudden gravitational collapse after an apparently
quiescent period in which little seems to happen macroscopically.Poon, Pirie,
and Pusey (1995); Harich  _et al._ (2016); de Graaf  _et al._ (2019);
Cipelletti, Martens, and Ramos (2020) This and other aging phenomena are
expected on thermodynamic grounds. A depletion gel is metastable. There is
therefore a driving force for evolution towards the lowest free energy state,
which is F-X coexistence.

What we have summarised so far pertains to bulk colloid-polymer mixtures. At
first sight, there is little incentive to study two-dimensional (2D) systems:
presumably, any difference to bulk behavior would be merely quantitative.
However, this intuition turns out to be incorrect.

An early 2D study extended the ‘primitive’ theory for bulk phase
behaviorLekkerkerker  _et al._ (1992) to calculate the phase diagram of a bulk
colloid-polymer mixture in the presence of a hard wall.Poon and Warren (1994)
The theory predicts that depletion-induced wall adsorption induces wall
freezing (= crystallization) at depletant concentrations below the bulk F-X
coexistence boundary. An interesting subtlety, Fig. 1(a), is that both
colloids and depletants at the surface are in osmotic equilibrium with the
bulk, with which they can exchange both species. Thus, unlike in the bulk,
which is a canonical ensemble where F-X coexistence is possible, surface
crystallization occurs in a grand-canonical ensemble, so that the crystal
fraction jumps from 0 to 1 at the critical bulk polymer concentration (at
fixed colloid volume fraction). Experiments using large
($R_{\mathrm{L}}=$0.23\text{\,}\mathrm{\SIUnitSymbolMicro m}$$) and small
($R_{\mathrm{S}}=$0.035\text{\,}\mathrm{\SIUnitSymbolMicro m}$$) charge-
stabilised polystyrene spheres (screening length
$\approx$5\text{\,}\mathrm{nm}$$) in which the small spheres act as depletants
Dinsmore  _et al._ (1997) confirmed this feature, and found reasonable
agreement with theory for the wall crystallization boundary.

The sedimentation height $z_{0}$ of the large spheres in these experiments,
defined such that the number of particles of radius $R$ in a dilute suspension
decreases with height according to $n(z)=n(0)e^{-(z-R)/z_{0}}$, is
$z_{0}=$160\text{\,}\mathrm{\SIUnitSymbolMicro m}$\gg R_{\mathrm{L}}$, as
assumed in the theory.Poon and Warren (1994) Two other experiments used much
lower $z_{0}/R$. Savage et al. Savage, Pei, and Dinsmore (2012) used
$R=$0.7\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ polystyrene spheres
($z_{0}\approx 8\,R$) with non-ionic surfactant micelles as depletants.
HobbieHobbie (1998) studied binary polystyrene colloids with
$R_{\mathrm{L}}=$1.45\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ ($z_{0}\approx
0.4\,R_{\mathrm{L}}$) and
$R_{\mathrm{S}}=$0.107\text{\,}\mathrm{\SIUnitSymbolMicro m}$$. In each case,
gravity and the particle-wall depletion attraction sufficed to confine all of
the large particles to an effectively 2D layer (at low enough bulk
concentration).

Now, the surface layer is a canonical ensemble of a fixed number of large
particles, but a grand-canonical ensemble of the smaller depletants, which
freely exchanges with the bulk, Fig. 1(b). In this ‘semi-grand canonical’
scenario,Lekkerkerker (1990) wall crystallization takes the form of F-X
coexistence, with an increasing crystal fraction as the bulk concentration of
depletants (which controls the surface depletant chemical potential)
increases. This is indeed seen in experiments,Hobbie (1998); Savage, Pei, and
Dinsmore (2012) where crystal nucleation is two-stepped, proceeding via an
intermediate gas-liquid phase separation, as predicted by theory.Evans, Poon,
and Cates (1997). Such kinetics is seldom seen in bulk, where the metastable
gas-liquid critical point nevertheless enhances crystal nucleation.ten Wolde
and Frenkel (1997)

Figure 1: Different ensembles for 2D colloid-depletant mixtures, with the
depletant represented as polymer coils, but they can also be smaller colloids
or micelles. (a) Colloids with $z_{0}\gg R$. Grand-canonical for both colloids
and depletants in the wall layer (bottom),Dinsmore  _et al._ (1997) with
exchanges of both species with the bulk (blue and red double arrows). (b)
Colloids with $z_{0}\lesssim R$. Grand-canonical for the depletants in the
surface layer, which exchange with the bulk (blue double arrow); but canonical
for the colloids, which are all sedimented, with no exchange with the
bulk.Hobbie (1998); Savage, Pei, and Dinsmore (2012)

Cerdà et al. have simulated a system mimicking a 2D colloid-polymer mixture:
nearly-hard disks interacting via an Asakura-Oosawa (AO) potential, which is
widely used to model polymer-induced depletion.Lekkerkerker and Tuinier (2011)
They worked at a size ratio $\xi=0.1$, and probed the behavior at a surface
colloid area fraction $\eta_{\rm c}=0.157$. Cerdà et al. found F-X coexistence
when the contact value of the AO potential, $U_{0}\geq
3.130\,k_{\mathrm{B}}T$. At the even higher $U_{0}=7\,k_{\mathrm{B}}T$, they
observed ramified clusters with fractal dimension $d_{\mathrm{f}}\approx 1.4$,
the exponent for diffusion-limited cluster aggregation (DLCA). However, the
local structure of the cluster strands is, unlike DLCA, clearly crystalline.
This finding has yet to be confirmed by experiments.

Here, we study a 2D colloid-polymer mixture in which we quantify how two-
dimensionality by comparing diffusivities next to a wall against Faxén’s
prediction.Faxén (1921); HB73; Chio and Tse (2020) The equilibrium phase
diagram is obtained and compared with theory.Poon and Warren (1994) We confirm
the locally-crystalline nature of the ramified clusters at the highest
attractions.Cerdà  _et al._ (2004) Quantitatively, the time dependence of the
average cluster size $\bar{n}\sim t^{z}$, with $0.5\lesssim z\lesssim 0.6$. We
discuss the possible origin of this exponent in the screened near-wall
hydrodynamics of clusters, and propose why the local structure of our ramified
clusters are crystalline, in striking contrast to the case of 3D depletion
gels.Griffiths, Turci, and Royall (2017)

##  II MATERIALS AND METHODS

We used $1.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$ diameter ($2R$) silica
spheres (Bangs Laboratories, density
$\rho\approx$2\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$$). The work we
report is part of a larger programGriffiths (2021) studying the effect of
self-propelled particles on colloidal gels, where the self-propelled particles
are motile Escherichia coli bacteria,Schwarz-Linek  _et al._ (2016) for which
we need to know the behavior of the cell-free system. We therefore dispersed
our colloids in a phosphate motility buffer (PMB) commonly used to study
motile E. coli ($6.2\,\mathrm{mM}$ K2HPO4, $3.8\,\mathrm{mM}$ KH2PO4 and
$0.1\,\mathrm{mM}$ EDTA at $p\ce{H}\approx 7.5$). We removed the NaCl usually
included in standard PMB to limit the ionic strength to $I=22.4\,\mathrm{mM}$
(screening length $2\text{\,}\mathrm{nm}$). Under these conditions, our
colloids do not visibly aggregate, so that any residual interparticle
attraction is $\lesssim k_{\mathrm{B}}T$. The colloids sediment rapidly as
single particles to the bottom of our sample chambers. Bulk volume fractions
of $10^{-4}\lesssim\phi_{\mathrm{c}}\lesssim 2\times 10^{-3}$ gave surface
area fractions of $0.04\lesssim\eta_{\rm c}\lesssim 0.8$.

To induce depletion attraction, we added sodium polystyrene sulphonate (NaPSS,
$M_{\rm w}\approx${10}^{6}\text{\,}\mathrm{Da}$$; Sigma Aldrich, used as
purchased). NaPSS behaves as a ideal neutral polymer in a medium with
$I=3.1\,\mathrm{M}$, while at $I=0.15\,\mathrm{M}$ we have a good solvent.Wang
and Yu (1988) We therefore expect swollen coils in our PMB. We estimated the
overlap concentration as the inverse intrinsic viscosity,Rubinstein and Colby
(2003) which came from extrapolating Kraemer and Huggins plots of the measured
viscosity as a function of polymer concentration,Rubinstein and Colby (2003);
Schwarz-Linek  _et al._ (2012) giving $c_{p}^{\ast}\approx
0.44\,\mathrm{wt}\,\%$. From this, we estimate a radius of gyration,
$r_{\mathrm{g}}$, of NaPSS in PMB using

$c_{p}^{\ast}=\frac{3M_{\rm w}}{4\pi r_{\mathrm{g}}^{3}N_{\rm A}},$ (1)

where $N_{\rm A}$ is Avogadro’s constant, giving
$r_{\mathrm{g}}\approx$45\text{\,}\mathrm{nm}$$, and a colloid:polymer size
ratio of $\xi\approx 0.06$.

The contact depletion attraction, $U_{0}$, is proportional to the polymer
activity, $a_{p}$, which in turn scales as the polymer concentration in a
reservoir in osmotic equilibrium with the colloid-polymer mixture.Lekkerkerker
_et al._ (1992); Lekkerkerker and Tuinier (2011) In our system, this is well
approximated by the bulk polymer concentration. This we report as a polymer
volume fraction, $\phi_{\mathrm{p}}$, estimated using a coil volume of
${4\pi}r_{\mathrm{g}}^{3}/3$, so that overlap, Eq. (1), corresponds to
$\phi_{\mathrm{p}}=1$. We work in the range $0\leq\phi_{\mathrm{p}}\leq 0.4$.

We sealed samples into $400\text{\,}\mathrm{\SIUnitSymbolMicro m}$-high glass
capillaries and aged them for $\approx$2\text{\,}\mathrm{h}$$ before video
recording in a Nikon Ti-Eclipse inverted microscope, typically using a $\times
50$ objective to resolve single particles. We tracked particles using a
Mikrotron high-speed camera (MC 1362) to determine bulk and near-wall
diffusivities. A Hamamatsu Orca 4.0 CMOS camera was used to identify clusters,
from which we obtained information on nearest neighbors, defined as all of a
particle’s neighbors whose centres are within
$2R+2r_{\mathrm{g}}+$0.1\text{\,}\mathrm{\SIUnitSymbolMicro
m}$\approx$1.7\text{\,}\mathrm{\SIUnitSymbolMicro m}$$, where the
$0.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ accommodates polydispersity. We
also used a $\times 10$ ($\mathrm{N.A.}=3$) objective to obtain images in
which single particles were not resolved. We thresholded these images,
identified clusters and calculated their areas using custom software. Results
for small clusters from the two methods agree up to a scaling constant between
cluster area in (pixel)2 and particle number. The calibrated low-resolution
method gives better statistics, especially for larger clusters.

##  III Theory

Previously, a mean-field van der Waals free-volume theory was used to predict
the wall freezing transition in a grand canonical system.Dinsmore  _et al._
(1997) In our semi-grand canonical ensemble, Fig. 1(b), the approach
simplifies considerably, as there is no need to equilibrate colloids in the
wall layer and the bulk. The semi-grand canonical free energy $F$ with ideal
depletants is

$F=F_{0}-\rho_{\mathrm{d}}k_{\mathrm{B}}T\langle{V_{\text{free}}}\rangle_{0},$
(2)

where $F_{0}$ is the free energy of a reference system of 2D-confined hard
spheres, and $\langle V_{\text{free}}\rangle$ is the volume available to
depletants of bulk number density $\rho_{\mathrm{d}}$ averaged over the
positions of the colloids in the reference system. The depletants are mutually
non-interacting (ideal), and exert an osmotic pressure
$\rho_{\mathrm{d}}k_{\mathrm{B}}T$ in the bulk, which acts as a reservoir.
They cannot approach $\leq r_{\mathrm{d}}$ to the surface of a colloid.

For $F_{0}$ we use published free energies of the hard disc fluid and
hexagonally-ordered crystalZeng and Oxtoby (1990) that account well for two-
dimensional freezing. For the former (cf. Eq. (18) in Ref. Zeng and Oxtoby,
1990):

$\frac{F_{0}^{\mathrm{(fluid)}}}{N_{\mathrm{c}}k_{\mathrm{B}}T}=\ln\left[\frac{4\eta_{\mathrm{c}}}{\pi}\right]-1-\ln(1-\eta_{\mathrm{c}})+\frac{\eta_{\mathrm{c}}}{1-\eta_{\mathrm{c}}},$
(3)

where $\eta_{\mathrm{c}}$ is the area fraction of colloids. After the first
two ideal terms, the next two terms give the excess free energy per particle
after Rosenfeld.Rosenfeld (1990) For the 2D crystal, we follow Hall’s
procedure for 3D hard-sphere crystalsHall (1972) and fitted

$\frac{F_{0}^{\mathrm{(crystal)}}}{N_{\mathrm{c}}k_{\mathrm{B}}T}=\ln\left[\frac{4\eta_{\mathrm{c}}}{\pi}\right]-1+c_{0}+c_{1}\beta_{\mathrm{c}}-2\ln\beta_{\mathrm{c}}\,$
(4)

to the 2D crystal curve in Fig. 2 of Ref. Zeng and Oxtoby, 1990 Here,
$\beta_{\mathrm{c}}=4(1-\eta_{\mathrm{c}}/\eta_{\mathrm{max}})$ where
$\eta_{\mathrm{max}}=\pi\sqrt{3}/6\approx 0.907$ is hexagonal close-packing.
Least-squares fitting gives $c_{0}\approx 3.08$ and $c_{1}\approx 0.30$. As
expected, as $\eta_{\mathrm{c}}\to\eta_{\mathrm{max}}$, $\beta_{\mathrm{c}}\to
0$ and $F_{0}^{\rm crystal}\to\infty$. With Eqs. (3) and (4) we find
coexistence between 2D F-X coexistence for $0.670<\eta_{\mathrm{c}}<0.732$,
agreeing with Ref. Zeng and Oxtoby, 1990.

For the free volume we use the standard expansion,Lekkerkerker and Tuinier
(2011)

$V_{\text{free}}=\mathrm{const}-N_{\mathrm{c}}V_{1}+{\textstyle{\displaystyle\sum}_{\>i>j}}\>V_{2}(r_{ij})+\dots$
(5)

where $V_{1}$ is the excluded volume for an isolated colloidal particle
allowing for the overlap with the depletion layer at the wall, $V_{2}(r)$ is
the overlap of the excluded volumes for a pair of colloids with centre-to-
centre distance $r$, and we sum over all pairs of particles $\\{i,j\\}$. We do
not require first two terms: they only contribute constants to the pressure
and chemical potential that do not affect phase behaviour in the semi-grand
ensemble. This is in contrast to the wall freezing transition in a bulk
colloidal suspension where the absolute value of the colloid chemical
potential is required, including the contribution from $V_{1}$, to solve for
the bulk-wall equilibria.Dinsmore  _et al._ (1997)

Figure 2: Phase behaviour of 2D-confined colloids in the presence of non-
adsorbing ideal depletants of exclusion radius $r_{\mathrm{d}}$ in the semi-
grand canonical ensemble for $\xi=r_{\mathrm{d}}/R=0.06$, calculated using the
approach described in the text. Horizontal tie lines span the F-X coexistence
region. The dashed line is the limit of thermodynamic stability of the fluid
phase (spinodal). The chained line is the fluid binodal estimated using the
simplified approach in Appendix A.

Despite 2D colloidal confinement, the volumes in Eq. (5) are 3D, which
precludes the use of scaled particle theory to estimate
$\langle{V_{\text{free}}}\rangle_{0}$. However, there will be no three-body
overlaps of the depletion layers of three spheres if
$r_{\mathrm{d}}/R=\xi<2/\sqrt{3}-1\approx 0.1547$,Reiss, Frisch, and Lebowitz
(1959). Moreover, there is no overlap of the depletion layers of two particles
and the wall if $\xi<1/4$. Both conditions are satisfied in our case. So, we
terminate Eq. (5) at the 2-body term and use the standard AO result

$V_{2}(r)=\frac{4\pi(R+r_{\mathrm{d}})^{3}}{3}\left(1-\frac{3r}{4(R+r_{\mathrm{d}})}+\frac{r^{3}}{16(R+r_{\mathrm{d}})^{3}}\right)\,.$
(6)

To average $V_{\text{free}}$ over the reference system, we write

$\frac{\langle{V_{\text{free}}}\rangle_{0}}{A}=\mathrm{const}-\frac{\eta_{\mathrm{c}}V_{1}}{\pi
R^{2}}+\frac{1}{2}\Bigl{(}\frac{\eta_{\mathrm{c}}}{\pi
R^{2}}\Bigr{)}^{2}\int_{2R}^{2(R+r_{\mathrm{d}})}\\!\\!2\pi
r\,{\mathrm{d}}r\,V_{2}(r)\,g(r)$ (7)

where $g(r)$ is the radial distribution function of hard discs in the
reference system, for which we we use the heuristic approximation

$g(r)=g_{c}+(1-g_{c})\frac{r-2R}{\lambda R},$ (8)

where $\lambda$ captures the decay of $g(r)$ from $g_{c}\equiv g(2R)$ at
contact. We use $g_{c}=(2-\eta_{\mathrm{c}})/2(1-\eta_{\mathrm{c}})^{2}$ from
scaled particle theory,Lebowitz, Helfand, and Praestgaard (1965) consistent
with Eq. (3) and the sum rule for the pressure in hard disc systems.Stones,
Dullens, and Aarts (2018) For simplicity we use Eq. (8) in both the fluid and
hexagonally-ordered crystal. For the fluid phase we use $\lambda=0.5$ (our
results are insensitive to $0.2\leq\lambda\leq 1$). This corresponds to the
typical contact peak in the radial distribution function in a dense
fluid.Barker and Henderson (1976) For the hexagonally-ordered crystal phase,
$\lambda=0.1$ ensures that the effective co-ordination number is around 6 at
$\eta_{\mathrm{c}}\approx 0.85$.

Substituting Eq. (8) into Eq. (7) we find a dimensionless specific (i.e. per
unit area) free volume

$\frac{\langle{V_{\text{free}}}\rangle_{0}}{A(2R)}=\mathrm{const}-\frac{V_{1}\eta_{\mathrm{c}}}{2\pi
R^{3}}+\xi^{3}(Pg_{c}+Q)\eta_{\mathrm{c}}^{2}\,.$ (9)

where $Q=2\xi(1+6\xi/5+\xi^{2}\\!/3)/3\lambda$ and
$P=4(1+\xi+\xi^{2}\\!/5)/3-Q$. Finally, combining Eqs. (2) and (9), we find a
dimensionless specific semi-grand free energy

$\frac{(2R)^{2}F}{Ak_{\mathrm{B}}T}=\frac{4\eta_{\mathrm{c}}F_{0}}{\pi
N_{\mathrm{c}}k_{\mathrm{B}}T}-\rho_{\mathrm{p}}(2r_{\mathrm{d}})^{3}(Pg_{c}+Q)\eta_{\mathrm{c}}^{2}\,,$
(10)

The constant and term linear in $\eta_{\mathrm{c}}$ have been dropped from Eq.
(9), and the multiplicative factor $\xi^{3}$ in the final term yields the
$(2r_{\mathrm{d}})^{3}$ which non-dimensionalises $\rho_{\mathrm{p}}$ in the
above.

Equation (10) with Eqs. (3) and (4) are solved numerically for phase
coexistence as $\rho_{\mathrm{p}}$ increases. A typical result is shown in
Fig. 2 for size ratio $\xi=0.06$. Ideal depletants broaden the two-phase
region as in the bulk,Lekkerkerker and Tuinier (2011); Lekkerkerker  _et al._
(1992) with the onset of the effect being at around
$\phi_{\mathrm{d}}=\pi\rho_{\mathrm{d}}(2r_{\mathrm{d}})^{3}\\!/6\lesssim
0.1$. We later compare this prediction to our experiments by identifying
$(r_{\mathrm{d}},\phi_{\mathrm{d}})$ with
$(r_{\mathrm{g}},\phi_{\mathrm{p}})$.

The limit of thermodynamic stability in the fluid can be found as the locus of
points where ${\mathrm{d}}^{2}\\!F/\mathrm{d}\eta_{\mathrm{c}}^{2}$ vanishes.
This gives the dashed ‘spinodal’ line in Fig. 2 (the fluid is unstable to the
right of the line). Finally, the chained line shows the prediction of the
fluid binodal using a simplified approach in which we only take into account
the cohesive polymer-induced AO ‘bond’ energies in the crystal (see Appendix
A).

##  IV Experimental results

###  IV.1 2D Confinement

At a density difference of
$\Delta\rho\approx${10}^{3}\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$$
with PMB, our particles sediment at
$\approx$1\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{\,}{\mathrm{s}}^{-1}$$.
A dilute suspension takes $\gtrsim$10\text{\,}\mathrm{min}$$ to establish
steady state in a $400\text{\,}\mathrm{\SIUnitSymbolMicro m}$-high
capillary,Weaver (1926) to give a height-dependent particle density (at
$\phi_{\mathrm{p}}=0$)Poon (2015) of $n(z)=n(0)e^{-(z-R)/z_{0}}$, with a
calculated $z_{0}=$0.23\text{\,}\mathrm{\SIUnitSymbolMicro m}$$. The
confinement is further increased when polymer is added to induce a depletion
attraction between the particles and the bottom capillary surface. So, we
study an essentially two dimensional (2D) layer of colloids at the bottom
capillary surface.

Figure 3: Diffusivity of our silica colloids next to a wall, $D$, as a
function of polymer volume fraction, $\phi_{\mathrm{p}}$. Experimental data
from particle tracking ( $\bullet$) and prediction using Eq. (17) (
$\blacksquare$). Inset: predicted average gap between particle (in nm) and
wall as a function of $\phi_{\mathrm{p}}$.

To quantify the confinement, we determined low $\phi_{\mathrm{c}}$
diffusivities by tracking. The bulk diffusivity
$D_{0}=$0.28\text{\,}{\mathrm{\SIUnitSymbolMicro
m}}^{2}\text{\,}{\mathrm{s}}^{-1}$$ calculated from the Stokes-Einstein
relation is reduced by near-wall hydrodynamics to

$D=\alpha D_{0},$ (11)

with $\alpha<1$, Fig. 3. Stronger attraction reduces $D$ by lowering the
particle-wall distance and so increases hydrodynamic hindrance. Faxén’s
approximate calculation predictsFaxén (1921); HB73; Chio and Tse (2020)

$\displaystyle\alpha(z)=\frac{D}{D_{0}}=1$ $\displaystyle-$
$\displaystyle\frac{9}{16}\frac{R}{z}+\frac{1}{8}\left(\frac{R}{z}\right)^{3}$
(12) $\displaystyle-$
$\displaystyle\frac{45}{256}\left(\frac{R}{z}\right)^{4}-\frac{1}{16}\left(\frac{R}{z}\right)^{5}.$

In the dilute limit, the probability of finding a particle at $z$ above the
wall scales as $e^{-U_{\rm tot}(z)/k_{\mathrm{B}}T}$, where the total
potential experienced by a near-wall particle is

$\frac{U_{\rm tot}^{\rm
w}(z)}{k_{\mathrm{B}}T}=\left\\{\begin{array}[]{ll}\infty&z\leq R\\\ U_{\rm
d}^{\rm w}(z)+U_{\rm g}(z)&R<z<r+2r_{\mathrm{g}}\\\ U_{\rm
g}(z)&R>z+2r_{\mathrm{g}}\\\ \end{array}\right..$ (13)

Here,

$\frac{U_{\rm g}(z)}{k_{\rm B}T}=\frac{z-R}{z_{0}}$ (14)

is the gravitational potential, and

$\frac{U_{\rm d}^{(w)}(z)}{k_{\rm
B}T}=-\left(\frac{3\phi_{\mathrm{p}}}{4r_{\mathrm{g}}^{3}}\right)R\left(2r_{\mathrm{g}}-z+R\right)^{2}$
(15)

is the particle-wall depletion interaction.Lekkerkerker and Tuinier (2011) We
define a dimensionless magnitude of the depletion potential at contact

$u_{0}^{(w)}=\left|\frac{U_{\rm d}^{(w)}(z=R)}{k_{\rm
B}T}\right|=\frac{3\phi_{\mathrm{p}}}{\xi}.$ (16)

The near-wall average diffusivity should be given by

$\displaystyle\bar{D}=\frac{D_{0}\int_{0}^{\infty}\alpha(z)e^{-U_{\rm
tot}^{\rm w}/k_{\mathrm{B}}T}\,{\rm d}z}{\int_{0}^{\infty}e^{-U_{\rm tot}^{\rm
w}/k_{\mathrm{B}}T}\,{\rm d}z}\times\frac{\mu_{0}}{\mu_{\rm p}},$ (17)

taking into account the increased viscosity of the polymer solution, $\mu_{\rm
p}$, relative to that of PMB, $\mu_{0}$, both of which we measured using
standard rheometry.

Equation (17) gives a good account of our data, Fig. 3. The calculated average
gap $\bar{z}-R$ between a particle and the surface, where $\bar{z}$ is
obtained from Eq. (17) with $\beta(z)$ replaced by $z$, drops from
$0.23\text{\,}\mathrm{\SIUnitSymbolMicro m}$ at $\phi_{\mathrm{p}}=0$ to
$2.5\text{\,}\mathrm{nm}$ at $\phi_{\mathrm{p}}=0.836$, Fig. 3 (inset). Our
system, especially at $\phi_{\mathrm{p}}>0$, is indeed highly confined to a
thin 2D layer, with the smallest gap being comparable to the expected
roughness of our particles.Rentsch  _et al._ (2006); Valmacco  _et al._ (2016)
This exercise also validates the use of the AO expression under our conditions
for the polymer-induced attraction between two spherical particles, for which
Eq. (15) is a special case (taking one sphere to have infinite radius).

Figure 4: (a) Equilibrium phase diagram of our system. (b) Single-phase fluid,
where the scale bar represents 50 $\rm{\mu m}$ and the inset scale bar
represents 10 $\rm{\mu m}$. (c) Fluid-crystal coexistence. (d) Ramified
cluster with crystalline local structure. (e) The bond orientation order
parameter, $\Psi_{6}$, as a function of polymer volume fraction,
$\phi_{\mathrm{p}}$ at three different values of $\eta_{\rm c}$ (see legend).
At each $\eta_{\mathrm{c}}$, three kinds of behavior are delineated,
demarcated by where $\Psi_{6}$ rises sharply, e.g. as quantified by where it
reaches 0.2 (dashed line), and by where $\Psi_{6}$ peaks (arrows).

###  IV.2 Phase diagram and cluster statistics

Visually, our system shows four regimes, Fig. 4(a). At
$(\phi_{\mathrm{p}}\lesssim 0.15,\eta_{\rm c}\lesssim 0.6)$, we find single
particles and transient cluster with $n\lesssim 4$ particles, Fig. 4(b). This
is a colloidal fluid. At higher $\eta_{\rm c}$ and $\phi_{\mathrm{p}}$. we
observe F-X coexistence of single particles and colloidal crystallites, Fig.
4(c). At the highest $\phi_{\mathrm{p}}$ and $\eta_{\rm c}\lesssim 0.75$, we
observe ramified clusters whose strands are crystalline, Fig. 4(d). Finally,
at $\eta_{\rm c}\gtrsim 0.75$, we observe a polycrystalline monolayer over a
small range of $\phi_{\mathrm{p}}$.

We quantified the degree of crystallinity of our system via the bond
orientation order parameter

$\Psi_{6}\equiv\bigl{\langle}|q_{6}(i)|^{2}\bigr{\rangle}.$ (18)

This is the all-particle average of the squared single-particle bond
orientation parameter, which, for particle $i$ and its set of $N_{i}$ nearest
neighbours, is given by

$q_{6}(i)=\frac{1}{N_{i}}\sum_{j=1}^{N_{i}}e^{i6\theta_{ij}},$ (19)

where $\theta_{ij}$ is the angle between the centre-to-centre vector from
particle $i$ to $j$ and an arbitrary fixed axis. Perfect crystallinity gives
$\Psi_{6}=1$.

Consider first the data for $\eta_{\rm c}=0.19$, Fig. 4(e). At
$\phi_{\mathrm{p}}=0$, $\Psi_{6}=0$. Increasing the polymer concentration, we
find $\Psi_{6}$ remaining low until $\phi_{\mathrm{p}}\approx 0.14$, whereupon
$\Psi_{6}$ increases sharply, evidencing transition to F-X coexistence. The
same behavior occurs at $\eta_{\mathrm{c}}=0.30$ and 0.39. We take the
transition to F-X coexistence in each case to be the first data point where
$\Psi_{6}\geq 0.2$. These points, open squares in Fig. 4(a), agree well with
the phase boundary determined by inspecting micrographs. At each
$\eta_{\mathrm{c}}$ studied, $\Psi_{6}$ reaches at sharp peak of $\lesssim
0.8$, and then falls. This is the onset of progressively more ramified
crystalline clusters with increasing fraction of edges. We argue below that at
long times, these clusters will percolation to form a gel. We therefore take
the peak position in $\Psi_{6}$ as the gel boundary, open triangles in Fig.
4(a), which again agrees with the visually-demarcated onset of ramified
clustering.

In the ramified cluster regime, we measured $c_{n}$, the number of clusters of
size $n$. Figure 5(a) shows how $c_{n}$ normalised by the total number of
cluster $c$ (so that $\sum_{n}c_{n}/c=1$) evolves with time for a sample with
$\eta_{\mathrm{c}}=0.2$ and $\phi_{\mathrm{p}}=0.326$. From these cluster size
distributions (CSDs), we extracted the average cluster size $\bar{n}(t)$, Fig.
5(c). It is clear that the clusters are growing throughout our experimental
time window, albeit with a reducing growth rate, suggesting that the system
will percolate at even longer times to form a gel. Our experimentally-
identified transition from F-X to ramified clusters is therefore the system’s
gel boundary.

Figure 5: Cluster statistics at $\eta_{\rm c}=0.2$ and
$\phi_{\mathrm{p}}=0.326$, which is in the ramified clusters regime. (a) Raw
normalised cluster size distribution, $c_{n}(t)$, at different times as in
legend (in seconds). (b) Scaled normalised cluster size distribution according
to Eq. (20). Blue = data sets for
$$10\text{\,}\mathrm{min}$<t<$100\text{\,}\mathrm{min}$$ Red = later data
sets. Continuous curve = fitting of blue data sets to Eq. (21), giving
$\lambda=-0.76$. (c) Evolution of average cluster size with time. The line has
slope $z=(1-\lambda)^{-1}=0.568$. The blue data in (b) are from the interval
between the vertical lines.

##  V Discussion

###  V.1 Equilibrium phase behavior

Apart from a difference in size ratio, our system should be directly
comparable to the simulations by Cerdà et al.,Cerdà  _et al._ (2004) who
studied nearly-hard particles at $\eta_{\rm c}=0.157$ interacting via an AO
attractive potential with dimensional range $\xi=0.1$ (we have $\xi=0.06$). At
$U_{0}=3.130\,k_{\mathrm{B}}T$, they observed a transition from ‘small
fluctuating clusters’ – our fluid phase – to large hexagonal close packed
clusters in a background of single particles – our F-X coexistence. Our
transition from fluid to F-X coexistence occurs at $\phi_{\mathrm{p}}\lesssim
0.16$, corresponding to $U_{0}^{\rm(c)}\lesssim 4k_{\mathrm{B}}T$, where
$U_{0}^{\rm(c)}={U_{0}^{\rm(w)}}/{2}$ (cf. Eq. (16)) is the dimensionless
colloid-colloid depletion attraction at contact.Lekkerkerker and Tuinier
(2011)

Next, we compare our experimental phase diagram with the theory outlined in
Section III. To do so, we need to relate the exclusion radius of the
depletant, $r_{\mathrm{d}}$, to a property of the polymer in our system.
Clearly, $r_{\mathrm{d}}=cr_{\mathrm{g}}$ for some dimensionless constant $c$.
The result for taking $c=1$, so that $\phi_{\mathrm{d}}=\phi_{\mathrm{p}}$, is
plotted in Fig. 4(a). (Taking $c=1$ means we reproduce Fig. 2 in Fig. 4(a).)
The theory gives a credible account of the shape of the F-X boundary for
$\eta_{\mathrm{c}}<0.67$ (where hard disks freeze). Quantitatively, the
agreement is much better than order of magnitude. Such agreement is
significant. Our theory predicts the onset of significant depletion effects on
the F-X transition occurs at $\rho_{\mathrm{p}}(2d)^{3}\sim O(1)$, because it
is this particular dimensionless combination that enters the free energy in
Eq. (10). If our model incorporates wrong physics, then another length, the
colloid radius $R$, may also enter into the non-dimensionalisation of
$\phi_{\mathrm{p}}$, potentially altering the predicted phase boundary by a
factor of $\xi$, $\xi^{2}$ or $\xi^{3}$. Given that our $\xi\sim O(10^{-1})$,
the good agreement we find without fine tuning $c$ confirms that depletion
correctly captures the essential physics of the phase behavior in our system.

It is interesting to compare our phase diagram with the 3D case. We do so via
the second virial coefficient, $B_{2}$, which is often used to compare results
for potentials of different shapes in the same spatial dimension.Noro and
Frenkel (2000) In our case, we normalise $B_{2}$ by the hard-sphere or hard-
disk values in 3D and 2D to give a dimensionless $b_{2}$.222We calculated
$b_{2}$ by numerical integration using an approximation to the AO potential
valid for $\xi\lesssim 0.1$ taken from Ref. Bergenholtz, Poon, and Fuchs, 2003
The F-X coexistence boundary for a 3D colloid-polymer mixture at $\phi_{\rm
c}<0.494$ (the onset of bulk crystallization) occurs at
$b_{2}^{\rm(3D)}\approx-0.67$, which is also the crystallization threshold for
many globular protein solutions. Poon (1997) Interestingly, our F-X boundary
at $\phi_{\mathrm{p}}\approx 0.16$, Fig. 4(a), corresponds to
$b_{2}^{\rm(2D)}\approx 0.72$. Backing out an equivalent
$\phi_{\mathrm{p}}=0.21$ for the $\xi=0.1$ system of Cerdà et al. via
$U_{0}^{\rm(c)}=\frac{3\phi_{\mathrm{p}}}{2\xi}$, we find that their F-X
transition occurs at $b_{2}^{\rm(2D)}=0.74$. This agrees with our experimental
value, but differs significantly from the 3D value of $\approx-0.67$. Thus, in
2D, the depletion attraction does not have to be as strong to bring the system
into F-X coexistence. In this sense, crystallization is easier in 2D.

Finally, we note that Li et al. reported an experimental study of a system of
sedimented PMMA colloids ($2R_{\rm L}=$3.27\text{\,}\mathrm{\SIUnitSymbolMicro
m}$$) with depletion attraction induced by smaller ($2R_{\rm
S}=$0.192\text{\,}\mathrm{\SIUnitSymbolMicro m}$$) PMMA colloids (giving
$\xi=0.05$).Li2019 In contrast to our experiments and theoretical prediction,
where the F-X coexistence region expands rather suddenly at
$\phi_{\mathrm{p}}\sim 10^{-1}$, Fig. 4(a), their F-X coexistence expands
gradually from the $\phi_{\mathrm{p}}=0$ region. We do not at present
understand the origins of this difference, and will not discuss their findings
further.

###  V.2 Gelation

Strikingly, our ramified clusters are locally crystalline, Fig. 4(d). This
contrasts with 3D gelation, in which the ramified clusters are locally
amorphous, displaying a variety of ‘locally-favored structures’ such as
icosahedra that are denser than crystalline packing but cannot tessellate
space.Patrick Royall  _et al._ (2008); Griffiths, Turci, and Royall (2017) The
absence of such structures in 2D is the origin of locally-crystalline ramified
clusters in our system; it is presumably also related to why crystallization
is ‘easier’ in 2D.

In 3D, gelation in colloid-polymer mixtures is triggered by a coarsening
spinodal gas-liquid phase separation texture that kinetically arrests. In our
system, the theoretical fluid spinodal (dashed line in Fig. 2) occurs far
above the experimental gel boundary, Fig. 4(a). However, the spinodal is where
density fluctuations diverge. At this point, the the mean field-van der Waals
approximation implied by Eq. (2) that underpins the free volume theory is
questionable since the two terms in Eq. (2) are likely to be of the same order
of magnitude. We therefore cannot rule out that the onset of ramified
structures is still co-incident with an underlying spinodal-like thermodynamic
transition in the fluid.333We note though that the free volume approach _can_
be justified for the F-X binodal calculation since there the effect of free
volume overlaps in the fluid phase is still small, and density fluctuations
are unimportant in the near-close-packed hexagonally-ordered crystal; indeed
this motivates the approximate approach described in Appendix A which already
yields a semi-quantitatively accurate prediction for the fluid binodal.

Cerdà et al.Cerdà  _et al._ (2004) found ramified, locally-crystalline
clusters at $U_{0}=7\,k_{\mathrm{B}}T$ at $\eta_{\rm c}=0.2$, which is
consistent with our data, Fig. 4(a), where $\phi_{\mathrm{p}}=0.25$
corresponds to a contact attraction of $6.25k_{\mathrm{B}}T$. Their CSD data
collapsed according to the ansatz

$c_{n}(t)=M_{1}[s(t)]^{-2}f\left[n/s(t)\right],$ (20)

where $M_{k}=\sum_{n}n^{k}c_{n}$ is the $k$-th moment of the CSD, and
$s=M_{2}/M_{1}$ is a measure of the (time-dependent) average cluster size. We
attempt this scaling for our data, Fig. 5(b). The resulting $f$ for the
earliest time, $t=$6\text{\,}\mathrm{min}$$ ($\approx 36\times$ the time for a
single particle to diffuse its own diameter), does not show the same peaked
behavior as data from all other times. We exclude these data from further
consideration. Equation (20) collapses data from
$$10\text{\,}\mathrm{min}$\lesssim t\lesssim$100\text{\,}\mathrm{min}$$, into
a universal $f$ that is peaked, Fig. 5(b) [blue]. Data from
$t\gtrsim$100\text{\,}\mathrm{min}$$, Fig. 5(b) [red], are increasingly noisy,
and show systematic deviations from a single universal curve, especially at
$n/s\gtrsim 2$.

For DLCA, the function $f$ in Eq. (21) takes the formvan Dongen and Ernst
(1985)

$f(x)=Ax^{-\lambda}e^{-ax},$ (21)

where $\lambda$ is the homogeneity exponent, defined such that if
$dc_{n}/dt={\textstyle\frac{1}{2}}\sum_{i+j=n}A(i,j)c_{j}c_{i}-\sum_{j=1}^{\infty}A(n,j)c_{n}c_{j}$
(a Smolouchowski equation), then $A(mi,mj)=m^{\lambda}A(i,j)$. We fitted this
form to the red data in Fig. Fig. 5(b)to obtain $\lambda\approx-0.76\pm 0.05$
(and $A=2.65\pm 0.2\times 10^{6}$, $a=2.58\pm 0.05$).444Our fitting was done
using the non-linear fitting function of Mathematica, which returns fitted
parameter values and their estimated errors.

Kinetic scaling arguments by KolbKolb (1984) predict that the mean aggregation
number in $d$-dimensional DLCA grows with time as $\bar{n}(t)\sim t^{z}$,
where the dynamic critical exponent $z=d_{\mathrm{f}}/(d_{\mathrm{f}}-(d-3))$,
and $d_{\mathrm{f}}$ is the fractal dimension of the ramified clusters.
Importantly, it can be shown that

$z=\frac{1}{1-\lambda},$ (22)

so that $\lambda=(d-3)/d_{\mathrm{f}}$. For two-dimensional DLCA,Meakin
(1983); Kolb, Botet, and Jullien (1983) $d_{\mathrm{f}}\approx 1.4$ implies
$\lambda\approx-0.71$, which is close to our fitted value.

The above assumes that the cluster mobility is inversely proportional to the
cluster radius (i.e. $\alpha=-d_{\mathrm{f}}^{-1}$ in Eq. (1) in Kolb) but
there are reasons to believe that the far-field hydrodynamic interactions may
be screened for wall-bound clusters (see Appendix B) which would instead make
the cluster mobility ultimately inversely proportional to the aggregation
number ($\alpha=-1$ in Kolb) yielding $z={1}/{2}$ and $\lambda=-1$. The best
fit of Eq. (20) with the constraint $\lambda=-1$ is almost visually
indistinguishable from the unconstrained best-fit $\lambda=-0.76$ on the scale
of Fig. 5(b). Our data therefore cannot distinguish between these two models.

Our best-fit value $\lambda=-0.76$ differs from Cerdà et al.Cerdà  _et al._
(2004) who find $\lambda=-0.35$, but this exponent is sensitive to details
such as the system concentration. More importantly, we should seek internal
consistency in the form of Eq. (22). For $\lambda=-0.76$, we expect $z=0.568$.
This dynamical exponent gives a reasonable account of our intermediate-time
$\bar{n}(t)$ data, Fig. 5(c). On the other hand, it is clear that our data
will not be able to decide between this exponent and the value of 0.5 expected
with near-wall hydrodynamic screening.

Systematic deviations from a pure power-law behavior occurs at the end of the
intermediate time window and beyond, Fig. 5(c). As time went on, particles
increasingly adhered to the capillary surface, especially at higher
$\phi_{\mathrm{p}}$. This is not surprising, considering the small particle-
wall gaps inferred from diffusivity, Fig. 3 (inset). The probability of
adesion increases with cluster size $n$. An adhered cluster can no longer
diffuse translationally, and has, at best, restricted rotational diffusivity.
Such adhesion will cause deviations from either of the predicted modes of
dynamical scaling.

##  VI Summary and concluding remarks

We have studied experimentally a layer of colloids at the bottom of a glass
capillary in the presence of smaller polymers. The combination of
gravitational sedimentation and the depletion attraction induced between the
spheres and the wall tightly confined the spheres to a 2D layer, in which we
have been able to deduce the sphere-wall gap from fitting measured
diffusivities to a hydrodynamic theory.Faxén (1921); HB73; Chio and Tse (2020)

The depletion attraction between the spheres induced by the polymers drives
them into F-X coexistence at colloid concentrations $\eta_{\mathrm{c}}$ very
much lower than that needed for 2D crystallization ($\eta_{\mathrm{c}}=0.67$)
in the absence of polymers. A free-volume theoryPoon and Warren (1994) adapted
to our semi-grand canonical system gives a good account of the F-X coexistence
boundary.

At high polymer concentration, the depletion attraction drives the formation
of ramified clusters that are locally crystalline, confirming a previous
simulation.Cerdà  _et al._ (2004) However, the cluster size distribution and
cluster growth dynamics show quantitative differences with these simulations,
which we speculate as due to particle adhesion and/or the screening
hydrodynamic interactions by the wall. Future simulations may explore the
validity of these proposals.

The 500-word 1954 letter by Asakura and Oosawa published in this
journalCipelletti, Martens, and Ramos (2020) together with Vrij’s later
detailed treatment based on the AO pictureVrij (1976) marked the start of
modern research into depletion-driven phenomena. A first theoretical account
of the phase behavior in colloid-polymer mixtures integrates out the polymeric
degrees of freedom and uses an AO potential for the inter-particle
interaction.Gast, Hall, and Russel (1983) Later, a ‘primitive model’ that
takes explicit account of the polymers’ centre-of-mass degrees of freedom
successively predicts polymer partitioning in coexisting phases.Lekkerkerker
_et al._ (1992) Its applicability was confirmed by bulk experiments.Ilett  _et
al._ (1995) Interestingly, this primitive model uses a semi-grand canonical
ensembleLekkerkerker (1990) as a calculational device.

It is gratifying that seven decades from ‘AO’ and three decades from the semi-
grand canonical model, we are able to perform experiments in a well-
characterized ensemble of this kind, fit our dynamical (diffusion) data by
appealing to an AO form of the interaction between particles and wall, Fig. 3,
and account for the equilibrium phase behavior using a modified version of the
original primitive model, Fig. 4(a).

###### Acknowledgements.

SEG was funded by an EPSRC studentship. NK, TV and WCKP were funded by ERC
Advanced Grant ERC-2013-AdG 340877-PHYSAP. The data that support the findings
of this study are openly available in [to be inserted].

##  Appendix A Approximate estimate of fluid binodal

In the presence of non-adsorbing polymers the cohesive free energy of the
hexagonally-ordered crystal can be approximated by calculating the energy
required to break the AO ‘bonds’, as in a solid-state physics problem.Kittel
(2004) By matching this to the colloid chemical potential in the fluid phase,
one can estimate of the location of the fluid binodal. Let the crystal co-
ordination number be $z$. Then

$F^{\mathrm{(crystal)}}\approx
F_{0}^{\mathrm{(crystal)}}-\frac{zN_{\mathrm{c}}}{2}\times\rho_{\mathrm{p}}k_{\mathrm{B}}T\,V_{2}(2R)$
(23)

where the second factor in the second term is the AO bond free energy, and for
simplicity we use the contact value of $V_{2}(r)$. Equating the resulting
colloid chemical potential to that in the two-dimensional fluid derived from
Eq. (3) gives

$\ln\left(\frac{\eta_{\mathrm{c}}}{1-\eta_{\mathrm{c}}}\right)+\frac{\eta_{\mathrm{c}}(3-2\eta_{\mathrm{c}})}{(1-\eta_{\mathrm{c}})^{2}}=\frac{\mu_{0}}{k_{\mathrm{B}}T}-\frac{z}{2}\times\rho_{\mathrm{p}}V_{2}(2R)$
(24)

where $\mu_{0}$ is the chemical potential in the unperturbed crystal, which in
the spirit of the approach we shall suppose constant.

For simplicity we neglect the excluded volume overlaps in the fluid (one can
show that they are small) and tacitly omit the $V_{1}$ terms which contribute
only a common constant to the chemical potentials. Since at coexistence the
chemical potentials of the fluid and crystal are the same, $\mu_{0}$ can be
obtained from the known fluid coexistence composition in the absence of added
polymer, viz. Eq. (24) should be verified by $\eta_{\mathrm{c}}\approx 0.670$
at $\rho_{\mathrm{p}}=0$. Solving Eq. (24) (with $z=6$) for
$\eta_{\mathrm{c}}$ as a function of $\rho_{\mathrm{p}}$ then provides an
estimate of the fluid binodal, shown for the present system as the chained
line in Fig. 2.

##  Appendix B Hydrodynamic interactions

We present heuristic arguments that in a wall-bound cluster containing $N$
particles, the wall ‘screens’ the hydrodynamic interactions such that the
cluster mobility $\sim 1/N$, at least in the scaling limit. We start with the
familiar result that a point force ${\mathbf{f}}$ in an unbounded fluid
generates a velocity field ${\mathbf{v}}$ at a distance ${\mathbf{r}}$ with
(Oseen tensor)Doi86

${\mathbf{v}}=\frac{1}{8\pi\mu
r}({\mathbf{I}}+\hat{{\mathbf{r}}}\,\hat{{\mathbf{r}}})\cdot{\mathbf{f}}\,.$
(25)

In this ${\mathbf{I}}$ is the unit tensor, $\hat{{\mathbf{r}}}={\mathbf{r}}/r$
with $r=|{\mathbf{r}}|$, and $\mu$ is the fluid viscosity. Similarly, Blake
and Chwang showed that a point force at a height $h$ above the wall generates
a flow field which behaves to leading-order in far-field asBlake74

${\mathbf{v}}=\frac{12hz}{8\pi\mu
r^{3}}\>\hat{{\mathbf{r}}}\,\hat{{\mathbf{r}}}\cdot{\mathbf{f}}$ (26)

Here, the no-slip boundary condition is coincident with the $z=0$ plane, and
${\mathbf{r}}=(x,y)$ is now the in-plane distance between the point $(x,y,z)$
where the velocity is measured and the point $(0,0,h)$ where the force is
applied. Crucially, according to Eq. (26), for $z\sim h$ the far-field decays
as $1/r^{3}$ rather than $1/r$ as in Eq. (25) for an unbounded fluid.not

At this point we recall that the mobility of a fractal cluster is essentially
determined by the behavior of $\langle{1/r_{ij}^{\>n}}\rangle$ where the
average is taken over all pairs of particles in the cluster,Wiltzius1987 and
according to the above we should take $n=1$ for freely-suspended clusters or
$n=3$ for wall-bound clusters. In terms of the pair distribution function
$g(r)\sim r^{\,d_{\mathrm{f}}-d}$,

$\langle{1/r_{ij}^{\>n}}\rangle\sim\frac{{\textstyle\int}\,r^{-n}\,g(r)\,r^{\,d-1}\,{\mathrm{d}}r}{{\textstyle\int}\,g(r)\,r^{\,d-1}\,{\mathrm{d}}r}\sim
R^{-n}$ (27)

where $R\sim N^{1/d_{\mathrm{f}}}$ is the cut-off in $g(r)$. But this _only_
holds when the integral in the numerator is dominated by this cut-off, which
requires $d_{\mathrm{f}}>n$. Plainly this is the case for freely-suspended
clusters (unless they happen to be fractal dust with $d_{\mathrm{f}}<1$), and
so one expects that the cluster mobility $\sim N^{-1/d_{\mathrm{f}}}$. This
scaling behavior has been widely confirmed,VanSaarloos1987; Warren1994;
Sorensen2011 and corresponds to the fact that the flow field is screened from
the interior of the cluster. For wall-bound clusters though, it is _not_ the
case that the integral in the numerator in Eq. (27) is dominated by the upper
limit since that would require $d_{\mathrm{f}}>3$ which is impossible in $d=2$
dimensions. Hence one concludes (perhaps a little tentatively!) that
hydrodynamic interactions should be negligible in far-field for wall-bound
clusters, or in other words the flow field, already screened by the wall, is
not further significantly reduced in the interior of the cluster. It follows
that the frictional drag should be extensive in the number of particles, and
the cluster mobility should scale as $\sim 1/N$ as claimed.

This conclusion obviously demands numerical verification, using Stokesian
dynamics or similar methods.Cichocki2000; Swan2007; Gauger2009 For now though,
we close with a couple of remarks higlighting the subtleties of this
hydrodynamic problem. First, it is interesting to note that the mobilities of
individual particles are in a sense _decoupled_ from the hydrodynamic
interactions, since the former are sensitive to the gap between the particles
and the wall (see section IV.1), which can be made arbitrarily small, whereas
the latter are essentially controlled by the heights of the particle centres
above the wall which are limited by the particle radii. Second, whilst the
presence of the wall couples the rotational and translational modes, in a
cluster the rotational modes are partly supressed if the particles are
mutually hindered from all rolling in the same direction (essentially as a
non-trivial consequence of the near-field hydrodynamics). So, it seems
doubtful that the individual particle friction coefficients are simply
additive, but this may not necessarily change the extensivity of the overall
cluster drag coefficient.

## References

  * Asakura and Oosawa (1954) S. Asakura and F. Oosawa, “On interaction between two bodies immersed in a solution of macromolecules,” J. Chem. Phys. 22, 1255–1256 (1954).
  * Vrij (1976) A. Vrij, “Polymers at interfaces and the interactions in colloidal dispersions,” Pure Appl. Chem. 48, 471–483 (1976).
  * Lekkerkerker and Tuinier (2011) H. N. W. Lekkerkerker and R. Tuinier,  _Colloids and the Depletion Interaction_, Lecture Notes in Physics (Springer Netherlands, 2011).
  * Lekkerkerker  _et al._ (1992) H. N. W. Lekkerkerker, W. C. K. Poon, P. N. Pusey, A. Stroobants, and P. B. Warren, “Phase-behavior of colloid plus polymer mixtures,” Europhys. Lett. 20, 559–564 (1992).
  * Ilett  _et al._ (1995) S. M. Ilett, A. Orrock, W. C. K. Poon, and P. N. Pusey, “Phase-behavior of a model colloid-polymer mixture,” Phys. Rev. E 51, 1344–1352 (1995).
  * Poon (2002) W. C. K. Poon, “The physics of a model colloid-polymer mixture,” J. Phys. Condens. Matter 14, R859–R880 (2002).
  * Evans, Poon, and Cates (1997) R. M. L. Evans, W. C. K. Poon, and M. E. Cates, “Role of metastable states in phase ordering dynamics,” Europhys. Lett. (EPL) 38, 595–600 (1997).
  * Pham  _et al._ (2002) K. N. Pham, A. M. Puertas, J. Bergenholtz, S. U. Egelhaaf, A. Moussaïd, P. N. Pusey, A. B. Schofield, M. E. Cates, H. Fuchs, and W. C. K. Poon, “Multiple glassy states in a simple model system,” Science 296, 104–106 (2002).
  * Zaccarelli (2007) E. Zaccarelli, “Colloidal gels: equilibrium and non-equilibrium routes,” J. Phys.: Condens. Matter 19, 323101 (2007).
  * Note (1) Note that if the colloids are too polydisperse to crystallise, this scenario still holds, except that, now, G-L coexistence is no longer metastable, but is the sole equilibrium thermodynamic phase transition in the system.
  * Poon, Pirie, and Pusey (1995) W. C. K. Poon, A. D. Pirie, and P. N. Pusey, “Gelation in colloid-polymer mixtures,” Faraday Discuss. 101, 65–76 (1995).
  * Harich  _et al._ (2016) R. Harich, T. W. Blythe, M. Hermes, E. Zaccarelli, A. J. Sederman, L. F. Gladden, and W. C. K. Poon, “Gravitational collapse of depletion-induced colloidal gels,” Soft Matter 12, 4300–4308 (2016).
  * de Graaf  _et al._ (2019) J. de Graaf, W. C. K. Poon, M. J. Haughey, and M. Hermes, “Hydrodynamics strongly affect the dynamics of colloidal gelation but not gel structure,” Soft Matter 15, 10–16 (2019).
  * Cipelletti, Martens, and Ramos (2020) L. Cipelletti, K. Martens, and L. Ramos, “Microscopic precursors of failure in soft matter,” Soft Matter 16, 82–93 (2020).
  * Poon and Warren (1994) W. C. K. Poon and P. B. Warren, “Phase behaviour of hard-sphere mixtures,” Europhysics Letters (EPL) 28, 513–518 (1994).
  * Dinsmore  _et al._ (1997) A. D. Dinsmore, P. B. Warren, W. C. K. Poon, and A. G. Yodh, “Fluid-solid transitions on walls in binary hard-sphere mixtures,” Europhysics Letters (EPL) 40, 337–342 (1997).
  * Savage, Pei, and Dinsmore (2012) J. R. Savage, L. Pei, and A. D. Dinsmore, “Experimental studies of two-step nucleation during two-dimensional crystallization of colloidal particles with short-range attraction,” in  _Kinetics and Thermodynamics of Multistep Nucleation and Self-Assembly in Nanoscale Materials_, edited by G. Nicolis and D. Maes (John Wiley & Sons, Ltd, 2012) Chap. 5, pp. 111–135.
  * Hobbie (1998) E. K. Hobbie, “Metastability and depletion-driven aggregation,” Phys. Rev. Lett. 81, 3996–3999 (1998).
  * Lekkerkerker (1990) H. Lekkerkerker, “Osmotic equilibrium treatment of the phase separation in colloidal dispersions containing non-adsorbing polymer molecules,” Colloids Surf. 51, 419–426 (1990).
  * ten Wolde and Frenkel (1997) P. R. ten Wolde and D. Frenkel, “Enhancement of protein crystal nucleation by critical density fluctuations,” Science 277, 1975–1978 (1997).
  * Faxén (1921) H. Faxén,  _dissertation_ (Uppsala Univ., 1921).
  * Chio and Tse (2020) C. C. Chio and Y.-L. S. Tse, “Hindered diffusion near fluid-solid interfaces: Comparison of molecular dynamics to continuum hydrodynamics,” Langmuir 36, 9412–9423 (2020).
  * Cerdà  _et al._ (2004) J. J. Cerdà, T. Sintes, C. M. Sorensen, and A. Chakrabarti, “Kinetics of phase transformations in depletion-driven colloids,” Phys. Rev. E 70, 011405 (2004).
  * Griffiths, Turci, and Royall (2017) S. Griffiths, F. Turci, and C. P. Royall, “Local structure of percolating gels at very low volume fractions,” J. Chem. Phys. 146, 014905 (2017), https://doi.org/10.1063/1.4973351 .
  * Griffiths (2021) S. E. Griffiths,  _The Effect of Attractive Forces on Active-Passive Interactions_ , Ph.D. thesis, The University of Edinburgh, Edinburgh (2021).
  * Schwarz-Linek  _et al._ (2016) J. Schwarz-Linek, J. Arlt, A. Jepson, A. Dawson, T. Vissers, D. Miroli, T. Pilizota, V. A. Martinez, and W. C. Poon, “ _Escherichia coli_ as a model active colloid: A practical introduction,” Colloids Surf. B 137, 2–16 (2016).
  * Wang and Yu (1988) L. Wang and H. Yu, “Chain conformation of linear polyelectrolyte in salt solutions: sodium poly(styrenesulfonate) in potassium chloride and sodium chloride,” Macromol. 21, 3498–3501 (1988).
  * Rubinstein and Colby (2003) M. Rubinstein and R. Colby,  _Polymer Physics_ (Oxford University Press, 2003).
  * Schwarz-Linek  _et al._ (2012) J. Schwarz-Linek, C. Valeriani, A. Cacciuto, M. E. Cates, D. Marenduzzo, A. N. Morozov, and W. C. K. Poon, “Phase separation and rotor self-assembly in active particle suspensions,” Proc. Natl. Acad. Sci. (USA) 109, 4052–4057 (2012).
  * Schneider, Rasband, and Eliceiri (2012) C. A. Schneider, W. S. Rasband, and K. W. Eliceiri, “NIH Image to ImageJ: 25 years of image analysis,” Nature Methods 9, 671–675 (2012).
  * Zeng and Oxtoby (1990) X. C. Zeng and D. W. Oxtoby, “Applications of modified weighted density functional theory: Freezing of simple liquids,” J. Chem. Phys. 93, 2692–2700 (1990).
  * Rosenfeld (1990) Y. Rosenfeld, “Free-energy model for the inhomogeneous hard-sphere fluid in $D$ dimensions: Structure factors for the hard-disk ($D=2$) mixtures in simple explicit form,” Phys. Rev. A 42, 5978–5989 (1990).
  * Hall (1972) K. R. Hall, “Another hard-sphere equation of state,” J. Chem. Phys. 57, 2252–2254 (1972).
  * Reiss, Frisch, and Lebowitz (1959) H. Reiss, H. L. Frisch, and J. L. Lebowitz, “Statistical mechanics of rigid spheres,” J. Chem. Phys. 31, 369–380 (1959).
  * Lebowitz, Helfand, and Praestgaard (1965) J. L. Lebowitz, E. Helfand, and E. Praestgaard, “Scaled particle theory of fluid mixtures,” J. Chem. Phys. 43, 774–779 (1965).
  * Stones, Dullens, and Aarts (2018) A. E. Stones, R. P. A. Dullens, and D. G. A. L. Aarts, “Communication: Contact values of pair distribution functions in colloidal hard disks by test-particle insertion,” J. Chem. Phys. 148, 241102 (2018).
  * Barker and Henderson (1976) J. A. Barker and D. Henderson, “What is "liquid"? Understanding the states of matter,” Rev. Mod. Phys. 48, 587–671 (1976).
  * Weaver (1926) W. Weaver, “The duration of the transient state in the settling of small particles,” Phys. Rev. 27, 499–503 (1926).
  * Poon (2015) W. C. K. Poon, “Colloidal suspensions,” in  _Oxford Handbook of Soft Condensed Matter_ , edited by E. Terentjev and D. A. Weitz (OUP, Oxford, 2015) pp. 1–49.
  * Rentsch  _et al._ (2006) S. Rentsch, R. Pericet-Camara, G. Papastavrou, and M. Borkovec, “Probing the validity of the Derjaguin approximation for heterogeneous colloidal particles,” Phys. Chem. Chem. Phys. 8, 2531–2538 (2006).
  * Valmacco  _et al._ (2016) V. Valmacco, M. Elzbieciak-Wodka, C. Besnard, P. Maroni, G. Trefalt, and M. Borkovec, “Dispersion forces acting between silica particles across water: influence of nanoscale roughness,” Nanoscale Horiz. 1, 325–330 (2016).
  * Noro and Frenkel (2000) M. G. Noro and D. Frenkel, “Extended corresponding-states behavior for particles with variable range attractions,” The Journal of Chemical Physics 113, 2941–2944 (2000), https://doi.org/10.1063/1.1288684 .
  * Note (2) We calculated $b_{2}$ by numerical integration using an approximation to the AO potential valid for $\xi\lesssim 0.1$ taken from Ref. Bergenholtz2003.
  * Poon (1997) W. C. K. Poon, “Crystallization of globular proteins,” Physical Review E 55, 3762–3764 (1997).
  * Patrick Royall  _et al._ (2008) C. Patrick Royall, S. R. Williams, T. Ohtsuka, and H. Tanaka, “Direct observation of a local structural mechanism for dynamic arrest,” Nature Materials 7, 556–561 (2008).
  * van Dongen and Ernst (1985) P. G. J. van Dongen and M. H. Ernst, “Dynamic scaling in the kinetics of clustering,” Phys. Rev. Lett. 54, 1396–1399 (1985).
  * Note (3) Our fitting was done using the non-linear fitting function of Mathematica, which returns fitted parameter values and their estimated errors.
  * Kolb (1984) M. Kolb, “Unified description of static and dynamic scaling for kinetic cluster formation,” Phys. Rev. Lett. 53, 1653–1656 (1984).
  * Meakin (1983) P. Meakin, “Formation of fractal clusters and networks by irreversible diffusion-limited aggregation,” Phys. Rev. Lett. 51, 1119–1122 (1983).
  * Kolb, Botet, and Jullien (1983) M. Kolb, R. Botet, and R. Jullien, “Scaling of kinetically growing clusters,” Phys. Rev. Lett. 51, 1123–1126 (1983).
  * Gast, Hall, and Russel (1983) A. P. Gast, C. K. Hall, and W. B. Russel, “Polymer-induced phase separations in nonaqueous colloidal suspensions,” J. Colloid Interface Sci. 96, 251–267 (1983).
  * Kittel (2004) C. Kittel,  _Introduction to Solid State Physics_ (Wiley, Hoboken, NJ, 2004).
  * (53) Eq. (26) is derived from the Blake tensor in the limit $r\to\infty$ at finite $z$, and differs from the better-known asymptotics for the flow field _distant_ from the wall (i.e. $r,z\to\infty$ with $z/r$ held constant) for which $|{\mathbf{v}}|\sim 1/r^{2}$.
  * Bergenholtz, Poon, and Fuchs (2003) J. Bergenholtz, W. C. K. Poon, and M. Fuchs, “Gelation in model colloid-polymer mixtures,” Langmuir 19, 4493–4503 (2003).


# Raindrops on Windshield: Dataset and Lightweight Gradient-Based Detection
Algorithm*

Vera Soboleva1,2,3 and Oleg Shipitko1,2 *This work was supported by Evocargo
LLC1Vera Soboleva and Oleg Shipitko are with Evocargo LLC, Moscow, Russia
<EMAIL_ADDRESS><EMAIL_ADDRESS>Soboleva and Oleg
Shipitko are with the Vision Systems Laboratory, Institute for Information
Transmission Problems, Moscow, Russia3Vera Soboleva is with the Moscow
Institute of Physics and Technology (National Research University),
Dolgoprudny, Russia

###### Abstract

Autonomous vehicles use cameras as one of the primary sources of information
about the environment. Adverse weather conditions such as raindrops, snow,
mud, and others, can lead to various image artifacts. Such artifacts
significantly degrade the quality and reliability of the obtained visual data
and can lead to accidents if they are not detected in time. This paper
presents ongoing work on a new dataset for training and assessing vision
algorithms’ performance for different tasks of image artifacts detection on
either camera lens or windshield. At the moment, we present a publicly
available set of images containing $8190$ images, of which $3390$ contain
raindrops. Images are annotated with the binary mask representing areas with
raindrops. We demonstrate the applicability of the dataset in the problems of
raindrops presence detection and raindrop region segmentation. To augment the
data, we also propose an algorithm for data augmentation which allows the
generation of synthetic raindrops on images. Apart from the dataset, we
present a novel gradient-based algorithm for raindrop presence detection in a
video sequence. The experimental evaluation proves that the algorithm reliably
detects raindrops. Moreover, compared with the state-of-the-art cross-
correlation-based algorithm [1], the proposed algorithm showed a higher
quality of raindrop presence detection and image processing speed, making it
applicable for the self-check procedure of real autonomous systems. The
dataset is available at $github.com/EvoCargo/RaindropsOnWindshield$.

##  I Introduction

As a part of the Intelligent Transportation System, autonomous vehicles use
video cameras to perceive the environment. They are used to detect objects on
the road and navigate. However, adverse weather conditions can negatively
affect the performance of the system. Raindrops or snow could block the
camera’s view and make the further movement of autonomous vehicles unsafe for
both passengers and pedestrians. In this regard, it is essential that
autonomous transportation systems can run self-check procedures and evaluate
whether they can still perform their intended function.

The main difficulty in image artifacts detection is that drops and other
contaminants can have all kinds of shapes, shades, and structures, which
significantly complicates their parametric description. Besides, artifact
detection algorithms should meet high requirements on speed and computational
efficiency. On the one hand, they must work in real-time, and on the other
hand, they are supplemental procedures in the autonomous system functioning
and should not be a significant consumer of computing resources.

(a) original image

(b) mask

(c) original image

(d) mask

Figure 1: Example of dataset samples labeling. White color denotes an area
with an artifact.

In this paper, we present a new publicly available dataset for training and
evaluating algorithms for various problems of detecting raindrops on an image.
The dataset contains sequences of images captured by a camera attached to the
vehicle during its movement. At the moment, it consists of $8190$ images,
$3390$ of which contain raindrops. All images are annotated with the binary
mask representing areas with raindrops. Annotation examples can be seen in
Fig. 1. We demonstrate the applicability of the dataset in the problems of
raindrop presence detection and raindrop segmentation. In addition to the
dataset, we present an artificial raindrops generation algorithm, which allows
the creation of various realistic raindrops to augment the data.

Apart from the dataset, we propose a novel algorithm for detecting raindrops
on a sequence of images. The algorithm is based on two assumptions. Firstly,
as observations show, raindrops mostly remain stable in their position in the
images. Secondly, the area covered by the drop becomes blurred, so there are
no sharp boundaries inside the raindrop image area. Thus, the sequence of
images is analyzed for the presence of the static regions without sharp edges
inside. For this, a gradient map is calculated for each frame, which is then
averaged over the entire sequence. The problem of detecting raindrops was
considered as a problem of binary classification of image sequences. The
performance of the algorithm was evaluated using ROC curve and AUC-ROC metric.
The experimental evaluation showed that the algorithm reliably detects
raindrops and demonstrates better classification quality and higher image
processing speed than the correlation-based algorithm [1]. Thus, a gradient-
based algorithm can be applied as a real-time procedure for autonomous visual
systems self-test.

The proposed algorithm can also solve the raindrop segmentation problem.
However, modern neural network approaches will indeed outperform it. In this
paper, we trained the original U-Net architecture’s neural network [2] to
compare our algorithm with it and demonstrate the dataset’s applicability for
this task. We evaluated the segmentation quality via the Dice coefficient
calculated over raindrops pixels, intersection over union (IoU), and pixel-
wise accuracy. According to the experiment results, U-Net showed high metrics
values on the test samples, which means that the proposed dataset is suitable
for training neural networks.

##  II Related work

Most of the existing weather conditions datasets are not suitable for solving
the problem of vehicle camera lens pollution detection and do not provide the
ability to work with a sequence of images since they contain independent
frames. For instance, Multi-class Weather Image [3] and Image2Weather [4]
datasets are only intended for classifying images by weather but not for
determining the presence of artifacts on the camera lens. On the other hand,
there are visual navigation aimed datasets such as the Oxford RobotCar dataset
[5] or KITTY [6], containing image sequences with various weather conditions.
However, unlike our dataset, they are not suitable for the raindrop
segmentation problem since it requires additional ground-truth mask labeling.

There are many works related to image artifact detection. Some of them use a
sequence of images, while others perform detection on a single image.
Algorithms from [1, 7, 8, 9] are based on the assumption that most artifacts
are temporally stable in their position in the image. Therefore they use a
sequence of frames for detection.

In [1] a pixel-wise normalized cross-correlation between pairs of images is
calculated to analyze a sequence of images for static regions. The idea is
that the corresponding pixels of two consecutive frames belonging to static
objects will show high correlation values. At the same time, the moving parts
of the scene will be weakly correlated. The main disadvantage of this approach
is that, in linear motion, distant objects remain almost stable in their
positions in the image over time and give high correlation values, leading to
false artifacts detections. Also, if the artifacts are very opaque and the
intensity variance inside the image area corresponding to them is close to
zero, the algorithm becomes computationally unstable. Moreover, if the
artifacts, on the contrary, are too transparent, then the intensity inside
them changes significantly, which leads to low correlation values and,
therefore, false-negative detection results.

The authors of [7] also assume that intensity change in raindrop area is less
than in area without it. Therefore, the pixel intensities’ spatial and
temporal derivatives are analyzed to detect artifacts. For this, the optical
flow and the pixel-wise change in intensity over time are calculated. The
major drawback of this method is that the algorithm can not work in real-time
because computing the optical flow of all pixels has high computational
complexity. If only pixel intensities’ temporal derivatives are used, the
algorithm becomes too sensitive to lighting conditions changes, making it
unsuitable for outdoor use. It should be noted that in [7], as well as in this
work, the artificial raindrops generation method is used to expand the
dataset.

In [8] a gradient map averaged over several frames is investigated, to which
nonlinear filtering is applied. The main difference between the algorithm
proposed in this work and [8] is that the algorithm from [8] solves the
problem of detecting tiny drops similar to point noise in the image. Gradient
maps are also used in [9], where the appearance of a drop in the image is
detected by changing the intensities between two sequential frames. Authors
use an image gradient to check if the candidate region’s boundaries have
either disappeared or become more blurred over time. The assumption that the
image area overlapped by the drop becomes blurred and does not have sharp
boundaries inside is also used in our algorithm.

There are also single image artifact detection algorithms. In [10], a raindrop
pattern is artificially created based on a photometric raindrop model.
Detection is then performed by comparing the generated raindrops with the
original image regions. The main disadvantage of this method is that it
detects only drops of a specific shape.

The authors of [11] propose to detect raindrops by detecting blurred regions.
They prove that blurred image areas, which are raindrops, correspond to a
lower Kurtosis coefficient than clear ones. Also, they use Discrete Wavelet
Transform to validate the top of the image (the sky) and Singular Value
Decomposition to validate the bottom (the road). The downside is that the
algorithm can not work with trees at the top of the image and too dense drops
and dirt.

There is another group of approaches that solve the problem with neural
networks. In [12] simplified modifications of VGG, ResNet, and InceptionV3
architectures are used for image classification. The time required to classify
one image (540x300 in size) by the degree and type of pollution is $8-12$ ms
with more than 96% accuracy. The main disadvantage of this approach is the
significant consumption of computing power by the model. The same disadvantage
is possessed by the method proposed in [13], where VGG-16 is used to detect
raindrops. To solve this problem, the authors of [14] demonstrate a modified
VGG-16. They significantly reduce the number of convolutional layers, which
allows them to reduce the number of model parameters by about $250$ times.

In addition, there are several works devoted to removing raindrops from images
[15, 16, 17]. However, using such algorithms in autonomous transportation
systems can be dangerous since it is never fully known what the algorithm will
recreate in place of the former drop.

##  III Dataset Description

The images for the dataset were captured by a camera attached to the vehicle
during its movement. The vehicle’s movement took place in urban areas and
highways, making the dataset ideal for training and assessing vision
algorithms for autonomous vehicle camera lens pollution detection.

The dataset represents sequences of video frames containing $8190$ images, of
which $3390$ contain raindrops. Images were labeled by outlining artifacts
with polygons. Labeling results are stored in JSON format. Besides, binary
masks were generated from this markup, which are also presented in the dataset
for convenience(Fig. 1). White color denotes an artifact area.

Most weather-related artifacts on a windshield or a camera lens have the
property of being almost stable in position over time. Therefore, there is no
need to label each image in the sequence separately. It is sufficient to
manually label the artifact on the image when it first appears and then copy
labeling for all subsequent frames where the artifact is in the same position.
In this work, every $50-100$-th image was manually labeled from the moment the
artifact first appeared.

##  IV Artificial raindrops generation algorithm

Collecting images with a variety of raindrops is a challenging and time-
consuming task. Thus, an image raindrop simulator would greatly simplify the
process and provide data augmentation to efficiently train and evaluate
different raindrop detection algorithms. In this paper, we propose an
algorithm that can generate diverse and realistic artificial drops in images.
The following steps describe the algorithm for generating one drop on one
image $A$ of size $(W,H)$.

  * •

A raindrop is assigned a random radius $R\in[R_{min},R_{max}]$, coordinates of
the drop center $(x,y)\in(0..W,0..H)$ and a shape: 0 – circle , 1 – egg, 2 – a
combination of two Bezier curves (Fig. 2). The egg shape is created by a
combination of a circle and a semi-ellipse.

  * •

An $AlphaMap$ is created – an empty (all zeros) grayscale image of size
$(5R,4R)$. A figure corresponding to the drop shape is drawn in the center of
the image. Depending on the selected in advance brightness (0-255), the
generated raindrop will be either more transparent or opaque. Next, a Gaussian
blur is applied to $AlphaMap$ to make the resulting drop more realistically
blurred.

  * •

Then, an image $P$ of size $(5R,4R)$ is created by cropping a rectangle
corresponding to the drop center from image $A$. Gaussian blur and barrel-
shaped radial distortion (fish-eye effect) are applied to it.

  * •

An image $P^{{}^{\prime}}$ is created by reducing the brightness of the image
$P$ with a factor of $0.3$ (factor of $1.0$ corresponds to the original image,
$0.0$ \- to completely black). After which $P^{{}^{\prime}}$ is inserted at
the appropriate place in the original image $A$ using $AlphaMap$ as a
transparency map. Where the $AlphaMap$ is $255$, the given image
$P^{{}^{\prime}}$ is copied as-is. Where the $AlphaMap$ is 0, the current
value of $A$ is preserved. Intermediate values will mix the two images
together. It is done to make the edges of the generated drop darker, as
natural drops usually have this property.

  * •

Finally, using $AlphaMap$ as a transparency map, the $P$ image is pasted to
the appropriate place on top of figure $P^{{}^{\prime}}$ to the original image
$A$.

Fig. 3 shows an example of egg-shaped drop generation.

Figure 2:  Geometry of artificial drop with different shapes. From left to
right: circle, egg, Bezier curves.

(a) (b) (c) (d)

Figure 3: Example of generation of a drop with egg shape. (a) - $AlphaMap$;
(b) - result of pasting $P^{{}^{\prime}}$ into original image, using
$AlphaMap$ as a transparency mask; (c) - result of pasting $P^{{}^{\prime}}$
into original image, using $AlphaMap$ as a transperency mask, without pasting
$P$ before; (d) - result of of sequential pasting of images $P^{{}^{\prime}}$
and $P$ into the original image using $AlphaMap$ as a transparency mask.

##  V Gradient-based algorithm for raindrops detection

Gradient maps are widely used in image processing to highlight the edges of
objects in the image: the edges in the image have a high gradient magnitude,
while the homogeneous areas have a low gradient.

The area of the image covered by a raindrop is blurred and, as a result, does
not have sharp boundaries inside, which means that the corresponding pixels
will have a low gradient value and will be dark on the gradient map. Due to
the raindrops positions stability in the image, the corresponding areas will
remain dark for at least some sequence of images with raindrops. At the same
time, the areas in which there is movement will have moving boundaries, which
means that the gradient in them will be larger. The algorithm proposed in this
work is based on this observation and uses the Sobel operator to compute
gradient maps.

The input of the algorithm is a sequence of $N$ grayscale images. For each
individual image, the horizontal and vertical gradients are calculated -
$G_{x}^{(n)}(x,y)$ and $G_{y}^{(n)}(x,y)$, respectively, $n=\overline{1,N}$,
where $(x,y)$ are the coordinates of the pixel in the image. They are used to
build a gradient map
$G^{(n)}(x,y)=\sqrt{G_{x}^{(n)2}(x,y)+G_{y}^{(n)2}(x,y)}$. Next, the gradient
is averaged pixel-wise over the entire sequence:

$G^{*}(x,y)=\frac{1}{N}\sum_{n=1}^{N}G^{(n)}(x,y).$

Gaussian filter is applied to the averaged gradient map to reduce salt-and-
pepper noise produced by small regions without edges on an image sequence but
not containing raindrops (see Fig. 4e). For actually detecting an artifact and
obtain the binary segmentation mask, the accumulated gradient map is inversely
binarized (white become black and wise versa) by thresholding with a certain
threshold $T_{b}$. Since the Gaussian filter not only reduces noise but also
reduces the size of the drop itself(see Fig. 4e as the result of applying
Gaussian blur to Fig. 4c), after binarization, dilation is applied, which
returns the raindrop to its previous size.

For the raindrops presence detection, the percent of artifact pixels (white
pixels) is calculated in the resulting binary map. If it exceeds the detection
threshold $T_{d}$, then the algorithm signals raindrops in the sequence of
images. A step-by-step raindrops detection algorithm is demonstrated in Fig.
5.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 4: This figure shows an example of artifact detection for a camera
having an adherent raindrop near the lens center. Two images (a) and (b) are
captured at $10$ frames interval. Calculation of the averaged gradient map
over these $10$ frames leads to the gradient map (c). In this map, the
artifact is visible at low gradient values (dark pixels). Fig. (d) illustrates
the result of averaged gradient map binarization. There is a noticeable noise
present in the binarized image. To filter it, a Gaussian blur is applied to
the gradient map as shown in (e) before binarization (see Fig.(f)). Finally,
we get a binary segmentation mask (g) by applying dilatation to the image (f).
Fig. (h) shows a ground-truth manual segmentation for this sequence. Figure 5:
Raindrop detection algorithm.

##  VI Experimental evaluation

To demonstrate the applicability of the dataset for solving various problems
and show the quality of the proposed algorithm, we consider the problem of
detecting drops in images in two forms. The first one is the binary
classification of image sequences into sequences with drops and sequences
without ones. And the second one is the raindrop region segmentation.

The dataset was divided into a training set ($719$ sequences of $10$ images)
and a test set ($100$ sequences of $10$ images) to select gradient-based
algorithm parameters and evaluate its performance. In both training and test
set, the number of sequences with drops was approximately equal to the number
of sequences without ones.

For all experiments, we used a personal computer with a 6-core Intel Core i7 -
9750H, 2.6 GHz, 16.0 GB RAM.

###  VI-A Algorithm parameters selection

As described in the previous section, the gradient-based algorithm requires
selecting four parameters:

  * •

$d$ the size of the Gaussian filter kernel,

  * •

$T_{b}$ the binarization threshold,

  * •

$m$ the size of the dilation kernel,

  * •

$T_{d}$ the detection threshold.

The parameters used in all further experiments were selected based on the
raindrop segmentation quality and the image processing time. The detection
threshold – the area occupied by raindrops on the image to signal the
successful detection – was fixed as $T_{d}=0.1$. The choice of such a value
for the detection threshold is based on the fact that overlapping only ten
percent of the image can make the autonomous system further movement unsafe.

The area occupied by drops is small relative to the image size, which makes
segmentation classes unbalanced. For all standard metrics such as Intersection
over Union (IoU) or Dice coefficient, the highest values are obtained when the
algorithm signals the total absence of raindrops.

To solve this, we used the Dice coefficient accumulated over all sequences as
follows: if there are $L$ sequences and $I_{i}$ is the intersection of the
ground truth raindrop area and the binary mask found by the algorithm for
sequence number $i$, and $U_{i}$ is the sum of pixels belonging to the
raindrop area in both images, then the metric is computed as follows:

$\text{Accumulated-
Dice}(d,m,T_{b},T_{d})=\dfrac{\sum\limits_{i=1}^{L}I_{i}}{\sum\limits_{i=1}^{L}U_{i}}.$

To select the optimal parameters, we tried to maximize this accumulated Dice
coefficient and minimize image processing time. With this metrics we obtained
the following optimal parameter values: $d=271$, $m=91$, $T_{b}=0.18$.

###  VI-B Image sequence classification

The problem of detecting raindrops on a sequence of images was considered as a
binary classification of image sequences. The quality of the sequence
classification was assessed using ROC curves and the AUC-ROC metric. ROC
curves were generated by varying the free parameter $T_{b}$ – the threshold
for gradient map binarization. In all experiments, $T_{b}$ varied from $0.1$
to $0.9$ with a step of $0.05$.

Fig. 6(a) shows the ROC curve of the gradient raindrop detection algorithm
calculated on sequences of $10$ images, using the $5\times 5$ Sobel operator.
The value of the AUC-ROC metric in this experiment turned out to be $0.83$,
and the processing time for the whole sequence varies from $0.7$ to $0.8$
seconds. Thus, the algorithm provides a high-quality classification of
sequences with a short image processing time, making it potentially applicable
for autonomous visual systems self-check procedure.

We also tested our algorithm on sequences of $100$ images. For this, $100$
sequences of $10$ images were re-sorted into $10$ sequences of $100$ images,
which was possible since together all these images made up a complete video.
As expected, the results showed that an increase in the number of images in
sequences leads to increased classification quality. Although it also
increases the time required for the detection of raindrops by $10$ times. The
ROC-AUC metric in this experiment turned out to be $1.0$. However, since the
sample consisted of only $10$ sequences, this does not mean that the algorithm
works perfectly.

We compared the proposed algorithm with the existing one based on pixel-wise
normalized cross-correlation (NCC) [1]. According to the authors, the NCC’s
main drawback is its poor performance on sequences taken during a straight
motion. In this case, distant objects in the images become almost stationary
relative to the observer, which leads to high correlation values between
subsequent frames and false-positive raindrop detections. Therefore, to
compare our algorithm with NCC, among all sequences of $10$ images, we
selected only those that were captured during the turning maneuver. So, we got
$20$ sequences, among which a half were with raindrops.

The ROC-curves received as a result of testing two algorithms on such
sequences are shown in Fig. 6c. As it could be seen from the plots, both
algorithms demonstrate a high quality of classification. However, the
gradient-based algorithm shows the best metric value: AUC-ROC(Grad) = $0.92$,
AUC-ROC(NCC) = $0.87$. In addition, ROC-curves on a full set of sequences (
Fig. 6a-b) showed that in contrast to NCC, the proposed algorithm also works
well on sequences captured during linear motion. This property gives the
gradient-based algorithm a great advantage since autonomous vehicle movement
is mainly rectilinear.

Also, we compared the proposed algorithm and NCC in image processing time.
Although in the implementation with the box filter [18], NCC, like the
proposed algorithm, has a linear computational complexity of one image
processing, it turned out that the processing time for one sequence of $10$
images for NCC is $1.4-1.5$ seconds, which is almost two times longer than the
processing time for one sequence by the proposed algorithm, which is $0.7-0.8$
seconds.

(a)

(b)

(c)

Figure 6: ROC-curves of the proposed algorithm and the algorithm based on
normalized cross-correlation (NCC) from [1]. ROC-curves plotted over: (a) -
sequences consisting of $10$ images, (b) - sequences consisting of $100$
images, (c) - sequences consisting of $10$ images and taken during the turning
maneuver.

###  VI-C Raindrop region segmentation

To demonstrate the dataset applicability for solving the raindrop segmentation
problem, we trained a neural network of the original U-Net architecture [2].
As a training sample, we took the same set of $7190$ images on which the
gradient-based algorithm parameters were selected. Also, $500$ images with
artificial drops were added to the training set to augment the data. The test
set consisted of the same $1000$ images on which the gradient-based method was
tested in the previous section. All the images in this experiment were resized
to $216\times 216$.

Three metrics were used to assess the quality of segmentation: Dice
coefficient calculated over raindrops pixels, Intersection over Union, and
accuracy. The average values of all three metrics on the test sample turned
out to be very high: IoU = $0.86$, Dice = $0.82$, accuracy = $0.97$. It means
that the dataset suits well for solving the segmentation problem and training
neural networks.

Besides, the proposed gradient-based algorithm was also evaluated on solving
segmentation problem on a sequence of images. In this case, all steps of the
algorithm remain the same until obtaining a binary mask. Metric values for the
gradient in this experiment on sequences of $10$ images resized to $216\times
216$ were as follows: IoU = $0.72$, Dice = $0.5$, accuracy = $0.95$. As
expected, U-Net significantly outperforms the proposed method. On the other
hand, the gradient-based algorithm requires significantly less time to process
one image. The average processing time for one image for U-Net is $0.7$
second, while the same time for the proposed algorithm is $0.025$ seconds.
Thus, the gradient method loses in segmentation quality compared to the U-Net.
However, it can work in real-time, consuming less computational resources,
which is very important for autonomous vehicles. Detailed results of
experiments are presented in table I.

TABLE I: Results of comparing quality and speed of raindrop segmentation for the gradient-based algorithm and U-Net. Dice score calculated over raindrops pixels, Intersection over Union, pixel-wise accuracy and average one $216\times 216$ image processing time were used to evaluate the performances of the algorithm. | IoU | Dice | Accuracy | Time(s)
---|---|---|---|---
Gradient-based algorithm | 0.72 | 0.50 | 0.95 | 0.025
UNet | 0.86 | 0.82 | 0.97 | 0.7

##  VII Conclusion

In this paper, we presented a publicly available dataset for training and
evaluating algorithms for various problems of detecting raindrops on the
camera lens or windshield in front of the camera. The dataset contains
sequences of images captured by a camera attached to the vehicle during its
movement. At the moment, it consists of $8190$ images, 3390 of which contain
raindrops. All images are annotated with the binary mask representing areas
with raindrops.

We demonstrated the applicability of the dataset in raindrop presence
detection and raindrop segmentation problems. In addition to the dataset, we
presented an artificial raindrops generation algorithm, which allows the
creation of various realistic raindrops to augment the data.

At the moment, we are actively working on expanding the dataset. One of the
main directions of dataset development is collecting images with other
weather-related artifacts such as glare, dirt, snow, and others.

Apart from the dataset, we proposed a novel algorithm for detecting raindrops
on a camera lens from a sequence of images. The problem of detecting raindrops
was considered as a problem of binary classification of image sequences. The
performance of the algorithm was evaluated using ROC curves and AUC-ROC
metrics. The experiment results showed that the algorithm reliably detects
raindrops and demonstrates better classification quality and higher image
processing speed than the already existing algorithm [1].

Also, the gradient-based algorithm was tested in the problem of raindrop
segmentation. We trained the original U-Net architecture’s neural network [2]
to compare with our algorithm and demonstrate the dataset applicability for
this task. To evaluate the algorithms’ segmentation quality, we used Dice
score calculated over raindrops pixels, Intersection over Union, and pixel-
wise accuracy. According to the experiment results, U-Net showed higher
metrics values on the test sample than the gradient-based algorithm. However,
the proposed algorithm requires significantly less time to process one image.
Thus, it can work in real-time, consuming less computational resources, which
is essential for an autonomous vehicle.

Further improvement of the proposed algorithm can be aimed at developing
additional verification of the detected raindrops to avoid misinterpretation
by the algorithm of other homogeneous and static regions, such as the sky.
Also, this algorithm can potentially be used to detect other types of
artifacts besides raindrops, such as snow, dirt, scratches on the lens, and
others.

## References

  * [1] Nils Einecke, Harsh Gandhi, and J. Deigmöller.  Detection of camera artifacts from camera images.  17th International IEEE Conference on Intelligent Transportation Systems (ITSC), pages 603–610, 2014.
  * [2] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  CoRR, abs/1505.04597, 2015.
  * [3] Z. Zhang and H. Ma.  Multi-class weather classification on single images.  In 2015 IEEE International Conference on Image Processing (ICIP), pages 4396–4400, 2015.
  * [4] Wei-Ta Chu, Xiang-You Zheng, and Ding-Shiuan Ding.  Camera as weather sensor: Estimating weather information from single images.  Journal of Visual Communication and Image Representation, 46:233–249, 2017.
  * [5] Will Maddern, Geoff Pascoe, Chris Linegar, and Paul Newman.  1 Year, 1000km: The Oxford RobotCar Dataset.  The International Journal of Robotics Research (IJRR), 36(1):3–15, 2017.
  * [6] Andreas Geiger, Philip Lenz, Christoph Stiller, and Raquel Urtasun.  Vision meets robotics: The kitti dataset.  The International Journal of Robotics Research, 32(11):1231–1237, 2013.
  * [7] S. You, R. T. Tan, R. Kawakami, Y. Mukaigawa, and K. Ikeuchi.  Adherent raindrop modeling, detectionand removal in video.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 38(9):1721–1733, 2016.
  * [8] Aurelien CORD and Didier Aubert.  Towards Rain Detection through use of In-Vehicle Multipurpose Cameras.  In IEEE Intelligent Vehicles Symposium (IV’11), page 6p., Germany, June 2011.
  * [9] Fawzi Nashashibi, Raoul de Charette, and Alexandre Lia.  Detection of unfocused raindrops on a windscreen using low level image processing.  In ICARCV, pages 1410–1415. IEEE, 2010.
  * [10] M. Roser and A. Geiger.  Video-based raindrop detection for improved image registration.  In 2009 IEEE 12th International Conference on Computer Vision Workshops, ICCV Workshops, pages 570–577, 2009.
  * [11] Vivek Akkala, Parth Parikh, Basavanna Mahesh, Ajinkya S. Deshmukh, and Swarup Medasani.  Lens adhering contaminant detection using spatio-temporal blur.  2016 International Conference on Signal Processing and Communications (SPCOM), pages 1–5, 2016.
  * [12] Alexey Ivanov and Dmitry Yudin.  Visibility Loss Detection for Video Camera Using Deep Convolutional Neural Networks: Volume 1, pages 434–443.  01 2019.
  * [13] Liang Hu, Li Chen, and Jun Cheng.  Gray spot detection in surveillance video using convolutional neural network.  pages 2806–2810, 05 2018.
  * [14] In-Gyu Bae, Seongbae Bang, and Wonha Kim.  Convolutional neural network for raindrop detection.  pages 102–105, 06 2019.
  * [15] R. Qian, R. T. Tan, W. Yang, J. Su, and J. Liu.  Attentive generative adversarial network for raindrop removal from a single image.  In 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 2482–2491, 2018.
  * [16] Jinwei Gu, Ravi Ramamoorthi, Peter Belhumeur, and Shree Nayar.  Removing image artifacts due to dirty camera lenses and thin occluders.  ACM Trans. Graph., 28, 12 2009.
  * [17] Y. Hirohashi, K. Narioka, M. Suganuma, X. Liu, Y. Tamatsu, and T. Okatani.  Removal of image obstacles for vehicle-mounted surrounding monitoring cameras by real-time video inpainting.  In 2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW), pages 857–866, 2020.
  * [18] M.J. McDonnell and College Park. Computer Science Center University of Maryland.  Box Filtering Techniques.  Report (University of Maryland, College Park. Computer Vision Laboratory). University of Maryland, 1980.