WorksPhysics and Astronomy
Closing of density wake in concentrated colloidal suspensions: two regimes governed by compressibility, friction, and self-propulsion
Abstract
When an obstacle moves through a compressible medium, it creates a density perturbation. Theoretical arguments regarding the shape of this density wake in Brownian suspensions either suppose infinite dilution or describe only the far field. Here, we study via experiments and numerical simulations the closing of the density perturbation just behind a circular obstacle in dense 2D suspensions of particles that are thermal or self-propelled. In the thermal case, we show that the characteristic length of the wake is a non-monotonic function of density, decreasing then increasing. We rationalise the scalings of both regimes by analytical arguments that highlight the role of the compressibility factor of the suspension. Based on this new understanding of the passive case, we show that the influence of self-propulsion observed in experiments and simulations can be understood as if the obstacle is effectively larger in an active bath, with non-trivial and distinct roles of the self-propulsion velocity and of the persistence time.
Introduction
The flow around an obstacle is a classical experiment in fluid dynamics since Stokes{sup}1. It is also a test bed for complex fluids{citep}2–5. Its spatially heterogeneous geometry enables to probe and discriminate complex rheologies{citep}6 with more efficiency than classical rheometric flows such as Couette{citep}7. In particular viscoelastic and viscoplastic fluids can display asymmetric flow field around a symmetric obstacle, unsteady flow{citep}5,8, and complex recirculations{citep}9. Viscoelastoplastic flows display features qualitatively different from viscoelastic and viscoplastic ones{citep}6. In granular media{citep}10–13 and suspensions{citep}14–17, moving an obstacle through the fluid also perturbs the particle density field, with a lower density region downstream and a higher density region upstream. Understanding the shape of this perturbation is essential to work out the mechanical response of the compressible medium to the obstacle motion{citep}15`. Furthermore, how the density wake opens and closes is of interest by itself, as it reveals the dynamical processes at play in the system and how emerging properties like pressure, viscosity or particle mobility couple with density.
The existence of a density wake stems from finite compressibility of the medium – particle concentration can be heterogeneous – and excluded-volume interactions between the obstacle and the particles. However, the closing of the density perturbation can have various physical origins. In the absence
of thermal motion or external vibration, granular particles do not diffuse and interact only at contact. In the absence of gravity, as for instance in a 2D horizontal granular layer, the wake is an empty cavity behind the obstacle which area is limited only by the possible extent of the jammed zone in front of the obstacle12. A density wake can thus exist even in the case of quasistatic motion of the probe. By contrast, in a granular pile the density wake is a dynamic phenomenon that balances obstacle velocity with rearrangement processes due to gravity18. Whether gravity effects are just the local free fall of grains on the obstacle size or a more global phenomenon akin to hydrostatic pressure is still debated18 and might depend on the history of the pile19.
In a dilute thermal suspension of hard spheres or hard disks, particle diffusion balances the density perturbation caused by obstacle motion15. This balance is summed up by the Péclet number , where is the velocity of the obstacle, the translational diffusion coefficient of the particles, and is the radius of the zone from which particle centres are excluded around the obstacle, with and the respective radii of the obstacle and a particle. At low Péclet the wake is a dipolar perturbation of the density field, decaying as in 2D, whereas at high Péclet the density becomes null behind the obstacle while a one particle-thick dense layer forms in front. Extensions of this thermal picture have been proposed for dilute suspensions of self-propelled particles20. When the persistence length of the particle motion is similar to , the effect of self-propulsion cannot be described as an effective diffusion, and a ballistic closing of the wake is observed20,21. However, the effect of density on wake closing remains unclear, both for passive and active suspensions. Indeed, most studies on obstacles in suspension of thermal or self-propelled particles have a microrheology focus: what is the drag force given the velocity of the probe, or what is the velocity given the applied force 15,20–24. Beyond the dilute regime, a simple replacement of by the collective diffusion coefficient in the expression of describes well microrheological responses17. Density wakes have been observed in simulations14,25 and experiments16,17 of dense thermal suspensions but their shapes have not been quantified. Hydrodynamic models have predicted power law decay ( in 2D) of the density perturbation at long range in dense assemblies of soft26 and hard colloids27, but do not describe the short range perturbation.
In this article, we address the question of wake closing in the context of dense passive and active colloidal systems by means of experiments, simulations, and continuum models. Experimentally, we consider an obstacle displaced at constant velocity in a monolayer of Janus colloidal particles that can be made self-propelled. We show that the density wake behind the obstacle has a characteristic length at short range. In the passive case, we characterise the dependence of this length with density and its scaling with obstacle velocity and size. Next, we characterize further the dependencies of this length in using numerical simulations of passive Brownian particles. We discover two regimes in density which we explain by introducing a passive continuum model that brings out the physical ingredients underlying our experimental observations: compressibility and friction. We then study how self-propulsion affects this scenario, both in simulations and in experiments, distinguishing between the respective effect of velocity and persistence time. Finally, we discuss how our new understanding of the density wake behind an obstacle can help the characterisation of passive or active suspensions.
Experimental setup
As in Ref28, we study 0.8 m-radius gold particles grafted with octadecanethiol29 and half-coated with Platinum. After suspension in pure water, the particles immediately sediment onto the flat bottom of the experimental cell, forming a bidimensional monolayer. The system is observed through a 20x fluotar objective of a Leica DMI 4000B microscope equipped with a custom-made external dark-field illumination and a camera (Basler ac-A2040-90um) recording at 10 fps. As shown in Figure 1(a), a very low in-plane apparent gravity is obtained by tilting the whole set-up with a small angle . We thus obtain a 2D sediment within the monolayer, with a density gradient along , as shown in Figure 1(b).

Figure 1. (a) Sketch of the experimental setup. Janus colloids form a monolayer that is sedimented by tilting the experimental setup by a small angle . A micromanipulator holds the obstacle that dips in the monolayer and is dragged at constant velocity along . (b) Detail of a dark-field microscopy image of the sedimented monolayer in the passive case. In-plane gravity acts along . The orange arrow represents the velocity at which the obstacle is dragged (from left to right, in all figures). The obstacle in motion appears as a circular zone depleted in colloids. The density wake is visible as a low-density area behind the obstacle.
The obstacle was forged from a borosilicate rod (BR-100-15, 1 mm, 15 cm) using a micropipette puller (P-1000 Sutter Instrument). We used two different protocols to obtain two obstacle sizes. After detachment of the two halves of the rod, we performed 30 cycles of fire polishing of the tip in order to obtain a smooth hemisphere. The rod is held by a micromanipulator (MP-285A from Sutter Instruments), so that the obstacle can be lowered in the monolayer within 0.5 m from the surface.
face of the container and then moved along an isodensity line (within ), as shown in Figure 1(b). Each measurement consists of two back-and-forth motions of at respectively and , separated by a waiting time. We find that a steady state is reached after , and thus average our data over the second half of the test. As the obstacle appears on the image only as a featureless dark zone contiguous to the wake, we use Euclidean distance transform to locate the obstacle and measure the apparent exclusion radius on each image. We have kept only experiments where the relative variations of within the time interval of interest are less than 5%. However, from experiment to experiment varies between and for the small obstacle depending on the distance between the tip of the obstacle and the surface of the container. For the large obstacle, varies between and . We use the trackpy library30 to locate the position of each particle, except in the vicinity of the large obstacle where uneven light reflection prevents quantitative image analysis. An undisturbed density profile is measured far from the obstacle, from which we can extract the equation of state of the suspension31. Fitting the divergence of this equation of state with respect to density with an analytical expression for a solid of hard disks32, we extract the effective hard disk radius of the particles . We can thus define a packing fraction as . The density field is then binned on a -spaced grid in the moving frame of the obstacle and time-averaged.
Our experimental parameters are the position along the density gradient at which we lower the obstacle, which sets the density, and the concentration of we add to induce self-propulsion of the particles, which sets the activity. The self-propulsion is often quantified by comparing the persistence length to the radius of the particle , where is the self-propulsion velocity under dilute conditions, and its persistence time (inverse rotational diffusivity)33,34,35. As in34,35, we quantify the activity by measuring the effective sedimentation length from the sedimentation profile in the dilute regime. The ratio between active and passive (no ) sedimentation lengths is then converted to , where and have been measured in a dilute case with vanishing .
Experimental results in passive cases
Figure 2 displays typical snapshots of a moving obstacle in a passive sediment. At all considered experimental conditions, we observe a density wake with an almost empty zone behind the obstacle. However, the length of this wake changes by about two orders of magnitude in units of the particle size.

Figure 2. Experimental micrograph of the wake for two different obstacle radii: (a) and (b) , at roughly the same density ( and ), and the same velocity . The two pictures are scaled so that appears the same to highlight the nonlinear dependence of wake length on obstacle size.
For instance, in Figure 2(a), the wake is a few particle size long, less than the obstacle diameter, whereas in Figure 2(b) the wake extends over hundreds of particle size. Between the two panels, only the obstacle diameter changed significantly (by a factor of ) although we adjusted the scales of the pictures so that the two obstacles are displayed at the same size. This representation highlights qualitatively that the wake length responds non linearly to the size of the obstacle.
A more precise description can be obtained by time-averaging the density in the frame of the obstacle, as shown in Figure 3(a). Here, we observe both the low-density wake behind the obstacle and an accumulation of particles in front. Contrary to the well-studied dilute case15, in this dense regime several layers of particles stack in front of the obstacle. More generally, the density field is structured in layers that follow the obstacle curvature in front and on the sides, but that detach behind, are advected and eventually close well-behind the obstacle. This is qualitatively consistent with 3D experiments at high packing fractions17, but strikingly different from recent analytical predictions in which the wake would flare out behind the obstacle before closing27.

Figure 3. (a) Map of the packing fraction time-averaged in the frame of the obstacle in a passive case where , see colour bar below. A Gaussian smoothing of size has been applied. The excluded zone is shown as a red circle of radius . The white disk represents the obstacle of radius . (b) Profiles along at various positions behind the obstacle indicated on (a) and (c) by vertical dashed lines of corresponding colours. A Gaussian smoothing of size has been applied. (c) Profile along that is -averaged on , i.e. dotted black lines in (a) and (b). A baseline measured at rest has been removed. The horizontal black dotted line is the density at rest at the position of the moving obstacle. The cyan and magenta curves are the best fits by exponential decay respectively in front and behind the obstacle. The respective decay lengths are shown as double arrows below (a). The horizontal gray line in (b) and (c) highlights the hard disk close packing fraction .
To study more quantitatively the density perturbation around the obstacle, we show in Figure 3(b) transverse density profiles at various positions along the wake. The blue curve is one radius ahead of the obstacle where the density reaches close packing. The orange curve is tangent to the exclusion zone just behind the obstacle and shows a deep minimum where the density goes to zero. This minimum becomes shallower (green curve) further away behind the obstacle until it disappears.
To quantify the wake more precisely, we average the density vertically over two particle diameters and plot the longitudinal density profile shown in Figure 3(c). Immediately in front of the obstacle, we observe that the density can reach values above close packing . This is reminiscent of the jammed area observed in granular experiments12. In order to extract a length characteristic of the closing of the wake, we fit behind the obstacle with an exponential decay (magenta curve). We call the resulting length the closing length. A similar procedure performed ahead of the obstacle (cyan curve) yields a second length that we call the front length. Determining require quantitative particle tracking in the immediate vicinity of the obstacle and is thus only accessible for the small obstacle.
In Figure 4(a), we show that the closing length in passive cases decreases with increasing packing fraction. Here the packing fraction is measured while the obstacle is at rest and averaged over the positions that the obstacle will occupy during its motion. We also observe that the closing length depends on the velocity of the obstacle. We thus construct the closing time in order to collapse the data for the two obstacle velocities, as shown in Figure 4(b).

Figure 4. (a) Experimental closing length as a function of packing fraction for two obstacle velocities in the passive case for small obstacles (). (b) Same data for the closing time. (c) Closing time rescaled as a function of packing fraction for two obstacle velocities in the passive case. Solid black curve corresponds to . (d) Same data with the inverse scaling in . Dashed black line is an affine guide for the eye.
To characterise further the scaling of we will turn to numerical simulations. From this characterisation, we will sketch analytical arguments that will extract the important physical ingredients to make sense of both numerical and experimental measurements.
Numerical simulations
To extend the validity of our experimental observations to a wider range of parameters, we perform numerical simulations of active Brownian particles, which have proved to be a suitable model for Janus colloids28,36.
We thus consider active Brownian particles in a box of size whose positions and orientations evolve in a two-dimensional space as
where is the friction coefficient, is the self-propulsion force, is the rotational diffusivity, is the magnitude of the gravitational force acting along the direction , and and are centred unit-variance Gaussian white noises. In the absence of interactions the bare self-propulsion velocity is , and the passive dilute translational diffusivity is . and are respectively the inter-particle forces resulting from a Weeks-Chandler-Andersen (WCA) potential, and the repulsive force between particles and the moving obstacle, whose centre of mass evolves as
where is the constant obstacle velocity. As in the experiments, the obstacle moves at constant velocity rather than at constant force and the obstacle position is used to probe different densities. See the Appendix for definitions of the different inter-particle potentials and for simulation details. We stress that, in such simulations, the viscosity, the pressure, and the compressibility are emerging properties, whereas the friction coefficient is a parameter of the simulation.
The closing length is estimated by the same procedure as in experiments. In Figure 4(c) we show the -dependency of the closing time at fixed obstacle size and friction coefficient in the passive case (). As in experiments, we observe a decreasing regime where the closing times at various obstacle velocities collapse. However, at higher densities the closing time also displays an increasing regime where the different obstacle velocities do no collapse, indicating a different scaling in . In Figure 4(d) we plot instead that displays a good collapse within the increasing regime. We thus observe two distinct regimes in densities of possibly different physical origins. The minimum between the two regimes shifts towards higher densities when increases. We also notice that the density just behind the obstacle is null () in the whole decreasing regime as observed in experiments. However the density starts to rise significantly above zero at the crossover, up to far in the increasing regime (not shown).
In Figure 5(a) we show the -dependency of at fixed obstacle velocity and friction coefficient in the passive case (). For all radii, we observe the decreasing regime followed by an increasing one. The former shifts toward higher densities as increases, whereas the latter is independent of . When rescaled by , as shown in Figure 5(c), the various values of collapse in the decreasing regime. In Figure 5(c), we confirm that increasing the obstacle radius pushes the transition between the two regimes towards higher packing fractions. In Figure 5(b), we show that the increasing regime disappears from our accessible packing fraction window as we increase the friction. In Figure 5(d), we show that the curves collapse on the decreasing regime once scaled by .

Figure 5. Closing length as a function of packing fraction in passive simulations with gravity. (a) Different obstacle radii. (b) Different friction parameter . (c,d) Same data as (a,b) rescaled by . Solid black line is that is the inverse compressibility factor of a hard disk fluid. It corresponds to the solid coloured lines in (a). Dashed black line in (a) is an affine guide for the eye, it corresponds to the dashed coloured lines in (c).
In summary, in the decreasing regime the closing length scales as times a decreasing function of , whereas in the increasing regime the closing length scales as and is independent of . The decreasing regime is characterised by a cavity (quasi zero density) behind the obstacle, whereas in the increasing regime the density is significantly above zero behind the obstacle. In the following, we will use these observations to extract the physics underlying both regimes.
Theoretical arguments
Continuum description
We now rationalise the observed regimes for the closing length. Let us describe the passive colloidal bath by its density , its pressure , and the friction coefficient of the colloids with the substrate and the surrounding fluid. Particle number is conserved, , and in the overdamped regime the particle current is driven by pressure gradients against friction,
A pressure difference acting over a distance thus moves the colloids at a velocity . This single estimate is all we need below; the only question is which pressure difference and which distance are relevant.
Two lengths and a Péclet number. Two lengths enter the problem. The first is structural: it is the obstacle radius . The second is dynamical: it is the distance over which the largest possible depression, i.e. , creates a flow that is just fast enough to keep up with the obstacle, , that is
where the compressibility factor , which quantifies the deviation from ideal gas behaviour. Equivalently, is the usual advection length built on an effective collective diffusion coefficient . The ratio of the two lengths defines a Péclet number,
which can equivalently be read as the ratio of the obstacle velocity to the refilling velocity that the full bath pressure imposes over the obstacle size. For systems of passive particles interacting via repulsive interactions, increases with density, thus decreases with density. In the dilute limit, and , so that Eq. (5) reduces to the Péclet number introduced above; in dense suspensions it is smaller by the factor .
Large Péclet number: refilling a cavity. For the obstacle moves faster than the bath can follow, and it leaves behind it an empty region of width , as observed in Figure 3. The bath then pushes colloids into this void with the full pressure acting over the distance , i.e. at velocity , and the cavity closes after a time . The corresponding closing length is
It reproduces the signatures of the decreasing regime: is proportional to and decreases with density through . The solid black curves in Figure 4(c) and Figure 5(c,Figure 5d) represent Eq. (6) using the compressibility factor of a hard disk fluid ^{237}, which is in good agreement with the data for packing fractions below 0.69, where a fluid state is expected.
Small Péclet number: no cavity. For the bath refills faster than the obstacle opens a gap. The density perturbation does not reach zero behind the obstacle. Nevertheless, we will assume that the pressure just behind the obstacle is negligible with respect to the bath pressure , thus . This assumption is justified at high packing where is a steep non linear function of density. In this regime the extent of the wake is itself
The scaling of Eq. (7) is compatible with the observations of the numerical simulations in the increasing regime: is independent of , scales as and increases with density through .
b. Crossover. Equations (6) and (7) can be summarised as for and for : the closing length is always larger than the obstacle, and is minimal, , at , corresponding to a compressibility factor . As the density increases, decreases and the system goes from the first regime to the second, which explains the minimum of in Figure 4 and Figure 5.
c. Relation to previous works. When the obstacle only weakly perturbs the bath, Eq. (3) can be linearised about the mean density , giving a diffusion equation for the density perturbation with the collective diffusion coefficient , with ^{26,27,38}. The same two lengths then appear, with replaced by : a “diffusive wake” behind the obstacle, of length , and a perturbation in front of it, of extent ^{27}. Our argument differs in two respects. First, we consider regimes where the pressure behind the obstacle is either zero (large ) or negligible with respect to the bath pressure (low and high packing), so that the linearisation is not justified. Refilling is then driven by the full pressure rather than by the stiffness , which turns into in Eq. (6). In other words, the relevant collective diffusion coefficient is instead of as in^{17}.
Second, the length of Eq. (7) was obtained by Banerjee et al.^{27} for the perturbation in front of the obstacle, while the perturbation behind is always described by the diffusive scaling of Eq. (6). By contrast, we find that Eq. (7) also controls the wake behind the obstacle at small , where the diffusive wake, of length , is hidden by the obstacle itself. This length cannot be recovered from the dipolar density perturbation found by Squires and Brady^{15} at small Péclet number. Furthermore, we note that in Figure 4(d) and Figure 5(a), the dashed black lines hint toward a increasing regime, which would imply in Eq. (7). This is consistent with phenomenological observations of Banerjee et al.^{27} but not with the equilibrium equation of state of the fluid or of the crystal of hard disks.
Compressibility factor
Our study on passive systems, in particular Eq. (6) and Eq. (7), suggests that the compressibility factor is the natural variable to describe the density dependence. However, our data suggest that we cannot rely on analytical expressions for hard disks in the dense regime. Furthermore, in active systems, the function depends on activity and can display non-monotonic behaviour as a function of for high enough activities^{35,39,40}. Thus, in order to study the effect of activity, we ought to plot our data as a function of rather than as a function of . To this end, in both our simulations under gravity and in our experiments, we measure the compressibility factor directly at the altitude of the obstacle from the sedimentation profile. Following^{31}, we define the hydrostatic pressure from an integral of the density estimated at rest or far from the obstacle:
where is the buoyant mass of a particle, the in-plane gravity. In terms of packing fraction, the compressibility is then
We measure the sedimentation length at the top of the passive sedimentation profile where the dilute sediment behaves as an ideal gas and . In practice, we integrate only up to a threshold altitude where the density is long enough to be in the dilute regime but high enough to be estimated reliably. The integral in is estimated from the dilute exponential law and is negligible at the densities considered in the present study. Previous studies have shown that the compressibility factor of active Janus colloids or ABPs can be obtained from the same method^{35}.
Active systems
In Figure 6(a) we show the -dependence of the closing time at various activities but at the same friction, obstacle radius and velocity. In the passive case, we confirm the dependence in the decreasing regime (dark blue solid line) and the dependence in the increasing regime (dashed magenta line). With increasing activity, we observe that the increasing regime remains unchanged, but that the decreasing regime shifts up, and the transition between the two regimes shifts towards larger values of . Experimentally, we observe only the decreasing regime, as shown in Figure 6(b). We attribute the absence of increasing regime to a higher translational friction in our experiments due to the proximity of the colloidal particles to the bottom of the container. As seen in Figure 5(b,Figure 5d), larger friction can push the increasing regime to inaccessibly high densities. Nevertheless, we confirm experimentally that the decreasing regime shifts up with activity.

Figure 6. Effect of activity. Rescaled closing time as a function of the measured compressibility factor at different activities for (a) simulations and (b) experiments. The rescaling is according to the decreasing regime (solid lines), since , and are constant in (a) a collapse of the increasing regime is also expected (dashed magenta line). On the colour bar, persistence lengths of the simulations are shown as red horizontal lines, whereas persistence lengths of the experiments are shown as vertical error bars.
Referring to Eq. (6), a shift up in the decreasing regime can be due either to a decrease in effective temperature, or to an effective increase of the friction, the obstacle velocity or the obstacle radius. However, according to Eq. (7) any of these changes would also affect the increasing regime, except for the radius . Therefore the effects of activity are only compatible with an effectively larger obstacle. We thus fit the decreasing regime at each activity to extract the effective obstacle radius . Extending arguments of Burkholder and
Brady20 outside of the dilute regime, we expect that active effects become important when becomes comparable to or larger than . We thus plot the increment of with respect to the passive case as a function of in Figure 7. In order to collapse simulations performed at different persistence times, we have to scale the effective radius increment with , i.e. the distance moved by the obstacle during a persistence time of the particles. Empirically, we found that numerical data display a crossover between a linear regime for and a cubic regime for . Although the persistence length range explored in experiments is comparable to simulations, experimental data explore only the and display the cubic regime with a different prefactor compared to simulations.

Figure 7. Increase of effective obstacle radius with activity in both experiments (Janus) and simulations (ABP, ), scaled by the distance moved by the obstacle during the persistence time of the particles, as a function of the ratio between particle velocity and obstacle velocity. Dashed lines are power laws of exponent 3. Dotted line is a power law of exponent 1.
Quantitative disagreement could be due to an underestimation of either or by a factor , or to physical ingredients absent from the ABP simulations: 3D nature of the rotational diffusion41, hydrodynamic or phoretic alignment of particle polarity with the surface of the probe28, or hydrodynamic interactions within the suspension42.
Conclusion and perspectives
We have quantified the length of the density wake at short range behind a moving obstacle in a dense monolayer of active colloidal particles in experiments and numerical simulations. Analytical arguments from a passive hydrodynamic model taking into account friction and compressibility of the suspension enabled us to identify two regimes where the length of the wake successively decreases and then increases with density. Experiments with passive particles are consistent with the decreasing regime where the physics is governed by the presence of an empty region just behind the obstacle, in particular the non linear scaling in obstacle size.
We found that endowing the particles with persistent self-propulsion has little effect on the increasing regime (higher densities) but that it increases the length of the wake in the decreasing regime (lower densities). The effect of activity on the length of the wake in both decreasing and increasing regimes can be understood in terms of an effectively larger obstacle. When the obstacle is slow with respect to the self-propulsion of the particles, we found an increment in effective radius proportional to the distance moved by the obstacle during the persistence time of the particles. In the decreasing regime, the effect of activity leads to a longer wake than in the passive case at constant compressibility. This finding goes against previous analytical and numerical investigations in the dilute regime that found a faster ballistic closing of the wake for active systems at constant density20,21. We attribute this apparent discrepancy to the fact that the compressibility factor controls wake closing and that depends on both density and activity. At low density where increases with the propulsion velocity39, we expect a decrease in wake length that could mask a modest increment in effective radius directly due to activity.
To sum up, the density wake behind an obstacle is a sensitive measurement of the compressibility of a suspension of passive or active particles, and its dependence on density. In particular at densities where the suspension at rest would be crystalline, the movement of the obstacle continuously destroys the crystalline order, allowing a measurement of the compressibility of the amorphous system behind the obstacle. However, the rear of the obstacle is not the whole story. In front of the obstacle we observe large densities, even higher than close packing. Here, we predict a much higher viscosity, or even an elasto-plastic behaviour that might dominate the shape of the density perturbation in front. More generally, our study has characterised a part of the density perturbation around a moving obstacle, that is not as simple as in the dilute case. From the whole density field, one could infer the average force acting on the obstacle and progress toward an understanding of the microrheology of passive and active suspensions.
Acknowledgments
We thank Shashi Thutupalli, Sriram Ramaswamy, and Patrice Bacchin for useful discussions. ML and CCB acknowledge support from the ANR ActivSoliPolar (ANR-24-CE30-6841). GD, CCB and HDA acknowledge support from the MITI Rheodaco. GS acknowledges support from the ERC Advanced Grant ActBio (funded as UKRI Frontier Research Grant EP/Y033981/1). VD dedicates this work to the memory of Rudi Podgornik, a mentor and friend.
Author declarations
a. Conflict of interest None.
b. Author contributions Investigation (simulations) by GS (lead). Investigation (experimental) by GD (lead) and BB. Data curation by GD, GS, BB, and ML. Formal analysis and software by GD, GS, HDA, and ML. Methodology by VD, HDA, and ML. Conceptualization, funding acquisition, project administration, supervision by FG, HDA, CCB and ML. Validation by HDA, CCB, and ML. Visualisation: ML. Writing original draft : GS, FG, VD, CCB, HDA and ML (lead). Writing review & editing :
c. Data availability statement Data for this article, including density maps from experiments and simulations, tables of intermediate results as well as the Matplotlib43 scripts to generate the figures are available at Zenodo at https://doi.org/10.5281/zenodo.23480622
Numerical details
To integrate the dynamics corresponding to Eqs. (5)-(6), we use a Euler scheme with adaptive time steps. At each iteration, the time step is set to , with being the interaction range of the WCA interparticle potential, is the largest absolute Cartesian component of the interparticle force acting on any particle at the current time step, and for passive simulations and for active simulations with . Interparticle forces results from a Weeks-Chandler-Andersen (WCA)
where is the depth of the potential, is the particle diameter, and is the interaction range. We use to retain only the repulsive part. The obstacle is modeled as a diverging radial potential
with and . Furthermore, particles can experience a gravitational force to further match the experiments. Since particles experience a gravitational force , we add a confining potential along the axis to avoid vertical currents. Particle are thus confined with a confining potential
with , , , and . We initialize particles in a box of sizes and with an initial density , and we let the particle relax before setting the obstacle into motion. Other simulation parameters are: , , and .
References
- [1]Stokes, G. G. (1851). On the Effect of the Internal Friction of Fluids on Pendulums. Cambridge Phil. Soc. Trans, 9, 8.
- [2]²L. Leal, M. Denn, and R. Keunings, “Lake arrowhead workshop special issue papers—introduction,” Journal of Non-Newtonian Fluid Mechanics 29, 1–8 (1988).DOI
- [3]³R. A. Brown and G. H. McKinley, “Report on the VIIIth international workshop on numerical methods in viscoelastic flows,” Journal of Non-Newtonian Fluid Mechanics 52, 407–413 (1994).DOI
- [4]⁴M. Alves, P. Oliveira, and F. Pinho, “Numerical Methods for Viscoelastic Fluid Flows,” Annual Review of Fluid Mechanics 53, 509–541 (2021).DOI
- [5]⁵S. J. Haward, C. C. Hopkins, S. Varchanis, and A. Q. Shen, “Bifurcations in flows of complex fluids around microfluidic cylinders,” Lab on a Chip 21, 4041–4059 (2021).DOI
- [6]⁶I. Cheddadi, P. Saramito, B. Dollet, C. Raufaste, and F. Graner, “Understanding and predicting viscous, elastic, plastic flows,” The European Physical Journal E 34, 1 (2011).DOI
- [7]⁷I. Cheddadi, P. Saramito, and F. Graner, “Steady couette flows of elastoviscoplastic fluids are nonunique,” J. Rheol. 56, 213–239 (2012).DOI
- [8]⁸C. C. Hopkins, S. J. Haward, and A. Q. Shen, “Upstream wall vortices in viscoelastic flow past a cylinder,” Soft Matter 18, 4868–4880 (2022).DOI
- [9]⁹M. T. Arigo and G. H. McKinley, “An experimental investigation of negative wakes behind spheres settling in a shear-thinning viscoelastic fluid,” Rheologica Acta 37, 307–327 (1998).DOI
- [10]¹⁰D. Chehata, R. Zenit, and C. R. Wassgren, “Dense granular flow around an immersed cylinder,” Physics of Fluids 15, 1622–1631 (2003).DOI
- [11]¹¹R. Candelier and O. Dauchot, “Journey of an intruder through the fluidization and jamming transitions of a dense granular media,” Physical Review E 81, 011304 (2010).pubmed.ncbi.nlm.nih.gov/20365365
- [12]¹²E. Kolb, P. Cixous, N. Gaudouen, and T. Darnige, “Rigid intruder inside a two-dimensional dense granular flow: Drag force and cavity formation,” Physical Review E 87, 032207 (2013).DOI
- [13]¹³A. Seguin, C. Coulais, F. Martinez, Y. Bertho, and P. Gondret, “Local rheological measurements in the granular flow around an intruder,” Physical Review E 93, 012904 (2016).DOI
- [14]¹⁴I. C. Carpen and J. F. Brady, “Microrheology of colloidal dispersions by Brownian dynamics simulations,” Journal of Rheology 49, 1483–1502 (2005).DOI
- [15]¹⁵T. M. Squires and J. F. Brady, “A simple paradigm for active and nonlinear microrheology,” Physics of Fluids 17 (2005).DOI
- [16]¹⁶A. Meyer, A. Marshall, B. G. Bush, and E. M. Furst, “Laser tweezer microrheology of a colloidal suspension,” Journal of Rheology 50, 77–92 (2006).DOI
- [17]¹⁷I. Sriram, A. Meyer, and E. M. Furst, “Active microrheology of a colloidal suspension in the direct collision limit,” Physics of Fluids 22, 062003 (2010).DOI
- [18]¹⁸P. Rognon and E. Guillard, “Drag forces in granular materials,” Soft Matter 21, 7005–7013 (2025).DOI
- [19]¹⁹F. Guillard, Y. Forterre, and O. Pouliquen, “Depth-Independent Drag Force Induced by Stirring in Granular Media,” Physical Review Letters 110, 138303 (2013).pubmed.ncbi.nlm.nih.gov/23581384
- [20]²⁰E. W. Burkholder and J. F. Brady, “Nonlinear microrheology of active Brownian suspensions,” Soft Matter 16, 1034–1046 (2020).DOI
- [21]²¹M. Knežević, L. E. Avilés Podgurski, and H. Stark, “Oscillatory active microrheology of active suspensions,” Scientific Reports 11, 22706 (2021), arXiv:2109.11218.arXiv
- [22]²²A. M. Puertas and T. Voigtmann, “Microrheology of colloidal systems,” Journal of physics. Condensed matter : an Institute of Physics journal 26, 243101 (2014).DOI
- [23]²³C. Reichhardt and C. J. O. Reichhardt, “Active microrheology in active matter systems: Mobility, intermittency, and avalanches,” Physical Review E 91, 032313 (2015).pubmed.ncbi.nlm.nih.gov/25871116
- [24]Deng, E., & Brady, J. F. (2022). Forced microrheology of active colloids. Journal of Rheology, 66, 955–972.
- [25]Weber, R., & Thakur, S. (2012). Hydrodynamic interactions in active colloidal crystal microrheology. Physical Review E, 86, 057302.
- [26]²⁶V. Démery, O. Bénichou, and H. Jacquin, “Generalized langevin equations for a driven tracer in dense soft colloids: construction and applications,” New Journal of Physics 16, 053032 (2014).DOI
- [27]²⁷J. P. Banerjee, R. Mandal, D. S. Banerjee, S. Thutupalli, and M. Rao, “Unjamming and emergent nonreciprocity in active ploughing through a compressible viscoelastic fluid,” Nature Communications 13, 4533 (2022).DOI
- [28]Fins, A., Wysocki, A., Ybert, C., Leocmach, M., & Rieger, H. (2024). How to steer active colloids up a vertical wall. Nature Communications, 15, 1710.
- [29]Li, J., Schlya, T. L., Liu, W., Shum, H., Rozen, I., Balazs, A. C., & Wang, J. (2015). Self-Propelled Nanomotors Autonomously Seek and Repair Cracks. Nano Letters, 15, 7077–7085.
- [30]³⁰D. Allan, T. Caswell, N. Keim, and C. van der Wel, “trackpy: Trackpy v0.3.2,” (2016).soft-matter.github.io/trackpy/v0.3.2
- [31]Piazza, R., Bellini, T., & Di Gregorio, V. (1993). Equilibrium sedimentation profiles of screened charged colloids: A test of the hard-sphere equation of state. Physical Review Letters, 71, 4267–4270.
- [32]³²W. W. Salsburg and W. W. Wood, “Equation of State of Classical Hard Spheres at High Density,” The Journal of Chemical Physics 37, 798–804 (1962).DOI
- [33]³³A. P. Solon, M. E. Cates, and J. Tailleur, “Active brownian particles and run-and-tumble particles: A comparative study,” The European Physical Journal Special Topics 224, 1231–1262 (2015).DOI
- [34]³⁴I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert, and L. Bocquet, “Dynamic clustering in active colloidal suspensions with chemical signaling,” Physical review letters 108, 268303 (2012).pubmed.ncbi.nlm.nih.gov/23005020
- [35]³⁵F. Ginot, I. Theurkauff, D. Levis, C. Ybert, L. Bocquet, L. Berthier, and C. Cottin-Bizonne, “Nonequilibrium equation of state in suspensions of active colloids,” Physical Review X 5, 011004 (2015).DOI
- [36]Klonvessa, C., Ybert, C., Cottin-Bizonne, T., Kawasaki, M., & Lecomach, A. (2022). Aging or dead: Origin of the non-monotonic response to weak self-propulsion in active glasses. The Journal of Chemical Physics, 156.
- [37]³⁷E. Helfand, H. L. Frisch, and J. L. Lebowitz, “Theory of the Two- and One-Dimensional Rigid Sphere Fluids,” The Journal of Chemical Physics 34, 1037–1042 (1961).DOI
- [38]³⁸A. Poncet, O. Bénichou, V. Démery, and G. Oshanin, “Universal long ranged correlations in driven binary mixtures,” Physical review letters 118, 118002 (2017).DOI
- [39]³⁹M. E. Cates and J. Tailleur, “Motility-induced phase separation,” Annu. Rev. Condens. Matter Phys. 6, 219–244 (2015).DOI
- [40]⁴⁰A. P. Solon, Y. Fily, A. Baskaran, M. E. Cates, Y. Kafri, M. Kardar, and J. Tailleur, “Pressure is not a state function for generic active fluids,” Nature Physics 11, 673–678 (2015).DOI
- [41]Hagen, B. T., Taeleen, S. van ., & Löwen, H. (2011). Brownian motion of a self-propelled particle. Journal of Physics: Condensed Matter, 23, 194119. arXiv:1005.1343
- [42]⁴²T. Yanagishima, Y. Liu, H. Tanaka, and R. P. A. Dullens, “Particle-Level Visualization of Hydrodynamic and Frictional Couplings in Dense Suspensions of Spherical Colloids,” Physical Review X 11, 021056 (2021).DOI
- [43]⁴³J. D. Hunter, “Matplotlib: A 2d graphics environment,” Computing in Science & Engineering 9, 90–95 (2007).DOI