Structure and dynamics of a layer of sedimented particles
Abstract
We investigate experimentally and theoretically thin layers of colloid particles held adjacent to a solid substrate by gravity. Epifluorescence, confocal, and holographic microscopy, combined with Monte Carlo and hydrodynamic simulations, are applied to infer the height distribution function of particles above the surface, and their diffusion coefficient parallel to it. As the particle area fraction is increased, the height distribution becomes bimodal, indicating the formation of a distinct second layer. In our theory we treat the suspension as a series of weakly coupled quasitwodimensional layers in equilibrium with respect to particle exchange. We experimentally, numerically, and theoretically study the changing occupancies of the layers as the area fraction is increased. The decrease of the particle diffusion coefficient with concentration is found to be weakened by the layering. We demonstrate that particle polydispersity strongly affects the properties of the sedimented layer, because of particle size segregation due to gravity.
pacs:
Valid PACS appear hereI Introduction
Being relevant to a wide range of practical scenarios, the behavior of colloid suspensions near solid surfaces has been thoroughly studied over the years. This research effort consists of several bodies of work, for each of which we can give only a few representative references. The first category of papers concerns the disruption of the structural isotropy of a threedimensional (3D) fluid suspension by the surface, e.g., the formation of a layered structure decaying away from the surface under equilibrium Van Winkle and Murray (1988); GonzálezMozuelos et al. (1991) and nonequilibrium ZuritaGotor et al. (2012) conditions. Another category addresses the effect of the anisotropic geometry on particle dynamics near a single planar surface — for isolated particles Happel and Brenner (1983); Perkins and Jones (1992); Cichocki and Jones (1998); Walz and Suresh (1995); Prieve (1999); Sholl et al. (2000); CarbajalTinoco et al. (2007); Blawzdziewicz et al. (2010), particle pairs Happel and Brenner (1983); Dufresne et al. (2000); Cichocki et al. (2007); ZuritaGotor et al. (2007), and a 3D suspension adjacent to a surface Anekal and Bevan (2006); Michailidou et al. (2009); Cichocki et al. (2010); Loppinet et al. (2012); Michailidou et al. (2013).
Regarding quasitwodimensional (quasi2D) layers of particles, most studies have considered the confinement of suspensions between two rigid surfaces. This research addressed structural properties of such confined suspensions CarbajalTinoco et al. (1996); Schmidt and Löwen (1997); Zangi and Rice (2000); Frydel and Rice (2003); Han et al. (2008), and the dynamics of single particles Lin et al. (2000); Dufresne et al. (2001); EkielJeżewska et al. (2008), particle pairs Cui et al. (2004); Bhattacharya et al. (2005a), and concentrated quasi2D suspensions Diamant et al. (2005); Bławzdziewicz and Wajnryb (2012); Baron et al. (2008). Another type of quasi2D suspensions has also been studied, where a particle layer is confined to a fluid interface Lin et al. (1995); Cichocki et al. (2004); Peng et al. (2009); Zhang et al. (2014).
In cases where the surface attracts the particles and the suspension is sufficiently dilute, the system can contain a single layer of surfaceassociated particles in contact with a practically particlefree solvent GonzálezMozuelos et al. (1991). A single layer can also form as a result of gravitational settling of particles toward a horizontal wall. This scenario is studied in the present work.
Sedimented colloidal particles undergo random Brownian displacements, which results in diffusive broadening of the fluctuating particle layer. The width of the particle height distribution above the bottom surface is characterized by the sedimentation length , i.e., the height at which the gravitational energy of a particle equals its thermal energy. The dynamics and height distribution of individual sedimented particles above the bottom surface were studied in Refs. 7; 8; 9 using total internal reflection microscopy. Particle monolayers at higher densities were investigated experimentally for a system in which the sedimentation length is much smaller than the particle diameter Skinner et al. (2010). It was shown that at high area fractions the suspension can assemble into quasi2D colloidal crystals, but formation of a nonuniform vertical microstructure was not observed, because of the small sedimentation length.
Here we are interested in the structure and dynamics of a surfaceassociated layer for which the sedimentation length is comparable to the particle diameter. We focus on the effects of the suspension concentration on the statistical height distribution of particles and their diffusion coefficient. Unlike the quasi2D suspensions confined between two surfaces or adsorbed at a fluid interface (which restricts particle configurations and motions in two directions) in the present system no constraints are imposed on the distance between the particles and the single wall. Thus, at sufficiently high area fractions, particles form a nontrivial stratified microstructure. This microstructure and its effect on particle dynamics are analyzed in our paper.
The article is organized as follows. Section II describes the experimental methods used to prepare the system, image the particles, and analyze the extracted data. In Sec. III we describe the theoretical background and numerical methods used to perform the simulations. In Sec. IV we present the results concerning the equilibrium structure of the quasi2D suspension observed in planes parallel to the bottom surface (the quasi2D radial distribution) and in the direction perpendicular to it (the height distribution). Section V addresses the diffusion of particles parallel to the surface, as affected by the surface proximity. We discuss our findings in Sec. VI.
Ii Experimental methods
ii.1 Quasi2D system of sedimented Brownian spheres
Quasi2D colloidal layers are created by placing a suspension of colloidal silica spheres in a glass sample cell m high. The particles are then allowed to sediment and equilibrate for 30 minutes at a temperature of approximately 24C before measurements start (Fig. 1). We use green fluorescent monodisperse, negatively charged silica particles (Kisker Biotech, PSIG1.5 Lot #GK0090642T) with diameter m, and mass density g/cm. Monolayers of area fraction are prepared by diluting the original suspension with double distilled water (DDW, 18 M), without and with the addition of salt at a concentration . The sample walls are cleaned and slightly charged by plasma etching to avoid particle attachment to the bottom wall of the cell. We observe that the aqueous medium above the colloidal monolayer is free of colloids. Since the particles are floating right above the bottom wall, we can treat the upper wall as a distant boundary.
ii.2 Imaging techniques
Particle position and motion in the – plane, perpendicular to the optical axis, are observed using epifluorescence microscopy (Olympus IX71). Images are captured at a rate of fps by a CMOS camera (Gazelle, Point Grey Research). We use inline holographic microscopy to image the dynamics of particles in three dimensions in dilute samples Kapfenberger et al. (2013). This imaging technique uses a collimated coherent light source (DPSS, Coherent, nm) to illuminate a sample mounted on a microscope. The light scattered from the sample interferes with the light passing through it, to form a hologram in the image plane. We reconstruct the light field passing through the sample by RayleighSommerfeld backpropagation and extract from it the particle location in three dimensions Lee and Grier (2007); Cheong et al. (2010); Kapfenberger et al. (2013). For holographic imaging measurements we use nonfluorescent silica particles with the same diameter (m, Polysciences Inc.). Additional details of the setup and measurement methods can be found elsewhere Kapfenberger et al. (2013).
We use confocal imaging to monitor particle positions in a dense layer in three dimensions. Our spinning disc confocal imaging system (Andor, Revolution XD) includes a Yokogawa (CSUX1) spinning disc, and an Andor (iXon 897) EMCCD camera. An objective lens (Olympus, x60, NA=1.1, water immersion) mounted on a piezoelectric scanner (Physik Instrumente, Pifoc P721.LLQ) is used to scan the sample in the axis, with a step size of 100 nm.
ii.3 Height calibration
A suspended tracer particle is subject to electrostatic and gravitational forces in addition to thermal fluctuations, affecting its height distribution Prieve (1999). The particle potential energy can be described as
(1) 
where is the vertical position of the tracer, is the gravitational acceleration,
(2) 
is the buoyant mass of the tracer ( is the mass density difference between silica and water), is the Debye screening length, and the amplitude depends on and the surface charges of both particle and glass surfaces. The corresponding probability distribution of the particle height is
(3) 
where is the thermal energy and is the normalization constant.
The height distribution of a single particle above the sample’s bottom was obtained from very dilute suspensions, using inline holographic imaging Lee and Grier (2007); Cheong et al. (2010); Kapfenberger et al. (2013) (see Fig. 2). Our holographic measurements provide values of relative particle positions, but not the absolute particle heights with respect to the bottom wall. We thus set the peak position to and focus on the height relative to this reference plane. The exponential decay on the right side of the probabilitydensity peak is governed by a decay length,
(4) 
(the sedimentation length), resulting from the competition between gravity and thermal forces. The exponentialdecay length determined from the holographic measurements agrees well with the calculated sedimentation length (4), without any fitting parameters (see Fig. 2). The electrostatic term of the probabilitydensity, which controls the steep rise of the probability, affects mostly the peak position rather than its shape. Since we shifted the peak position to , the fitting of the entire probabilitydensity using Eqs. (1) and (3) was insensitive to the value of . Reasonable fits were obtained for in the range of nm. Better estimations of and are given in Sec. IV.1.1, using mobility measurements.
The applicability of the holographic imaging is limited to lowdensity suspensions, whereas the confocal imaging can be also used at higher concentrations. On the other hand, confocal height measurements suffer from spherical aberrations due to multiple changes in refractive index in the imaging path. This leads to a systematic error in measuring , which can be eliminated by proper calibration. We calibrate the confocal measurement of the relative vertical particle positions by requiring the exponential decay of the height distribution to agree with the known, and verified, value of .
In Fig. 3(a) we show the particleheight distribution at for two different particle sizes (m). The distributions are shifted so that the highest probability is located at . Scaling the logarithm of the distributions by [inset of Fig. 3(a)] shows that the normalized decay constants for the two particle sizes have approximately the same value, from which we calibrate the confocal microscope’s height measurements. In Fig. 3(b), the height distributions extracted by the two methods (holographic and confocal imaging) are overlaid. This figure emphasizes the higher accuracy of holographic imaging over confocal imaging, especially around , where the increase in distribution should be very steep Prieve (1999); Kapfenberger et al. (2013). The difference between the curves can also be attributed to polydispersity, since the holographic imaging is a singleparticle measurement while the confocal imaging is a multipleparticle measurement, and its corresponding curve represents an average over 40 particles.
Iii Numerical methods
iii.1 The system
iii.1.1 Particles and their interactions
Silica particles are modeled as Brownian hard spheres with or without electrostatic repulsion (depending on the salt concentration), immersed in a fluid of viscosity . The bottom wall is treated as an infinite hard planar surface. Creepingflow conditions and no slip boundary conditions at the particle surfaces and at the wall are assumed.
In a salt solution with , the Debye length is only about 5 nm, and therefore electrostatic interactions are screened out. The particles thus interact only via infinite hardcore particle–particle and particle–wall potentials and the gravity potential , and no other potential forces are involved. The strength of the gravity force is described by the sedimentation length (4).
In addition to the hardcore repulsion, in DDW with no added salt () particles are assumed to also interact via particle–wall and particle–particle Debye–Hückel potentials,
(5) 
and
(6) 
where is the Debye screening length, and are the potential amplitudes, and is the distance between the particle centers. The consideration of Debye–Hückel potentials in the saltfree case is based on our experimental measurement nm. A finite Debye screening length in DDW stems from the presence of residual ions in the solution Behrens and Grier (2001).
iii.1.2 Suspension polydispersity
To determine the effects of the suspension polydispersity on the nearwall microstructure and dynamics, we have performed numerical simulations for a hardsphere (HS) system with a Gaussian distribution of particle diameters,
(7) 
where and are the actual and average particle diameters, and is the standard deviation. All the particles have the same mass density ; hence, particles of different sizes have different buoyant masses and different sedimentation lengths (4). The dimensionless sedimentation length based on the average particle diameter , is defined as
(8) 
where
(9) 
The area fraction based on the average particle diameter is
(10) 
where is the number of particles per unit area. Since the particles are free to move in the direction, the area fraction can exceed 1.
iii.1.3 System parameters
The simulations were carried out for the following system parameters: For the dimensionless sedimentation length (8) we use the value
(11) 
calculated from the particle size and density. Based on the comparison between the calculated and measured values of the equilibrium average of the lateral selfdiffusion coefficient for isolated particles in DDW, we estimate that the Debye length and the amplitude of particle–wall electrostatic repulsion are
(12) 
These values are used for saltfree suspensions at all suspension concentrations. Assuming that the charge densities of the particle and wall surfaces are similar, we take
(13) 
for the interparticlepotential amplitude, as follows from the Derjaguin approximation Israelachvili (1992).
The simulations were performed in the range of area fractions . For polydisperse HS systems the calculations were carried out for , 0.15, 0.20, and 0.25 (we estimate that for the silica particles used in the experiments). For particles interacting via the Debye–Hückel potentials (5) and (6) only monodisperse suspensions were considered.
iii.2 Evaluation of the equilibrium distribution
iii.2.1 Low density limit
For monodisperse suspensions at low particle concentrations, the equilibrium particle distribution is given by the normalized Boltzmann factor (3). To determine the particle distribution for a dilute polydisperse suspension, the particlesizedependent Boltzmann factor for individual particles, , is convoluted with the particlesize distribution (7),
(14) 
For a HS system
(15) 
according to equations (1)–(4), where is the Heaviside step function, and the sedimentation length is particlesize dependent due to the variation of particle mass.
iii.2.2 Monte–Carlo simulations
To determine the equilibrium microstructure of a sedimented suspension at finite particle area fractions, equilibrium Monte–Carlo (MC) simulations were performed for 2Dperiodic arrays of spherical particles in 3D space (with periodicity in the horizontal directions and and the box size ). The particles interact via infinite hardcore repulsion and the pairadditive potential
(16) 
which includes the gravity term and particle–wall and particle–particle screened electrostatic potentials (5) and (6). Here is the particle configuration (with denoting the position of particle ), is the vertical coordinate of particle , and is the relative particle distance.
A purely HS system with was modeled for monodisperse particles and for polydisperse particles with the Gaussian size distribution (7). For systems with nonzero electrostatic repulsion only monodisperse particles were considered.
The initial configuration was prepared by placing particles randomly in a vertical cuboid box with the square base and the height . The size of the 2Dperiodic cell was determined to obtain the required area fraction of the sedimented particle layer. The suspension was allowed to sediment by following the MC randomwalk dynamics in the configurational space Frenkel and Smit (2002) (as described below). After the equilibrium state was reached, suspension properties were obtained by averaging the quantities of interest over at least 200 independent configurations.
Our adaptive simulation procedure was performed by repeating the MC steps defined as follows:

A randomly selected particle is given a small random displacement, , where is chosen from a 3D Gaussian distribution with the standard deviation adaptively adjusted to the current mean gap between particles. This displacement results in the change of the configuration from to .

According to the Metropolis detailed balance condition, the new configuration is accepted with the probability
(17) provided that there is no particle–particle or particle–wall overlap.
To let the system reach an equilibrium state , the MC step (a) and (b) is repeated times. The next independent equilibrium configuration is obtained from the previous configuration by performing MC steps. The particle height distribution and other equilibrium quantities are obtained by averaging over 200 independent configurations .
iii.3 Hydrodynamics and selfdiffusion
iii.3.1 Low density limit
In the absence of a wall, the selfdiffusion coefficient of an isolated solid sphere with diameter is given by the Stokes–Einstein expression
(18) 
The selfdiffusion coefficient of a sphere with diameter at a distance from the wall is smaller by a factor
(19) 
where the normalized mobility coefficient depends on the dimensionless particle position and no other parameters. Relation (19) refers to the lateral component of the selfdiffusion coefficient (parallel to the wall), which was measured in our experiments. However, an analogous expression also holds for the normal component.
For monodisperse particles in the dilutesuspension limit, the effective selfdiffusion coefficient averaged across the suspension layer is obtained by integrating (19) with the Boltzmann distribution (3),
(20) 
For a polydisperse suspension, an additional average over the particlesize distribution (7) is needed,
(21) 
iii.3.2 Computations for larger densities
The effective selfdiffusion coefficient for suspensions at higher concentrations was evaluated using a periodic version Bławzdziewicz and Wajnryb (2008) of the Cartesianrepresentation algorithm for a system of particles in a parallelwall channel Bhattacharya et al. (2005b, c). In our approach, periodic boundary conditions in the lateral directions are incorporated by splitting the flow reflected by the particles into a shortrange nearfield contribution and a longrange asymptotic Hele–Shaw component. The nearfield contribution is summed explicitly over neighboring periodic cells, and the Hele–Shaw component is evaluated using Ewald summation method for a 2D harmonic potential Cichocki and Felderhof (1989); Bławzdziewicz and Wajnryb (2008).
The onewall results were derived from the twowall calculations using an asymptotic procedure based on the observation that in the particlefree part of the channel the velocity field tends to a combination of a plug flow and a shear flow. All other flow components decay exponentially with the distance from the particle layer. The onewall results are obtained by eliminating the shear flow and retaining only the plug flow generated by hydrodynamic forces induced on the particles Sadlej et al. (2009). The calculations were performed for the distance to the upper virtual wall , which is sufficient to obtain highly accurate onewall results.
The selfdiffusion coefficient is determined by averaging the trace of the lateral translational–translational particle mobility, evaluated using Hydromultipole codes based on the above algorithm, with the multipole truncation order Abade et al. (2010). The averaging was performed over equilibrium configurations of particles in a 2Dperiodically replicated simulation cell. Independent equilibrium configurations were constructed using the MC technique described in Sec. III.2.2.
Iv Structure of the quasi2D suspensions
iv.1 Experimental results
A typical image of our quasi2D colloidal suspension is shown in Fig. 1. For each area fraction and salt concentration, the suspension can be characterized by the structure in the – plane (parallel to the cell floor) and the density profile in the direction (perpendicular to the floor). In this section we discuss results of our measurements of the microstructure of a sedimented particle layer.
iv.1.1 Mean particle height at low area fractions
As mentioned in Sec. II.3, our imaging techniques do not yield absolute particle heights. To estimate the mean particle distance from the bottom wall (the mean height) in a dilute suspension layer, we observe particle dynamics in the horizontal directions, and compare measurement results with theoretical calculations of the effect of the wall on the lateral particle diffusion. Using fluorescence imaging, we determine the projection of particle trajectories onto the – plane, , and extract the effective selfdiffusion coefficient,
(22) 
where is the time interval. The positiondependent diffusivity in the – plane of a single particle near a planar wall is given by the following expansion in the particle–wall distance Happel and Brenner (1983); Perkins and Jones (1992),
where is the wall position. The expansion (IV.1.1) is accurate within 5% to 1% as increases from up to , in the range where sedimented particles spend most of the time in a lowdensity suspension under equilibrium conditions. Here is the sedimentation length (4).
From expression (IV.1.1) and extracted according to (22), we can calculate the suspension’s mean distance from the wall (where in (IV.1.1) is replaced by a mean value ). This calculation holds in the limit , where there are no particleparticle interactions. We measure from the particle trajectories, , in extremely low area fraction solution, (in saltfree water), and obtain a mean distance from the wall m, corresponding to a mean gap of 0.3–0.4 m between the particle surface and the wall. We also extract for different salt concentrations by extrapolating (measured at various area fractions) to (see Sec. V.1), obtaining m for and m for . The latter matches the average height extracted from the diffusion of tracers in the extremely low density suspension.
For nm (added salt), the mean height calculated from the Boltzmann distribution (3) is dominated by the exponential decay due to gravity and is practically independent of in the particle–wall potential (1). For , using the particle mass as determined from Eq. (2) with no fitting parameters, we get m, in agreement with the diffusivitybased measurements of . This result confirms that in the addedsalt case we can neglect the electrostatic repulsion from the wall. For the saltfree case, taking nm, we obtain m for in the range 5–15 . These values are consistent with those obtained by fitting the measured height distribution to the theoretical expression (1) (see Sec. II.3).
Since Eq. (IV.1.1) does not include lubrication correction for small particle–wall gaps , it overpredicts for ; however, the accuracy of the approximation is sufficient for the purpose of the present estimates. In our calculations discussed in Secs. III and V.2, highly accurate Hydromultipole results were used instead of the farfield approximation (IV.1.1).
iv.1.2 Radial distribution in the horizontal plane
To verify that no crystalline or hexatic structures are formed at higher values of the area fraction, we evaluate from the experiment the radial distribution function and the full 2D pair distribution in the – plane, for both the saltfree and saltadded suspensions. No dependence on was found. The radial distribution for several values of the area fraction is shown in Fig. 4(a) for the saltfree system and in 4(b) for the saltadded system.
For monodisperse hard spheres the first peak of should correspond to the diameter of the sphere. Our measurements show that the first peak is at m for suspensions without salt and at m for suspensions with . The difference between these two numbers implies that the effective shell around the particles in the saltfree samples is around 4050 nm, which provides an estimation for the screening length in DDW without the addition of salt. This estimate of is consistent with the other two mentioned above.
iv.1.3 Vertical density profile
The height distributions of the silica particles at different area fractions of the sedimented particle layer were acquired using confocal imaging and conventional image analysis Crocker and Grier (1996). These distributions for saltadded suspensions with are plotted in Fig. 5(a) for several values of the areafraction . Since we cannot precisely measure the position of the wall, the distributions are shifted so that their first peak (close to the wall) is located at . These distributions indicate the formation of a second layer of particles for area fractions . The observed center of the second layer is located m above the center of the first layer. The layer separation is thus significantly smaller than the expected separation m (which is similar to the peak separation for the radial distribution). See further discussion in Secs. VI and Appendix A.
To highlight the onset of the formation of the second layer, we look at the subtraction of the height probability distribution of the lowest area fraction from the distribution of all area fractions, [Fig. 6(a)]. Two phenomena are expected when a second layer is formed: (i) negative values at m, corresponding to a reduction in the fraction of particles populating the first layer, (ii) positive and increasing values at m, corresponding to the formation and increasing population of the second layer. The values of at and m are plotted in Fig. 6(b). The two expected phenomena are observed at approximately , indicating the area fraction above which a second layer becomes occupied. At area fractions smaller than 0.3 we still obtain negative values of at m, and positive values at m, however these values are relatively low and can correspond to the broadening of the exponential distribution due to increase in .
iv.1.4 Particlelayer occupation fractions
For the area fractions at which a clear second peak in the particle distribution is seen in Fig. 5(a) (i.e., for and ), we fit the area around each peak to a Gaussian function and define the point of intersection between the two Gaussians as the effective boundary between the two layers. Figure 5(b) shows the two distributions with the Gaussian fits and our definition of that boundary, which turns out to be at a distance of m above the peak of the first layer in both densities.
Using this boundary, we evaluated the occupation fractions of the bottom () and top layer (), where is the area fraction of particles in layer . The results are shown in Fig. 7(a) for a suspension in solution as a function of the total area fraction . As expected, the fraction of particles populating the second layer grows as the total area fraction of the suspension is increased.
An additional independent measurement of the layer occupation fractions is done using epifluorescence microscopy, which enables us to image the different layers separately [see Fig. 1(a)]. The occupation fraction of each layer is determined by counting the number of particles observed therein. The occupation fractions measured using the epifluorescence imaging technique are plotted in Fig. 7(a) along with the results obtained from the confocal microscopy. The two methods yield similar results.
Alternatively, we can represent the layeroccupation results in terms of the area fractions and of the first and second layers [see Fig. 7(b)]. Both and increase as is increased, and seems to saturate at .
iv.2 Numerical simulations
Here we present results of MC simulations of the equilibrium microstructure of a HS suspension in the near wall region. The HS potential corresponds to the system with , for which the electrostatic repulsion is negligible. Since the suspension used in the experiments is polydisperse, we consider both monodisperse and polydisperse systems.
iv.2.1 Nearwall particle distribution
Figure 8 shows MC results for the suspension density profile for a monodisperse suspension and polydisperse suspensions at several area fractions. Similarly to the experimental results, the simulations show that there is a single layer of sedimented particles at low area fractions , and a twolayer microstructure at higher area fractions. (Development of a third layer for is also noticeable in the region .) Suspension polydispersity results in broadening of the peaks of the particle distribution.
A direct comparison between the experimental and simulation results is presented in Fig. 9 for two values of the area fraction . At low area fractions [Fig. 9(a)] the agreement between the experiments and simulations is good. (The standard deviation of the particlesize distribution for which the simulations match the experimental data, , is larger than the estimated standard deviation based on the manufacturer’s specifications; the additional spread of the experimentally observed peak can be attributed to random errors of the particle height evaluation from the confocalmicroscopy images.)
A comparison of the numerical and experimental results at a higher area fraction, as shown in Fig. 9(b) [also see Figs. 5 and 8], reveals that (i) the experimentally observed second maximum of the density distribution develops at lower area fractions than the corresponding maximum in the numerical simulations; (ii) the experimental second peak is narrower, and its position is shifted towards the wall. In contrast, the plots of the excess distribution with respect to the lowdensity limit, shown in Fig. 6(a) and the insets of Fig. 8, indicate that the onset of the formation of the second layer occurs at approximately the same area fraction according to the simulations and experiments. Moreover, the measured and calculated occupation fractions of the layers are similar for all area fractions, as depicted in Fig. 10(a). A possible source of the observed discrepancies between the experimental and numerical results for the particle distribution is described in Appendix A. It also provides a plausible explanation of the fact that the agreement between the experiments and MC simulations for the layer occupation fractions is quite good in spite of the discrepancies for .
iv.2.2 Polydispersity effects
The results in Fig. 10(a) show that the occupation fraction of the first two layers is relatively insensitive to the suspension polydispersity; in contrast, the occupation fraction of the third layer strongly increases with the standard deviation of particle diameter. This increase stems from the presence of smaller (lighter) particles in polydisperse systems: smaller particles tend to migrate into the top layer, as evident from Fig. 10(b). For dilute suspensions, the particlesize segregation results in variation of the slope of with the distance from the wall, as illustrated in Fig. 11. We estimate that this variation causes an approximately 20 % uncertainty of the calibration of the confocal height measurements described in Sec. II.3.
iv.3 A quasi2D model of the equilibrium layered microstructure
Here we present a semiquantitative theoretical model for evaluating the occupation fractions of the particle layers in a sedimented colloidal suspension. Our theory is based on the assumption that the suspension microstructure can be approximated as a collection of weakly coupled quasi2D layers in thermodynamic equilibrium with respect to particle exchange.
The equilibrium condition for layers and is
(24) 
where is the chemical potential of layer , and is its position. In our model, is approximated as the chemical potential of a 2D harddisk fluid of area fraction . All disk diameters are equal to the sphere diameter , which corresponds to a layer of spheres with the same vertical position .
In the low areafraction limit, the chemical potential of a harddisk fluid is
(25) 
where depends only on the temperature . According to the equilibrium condition (24) and equation of state (25), we thus have
(26) 
with the ratio given by the Boltzmann factor
(27) 
where is defined by Eq. (4) and . We assume that the layer separation is independent of .
For finite area fractions, relation (26) is replaced with
(28) 
where the layer occupation ratio depends on the area fraction in the adjacent layers. The factor is determined from the equilibrium condition (24) with the help of the Gibbs–Duhem relation
(29) 
where is the lateral 2D pressure within the layer. Combining (24), (28), and (29) yields
(30) 
where
(31) 
The differential equation (30) is solved for with the boundary condition (27) at . Occupation fractions are then determined by iteration, applying Eq. (28) and the relation
(32) 
We have solved Eq. (30) and determined the occupation fractions using the scaledparticletheory equation of state for hard disks Helfand et al. (1961),
(33) 
The results of our calculations are presented in Fig. 12 for a HS system with the same value of the sedimentation length (11) as in our MC simulations. Based on the separation between the first and second peak of the suspension density profile shown in Fig. 8(a), the calculations were performed for .
The theoretical results in Fig. 12 are compared with the MC simulations of a monodisperse HS suspension with the boundaries between the layers set to and , consistent with the peak positions. The agreement between our simple theory and simulations is quite good. A similar agreement was obtained for other values of the dimensional parameter (results not shown).
The layer boundaries used in Sec. IV.2 and V.2 to compare the MC results with experiments differ from the boundaries used in the above model by approximately 10 %. Due to the observed deviation between the measured and simulated particle distributions (see Fig. 9), it is not possible to define the layer boundaries in a unique, equivalent way for the experimental and simulated systems. Therefore, the layer boundaries and used in Sec. IV.2 and V.2 were chosen based on the comparison between the experimental and numerical results for the occupation fractions and selfdiffusivities in particle layers.
V Particle dynamics
v.1 Experimental results
The short time selfdiffusion coefficient in the – plane, , is determined for different total area fractions of the sedimented particles by extracting the mean square displacement (22) from 2D epifluorescent images of the first and second particle layer. The meansquare displacement is measured over a time interval that is small compared to the structural relaxation time of the suspension, to ensure that the measurements yield the shorttime selfdiffusion coefficient.
The results are shown in Fig. 13 for suspensions with salt concentration and saltfree suspensions with . The selfdiffusion coefficient is expected to decrease as the particle concentration increases; indeed, we observe this decrease for both salt concentrations and in both layers, for . In the case of , corresponding to nm, the particles can get much closer to the cell floor, which in turn results in lower values of the selfdiffusion coefficient compared to suspensions with .
Using a linear fit to the values of for the low area fractions, where there is no observable second layer, we can extrapolate to and extract the selfdiffusivity of a single particle. The extrapolated results agree well with the measurements at very low concentrations , as discussed in Sec. IV.1.1.
From the known occupation fractions and for each we can weigh the contribution of each layer to the total selfdiffusivity, and construct an effective of the whole suspension (Fig. 13). As expected, for the effective selfdiffusion coefficient decreases as is increased in both salt concentrations. For larger we observe a flattening of , which clearly indicates that the second layer becomes dominant in those area fractions. This observation is supported also by the saturation of at [Fig. 7(b)].
v.2 Numerical simulations
The results of our numerical simulations for the shorttime lateral selfdiffusion coefficient in a HS system are presented in Fig. 14 for a monodisperse suspension and for polydisperse suspensions with and . Figure 15 shows the corresponding results for a system of monodisperse hard spheres with particle–wall and particle–particle electrostatic repulsion (5) and (6).
The results depicted in Fig. 14 indicate that for moderately polydisperse suspensions (in the range corresponding to the polydispersity of silica particles used in the experiments), the selfdiffusion coefficient is only moderately dependent on . For larger values of the variance of particle diameters, the normalized selfdiffusion coefficient significantly increases with the degree of the polydispersity, because the mobility is dominated by small particles. This increase is illustrated in Fig. 16 for a suspension in the lowareafraction limit .
The results of our hydrodynamic calculations for a HS suspension and for a suspension with screened electrostatic repulsion are compared with experimental results for suspensions with (Fig. 14) and (Fig. 15). For the system with salt, the measured values are slightly closer to their numerical analogs than in the absence of salt. Summarizing, the experimental and numerical results agree well both for the overall selfdiffusivity and for the selfdiffusivity in individual particle layers.
Vi Discussion
In this paper we have studied in detail the structure and dynamics of quasi2D colloidal suspensions near a wall, comparing experiment and theory. Our central result is a rather sharp formation of a distinct second layer at an area fraction of . This value is much lower than the area fraction required for closepacking or other 2D structural changes such as the formation of hexatic or crystalline order. One important consequence of this result concerns the apparent selfdiffusion of the particles in the suspension and its dependence on particle density. Due to the higher mobility of the particles in the elevated layer, the effective diffusivity is higher and levels off as particle density increases. The experimentally observed behavior could be interpreted incorrectly if one is unaware of the layering (or stratifying) effect.
We find good agreement between experimental and simulation results for the occupation fractions of the first and second layers and for the lateral selfdiffusivity (both for the entire suspension and in the individual layers). However, we also find an unexpected discrepancy in the position and the height of the second peak in the nearwall particle distribution. While the source of this discrepancy is unknown, one possibility, related to optical aberrations, is suggested in Appendix A. On the other hand, the difference between theory and experiment might also be a result of an actual physical effect, such as more complicated electrostatic interactions setting in at higher layer densities.
Another new insight put forth in this study is the significant effect that polydispersity has on the occupation and composition of layers close to the bottom wall, even in the case of a relatively small dispersion of particle sizes. The effect of polydispersity is evident already at low densities, since the smaller and larger particles segregate into the upper and lower layers, respectively. We expect the phenomena described here to be quite general and to be manifested in any such system where the sedimentation length is of the order of the particle diameter. This conclusion is supported by the appearance of the phenomena both in experiments and in Monte–Carlo simulations.
An important outcome of this paper is the construction of a very simple theoretical quasi2D model of the layered microstructure in thermodynamic equilibrium. Such systems have been analyzed earlier using densityfunctional theory Chen and Ma (2006), but our theoretical model is much simpler and easier to apply. We have demonstrated that the model approximates well the experimental and numerical results for the system studied in this work.
We conclude with three open issues. Layering phenomena near a wall are well documented in 3D suspensions as well Van Winkle and Murray (1988); GonzálezMozuelos et al. (1991); ZuritaGotor et al. (2012). An interesting question is whether this perturbation to the 3D pair correlation function could be fundamentally related to the sequential layering reported here. The structural features near the wall should also affect two and manyparticle dynamics in the quasi2D suspensions, which can be characterized by twopoint microrheology. Finally, taking a more detailed account of interparticle forces such as strong electrostatic interactions may hopefully provide deeper understanding of the effects observed in this work.
Acknowledgements.
H.D. wishes to thank the Polish Academy of Sciences for its hospitality. This research has been supported by the Israel Science Foundation (Grants No. 8/10 and No. 164/14) and by the Marie Curie Reintegration Grant (PIRG04GA2008239378). A.S.–S acknowledges funding from the TelAviv University Center for Nanoscience and Nanotechnology. M.L.E.–J. and E.W. were supported in part by Narodowe Centrum Nauki (National Science Centre) under grant No. 2012/05/B/ST8/03010. J.B. would like to acknowledge the financial support from National Science Foundation (NSF) Grant No. CBET 1059745.Appendix A
We present a simple model to support a hypothesis that the discrepancy between the measured and calculated nearwall particle distributions stems from optical aberration caused by nonuniform optical properties of the suspension in the nearwall region. We assume that such aberration produces a nonlinear rescaling of the coordinate ,
(34) 
where is the actual and is the measured particle position. The rescaling (34) results in the corresponding transformation of the particle density
(35) 
To demonstrate that a distortion (34) can produce the observed shift and change of height of the features of the distribution , we consider an ad hoc distortion model with the transformation between the measured and actual vertical coordinates given by the equations
(36a)  
and  
(36b) 
where , , and are the model parameters. The transformation (36) describes positiondependent coordinate contraction with the amplitude gradually increasing in the region (the region where the second peak occurs according to the experimental data). The overall deviation of the Jacobian (36) from unity is proportional to the area fraction of the suspension layer.
Figure 17(a) compares the distorted distribution (35) with the corresponding untransformed distribution obtained from MC simulations of a HS suspension at the area fraction . Figure 17(b) presents the distorted distribution for the set of area fractions for which experimental results are depicted in Fig. 5(a). The parameter values of the transformation (36) are given in the figure caption.
The results show that the coordinate transformation (36a) shifts the position of the second particle layer to the left and produces a corresponding enhancement of the peak of particle distribution, similar to the experimentally observed features of the distributions depicted in Figs. 5(a) and 9(b). Thus our calculations provide indirect support to our opticaldistortion hypothesis. The distortion hypothesis can also explain why the measured and calculated occupation fractions and selfdiffusivities of the particles in the top layer agree well (see Figs. 10 and 14), in spite of the fact that the observed and calculated positions of the layer differ significantly.
It is an open question what the source of the distortion (34) might be. Since the suspension is imaged from above in our confocalmicroscopy system, we hypothesize that reflection of laser light from the first (bottom) particle layer results in stray illumination of the second layer, producing distorted particle height measurements. The optical distortion hypothesis can be verified by experiments using refractiveindex matched suspensions, but such investigations are beyond the scope of the present study.
References
 Van Winkle and Murray (1988) D. H. Van Winkle and C. A. Murray, J. Chem. Phys. 89, 3885 (1988).
 GonzálezMozuelos et al. (1991) P. GonzálezMozuelos, M. MedinaNoyola, B. D’Aguanno, J. M. MéndezAlcaraz, and R. Klein, J. Chem. Phys. 95, 2006 (1991).
 ZuritaGotor et al. (2012) M. ZuritaGotor, J. Bławzdziewicz, and E. Wajnryb, Phys. Rev. Lett. 108, 068301 (2012).
 Happel and Brenner (1983) J. Happel and H. Brenner, Low Reynolds Number Hydrodynamics (Martinus Nijhoff, The Hague, 1983).
 Perkins and Jones (1992) G. S. Perkins and R. B. Jones, Physica A 189, 447 (1992).
 Cichocki and Jones (1998) B. Cichocki and R. B. Jones, Physica A 258, 273 (1998).
 Walz and Suresh (1995) J. Y. Walz and L. Suresh, J. Chem. Phys. 103, 10714 (1995).
 Prieve (1999) D. C. Prieve, Adv. Colloid Interface Sci. 82, 93 (1999).
 Sholl et al. (2000) D. S. Sholl, M. K. Fenwick, E. Atman, and D. C. Prieve, J. Chem. Phys. 113, 9268 (2000).
 CarbajalTinoco et al. (2007) M. D. CarbajalTinoco, R. LopezFernandez, and J. L. ArauzLara, Phys. Rev. Lett. 99, 138303 (2007).
 Blawzdziewicz et al. (2010) J. Blawzdziewicz, M. L. EkielJeżewska, and E. Wajnryb, J. Chem. Phys. 133, 114703 (2010).
 Dufresne et al. (2000) E. R. Dufresne, T. M. Squires, M. P. Brenner, and D. G. Grier, Phys. Rev. Lett. 85, 3317 (2000).
 Cichocki et al. (2007) B. Cichocki, M. L. EkielJeżewska, and E. Wajnryb, J. Chem. Phys. 126, 184704 (2007).
 ZuritaGotor et al. (2007) M. ZuritaGotor, J. Bławzdziewicz, and E. Wajnryb, J. Fluid Mech. 592, 447 (2007).
 Anekal and Bevan (2006) S. G. Anekal and M. A. Bevan, J. Chem. Phys. 125, 034906 (2006).
 Michailidou et al. (2009) V. N. Michailidou, G. Petekidis, J. W. Swan, and J. F. Brady, Phys. Rev. Lett. 102, 068302 (2009).
 Cichocki et al. (2010) B. Cichocki, E. Wajnryb, J. Bławzdziewicz, J. K. G. Dhont, and P. Lang, J. Chem. Phys. 132, 074704 (2010).
 Loppinet et al. (2012) B. Loppinet, J. K. G. Dhont, and P. Lang, Eur. Phys. J. E 35, 62 (2012).
 Michailidou et al. (2013) V. N. Michailidou, J. W. Swan, J. F. Brady, and G. Petekidis, J. Chem. Phys. 139, 164905 (2013).
 CarbajalTinoco et al. (1996) M. D. CarbajalTinoco, F. CastroRoman, and J. L. ArauzLara, Phys. Rev. E 53, 3745 (1996).
 Schmidt and Löwen (1997) M. Schmidt and H. Löwen, Phys. Rev. E 55, 7228 (1997).
 Zangi and Rice (2000) R. Zangi and S. A. Rice, Phys. Rev. E 61, 660 (2000).
 Frydel and Rice (2003) D. Frydel and S. A. Rice, Phys. Rev. E 68, 061405 (2003).
 Han et al. (2008) Y. Han, Y. Shokef, A. M. Alsayed, P. Yunker, T. C. Lubensky, and A. G. Yodh, Nature 456, 898 (2008).
 Lin et al. (2000) B. H. Lin, J. Yu, and S. A. Rice, Phys. Rev. E 62, 3909 (2000).
 Dufresne et al. (2001) E. R. Dufresne, D. Altman, and D. G. Grier, Europhys. Lett. 53, 264 (2001).
 EkielJeżewska et al. (2008) M. EkielJeżewska, E. Wajnryb, J. Bławzdziewicz, and F. Feuillebois, J. Chem. Phys. 129, 181102 (2008).
 Cui et al. (2004) B. Cui, H. Diamant, B. Lin, and S. A. Rice, Phys. Rev. Lett. 92, 258301 (2004).
 Bhattacharya et al. (2005a) S. Bhattacharya, J. Bławzdziewicz, and E. Wajnryb, J. Fluid Mech. 541, 263 (2005a).
 Diamant et al. (2005) H. Diamant, B. Cui, B. Lin, and S. A. Rice, J. Phys.: Condens. Matter 17, S4047 (2005).
 Bławzdziewicz and Wajnryb (2012) J. Bławzdziewicz and E. Wajnryb, J. Phys. Conf. Ser. 392, 012008 (2012).
 Baron et al. (2008) M. Baron, J. Bławzdziewicz, and E. Wajnryb, Phys. Rev. Lett. 100, 174502 (2008).
 Lin et al. (1995) B. Lin, S. A. Rice, and D. A. Weitz, Phys. Rev. E 51, 423 (1995).
 Cichocki et al. (2004) B. Cichocki, M. L. EkielJeżewska, G. Nagele, and E. Wajnryb, J. Chem. Phys. 121, 2305 (2004).
 Peng et al. (2009) Y. Peng, W. Chen, T. M. Fischer, D. A. Weitz, and P. Tong, J. Fluid Mech. 618, 243 (2009).
 Zhang et al. (2014) W. Zhang, S. Chen, N. Li, J. Z. Zhang, and W. Chen, PLoS ONE 9 (2014).
 Skinner et al. (2010) T. O. E. Skinner, D. G. A. L. Aarts, and R. P. A. Dullens, Phys. Rev. Lett. 105, 168301 (2010).
 Kapfenberger et al. (2013) D. Kapfenberger, A. SonnSegev, and Y. Roichman, Opt. Express 21, 12228 (2013).
 Lee and Grier (2007) S.H. Lee and D. G. Grier, Opt. Express 15, 1505 (2007).
 Cheong et al. (2010) F. C. Cheong, B. J. Krishnatreya, and D. G. Grier, Optics Express 18, 13563 (2010).
 Behrens and Grier (2001) S. H. Behrens and D. G. Grier, The Journal of Chemical Physics 115, 6716 (2001).
 Israelachvili (1992) J. N. Israelachvili, Intermolecular and Surface Forces (Academic Press, London, 1992).
 Frenkel and Smit (2002) D. Frenkel and B. Smit, Understanding Molecular Simulation. From Algorithms to Simulations (Academic Press, New York, 2002).
 Bławzdziewicz and Wajnryb (2008) J. Bławzdziewicz and E. Wajnryb, Phys. Fluids. 20, 093303 (2008).
 Bhattacharya et al. (2005b) S. Bhattacharya, J. Bławzdziewicz, and E. Wajnryb, Physica A 356, 294 (2005b).
 Bhattacharya et al. (2005c) S. Bhattacharya, J. Bławzdziewicz, and E. Wajnryb, J. Fluid Mech. 541, 263 (2005c).
 Cichocki and Felderhof (1989) B. Cichocki and B. U. Felderhof, Physica A 159, 19 (1989).
 Sadlej et al. (2009) K. Sadlej, E. Wajnryb, J. Bławzdziewicz, M. EkielJeżewska, and Z. Adamczyk, J. Chem. Phys. 130, 144706 (2009).
 Abade et al. (2010) G. Abade, B. Cichocki, M. L. EkielJeżewska, G. Naegele, and E. Wajnryb, J. Chem. Phys. 132, 014503 (2010).
 Crocker and Grier (1996) J. C. Crocker and D. G. Grier, J. Colloid. Interf. Sci. 179, 298 (1996).
 Helfand et al. (1961) E. Helfand, H. L. Frisch, and J. L. Lebowitz, J. Chem. Phys. 34, 1037 (1961).
 Chen and Ma (2006) H. Chen and H. Ma, J. Chem. Phys. 125, 024510 (2006).