The rheology of confined colloidal hard discs

Colloids may be treated as “big atoms” so that they are good models for atomic and molecular systems. Colloidal hard disks are therefore good models for 2d materials and although their phase behavior is well characterized, rheology has received relatively little attention. Here we exploit a novel, particle-resolved, experimental set-up and complementary computer simulations to measure the shear rheology of quasi-hard-disc colloids in extreme conﬁnement. In particular, we conﬁne quasi–2d hard discs in a circular “corral” comprised of 27 particles held in optical traps. Conﬁnement and shear suppress hexagonal ordering that would occur in the bulk and create a layered ﬂuid. We measure the rheology of our system by balancing drag and driving forces on each layer. Given the extreme conﬁnement, it is remarkable that our system exhibits rheological behavior very similar to unconﬁned 2d and 3d hard particle systems, characterized by a dynamic yield stress and shear-thinning of comparable magnitude. By quantifying particle motion perpendicular to shear, we show that particles become more tightly conﬁned to their layers with no concomitant increase in density upon increasing shear rate. Shear thinning is therefore a consequence of a reduction in dissipation due to a weakening in interactions between layers as the shear rate increases. We reproduce our experiments with Brownian Dynamics simulations with Hydrodynamic Interactions (HI) included at the level of the Rotne–Prager tensor. That the inclusion of HI is necessary to reproduce our experiments is evidence of their importance in transmission of momentum through the system


I. INTRODUCTION
Confined fluids often exhibit modified flow behavior compared to the bulk due to coupling to the boundary.This may be static (i.e.structural alterations due to energetic or entropic interactions) or dynamic (e.g. the hydrodynamic influence of the wall).The shear viscosity of simple liquids increases by orders of magnitude in films less than ∼ 7 molecules thick [1][2][3] .This is a consequence of confinement-induced solidification near the boundary, driven by van der Waals interactions [4][5][6] .By contrast, the viscosity of water increases more modestly under similar conditions 7,8 as confinement-induced solidification is suppressed by the hydrogen bond network.
On the mesoscopic scale, rheological measurements of a) Electronic mail: paddy.royall@epsci.psl.eutwo-dimensional soft materials typically focus on interfacially adsorbed components including nanoparticles 9,10 , proteins [11][12][13][14] , lipid [15][16][17][18] , surfactant 19 and polymer 10,[20][21][22] monolayers or asphaltene films [23][24][25] , 2d foams 26,27 , Hele-Shaw emulsions 28 lipid bilayers 29 dusty plasmas 30 and colloids 11,[31][32][33][34] .The latter can be tailored to exhibit interactions similar to simple models of atomic and molecular systems, and since they explore phase space in much the same way colloids provide a means to probe the same phenomena, for example the effect of confinement, and other external fields such as shear, upon phase behavior [35][36][37] .Due to their mesoscopic size, colloids may be confined by walls that are smooth on the particle scale or by walls that are rough (like atomic and molecular systems).Another class of materials at a somewhat larger lengthscale is (athermal) granular matter and here too individual particles may be studied, and even the force networks between them 38 .Under vibration, granular mat-ter can mimic the phase behavior of thermal systems 39 .The rheology of granular matter has received considerable attention [40][41][42][43][44][45][46] .Like suspensions of colloids, granular materials become very much more viscous at high packing fraction, however the limiting case of divergence viscosity is jamming in the case of grains, rather than the thermal glass transition (which occurs at a lower packing fraction) in the case of colloids 47 .At the microscopic level, the shear response of colloidal and granular systems is distinguished in that colloidal particles usually do not come into physical contact with each other or with the walls of the sample cell.Therefore, dissipation occurs through hydrodynamic coupling and Brownian motion rather than friction due to direct contacts as in granular matter.Very occasionally, under extreme shear rates for example, contact can occur, leading to very sudden shear thickening 48,49 .
In colloidal systems at volume or area fractions where glassy dynamics are encountered, researchers report both increases and decreases in viscosity with respect to the bulk depending on the degree of confinement and the boundary details.Relaxation may be accelerated (decelerated) near smooth (rough) walls [50][51][52] .These effects are attributed to boundary-induced structural modification dependent on the shape, roughness, wetting characteristics or interaction potential of the wall (to name but a few possibilities) [53][54][55] .As a general rule, the formation of well-defined particle layers is associated with faster relaxation or reduced viscosity.
The bulk rheology of a hard-sphere-like (or hard-disclike) colloidal suspension depends on its volume (or area) fraction, φ 56 , and the shear rate, γ.At low shear rates, viscosity decreases with γ57-60 .For volume (area) fractions approaching the hard sphere (disc) glass transition, shear thinning occurs after the applied stress exceeds a yield stress on the scale of k B T /a 3 , where k B T is the thermal energy and a is the particle radius [59][60][61][62][63][64] .Upon increasing γ, hard sphere suspensions exhibit a Newtonian plateau followed by shear thickening at very high shear rates 57,58,65,66 .Shear thinning in these systems is attributed to stratification of particles decreasing resistance to flow, while shear thickening at high γ is a consequence of frictional particle surfaces 67 .
While the rheology of (3d) colloidal hard spheres has been extensively studied [68][69][70] attention to their 2d analogue, hard discs has focussed on quiescent (unsheared) systems [71][72][73][74][75][76][77][78][79][80] .Here, we focus on the rheology of confined colloidal hard discs.We present experiments on a quasi-hard-disc system of spherical colloids adjacent to a solid substrate confined by 27 identical particles optically trapped along a circular boundary and subjected to shear by boundary rotation.These are accompanied by computer simulation of a 2d system using Brownian Dynamics with hydrodynamic interactions (HI) included at the level of the Rotne-Prager tensor.Under quiescent conditions, our system has a bistable state between a hexagonal configuration that distorts the flexible walls [Fig. 1  (b)] and a layered fluid where hexagonal ordering is sup-pressed [Fig. 1 (a)].Hexagonal configurations exhibit voids adjacent to the wall, enabling relaxation mechanisms that are absent in the bulk 78 .Under shear, hexagonal configurations rotate as a rigid body, slipping only at the wall interface, while layered fluid configurations slip between each layer.At higher shear rates, the system is shear melted and only the layered fluid is observed 81 .It has subsequently been shown that simulations of a similar system exhibit distinct regions of shear thinning and thickening as a function of shear rate 82,83 .Here, we develop a layer-by-layer approach to determine local rheological properties directly from experiments and simulations via drag forces, and we characterise motion perpendicular to shear via mean squared displacements.Our approach opens a route to particle-level analysis of rheological properties in experimental systems.
The key finding of this work is that, in contrast with many molecular liquids [1][2][3] , strong confinement has little effect on the shear rheology of colloidal quasi-harddiscs.Our measurements are qualitatively and quantitatively aligned with simulations of hard discs under steady shear 62 .We measure a flow curve indicative of a yield stress, and shear thinning towards a Newtonian plateau.Motion perpendicular to shear is suppressed with no change in local density as shear rate is increased.We hypothesize that this suppression of perpendicular motion reduces the coupling between adjacent particle layers, reducing drag and dissipation, and consequently we measure a reduction in viscosity.This interpretation is illustrated in Fig. 1 (h-i).Our simulations reveal that the mechanism of momentum transfer in this system is dominated by hydrodynamic interactions, and the excluded volume interactions between the particles play a rather minor role.
This paper is organised as follows.Section II describes the experimental and simulation procedures, and details of our stress calculation.Full details of the drag coefficient determination, simulations and interlayer hopping are provided in the Appendices.Section III A presents key structural and dynamic measurements required for the measurements of local shear rate, stress and viscosity in Section III B. Section III C quantifies the motion perpendicular to shear.Finally, Section IV combines these measurements to develop our explanation of shear thinning.We conclude this work in section V. Coupling between layers depends on the range of radial motion within adjacent layers.In a slowly driven system, (h), particles explore a larger radial extent than in a quickly driven system, (i).For larger local shear rate, particles in adjacent layers are further apart, interlayer interactions are weaker and the effective viscosity is lower.This manifests as shear thinning.
vent, the particles quickly sediment, forming a quasitwo-dimensional layer adjacent to the lower glass coverslip, which is treated with Gelest Glassclad 18 to prevent particle adhesion.The gravitational length is l g /a = 0.030 ± 0.002, resulting in negligible out-of-plane particle motion in the z direction, while still being far enough from the substrate that any direct interactions with the substrate through contacts can be safely neglected.Interparticle interactions are of Yukawa form with a Debye length λ D ≈ 25 nm, which is sufficiently short that they may be considered quasi-hard-discs 84 .In the dilute limit, the time for a particle to diffuse its radius (the Brownian time) in the substrate-adjacent, quasi-two-dimensional layer is measured to be τ B ≈ 70 s.
Computer-controlled holographic optical tweezers in an inverted microscope are used to manually gather N = 75 particles, 27 of which are optically trapped to form a circular boundary of radius R which confines the remaining n conf = 48.The spring constant of the optical traps maintaining this boundary is extracted from the trapped particle trajectories and is determined to be k = 105(1) k B T a −2 .The colloids sediment to the bottom of a sample cell, and thus there is a (relatively) large amount of solvent above.We presume that any lo-cal heating effect is swiftly dissipated to the surrounding solvent.When we tested the system, for example at different strengths of the laser tweezers, we saw no signs of any such heating.
Using the optical tweezers, the boundary particles are translated along a circular path through the periodic displacement of the optical tweezers array in discrete steps of length a/4.The frequency with which the array of traps is updated defines the boundary rotation speed, which we characterise using the boundary Péclet number, Pe = τ B /τ D , where τ D is the time taken to drive a boundary particle a distance a.We consider experiments with Pe in the range 1.75 ≤ Pe ≤ 19.25.
Following the initiation of boundary rotation, the system is left to complete 5 full rotations before micrographs are acquired at a rate of 2 frames per second for up to 3 hours.No indication of anything other than a steadystate behavior was found.Particle trajectories are extracted from micrographs using standard algorithms 85 .Circular polar co-ordinates, r and θ, are defined from the system centre.

B. Simulation
We perform 2d Brownian Dynamics simulations of strongly screened charged colloids interacting via a Yukawa pair potential under highly ionic conditions.
with r denoting the inter-particle separation.The inverse screening length κ is chosen as κa = 14.85 and the contact potential V (r = 2a) ≈ V 0 k B T where V 0 = 0.85, informed by the experimental parameters, ensuring quasihard-disc behavior 81 .In particular, we found empirically that κa = 14.85 was sufficiently hard to closely reproduce key dynamical properties of the experimental system such as the way in which layers of particles slip past one another.Further details are available in Appendix A and Ref. 81 .We neglect polydispersity.Boundary particles experience harmonic potentials mimicking optical traps and are translated along a circular path at a prescribed rate.Hydrodynamic effects are accounted for at the Rotne-Prager level 86,87 in the presence of a planar substrate (Blake's solution 88 ), giving the overdamped equation of motion for a particle trajectory r i in a time step δt The mobility tensor µ ↔ ij comprises both the self-mobility and the entrainment of particle i by the hydrodynamic flow field created by conservative forces F j on particle j.This force stems from pair interactions and, for the boundary particles, the harmonic potentials.The random displacement, δW i , is sampled from a Gaussian distribution with zero mean and variance 2D 0 δt (for each Cartesian component) fixed by the fluctuationdissipation relation, where D 0 is the diffusion coefficient.Further details are provided in Appendix A. We have previously shown that these simulations faithfully reproduce experiments both qualitatively and quantitatively 81 .

C. Stress Calculation
Figure 1(a) shows that the particles organize into concentrated layers.This motivates our determination of the stress, by considering each layer, i.We presume that the angular velocity and radial position are constant within each layer, which, as we shall show in Section III is reasonable for our purposes.We thus perform a series of force balances to obtain the stress across each layer, σ i .Within each layer there are three forces which must sum to zero at steady state.There is some driving force (either optical forces or the transmitted force due to the motion of an external layer), F ext .This is balanced by some self-drag force representing dissipation within the layer, F self , and some additional force that drives the next internal layer, F int .In the case of layer 1 (the central layer), there is only self-drag.
The dilute limit single particle drag coefficient near the substrate, ζ emp = 1.9 × 10 −7 kg s −1 , is extracted from measurements of Brownian motion in the dilute limit, without any optical traps.However, multiple particles moving along a circular path of radius r near a substrate each experience reduced hydrodynamic drag due to the presence of the other particles 89 .This is the drafting effect.This drag reduction depends on the number of particles and r, which suggests that particles in different layers experience different drag coefficients.Based on the discussion in Appendix B we assume a constant drag coefficient throughout the system, ζ = 0.34ζ emp , which is reduced compared to the dilute limit drag coefficient due to many-body hydrodynamic effects, but is independent of radius.This effect is discussed in greater detail and this approximation is justified in the Appendix B. We therefore assume that the self-drag force has the same form for all particles and is proportional to velocity, Consider layer 5, the boundary, which consists of n 5 = 27 particles, each of which is located at radial position r 5 and subjected to an optical force F opt .All n 5 particles move with angular velocity ω 5 and each experiences a self-drag proportional to its velocity.Balancing the forces on layer 5 yields: where F 4 is the as of yet unknown force transmitted inwards to drive the motion of all n 4 = 21 particles forming layer 4. The force balance for layer 4 is then: where, once again, F 3 is the unknown force required to drive layer 3. Propagating the layer-by-layer force balance inwards to layer 1 reveals the retrospectively trivial result that the optical driving forces are exactly balanced by the self-drags experienced by all of the particles: Replacing the optical forces by the sum of the drag forces gives the unknown force required to drive layer i, F i , as the sum of the drag forces on layer i and all layers that are internal to i: Assuming that the inter-layer forces act at circular contact lines between layers at radial locations intermediate between the two layer centres, we calculate the twodimensional tangential stress on layer i:

III. RESULTS
We begin by recapitulating the structural and dynamic features necessary for our rheological analysis.We subsequently implement the analysis described in section II C to measure the shear rheology combining particle-level structural and dynamic information to extract stresses from drag forces.Finally, we quantify particle motion in the radial direction, perpendicular to shear.

A. Angular Velocity Profiles and Structure
Without shear, the system can adopt either layered fluid or locally hexagonal structures, shown in Fig. 1 (a) and (b) respectively 78 .When sheared at Pe 3, hexagonal structures can persist and rotate as a rigid body 81 .The hexagonal structure is characterised by multiple peaks in the radial density profile [black line in Fig. 1 (d)] and a flat angular velocity profile [black line in Fig. 1 (e)].For Pe 3, the system is shear melted and only layered structures are observed [coloured data in Fig. 1 (d)].The location of the ith peak in the layered density profile is labelled r i .Layer populations are fixed at n 1 = 3, n 2 = 9, n 3 = 15, n 4 = 21 and n 5 = 27.In these shear melted experiments, angular velocity decreases in a step-like manner from the boundary to the centre coloured data in Fig. 1 (e)], representing slipping between adjacent layers.The average angular velocity within layer i is denoted ω i .
The area fraction is estimated for layer i as φ i = πa 2 / A V i where A V i is the average Voronoi cell area for particles in layer i 79 .Figure 1 (e) shows that φ varies from φ 1 ≈ 0.8 at the centre, to φ 4 ≈ 0.76 adjacent to the boundary.At these densities, bulk hard discs are crystalline 90 .However, in shear-melted experiments, hexagonal ordering is inhibited by the curved boundary and the application of shear and our system is liquidlike 78,81 .
Shearing the system modifies the area fraction profile compared to the unsheared system [Fig. 1 (e)], enhancing φ 4 and suppressing φ 3 and φ 2 .However, once the system is shear melted, the area fraction profile is insensitive to Pe in the range of Pe investigated, as shown in Fig. 1 (f).

B. Shear Rheology
To characterize the shear rheology, we follow the analysis in section II C and treat the shear-melted system as a series of coupled layers, following reference 26 .Layers have fixed populations, n i .We assume all particles are located at radial position r i corresponding to the peaks in the density profile [Fig. 1 (d)] and move with angular velocity ω i [Fig. 1 (e)].Viscosity is η = σ/ γ where σ is the stress and γ is the shear rate.The shear rate experienced by layer i is obtained from the angular velocity profile as where ∆ω i+1,i = ω i+1 − ω i is the step down in ω between layer i + 1 and layer i, and the radial separation between the layers ∆r i+1,i = r i+1 − r i .We determine the stress throughout the system as discussed in section II C. In particular, the corresponding stresses, σ i , are obtained via a force balance on each layer, the key result is given in Eq. 7, that the force required to drive layer i is the sum of the drag forces on layer i and all layers internal to i.This force acts along a circular interlayer contact line, giving the stress on layer i, σ i , from which we define an effective viscosity, η i = σ i / γi , for each layer.
Figure 2 shows the results of this analysis for experiments (squares) and simulations (crosses).Stress is scaled by ζ −1 τ B , shear rate by τ B and viscosity by ζ −1 to facilitate the comparison.Points in the main panels are coloured according to the driving Pe, while the insets present the same data coloured by the layer number.Error bars are determined from angular velocity fluctuations in each layer determined from particle tracking.We see in Fig. 2 that the simulations do not reach such low shear rates low rates as the experiments.Now the only parameter that is directly set in both the experiments and the simulations is the rotation speed of the outer particle layer.Shear rates experienced by internal layers are determined by the physics of the system, ie the coupling between adjacent layers.So, the shear rates plotted are measurements, not directly controlled quantities.The viscosity of hard particle suspensions depends on shear rate and volume/area fraction [56][57][58][59][60][62][63][64] . The shar rate is largest in layer 4 and smallest in layer 1 and lower viscosity is measured for greater shear rate.However, Fig. 1 (f) and (g) show that local area fraction is largest in layer 1 and decreases towards the boundary.Therefore, measuring larger viscosity nearer the centre of the system could be a consequence of increased density, and unrelated to the local shear rate.However, Fig. 1 (g) shows that the area fraction of each layer is independent of Pe, and the inset to Fig. 2 (b) shows that the viscosity of a given layer does decrease as the shear rate increases.We expect the spatial variation in area fraction to make some contribution to the viscosity, but the shear rate dependence is unambiguous.
The data in Fig. 2 are remarkably consistent with theoretical predictions and simulations of unconfined binary hard discs under steady shear 62 .The flow curve suggests a dynamic yield stress 63,91,92 , beyond which, shear thinning is found.At high γ, the response approaches a Newtonian region.The solid (dashed) black lines show Herschel-Bulkley fits to the experimental (simulated) data of the form σ = σ y + k γν , where σ y is the dynamic yield stress and ν is the high shear rate exponent, which is unity for Newtonian flow.The fits yield ν = 0.82 ± 0.03 in experiment and ν = 0.97 ± 0.06 in simulation, consistent with weakly shear thinning (experiment) and Newtonian (simulation) flow at the largest shear rates studied.We relate the yield stress to the change in structure upon the application of shear, from a hexagonal to layered fluid configuration.It is possible to fit the data in Fig. 2 with a power-law σ = k γν (rather than the Herschel-Bulkeley fit).While the error bars in the case of the data points corresponding to low shear rates are large enough that the this regime can be fitted with a power law, in fact the fit at high shear rate is much worse than the Herschel-Bulkeley fit shown.Further dependencies are of course possible, such as a power-law with a shear-rate dependent exponent.While this would be most interesting to explore in the future, the quality of our existing data means that it may be challenging to accurately discriminate between such more complex dependencies.
The dynamic yield stress is ∼ 0.1k B T /a 2 .A dynamic yield stress is the stress required to maintain flow, and is smaller than the static yield stress, which must be exceeded during flow start-up.Although we focus on yielded systems, at very low Pe the system is unyielded and rotates as a rigid body 81 .Rigid-body rotation requires the local stress to be less than the static yield stress.Therefore, the stress measured in unyielded systems is a lower limit estimate of the static yield stress.In experiment, the largest stress for which the system does not yield, is ∼ 1.9 k B T /a 2 , which is comparable to yield stresses at the scale k B T /a 3 in hard spheres 59,60,63,64 .
The crossover from shear thinning to Newtonian flow occurs in the range γτ B ∼ 10 −1 to 10 0 , which is coincident with this crossover in simulations and theory of hard discs 62 , hard spheres 58 , and charged colloids 93,94 .That the rheology under such strong confinement bears even a qualitative resemblance to bulk behaviour is a remarkable finding.Many molecular liquids and complex fluids exhibit increases in viscosity of many orders of magnitude when confined to a few particle layers 1-3,50-52 , but this is not the case for colloidal hard discs, which reproduce their predicted bulk rheology when confined to a Radial MSD / a 2 (j) All layers in a single free centre simulation driven at Pe = 14.Points in (a-h) coloured according to driving Pe and in (i) and (j) according to layer number as indicated in the legends.Lines in (i) and (j) show fits to Eq. 9 as described in the text.
system only 8 particles across.We emphasize that the yield stress is likely related to the shear-induced melting of hexagonal order in this system.This ordering we have investigated previously 81,84 .

C. Motion Perpendicular to Shear
At the particle population of interest, and on the timescale of the experiment, particles are confined to layers (Appendix C).However, within the layers, they exhibit positional fluctuations in the radial direction as indicated by the width of the peaks in Fig. 1 (d).The mean squared displacement (MSD) in the radial co-ordinate is shown in Fig. 3 for all experiments (a-d) and simulations (e-h).Black points in the experiment panels show the radial MSD for Pe = 0, averaged over 5 experiments.These MSDs grow with time up to a plateau which represents the confinement within layers.In the case of the experiments, the increase to a (somewhat noisy) plateau is continuous, in some simulation data, notably for layer 2, there is some evidence of two timescales in the radial MSD [Fig.3(g)].Given that this two timescale behavior is only evidenced in layer 2 of the simulations, and not at all in the experiments, we focus on a single timescale, τ rad as determined below in Eq. 9.
Within each layer, the radial MSD approaches its plateau more quickly as Pe is increased, indicating a coupling between radial and tangential motion.The shear rate is largest in layer 4 and decreases towards the system centre.The long-time plateau is reached most quickly in layer 4, and progressively more slowly in layers 3 and 2. Therefore, larger γ causes faster radial motion.Additionally, in all but the central layer, shearing increases the plateau above that measured in the unsheared system (black data) indicating that the amplitude of radial motion is increased under shear.
Radial MSDs are fit with a function of the form which captures their growth with A the plateau height.
Example fits for all four layers in an experiment and a simulation at Pe = 14 are shown in Fig. 3 (i) and (j).We note that not all of the fits in (i) and (j) are of high quality.While we did find an improved fitting with the use of two timescales (i.e. two independent contributions to the right hand side of Eq. 9), in fact the small timescale proved to be very scattered, and its physical significance was unclear.A more sophisticated treatment than we have performed here with Eq. 9 would be interesting to explore in the future.
Figure 4 shows the radial MSD plateau height A and timescale τ rad as a function of the local shear rate for all experiments and simulations.The fit to the plateau height [Fig.4 (a)] reveals a downward trend with increasing γ in both experiment and simulation for all layers except layer 1. Figures 3 (a-d) show that shear enhances radial motion compared to the unsheared case (with the exception of layer 1).But, fitting reveals that this enhancement is greatest for lower shear rates.Particles are maximally radially mobile at the onset of shear melting to a layered fluid, and become increasingly confined in their layers at higher γ.
The timescale τ shows a similar dependence on γ [Fig.4 (b)], decreasing as γ increases, indicating a faster approach to the long-time plateau.When subjected to faster shear, particles more quickly explore their full range of motion perpendicular to shear.This is very clearly evident for layers 2 to 4 in experiments.Once again, layer 1 does not follow the trend evident in layers 2 to 4.

IV. DISCUSSION
The rheology of our system is similar to that of bulk colloidal hard discs, with evidence for a yield stress and shear thinning 62 .Our particle-resolved approach allows us to address the origins of this yield stress and shear thinning.The former we believe is related to the breakdown of hexagonal configurations found in the absence of shear.The latter seems to be similar to the 3d case of shear-induced stratification 58 .We have shown that motion parallel and perpendicular to shear are coupled.At larger local shear rate, particles more quickly explore their layer radially.Simultaneously, the extent of radial motion within the layer is reduced and they are more tightly confined.This is concurrent with shear thinning.
In colloidal hard sphere (and disc) systems, shear thinning is attributed to increased stratification parallel to shear as shear rate is increased 94,95 .When particles are organised into layers, interactions between them become weaker and consequently the resistance to flow is reduced.Concentric particle layers are enforced by the boundary of our system, and therefore the flow resistance is inherently lower than in a bulk suspension at comparable density.However, Fig. 2 shows that shear thinning is observed.We explain this fact using the data of Section III C. Radial motion is suppressed as shear rate is increased, but Fig. 1 (f) shows that local area fraction is independent of Pe.Faster shearing does not increase the packing density of particles, yet they exhibit more tightly confined dynamics.Radial motion brings particles to interact with particles in adjacent layers, which dissipates energy and resists flow.Suppressed radial motion with no concomitant increase in area fraction means that particles in adjacent layers are, on average, further apart and their interactions with one another are weakened.This results in a reduction in drag between layers, and therefore flow resistance, or viscosity.This interpretation is illustrated in Fig. 1 (h) and (i).Thus, the radial dynamics suggest the rheology -at larger shear rates, particles are more tightly confined to their layers at constant density.This reduces interlayer drag and leads to shear thinning.
What is the nature of the interactions that dominate the rheology of our system?Our simulations include hydrodynamic interactions at the Rotne-Prager level for the particle-particle interactions and also the Blake tensor for the coupling to the substrate.We found it necessary to include HI at this level to reproduce the behavior of the experiments.Pure Brownian Dynamics (without HI between the particles) leads to weak momentum transfer through the system, due to a far higher level of slip between the layers than is encountered in the experiments.Therefore we conclude that momentum transfer is dominated by HI and that steric effects due to the excluded volume of the particles play a secondary role.Thus, compared to molecular systems, where van der Waals interactions can lead to solidification near the boundary, here instead there are two important differences: (i ) there is no equivalent of the long-ranged van der Waals interactions in our hard discs, with no mechanism for solidification.(ii ) The momentum transfer is dominated by solvent-mediated hydrodynamic coupling, which is absent in molecular systems.The importance of hydrodynamic coupling between the particles suggests that similar behavior might be encountered in wet granular matter 48,49 .However in our system, the lack of direct particle-particle contacts (and increased ordering) enables a purely shear-thinning regime.

V. CONCLUSION
We have investigated the rheology of colloidal hard discs confined in a layered fluid configuration under shear and hexagonal configuration with very weak or no shear in experiment and simulation.Using a particle-level analysis, we infer the inter-layer forces from drag forces and the driving force exerted by optical tweezers.Since this system is dissipative and in steady state, balancing the forces on each layer allows the measurement of the local viscosity.This is the first experimental measurement of the viscosity of a hard-disc-like colloidal system under steady shear.
The flow curve is indicative of a dynamic yield stress and shows that the confined hard disc system exhibits shear thinning at low shear rates and approximately Newtonian behavior for γτ B 0.1.The dynamic yield stress is ∼ 0.1 k B T /a 2 .We find evidence for a static yield stress with a lower bound of 1.9 k B T /a 2 .This is remarkably similar to sheared, unconfined, bidisperse hard discs and bulk hard spheres.This is by no means an anticipated result as strongly confined systems regularly exhibit very different rheological responses to their bulk counterparts 1,3,51,52 .In our system, there is a change in structure, in that the system undergoes shear melting, which we relate to the yield stress.In the future, it would be intriguing to explore the response of this system to oscillatory shear and to investigate any yield stress in more detail.
Shear thinning in colloidal hard particle systems is due to shear-induced layering progressively reducing off-axis interparticle collisions as shear rate is increased, reducing the viscosity.By measuring the particle motion perpendicular to the direction of shear, we show that this is also the case in our system.At higher shear rate, radial motion is suppressed without a change in local density and therefore interlayer particle collisions are reduced, reducing the coupling and dissipation between layers, and the therefore the viscosity.
Colloidal hard discs and spheres under extreme confinement behave remarkably similarly to their bulk counterparts in both 2 and 3 dimensions and are qualitatively different to systems dominated by van der Waals interactions, for which viscosity increases massively on increasing confinement [1][2][3] .We attribute this to an absence of long-ranged vdW interactions in our system.Furthermore, we infer from our computer simulations that hydrodynamic coupling between the particles is the dominant mechanism of momentum transfer with excluded volume interactions playing a secondary role.This latter observation suggests that similar behavior might be observed in wet granular matter, in the case that interactions due to contacts between particles are not dominant 48,49 .
We have considered a particular geometry here, where a population of quasi-hard discs are effectively "corralled" by 27 tweezer particles arranged in a circle.Of course other geometries are possible.In the case of planar shear.we expect that behavior we observe would also be found, as one would expect particles to form layers parallel to the confinement, as is the case here.Depending on the particle spacing one might expect coupling between the packing of free particles and of wall particles.
In the future, it would be attractive to carry out a more complete inclusion of the HI than we have done here, for example with Lattice-Boltzmann dynamics or Stochastic Rotation dynamics.In particular, it would be useful to enquire whether such a description would exhibit the shear-thinning behavior seen in the experiments, and furthermore to explicitly probe lubrication phenomena neglected in the simulations we have performed here.It would also be interesting to develop a better description than the fit we have used to describe the radial MSDs in Eq. 9.
Jens Eggers, Yael Roichman, Thomas Speck and Todd Squires.HL was supported by the German Research Foundation (DFG) within the project LO 418/20-2.IW and CPR acknowledge support from the European Research Council (ERC Consolidator Grant NANOPRS, project number 617266). the 27 driven wall particles also from the harmonic trap potential V t .The random displacement δW i is sampled from a Gaussian distribution with zero mean and variance 2D 0 δt (for each Cartesian component) fixed by the fluctuation-dissipation relation, where D 0 is the diffusion coefficient.
In the simulations, the length scale is set by κ, the energy scale by k B T , and the time scale by τ = 1/(κ 2 D 0 ).The inverse screening length κ has been chosen as κa = 14.85,where the experimental value of the radius served as a reference.Consequently, the corral radius has been set to κR 0 = 128 yielding the experimental ratio of R 0 /(2a) ≈ 4.31.The high screening at κa = 14.85 together with the contact potential chosen as V (r = 2a) ≈ 0.85k B T ensures the quasi hard-disc-behaviour.Another crucial parameter in our system is the trap strength which has been set to λ = 0.42κ 2 k B T in order to mimic the laser trap strength in the experiments.The time step is chosen as δt = 10 −4 τ .Our simulations run for up to 8×10 3 τ , corresponding to approximately 35.5τ B with τ B being the experimental Brownian time, ie., the time one colloid needs to diffuse a length equal to its radius.

Appendix B: Determining the Drag Coefficient
The analysis described in the next section relies on knowledge of the drag forces acting on each each particle to extract dimensional forces, and therefore stresses and viscosities.Thus, it is necessary to know the single particle drag coefficient for our a = 2.5 µm radius polystyrene spheres undergoing quasi-two-dimensional motion near the substrate.To this end, we track the diffusive motion of these particles in the dilute limit, without any optical traps, and show the resulting mean squared displacement in Fig. 5 (a).The Stokes-Einstein relation says that this is linear in time, with gradient 4k B T /ζ.Thus, we perform a linear fit to the experimental data in Fig. 5 (a) and obtain an empirical measurement of the single particle drag coefficient to be ζ emp = 1.9 × 10 −7 kg s −1 .This represents an approximate four-fold increase over the Stokes drag coefficient ζ 0 = 6πηa = 4.7 × 10 −8 kg s −1 for spheres of 2.5 µm radius in water.
However, it is known that multiple particles moving along a circular path of radius r each experience reduced hydrodynamic drag due to the presence of the other particles 89 .This is the drafting effect.Ladavac & Grier give an approximate result at the level of the Stokeslet approximation for the single-particle drag coefficient when n particles of radius a are equally spaced on a ring of radius r, where r a, and near a planar substrate 89 .
This result is truncated at order (h/r) 3 .Here θ ij = (2π/n)(j − i) is the angular separation between particles i and j and h is the distance between the substrate and the particle centres.
Using the radial location of each layer's peak, r i , and the layer populations n i , we calculate this modification to the drag coefficient for our system, assuming h = a + l g , that is that the particles are located one gravitational length above the substrate.The dependance on the radius of the circular path and the number of particles suggests that particles in different layers of our system will experience different drag coefficients.The results of this calculation are shown in Fig. 5 (b).The drag correction calculated using equation B1 for layer 1 is negative, an unphysical result that is a consequence of the fact that equation B1 is valid only for r a, which is not true for layer 1 as r 1 ≈ a.Therefore, ignoring the invalid layer 1 result, it is evident that the drafting effect results in a reduced drag coefficient compared to the dilute limit singleparticle wall-corrected drag coefficient ζ w 0 .Furthermore, the anticipated dependence of ζ w N on radial position is ev-ident, though weak.If we identify our empirically measured wall-corrected drag coefficient, ζ emp , with ζ w 0 , then we anticipate that particles in our system will experience a drag coefficient of approximately 0.34 ζ emp due to the drafting effect, where 0.34 is the average of ζ w N /ζ w 0 over layers 2 to 5, indicated by the dashed line in Fig. 5 (b).
Estimating the true drag correction is further complicated by the fact that our system consists of a series of concentric and closely interacting layers.Equation B1 considers a single circular particle layer in isolation, and therefore cannot be a true description of our system.The effective drag coefficient experienced by a particle in layer i likely depends on the motion of particles in layers i + 1 and i − 1 in addition to the other particles in layer i.Therefore, the calculated reduction in drag is really only an order-of-magnitude estimate of the drafting effect.The sum in equation B1 is dominated by contributions from particles j = 2 and j = n, which are the neighbouring particles of particle 1.Since the separation between neighbouring particles is approximately the same in all layers, and since the radial dependence of of ζ w N is predicted to be weak, we treat the drag coefficient as being independent of radial position and equal to ζ = 0.34 ζ emp .This assumption is necessary to estimate the drag coefficient in layer 1, for which equation B1 is invalid, as indicated its prediction of a negative drag coefficient in this region.In the more densely packed system at n conf = 48 the populations of all layers are constant in time, as required for the rheological analysis described in the main manuscript.When the population is reduced to n conf = 44, however, particles occasionally move between layers.When layer populations are not fixed, the rheological analysis as described in the main manuscript cannot be applied, and so we have focused our attention on the n conf = 48 system.

Figures 1 (
Figures1 (a) and (b) show micrographs and (c) shows a side-view schematic of the system.Polystyrene spheres of radius a = 2.5 µm and polydispersity 2% are suspended in a 3 : 1 mass ratio mixture of deionised water and ethanol, loaded into a cell constructed using three glass coverslips and a microscope slide and sealed with epoxy.Owing to their density mismatch with the sol- FIG. 2.(a) Dimensionless stress and (b) viscosity as a function of dimensionless strain rate from experiments (squares) and simulations (crosses).Solid (dashed) line shows Herschel-Bulkley fit to experimental (simulated) data.Points are colored according to the driving Pe.Error bars are determined from angular velocity fluctuations in each layer.Insets show same data colored according to particle layer as indicated in the legend.

4 FIG. 3 .
FIG. 3. Radial mean squared displacements.All experimental data in (a) layer 4, (b) layer 3, (c) layer 2, and (d) layer 1. Black points and grey shaded regions represent average and standard deviation of behaviour in 5 unsheared experiments.All simulation data in (e) layer 4, (f) layer 3, (g) layer 2, and (h) layer 1. (i) All layers in a single experiment driven at Pe = 14.(j)All layers in a single free centre simulation driven at Pe = 14.Points in (a-h) coloured according to driving Pe and in (i) and (j) according to layer number as indicated in the legends.Lines in (i) and (j) show fits to Eq. 9 as described in the text.

FIG. 4 .
FIG.4.Parameters extracted from fits to radial MSDs as a function of shear rate.(a) Plateau height A (Eq. 9).(b) Timescale τ rad (Eq.9).Solid (dashed) lines guide the eye to trends in experimental (simulated) data.Points are colored according to layer number indicated in the legend.

FIG. 5 .
FIG. 5. (a) Mean squared displacement in the dilute limit in the absence of optical traps.Line is linear fit used to extract the empirical drag coefficient.(b) Drag coefficient correction term of the form in 89 as a function of radial location, calculated using experimentally measured layered structure.Horizontal dashed line shows the average over layers 2 to 5.
Appendix C: Interlayer Hopping at Lower Population

Figure 6
Figure6shows the time evolution of layer populations in experiments driven at Pe = 17.5 for confined populations (a) n conf = 44 and (b) n conf = 48.In the more densely packed system at n conf = 48 the populations of all layers are constant in time, as required for the rheological analysis described in the main manuscript.When the population is reduced to n conf = 44, however, particles occasionally move between layers.When layer populations are not fixed, the rheological analysis as described in the main manuscript cannot be applied, and so we have focused our attention on the n conf = 48 system.

1 FIG. 6 .
FIG. 6.Time evolution of layer populations ni measured in experiments driven at Pe = 17.5 with confined populations (a) n conf = 44 and (b) n conf = 48.