Toward a flow-dependent phase-stability criterion: Osmotic pressure in sticky flowing suspensions

Sticky colloidal suspensions are prevalent in biological, pharmaceutical, industrial, and personal-care fluids, in materials as diverse as amyloid aggregates, cells, monoclonal antibodies, adhesives, detergents, and many others. Common among these examples is the suspension of microscopic particles, micrometers or less in size, in a continuum liquid, where Brownian diffusion competes with attractive interparticle forces induced by various physico-chemical phenomena that can drive the formation of aggregates or clusters, as well as phase separation. Prediction of such aggregation, phase separation, and crystallization via equilibrium theories is routine; connections to nonequilibrium phase behavior such as gelation and vitrification are also well-known.^1–16 Colloidal systems can be engineered to enhance or avoid such stability—for example, resisting phase separation is a key goal in the design of monoclonal antibodies.^17,18

The presence of flow induced by imposed shear or gravitational forces can lead to further changes in suspension stability: for example, the creaming action of activated particles in shaving foam,¹⁹ the interplay between sedimentation and gelation,^20–23 and the injection of antibodies^17,18 that involves motion relative to a background suspension that can help mix the product or trigger aggregation or gelation. Mechanistically, the suspension structure and its stability depend on an interplay between colloidal-scale attractive forces, Brownian motion, and the flow or external body forces that cause relative particle motion. Many prior experimental and computational studies of attractive suspensions have focused on equilibrium behavior such as phase separation, gelation in the absence of flow,^7,8,10,24,25 or the mechanical response to flow or forces of large fully formed gel networks.^26–28 However, in many cases, flow and attractive forces emerge simultaneously, such as in gel collapse under gravitational forces^{20–23,26,29–31} or phase separation during flow,^32–35 but there are few studies on how flow itself moves the region of stability. A recent study of the yield in reversibly bonded colloidal gels recasts the mechanical yield under weak stress as nonequilibrium phase separation re-initiated by a release from kinetic arrest,^27,35 where negative osmotic pressure plays a key role.^26,36 This connection between phase and flow in attractive suspensions motivates the present theoretical study of the microstructure and rheology of attractive colloidal suspensions with a focus on the interplay between externally imposed forces, attractions, and Brownian motion.

By analogy to molecular systems, the phase envelope of solvent-suspended colloids has been well-established by the second virial coefficient^7,8,37,38 B₂, which describes the effect of pair-level interactions on the osmotic pressure in a virial expansion.³⁹ Attractive forces decrease B₂, reflecting a negative contribution to osmotic pressure. One key outcome is that the connection between a macroscopic thermodynamic variable, the osmotic pressure, and a pair-level interaction strength, B₂, sets the onset of equilibrium phase separation. The threshold value of B₂ is given by the Vliegenthart–Lekkerkerker (VL) criterion,⁸ which establishes the colloidal gas–liquid critical point at B₂ = −6v_p (where v_p is the particle volume). The semi-empirical VL criterion serves as a remarkably consistent predictor of colloidal gas–liquid phase separation across a wide variety of short-ranged interaction potentials.^7,40,41 The second virial coefficient has also shown utility as a model-agnostic predictor of linear-response rheology when hydrodynamic interactions (HI) are weak.⁴² The structural distortion induced by weak flow and the viscosity in the Newtonian plateau are quantitatively predictable through B₂ irrespective of the detailed functional form of the interparticle attraction.⁴² 

However, nonlinear rheology arising from stronger flows cannot be predicted by the second virial coefficient because B₂ encodes only the equilibrium structure, while strong flow shifts the balance of microscopic forces and displaces the structure from equilibrium.⁴³ Structural instability in suspensions is known to arise in the presence of flow, where relative particle motion can tip the balance of attractive, repulsive, Brownian, and other microscopic forces. For example, flow can induce flocculation when shear or extensional flow brings particles into contact, where attractive forces form durably bonded aggregates.^44,45 More recently, gravity has been shown to trigger the re-emergence of phase separation in reversibly bonded colloidal gels, not by driving sedimentation or network rupture but by rupturing a tiny fraction of interparticle bonds that then permits a cascade of bond relaxation throughout the gel, driven by negative osmotic pressure that, in turn, drives condensation.²⁶ In addition to attractive and Brownian forces, hydrodynamic interactions influence relative particle motion and, by extension, the interplay between flow, microstructure, and rheology.^46–54 The tendency of hydrodynamic interactions to couple particle motion and stabilize particle trajectories can also affect the stability of flowing suspensions.^44,55–57

The phase stability of attractive suspensions is influenced by flow strength and hydrodynamic interactions; although the second virial coefficient B₂ can predict the onset of thermodynamic colloidal phase instability via the VL criterion,⁸ it cannot predict phase instability triggered by flow. We recognize the potential of expanding this criterion to flowing, nonequilibrium systems, and in the present study, we aim to determine flow-induced changes in the osmotic pressure that could trigger particle aggregation and possible phase separation. We hypothesize that such an expanded criterion could be useful in predicting phase instability in flowing, attractive suspensions. The connection of B₂ to osmotic pressure provides a guide to a new flow-dependent phase-stability criterion: whereas the second virial coefficient accounts for the equilibrium effects of attractions on osmotic pressure, the nonequilibrium osmotic pressure can be expressed via statistical mechanics as a function of the detailed particle microstructure,^51,53,54,58 automatically encoding the effects of flow and any microscopic forces represented in the Smoluchowski equation governing microstructural evolution. The microstructural perspective is a composite view of many particle trajectories throughout the suspension; indeed, suspension stability is underpinned by the stability of particle trajectories. For example, sedimentation can lead to channel formation²² and disruption of gelation,^20,21 and the relative motion driven by motor proteins in biological cells can stir cytoplasm.^59,60

In the present study, we aim to develop a nonequilibrium criterion of phase instability analogous to the second virial coefficient by identifying the combination of flow strength and interparticle attraction strength that leads to negative osmotic pressure, the hallmark of condensation. To do so, we build a simple model that directly connects particle trajectories to suspension-averaged osmotic pressure in which we can theoretically analyze the individual and cooperative influences of interparticle attractions, Brownian motion, and hydrodynamic interactions on the average osmotic pressure in a colloidal suspension, as it evolves with changes in flow strength. The microrheological model system is an excellent framework for studying detailed particle motion and its impact on rheology; thus, we utilize and expand active microrheology theory of hydrodynamically interacting colloids to include attractive forces and study relative particle motion, structure, and rheology that pertain to phase behavior.

In active microrheology, a colloidal probe or set of probes is driven by an external force such as a gravitational or magnetic field through a complex fluid. The relative motion between the probe and background bath particles is monitored to infer material and flow properties. Probe motion fluctuates as it encounters microstructural features of the bath, where its mean motion reveals suspension viscosity,^47,61,62 its fluctuations reveal flow-induced diffusion,^63–66 and the combination leads to suspension stress.^53,54,67 Batchelor established the fundamental theory that connects the statistical mechanics of colloidal particle motion to Stokes-flow hydrodynamics and predicted that Brownian motion destroys Stokes flow reversibility to produce non-Newtonian rheology,⁴⁶ which later studies confirmed.^51,62 In combination with nonequilibrium statistical mechanics, Batchelor’s theory connects probe motion to material and flow properties such as the viscosity and osmotic pressure. This technique is an important complement to traditional macroscopic shear rheology, where microrheology permits for the interrogation of nanoscopic quantities of fluid such as droplets, rare biological fluids, or the interior of biological cells.^68–73 Prior theoretical models have established viscosity, stress, diffusion, and osmotic pressure over a wide range of probe forcing strengths for purely repulsive hard spheres, both with and without hydrodynamic interactions.^{53,54,61–63,67,74–76}

However, in many material systems, there are physico-chemical conditions that produce attractive forces between the background particles or between a probe and the bath particles, which influence both phase behavior and rheology. The influence of these microscopic forces on equilibrium phase behavior is well-established in theoretical models,^8,24,25 and statistical–mechanical models of flow rheology are robust for flowing suspensions with hydrodynamic, frictional, or repulsive forces.^{51,53,54,61,62,75–77} By contrast, the rheology of attractive, flowing suspensions, both for shear flow and for microrheology, is less well understood, especially in the presence of hydrodynamic interactions.^{42,43,78–82}

In our recent work, we constructed the basic framework for how attractive forces alter the interaction between the probe and bath particles in a freely draining suspension, and we identified the impact of such attractions on the particle configuration and the viscosity for a range of external forces.^42,43 The consideration of both attractive forces and hydrodynamic interactions on suspension viscosity under stronger flows is a natural extension, which we consider in the present work. Attractive forces also directly affect osmotic pressure and can thus trigger phase instability, but to our knowledge, there are no studies of the impact of attractions on hydrodynamically interacting colloids on nonequilibrium osmotic pressure in active microrheology. In the present study, we use active microrheology to investigate this interplay between attractive forces, hydrodynamic interactions, Brownian motion, and external forces on suspension structure and how these influence the nonequilibrium osmotic pressure as it evolves with flow strength.

The remainder of this paper is organized as follows: In Sec. II, we describe the model colloidal system and our method for tuning the strength of external forcing and interparticle forces, both attractive and repulsive. We derive an equation for equilibrium and nonequilibrium microstructure, relating it to mean probe motion, microviscosity, and nonequilibrium osmotic pressure in Sec. III. We present our results for microstructure, viscosity, and nonequilibrium osmotic pressure in Sec. IV. Finally, in Sec. V, we discuss the applications of our results and possible extensions to our model.

We consider a set of N rigid, neutrally buoyant, spherical Brownian particles, all of hydrodynamic size a_h, suspended in a volume V_s of a Newtonian fluid of density ρ and viscosity η. One of the bath particles (the “probe”) is driven by a fixed external force $Fext$ through the suspension, setting the suspension into motion. The strength of inertial forces relative to viscous stresses is small, as defined by the Reynolds number, Re ≡ ρUa_h/η, where U is the characteristic probe velocity and the particles are small, Re ≪ 1. Fluid mechanics are thus governed by the Stokes equations and particle motion by low-Reynolds number hydrodynamics. The number density of probes, n_a, is small relative to that of the bath particles, n_b. Probe motion distorts the configuration of surrounding bath particles from their equilibrium configuration, producing an entropic force that resists probe motion as the Brownian motion of the bath particles acts to recover their equilibrium microstructure. Configuration-dependent colloidal-scale interparticle forces also influence probe motion, including repulsive forces arising from surface roughness, electrostatic repulsion, and surface-adherent polymer brushes, as well as attractive forces from, e.g., van der Waals forces, bridging polymer chains, and depletant particles. Repulsive forces can alter the strength of hydrodynamic forces, and a concomitant interplay with Brownian motion and flow influences rheology.

The excluded-shell model is a simplified means of representing a no-flux thermodynamic surface surrounding the no-slip hydrodynamic surface of a particle⁸³ (Fig. 1). The thermodynamic radius a represents repulsion that extends beyond the no-slip surface of the particle and can weaken the influence of hydrodynamic interactions by keeping particles far apart. Similarly, repulsions can also alter range-dependent attractive forces. We represent repulsion between particles via a conservative interparticle force derivable from a spherically symmetric potential that varies only with the relative position of an interacting pair. Since the absolute position does not matter in the absence of external fields or system boundaries, we place the probe at a position x₁ and a bath particle (particle 2) at x₂, giving a relative configuration r = x₂ − x₁ (Fig. 1). The interparticle force exerted on a bath particle by the probe is then given by $F2P=−∇rV(r)$⁠, where r points from the center of particle 1 toward the center of particle 2. A force of equal strength acts in the opposite direction on particle 1. For hard spheres with no other repulsion other than steric exclusion, the interparticle force is zero away from contact; it acts only repulsively and only at contact at the thermodynamic surface, where an infinitely strong repulsion bars overlap; the repulsive potential V^rep is then given by

Here, k is Boltzmann’s constant, T is the absolute temperature, and the center-to-center distance r is the magnitude of r.

Repulsions affect the strength of hydrodynamic forces between two particles, which scale as 1/r at wide separations and inversely as the gap thickness between a closely interacting pair. This influence is parameterized by the dimensionless repulsion range, κ,

When κ → 0, particles can approach one another closely enough to experience lubrication interactions and strong hydrodynamic coupling, and when κ → ∞, long-ranged interparticle repulsions prevent the no-slip surfaces from approaching closely, and hydrodynamic interactions are negligible. Although the excluded-shell model is a simple approximation of interparticle repulsion, it has been employed to great effect in recovering the important qualitative behavior of colloidal dispersions in both shear macrorheology and probe microrheology.

Colloidal particles may also exert attractive forces on one another, ranging from van der Waals forces and Coulombic interactions between charged colloids, which can be modeled by the DLVO potential, and colloidal-scale reversible attractions arising from adhering or non-adhering polymer chains. Depletion interactions caused by suspended polymers are frequently encountered in biophysical systems and industrial applications; the simplest model for such phenomena, the Asakura–Oosawa potential, is often employed to model depletion attractions.⁸⁴ 

In the present study, we model a radially symmetric attractive potential, where we make the assumption that strength and range can be tuned separately. The attractive force can be made dimensionless by scaling it relative to thermal fluctuations: F^attr ∼ 2kT/a. We model an AO potential V^attr with attraction range 2a ≤ r ≤ 2a(1 + Δ) (Fig. 2), where Δ is the range of the attractive force beyond the thermodynamic surface relative to the particle’s thermodynamic size a,

where the attraction strength at contact is made dimensionless as A ≡ F^attr(r = 2a)/(2kT/a). The interparticle force arising from this potential exerts an equal and opposite force on an interacting pair. The force exerted on particle 2 (a bath particle) by particle 1 (the probe) is $F2attr=−∇rVattr$⁠, and the force exerted on particle 1 by particle 2 is $F1attr=+∇rVattr$⁠. Given Eq. (3), the attractive force $F1attr$ within the AO attraction range is

Here, $r̂$ is the unit vector pointing from the probe toward the bath particle, along the line of centers. The AO potential used in our study is shown in Fig. 2 alongside the square-well and excluded-shell potentials for comparison. The range of attraction is arbitrary; in the present study, we focus on short-ranged attractions, where the attraction range is smaller than the particle size, Δ < 1.

As the probe travels through the bath, it interacts with other particles that exert short-range attractions and hard-sphere repulsion. These encounters alter the probe velocity, from which one can infer suspension viscosity via application of Stokes’s drag law.⁶¹ The probe attains its mean speed over many interactions with other particles, i.e., after sufficient sampling of bath configurations. The probe, in turn, distorts the bath particle microstructure, which induces nonequilibrium viscosity and stress. In this section, we derive a pair-level Smoluchowski equation governing microstructural evolution followed by a formulation of the microviscosity and the particle-phase stress from which we obtain the nonequilibrium osmotic pressure.

We begin by restricting our scope to the semi-dilute limit, where pair-level interactions influence the microstructure—the microscopic arrangement of particles—but three-body and higher order interactions may be neglected. In the semi-dilute limit, the steady-state particle configuration can be described by the probability function P₂(x₁, x₂) of finding the probe at a position x₁ and any other bath particle at a position x₂. This probability function obeys a Smoluchowski equation, which in the pair limit reads

where the temporal evolution of particle distribution ∂P₂/∂t is balanced by the sum of all particle fluxes, $∑α=12jα=∑α=12UαP2$⁠, where U_α is the velocity vector of particle α. Here, j₁ = U₁P₂ denotes the probability flux of the probe and j₂ = U₂P₂ denotes that of a bath particle, and gradients are taken with respect to the absolute positions of the particles. The particle flux comprises contributions to the particle velocity from external (⁠$F1ext$⁠, acting only on the probe), interparticle (⁠$FαP$⁠), and entropic (Brownian) forces [$FαB≡−kT∇αlnP2(x1,x2)$], where particle motion and forces are coupled by the hydrodynamic mobility $MαβUF$⁠, which couples the velocity of particle α to the force exerted on particle β. The particle fluxes can thus be expressed as

Here, we have made use of the fact that the mobility tensor is related to the translational diffusion tensor by $Dαβ≡kTMαβUF$ as required by the fluctuation–dissipation theorem.⁸⁵ The mobility tensor depends on the particles’ hydrodynamic size and the configuration of particle pairs, namely, their center-to-center separation, and is given by

where $xαβA(r)$ and $yαβA(r)$ are the configuration-dependent hydrodynamic mobility functions and I is the identity tensor. Because the absolute position does not matter, we adopt a relative coordinate system with the probe centered at z = x₁ and a bath particle at r = x₂ − x₁. The probability density function can then be related to the pair-distribution function g(r),

Here, n_b is the bath particle number density, P₁(z) denotes the probability of finding a probe at a given absolute position z, and P_1/1(r|z) denotes the probability of finding a bath particle at relative position r, given a probe position z. The interparticle forces may be expressed as gradients of an interparticle potential, $FαP=−∇αV$⁠. This potential includes both attractive and repulsive contributions, V ≡ V^rep + V^attr, both of which we described in Sec. II. At the steady state, ∂g/∂t = 0, we may integrate over all positions of the probe, z, to obtain

where all gradients are taken with respect to relative position, r, and the relative Brownian diffusivity is defined as

The hydrodynamic functions $G(r/ah)=x11a−x12a$ and $H(r/ah)=y11a−y12a$ describe relative mobility along and transverse to the line of centers, respectively, and have analytical forms.^48,86,87 The Smoluchowski equation can be made dimensionless by scaling quantities as

When the probe and bath particles are of the same size, the relative mobility can be expressed in terms of the relative diffusivity as M₁₁ − M₁₂ = −D_r/2kT. With these scalings, the governing equation and the boundary conditions—rigid spheres permit no relative flux at contact, j₁ − j₂ = 0, and in a semi-dilute suspension, there is no long-range order—can be expressed in the dimensionless form as

where the subscripts on the gradients with respect to r have been omitted for brevity and the Péclet number, Pe = F^ext/(2kT/a), emerges naturally from the scalings in Eq. (11). The unit vector $F̂ext≡Fext/Fext$ points along the applied force in $−êz$⁠.

External forcing breaks the symmetry of the surrounding bath, inducing structural distortion that varies with relative position r, attraction strength A, and forcing strength Pe. We express the total, distorted microstructure, g(r; Pe, A), in terms of its displacement from the equilibrium microstructure g^eq(r; A), which is radially symmetric and Boltzmann-distributed: g^eq = exp(−V/kT). The total distorted microstructure can be expressed as a sum of the equilibrium microstructure plus the distortion f(r; Pe, A) as

Equations (12) and (13) can be solved to obtain the microstructure resulting from any strength of external forcing and attractive interactions and for conditions where the strength of hydrodynamic interactions can be very strong, very weak, or anywhere in between. Analytical solutions can be obtained in limiting conditions for each of these, and numerical solutions are employed here to obtain solutions over the full range of attractions, hydrodynamic interactions, and flow strength (the details are provided in  Appendix D). A centered finite-difference scheme is employed, where gradients are discretized in spherical coordinates to obtain a sparse banded matrix, which was solved in Matlab using a LaPack banded direct solver to obtain the distorted microstructure. When flow is strong, an upstream boundary layer forms, which produces very steep density gradients; these are faithfully represented by a dense placement of radial grid points within the boundary layer that scales with Pe. Similarly, strong attractions induce density gradients within the attraction range, and the density of grid points may be further increased with the scaling in dimensionless attraction strength, A.

Stokes’s drag law sets the velocity U⁰ of an isolated spherical probe of size a_h forced at the low-Reynolds number through a Newtonian solvent of viscosity η; for a force $Fext$⁠,

from which the solvent viscosity can be deduced. The effective viscosity of a suspension of colloids can similarly be deduced by application of Stokes’s drag law, where the force probe encounters many configurations of colloidal particles as it travels through the bath, leading to a well-defined mean probe velocity,⁶¹ ⟨U₁⟩,

where the angle brackets around $U1$ denote an average of U₁ over all configurations weighted by the probability distribution function. As mentioned above, the changes in the probe velocity that cause its mean speed to differ from its Stokes velocity arise from interparticle forces, Brownian forces, and hydrodynamic interactions, as previously described by Khair and Brady⁶² and, subsequently, Mohanty and Zia,⁷⁵ 

where $⟨U1H⟩$⁠, $⟨U1B̄⟩$⁠, and $⟨U1P⟩$ signify the volume-averaged contributions of hydrodynamic, Brownian, and interparticle forces to the mean probe motion, respectively, averaged over all configurations. The overbar on $⟨U1B̄⟩$ signifies that the Brownian velocity is a time-averaged velocity; it indicates that the configuration-dependent Brownian velocity is an average over timescales that are much longer than the solvent timescale. The angle brackets signify an ensemble average velocity, where the velocity of the probe through a given configuration is weighted by the probability of occurrence of that configuration as

The hydrodynamic mobility functions $x11A(κr+r)$⁠, $y11A(κr+r)$⁠, G(κr + r), and H(κr + r) were described previously in Sec. III A. The interparticle force F^P can comprise both repulsions and attractions, as described in Sec. II.

The effective viscosity is inferred as the change in the mean probe speed, obtained in experiments by monitoring speed over long times as the probe encounters many representative surrounding particle configurations; it can be predicted theoretically from Eq. (15) from the ensemble-averaged velocity,

The effective viscosity comprises two contributions: solvent drag plus the microviscosity—the drag exerted by the microstructure. As with the velocity in Eq. (16), the viscous drag of the microstructure comprises hydrodynamic, Brownian, and interparticle contributions,

where we will drop the superscript micro from the individual components for clarity. To leading order in bath particle volume fraction ϕ_b ≪ 1, the hydrodynamic, Brownian, and interparticle contributions to the viscosity (scaled on the solvent viscosity) are given by

Equation (22) is valid for any spherically symmetric interaction potential, whether attractive or repulsive.

The osmotic pressure in colloidal suspensions is very well described, at equilibrium, by the same framework used to describe osmotic pressure in molecular fluids.^53,58,67 At infinite dilution, the ideal osmotic pressure is nkT because particle interactions are negligible. Interactions between colloidal particles give rise to departures from ideal osmotic pressure, which, at equilibrium, can be expressed as a virial expansion in the number density,³⁹ 

where Π^eq denotes the equilibrium osmotic pressure and n is the particle number density. The sign of B₂ reflects the tendency of a system to expand or contract, reflecting the influence of repulsive or attractive interparticle forces, respectively. This sets equilibrium phase behavior such as the gas–liquid critical point.^7,8,88 The second virial coefficient is a function of the interparticle pair potential V(r),

For hard spheres, an analytical expression may be obtained: $B2HS=16πa3/3$⁠, reflecting the positive osmotic pressure that is entropic in origin, arising from the particles’ excluded volume. The hard-sphere exclusion provides a natural scaling that produces the reduced second virial coefficient,

Physically, because this scaled quantity can be positive or negative, it reveals how attractions and repulsions decrease or increase the effective volume of particles in comparison to purely repulsive hard spheres. Interparticle potentials increase $B2*$ if the potential is positive (repulsions) and decrease it if the potential is negative (attractions); the sign of $B2*$ indicates the tendency of a system to expand or contract at equilibrium.

At equilibrium, the osmotic pressure Π^eq is fully described by Eq. (23), where the second virial coefficient [Eq. (24)] encapsulates the effects of pair interactions in the equilibrium structure, g^eq = exp[−V(r)/kT]. However, imposed flow alters the relative strengths of microscopic forces, which induces structural asymmetry that, in turn, alters the osmotic pressure. This nonequilibrium osmotic pressure and asymmetry cannot be described by Eqs. (23) and (24) and instead requires nonequilibrium statistical mechanics, where the osmotic pressure is first obtained mechanically from the trace of the suspension stress,

The particle-phase stress, ⟨Σ⟩, can be divided into nonhydrodynamic and hydrodynamic contributions as follows:

where I is the isotropic unit tensor and the angle brackets denote an ensemble average over all positions of the bath particles. We follow the approach of Zia and co-workers,^53,54,67 considering a material through which there is a large number of probe particles moving in response to an external force. The probes are taken to be so dilute that they do not interact with one another (ϕ_a ≪ 1), and they are dilute relative to the bath particles, n_a ≪ n_b. The bath particles are also dilute (ϕ_b ≪ 1); moreover, they are force-free. As a result, bath–bath interactions are also negligible, and the bath particles only experience stresses when encountered locally by the probe particles. As a result, the particle-phase suspension stress arises entirely from the probe phase and its interactions with nearby bath particles. A detailed discussion of the averaging, following Batchelor’s method,^46,89 is given in Appendix A of the work of Chu and Zia.⁶⁷ The first term on the right-hand side of Eq. (27) is the ideal osmotic pressure arising from the thermal motion of bath particles. The second term $narFP$ is the elastic interparticle stress.^26,53,67 It comprises both repulsive and attractive interparticle forces experienced between the probe and nearby bath particles. The third term on the right-hand side of Eq. (27), $ΣH$⁠, accounts for hydrodynamic stresses arising from particle motion. Such motion can be induced by externally imposed flow, by Brownian disturbance flows, and by interparticle forces, the last of these that produces a dissipative interparticle stress,

Each of these average stress contributions—external force-induced hydrodynamic stress $ΣH,ext$⁠, Brownian drift-induced hydrodynamic stress $ΣB$⁠, and dissipative interparticle stress $ΣP,dis$—is given by the probe number density times the volume-averaged stresslet: $Σ=naS$⁠. The stresslet is the symmetric part of the stress tensor, and it may also be expressed as a sum of these contributions,

The linearity of Stokes flow connects the stresslet to particle motion via the hydrodynamic resistance and mobility tensors,^47,86,90

where $Uext$ and Ω^ext, respectively, denote particle translation and rotation due to external forces and U^P and Ω^ext denote translation and rotation due to interparticle forces. The hydrodynamic mobility tensors M_UF and M_ΩF couple forces (whether external or interparticle) to translational and rotational velocities, respectively, and the hydrodynamic resistance tensors R_SU and R_SΩ couple particle motion to surface tractions. A brief discussion of the hydrodynamic resistance and mobility functions and their contributions to the stress tensor $Σ$ is given in  Appendix A. In Eq. (30), the external force-induced hydrodynamic stresslet S^H,ext and the dissipative interparticle stresslet S^P,dis are expressed in terms of particle motion induced by external forces (U^H,ext) and interparticle forces (U^P,dis); the Brownian drift stresslet S^B arises due to disturbance flows caused by the Brownian motion of particles diffusing toward regions of higher particle mobility. This approach to defining the stress follows the work of Chu and Zia;^54,67 it disambiguates prior definitions in the literature by recognizing that the hydrodynamic stress comprises not only external force-induced stress but also contributions from Brownian drift and viscous dissipation from interparticle forces. Additionally, the inclusion of both elastic and dissipative interparticle stresses within the total interparticle stress recovers the net zero contribution to stress from hard-sphere interparticle forces in the pure-hydrodynamic limit of κ → 0, where lubrication forces replace the role of hard-sphere repulsions.⁶⁷ 

The average stresses denoted above encode the interplay between flow and microscopic forces, which influences the osmotic pressure both directly and indirectly. Flow and particle interactions set the osmotic pressure directly through the configuration-dependent microscopic forces (mobility, resistance, and interparticle potential); flow and microscopic forces determine the probability of a given microstructural configuration g(r; Pe, A) that produces such interactions, and this weighting gives the indirect effects of flow and interactions. We derive the statistical-mechanics expressions for each of the stresses as the weighted influence of each configuration on the relevant forces,

The scalar hydrodynamic functions $XαβP$⁠, $XαβG$⁠, $YαβG$⁠, and $YαβH$ are elements of the resistance tensor, relating hydrodynamic stress on the surface of particle α induced by the motion of particle β. The scalar mobility functions $xαβA$⁠, $yαβA$⁠, and $yαβB$ are elements of the hydrodynamic mobility tensor,^48,87,91 relating the velocity of particle α induced by a force on particle β. In Eqs. (31)–(33), these hydrodynamic functions are evaluated at r(κ + 1); a brief description of these functions is given in  Appendix A. The flow-induced osmotic pressure contributions derive directly from the trace of stress, following Eq. (26),

At equilibrium, these contributions are exactly zero; the equilibrium contribution of $4B2*$ has been subtracted from the trace of Eq. (34) to give the nonequilibrium (elastic) interparticle osmotic pressure $ΠneqP,el$⁠. The resulting nonequilibrium contribution to the osmotic pressure can be positive—a tendency for the particle phase to expand in space—or negative—a tendency for the particle phase to condense.

We solved the equations formulated in Sec. III governing the distribution of colloids in space in the presence of attractive forces of varying strengths and strong hydrodynamic interactions. The solutions for these equations are presented and discussed here with a focus on the interplay between hydrodynamics and attractive forces on structural asymmetry. This is followed by a section on the non-Newtonian viscosity of the flowing, sticky suspension. In Sec. IV C, we examine the nonequilibrium osmotic pressure and place it in the context of stability, where we construct a phase map centered on the osmotic pressure and how it can lead to phase separation.

The microscopic balance of forces between particles in suspension changes with flow, altering the microstructural configuration that influences viscosity. In a purely repulsive suspension under weak probe forcing (the linear-response regime), the structural distortion in active microrheology is a diffusive dipole,⁴⁷ regardless of the strength of hydrodynamic interactions,^61,62 with symmetric particle accumulation and depletion regions on the upstream and downstream faces of the probe, respectively. In the current study, we extend our prior work in freely draining, attractive suspensions^42,43 to solve the Smoluchowski equation governing this dipolar structure at arbitrary attraction strengths, now generalized to arbitrary strengths of hydrodynamic interactions (the details are given in  Appendix B), to obtain distortion, f(r), which is plotted in Fig. 3. The purely repulsive structure is shown for the freely draining and strong-hydrodynamic limits in Figs. 3(a) and 3(f), respectively. Brownian motion tends to separate particles, either randomly in a freely draining suspension or via Brownian drift with strong hydrodynamic interactions. Meanwhile, advection acts to sweep particles downstream as they pass the probe. Attractive interparticle forces hinder diffusion and thus condense the dipole closer to the probe as attraction strength grows [Figs. 3(b)–3(d)]. Beyond a threshold attraction strength (⁠$B2*=−0.5$ for freely draining systems), structural condensation gives way to a structural reversal regime: strong attractions [Fig. 3(e)] reverse the accumulation and depletion regions because passing particles bond to the probe, and they rotate as doublets under the advective flow. Doublet rotation transfers particle density to the downstream region.⁷⁸ At a short range, hydrodynamic interactions enhance this downstream transfer because lubrication forces make pair couplings more durable [Figs. 3(f)–3(j)]. Overall, lubrication forces condense the distorted structure and shift structural reversal to lower attraction strengths. Brownian motion acts in opposition to hydrodynamic forces by destroying structural symmetry; we turn next to the effects of stronger forcing and greater structural asymmetry.

When flow is stronger than Brownian motion (Pe ≥ 1), structural asymmetry becomes more pronounced. We can examine dual limits: weak vs strong hydrodynamic interactions and repulsive vs strongly attractive interactions, represented in Figs. 4(a), 4(f), 4(e), and 4(j). Looking only at these extremes, one observes qualitatively identical behavior: when attractions are negligibly weak [Figs. 4(a) and 4(f)], the particles accumulate upstream of the probe and a trailing depletion region is formed because Brownian motion is too weak to heal the structural disturbance of the probe.^61,62 When attractions are strong [Figs. 4(e) and 4(j)], doublet formation leads to complete structural reversal: particle density is swept from the upstream region into a trailing accumulation zone downstream. In between these two limits, the onset of structural reversal manifests first in the far-field region [Fig. 4(d)] and then, at higher attractions, also close to contact in the freely draining limit [Fig. 4(e)].

The opposite is true when hydrodynamics are strong: the structure is first drained close to contact as attractions grow but remain moderate [Fig. 4(h)]. Once attractions are strong, the structural asymmetry propagates to the far field [Fig. 4(j)]. Interestingly, when hydrodynamics are strong, there is an interim structure of moderate attractions [Fig. 4(g)] where particles accumulate all over the entire probe surface, posing an interesting possibility for macroscopic condensation and phase destabilization. Mechanistically, hydrodynamic interactions hinder relative motion; as Brownian motion weakens under strengthened flow, the probe encounters a greater density of bath particles, which are swept downstream by advection and retained by hydrodynamically enhanced bond durability, shifting the reversal point to lower attraction strength. Overall, moderate advection shifts structural reversal to lower attraction strength by assisting doublet rotation, and hydrodynamics increase bond durability to enhance this effect. The natural next step is to identify the flow strength that can disrupt these bond dynamics.

We now turn our attention to the advection-induced changes in the balance of microscopic forces. Under strong forcing (Pe > 10), an upstream boundary layer and trailing wake form in purely repulsive colloidal suspensions^61,62 [Figs. 5(a)–5(c)]. Stronger flow compresses the boundary layer close to the probe surface, and hydrodynamic interactions delay downstream detachment of the boundary layer, narrowing the wake [Figs. 5(b) and 5(c)]. Strong attractions combine with strong advection to drain particle density away from the upstream region (forming a depletion boundary layer) and shift them downstream into the wake [Figs. 5(d)–5(f)]. The shape of this wake is set by the competition between advection and attractions: if attractions are strong, the wake becomes an accumulation wake, but if advection can overpower attractive bonds, the accumulation wake peels away into a trailing accumulation annulus, surrounding the recovered depletion wake.⁴³ Strong hydrodynamic interactions tip the balance in favor of attractions: by slowing relative motion, strong hydrodynamics enable attractive bonds to resist advective breakage at flow strengths more than double that required to form a depletion wake in freely draining suspensions [Fig. 5(f)]. This can also be understood from a trajectory-analysis perspective, shown in  Appendix F: strong hydrodynamic interactions dramatically lengthen the average time of particle encounters, and as bath particles are transferred downstream by advection, the increased duration of particle encounters, averaged over time, increases particle density in the wake. Additional contour plots of microstructure f(r) are shown in  Appendix E over a broad range of attraction and forcing strengths, illustrating the transitions between accumulation and depletion wakes.

When forcing grows even stronger (Pe ≥ 100), it compresses the upstream boundary layer to a small region within the attraction range. In this regime, three distinct distortion regions arise upstream. Far upstream, advection is dominant, and attractions do not affect the microstructure; particles accumulate because hydrodynamic entrainment slows their relative motion, increasing particle density within this region. Figure 6 shows the structural distortion along the upstream stagnation streamline f(θ = π) as it varies with surface-to-surface separation r − 2 at several different strengths of flow and attraction. At low Pe [Fig. 6(a)], structural distortion decays smoothly as surface-to-surface separation r − 2 increases, but stronger forcing [Figs. 6(b)–6(d)] pushes the distorted structure closer to the probe surface. In Fig. 6(d), Pe = 1000, in the green-shaded region beyond the attraction range (r − 2 > 0.2), all the curves collapse together, indicating that the upstream structural distortion f(θ = π) in the outer region is unaffected by attraction strength. Another region forms near the probe surface: an O(Pe⁻¹)-thin boundary layer develops where Brownian motion once again becomes important, competing with attractions and advection. Within this region, lubrication interactions dramatically slow relative motion and Brownian motion acts to disrupt doublet formation. As a result, particles spend more time in this region, and average particle density increases. In Fig. 6(b), the innermost region is shaded red; here, we find that upstream particle accumulation can occur even for suspensions with attraction strengths higher than the low-Pe reversal threshold, such as A = 20 (⁠$B2*=−1.37$⁠), meaning that if Brownian motion is significant relative to attractions, a very thin accumulation region will always form. The third region forms between the boundary layer and the outer region (i.e., within the range of attractions, but outside the boundary layer), where Brownian motion is insignificant, but attractions and advection are important, acting in tandem to rapidly bring particles closer to the probe surface. This rapid motion means particles spend less time within this region compared to purely repulsive particles; as a result, this middle region is drained of particle density, as seen in Fig. 6(d): structural distortion in the middle region along the upstream streamline is negative, f(θ = π) < 0, even for A = 10 (⁠$B2*=0.56$⁠), where attractions are still far too weak to reverse the linear-response distortion. That is, when advection is dominant but attractions are present, two accumulation regions envelop a depletion zone; growing particle accumulation within the boundary layer is expected to dominate the rheology.

The structure resulting from the interplay of advection, Brownian motion, and attractions sets rheology and influences suspension stability. We examine the consequences of hydrodynamically enhanced bond durability for viscosity in Sec. IV B.

The evolving microstructure studied in Sec. IV A sets rheology. In the linear-response regime, attractions increase the equilibrium particle accumulation, g^eq, which, in turn, is expected to increase the equilibrium hydrodynamic drag⁴² [Eq. (22)]. External forcing produces flow and distorts the microstructure, g(r), and, in tandem with a changing balance of microscopic forces during flow, is expected to affect the microviscosity. In this subsection, we present the impact of the combined effects of hydrodynamic interactions, attractive forces, Brownian motion, and flow on the non-Newtonian microviscosity. We will first examine how attractions alter the viscosity trends such as flow-thinning in the linear-response and low-Pe regimes familiar in freely draining suspensions. Then, we will examine the opposite limit of very strong flow, showing how attractive forces alter familiar hydrodynamic effects such as flow-thickening. Finally, we will study the intermediate regime where the interplay between attractions and hydrodynamic interactions gives rise to entirely new non-Newtonian rheology.

In our previous study of the linear-response microrheology of freely draining, attractive suspensions, we showed that the low-Pe Newtonian plateau is decreased by weak attractive forces owing to a shift from energy dissipation to energy storage in stretched bonds,⁴² until strong attractions drive the plateau back up by flipping the microstructure. In Fig. 7, which plots the microviscosity η^micro at various attraction strengths as it evolves with Pe, the combined influence of attractions and hydrodynamics on the linear-response microviscosity is evidenced by the set of Newtonian plateaus in the low-Péclet region of the plot. The monotonic increase in the plateau with attraction strength reveals that linear-response attraction-thinning is absent when hydrodynamic interactions are strong. Both structure and hydrodynamic mobility are responsible for this behavior, as evident from Eq. (22). The disappearance of attraction-thinning is explained as follows: although the reversal in distorted structure f₁(r) that we observe for freely draining, attractive suspensions also occurs with strong hydrodynamic interactions (cf. Sec. IV A, Fig. 3), the key difference is that hydrodynamics produce equilibrium drag, η^H,eq, on the equilibrium microstructure g^eq, driven by changes in the relative mobility. Because g^eq increases monotonically with attraction strength, so does η^H,eq and, consequently, the low-Pe Newtonian plateau in η^micro. In the weakly nonlinear regime, Pe ≲ 1 in Fig. 7, the familiar flow-thinning emerges for purely repulsive suspensions, but in the presence of weak attractions, this flow-thinning regime is suppressed and, with moderate attractions, is ultimately replaced by low-Pe flow-thickening, a mechanism distinct from lubrication-based flow-thickening.

Recognizing that stronger flow changes the balance between attractions, Brownian motion, and advection, we separate the individual contributions to the microviscosity and plot them in Fig. 8. The nonequilibrium hydrodynamic, Brownian, and attractive contributions to microviscosity, as defined in Eq. (22), at attraction strength A = 12 (⁠$B2*=0.37$⁠), are plotted in Fig. 8 as they evolve with flow strength. The corresponding contour plots of the microstructure at Pe = 0.1, 1, 10, 100, and 1000 are also shown in Fig. 8. The attractive viscosity contribution, η^attr, is initially negative and of magnitude nearly equal to the Brownian contribution η^B, showing how attractions weaken the Brownian drift and thus suppress flow-thinning. The curve for the hydrodynamic viscosity η^H,neq rises with increasing Pe (until Pe ∼ 10) because weakly nonlinear flow increases nonequilibrium particle density around the probe due to hydrodynamic strengthening of attractive bonds. As a result, hydrodynamic drag increases. This increase dominates non-Newtonian behavior up to Pe ≈ 10. Overall, when external forcing is weak to moderate, i.e., too weak to break attractive bonds, hydrodynamic drag tends to thicken attractive suspensions by both retaining particles through increased bond durability (a structural effect) and slowing relative motion (a direct effect from reduced mobility). Although attractions and hydrodynamic interactions can combine to thicken a suspension under weaker flows, attractive forces influence the high-Pe flow-thickening regime differently.

We return to Fig. 7(a) to examine the opposite limit of strong forcing, Pe ≫ 1. Without attractions, the familiar flow-thickening behavior in hydrodynamically interacting suspensions emerges near Pe ∼ 1: the total microviscosity increases as flow strength increases. For smooth particles, this behavior arises from lubrication forces between the probe and bath particles, which emerge when strong flow compresses the upstream boundary layer enough to produce many lubrication interactions between particles.^51,62 However, attractions can delay this onset by two orders of magnitude in Pe: when attraction strength is increased to A = 15, Fig. 7(a) shows that flow-thickening is suppressed until Pe > 300. Beyond this value, a small increase in η^H,neq is observable in Fig. 8(b), corresponding to flow-thickening in its initial stages. The microscopic origins of this delay are examined in Fig. 8(b). The curve for the nonequilibrium hydrodynamic viscosity actually becomes negative at Pe ≈ 10, continues to decrease until Pe ≈ 100, and then slowly grows less negative. This surprising negative contribution to the viscosity is driven by the boundary-layer dynamics reported in Fig. 6, where the red region nearest the probe surface is drained of particle density because attractions strongly suppress Brownian drift. Equation (22) shows that the hydrodynamic contribution to the microviscosity, η^H,neq, directly depends on the particle density in this region. Thus, as a result of the boundary-layer drainage, the hydrodynamic drag decreases. Indeed, particle accumulation is absent for very strong attractions, $f(r−2<2Δ;B2*≪0;θ=π)<0$⁠, permitting no lubrication interactions, and therefore, none of the dramatic slowing of relative motion that produces flow-thickening. When flow strength increases to Pe = 100, particle accumulation within the boundary layer reappears for A = 15 [Fig. 6(c)], but it is only when flow strength has compressed the boundary layer enough for diffusion to outcompete attractions within the boundary layer, as in Fig. 6(d), Pe = 1000, that lubrication interactions begin to slow relative motion and recover high-Pe flow-thickening. This flow-thickening regime is further delayed as attraction strength increases, e.g., A = 20 in Fig. 8(c), where little to no flow-thickening occurs even at Pe ∼ 1000. The rightward shift of the flow-thickening regime shown in Fig. 7(a) with increasing attraction strength reflects the increasingly strong boundary-layer compression required to accumulate particles near contact upstream—and counteract advective drainage of particle density in both the downstream wake and the attractive well beyond the diffusive boundary layer; cf. the drained “middle region” of Figs. 6(c) and 6(d). Only when boundary-layer accumulation is recovered can lubrication interactions slow relative motion to flow-thicken the suspension. In attractive suspensions, low-Pe and high-Pe flow-thickening regimes arise from entirely different structural phenomena; they are separated by an intermediate-flow-thinning regime, but the micromechanical origins of this flow-thinning regime are different from those in hard-sphere dispersions.

We now turn our attention to the intermediate-flow-thinning regime in the middle region of Fig. 7(b) that occurs in strongly attractive suspensions. The individual microscopic forces that produce this behavior are examined in Figs. 8(b) and 8(c), where the height of the low-Pe plateau is set by bond strength and as such is higher in Fig. 8(c) than in Fig. 8(b); bonds between the probe and nearby bath particles are too strong for Brownian motion to break. However, when flow becomes strong enough to break bonds, the downstream region splits into an accumulation annulus around a trailing depletion wake (contour plots beneath each viscosity plot), reducing the hydrodynamic viscosity below its equilibrium value, η^H,neq < 0, and producing hydrodynamic flow-thinning. Mechanistically, flow becomes strong enough to tear bath particles away from the probe’s downstream surface, draining away the particle density that had accumulated near equilibrium, g^eq. However, bonds stretch before they break, and this leads to an attractive flow-thinning regime that precedes bond breakage: the attractive viscosity η^attr decays under moderate forcing, 1 < Pe < 10, even without wake-splitting because flow stretches bonds downstream and thins out the accumulation wake, leading to overall flow-thinning in η^micro. For strong attractions such as A = 20 [Fig. 8(c)], this bond-stretching regime flow-thins the suspension as much as the true bond breakage regime that follows at higher Pe. Overall, in strongly attractive, hydrodynamically interacting suspensions, flow-thinning is a two-stage process consisting of initial bond stretching, where increasing asymmetry reduces attractive drag η^attr, followed by subsequent bond breakage, where advection becomes strong enough to separate particles downstream and drain away the equilibrium structure, g^eq.

Motivated by the detailed evidence of evolving bond dynamics, we next look at a rheological quantity that emerges directly from particle bonds and interparticle forces: the osmotic pressure. The nonequilibrium osmotic pressure $Πneq$ in a flowing attractive suspension cannot be described simply by the second virial coefficient $B2*$ because, first, hydrodynamic interactions influence $Πneq$—and hydrodynamic interactions are not encoded in $B2*$—and second, $B2*$ does not account for flow-dependent changes in bond dynamics. This flow-induced change in osmotic pressure is a rheological quantity obtainable from the suspension stress [Eq. (26)].

Interparticle forces set the osmotic pressure: repulsive forces contribute positive osmotic pressure and attractive forces contribute negative osmotic pressure. Interparticle bonds allow a pair to experience both attractive and repulsive forces when flow and Brownian motion change pair separation: a bond can be compressed, relaxed, stretched, or broken, which contribute positive, zero, negative, and zero osmotic pressure, respectively. Furthermore, hydrodynamic stresses arising from particle motion also change the flow-induced osmotic pressure [Eq. (35)]. In this section, we will examine the tandem influence of attractions and hydrodynamic interactions on osmotic pressure and explain the behavior from the perspective of bond dynamics influenced by flow.

The first effects of flow (and hydrodynamic interactions) on the osmotic pressure appear under weakly nonlinear forcing at O(Pe²). Previous theoretical work established that for hard-sphere suspensions, $Πneq$ is always positive, but hydrodynamic interactions suppress flow-induced osmotic pressure.^53,54 Although attractions decrease the equilibrium osmotic pressure, they also change the influence of particle motion on nonequilibrium, flow-induced osmotic pressure. We solved the Smoluchowski equation for linear and weakly nonlinear flow to obtain the O(Pe²) structural distortion, presented in  Appendix B, from which we then obtained the effects of weak flow on osmotic pressure, $Πneq$⁠. The evolution of this weak flow-induced osmotic pressure, $Πneq/Pe2$⁠, with attraction strength is plotted in Fig. 9 for various strengths of hydrodynamics, showing that $Πneq$ evolves non-monotonically with attraction strength, $1−B2*$⁠. Weak-to-moderate attractions suppress nonequilibrium osmotic pressure [Fig. 9(a)], but this trend reverses near $B2*≈−0.5$⁠. When attraction strength is increased such that $B2*<−1$⁠, $Πneq$ becomes linear in the second virial coefficient, $Πneq∼−B2*$⁠, as shown in Fig. 9(b). This behavior is reminiscent of our previous results for microviscosity in the linear-response regime.⁴² However, we show in  Appendix C that unlike the linear-response structure and viscosity, the strong-attraction scaling of $Πneq$ with $B2*$ is not model agnostic in the absence of hydrodynamics.

Hydrodynamic interactions suppress flow-induced osmotic pressure at all attraction strengths, but most notably, strong hydrodynamic interactions can give rise to negative flow-induced osmotic pressure. In the weak-forcing limit of Pe ≪ 1, this negative osmotic pressure is small in magnitude relative to the equilibrium osmotic pressure, $Πneq∼Pe2≪|4B2*|$⁠, but it reveals a tendency for flow to aggregate particles, which may induce phase separation in denser suspensions. Negative $Πneq$ appears at $B2*$ close to the Vliegenthart–Lekkerkerker critical-point criterion,⁸ suggesting that weak flow could precipitate phase separation in an equilibrium-stable suspension.

Overall, under weakly nonlinear forcing, attractions first reduce and then increase the osmotic pressure as they grow stronger, with a well-defined minimum flow-induced osmotic pressure, $Πneq$⁠. Hydrodynamic interactions suppress flow-induced osmotic pressure over all attraction strengths, but most notably, they induce negative flow-induced osmotic pressure for a specific range of attraction strengths between $0.7<1−B2*<3$⁠. In the weakly nonlinear regime, negative osmotic pressure is only O(Pe²) but could potentially condense a suspension when hydrodynamics are strong. This non-monotonic behavior suggests a competition between Brownian, attractive, and hydrodynamic forces that changes with $B2*$⁠, arising from transitions in the distorted microstructure.

For a micromechanical understanding of the minimum in flow-induced osmotic pressure at moderate attraction strengths, we examine the contributions of different microscopic forces to nonequilibrium osmotic pressure. In freely draining suspensions, repulsive and attractive interparticle forces are the sole contributions to osmotic pressure, Π^rep and Π^attr, respectively. Their nonequilibrium contributions are plotted in Fig. 10(a), which shows that as $B2*→0$⁠, attractions grow stronger and almost completely negate the contributions of repulsive forces to $Πneq$⁠. Beyond this threshold attraction strength, attractive forces become the dominant contribution to osmotic pressure. Near $B2*=−0.5$⁠, $Πneqrep$ and $Πneqattr$ cross over and change sign; this corresponds to a reversal in the O(Pe²) distorted structure, as we show in  Appendix B.

When hydrodynamic interactions are strong, they produce dissipative osmotic pressure contributions $ΠneqP,dis$ that counteract elastic interparticle forces [Eq. (35)], weakening $Πneqattr$ and eliminating $Πneqrep$⁠, which is replaced by an entropic contribution from the Brownian drift,⁵⁴ $ΠneqB$⁠. In this strong-hydrodynamic limit, the weakened attractive and the Brownian osmotic pressure contributions behave qualitatively similar to the freely draining $Πneqattr$ and $Πneqrep$⁠, respectively [Fig. 10(a) vs Fig. 10(b)]. However, now, a direct hydrodynamic osmotic pressure arises from probe forcing, $ΠneqH,ext$⁠, which is not counteracted by the other microscopic forces. These hydrodynamic forces increase $Πneq$ when particles accumulate upstream and decrease $Πneq$ when particles accumulate downstream. Hydrodynamic interactions enhance downstream particle retention, and therefore, they reverse the sign of $ΠneqH,ext$⁠, reducing the total flow-induced osmotic pressure $Πneq$ enough to reverse its sign as it reaches the minimum value in Fig. 10.

Summarizing Figs. 9 and 10, the non-monotonic evolution of $Πneq$ is explained by the separate contributions from attractions and entropic forces, $Πneqattr$ and $Πneqrep$ or $ΠneqB$⁠, which are nearly equal in magnitude and opposite in sign and eventually cross over and reverse sign because attractions flip the O(Pe²) distorted structure. Strong hydrodynamic interactions add an external force-induced hydrodynamic contribution $ΠneqH,ext$ that is not balanced by the other forces and becomes negative near the minimum $Πneq$⁠, producing negative flow-induced osmotic pressure, $Πneq<0$⁠. Hydrodynamic interactions and attractions reinforce one another: attractions create bonds, hydrodynamics strengthen those bonds, and together, these reverse the dipolar structure that, in turn, gives rise to negative hydrodynamic and flow-induced osmotic pressure. The nonequilibrium osmotic pressure is small compared to the equilibrium osmotic pressure but may persist at higher Pe, which could trigger flow-induced phase separation.

For strongly nonlinear forcing, hydrodynamic interactions dramatically influence the microstructure and $Πneq$⁠. Hydrodynamic effects cannot be captured by $B2*$⁠, and they become much more pronounced when advection becomes stronger than diffusion, Pe ≳ 1. In Fig. 11, flow-induced osmotic pressure (scaled by Pe) is plotted as a function of flow strength for several different attraction strengths. Figure 11(a) reveals that $Πneq$ is strictly positive for hard spheres and weakly attractive suspensions (0 ≤ A ≤ 17) for all Pe and is monotonically suppressed by attractions over all Pe. However, stronger attractions produce a low-Pe pressurized peak followed at higher Pe by depressurization and then reversal again to pressure increase. This pressurize/decompress/pressurize sequence becomes more pronounced as attractions get stronger, leading to a dramatic overpressure condition followed by negative osmotic pressure at Pe = O(1). The lower-Pe osmotic pressure increase was discussed above in the weakly nonlinear limit and results from a decrease in the number of nearby particles with increasing flow strength, decreasing the density of stretched bonds and increasing the density of unbound pairs. The reversal near Pe ∼ 1 occurs because moderate flow strength increases the flow of particles to the probe but quickly transfers those particles downstream—where flow is too weak to break bonds or sweep them away. The resulting accumulation wake exerts attractive forces on the probe that make a negative contribution to $Πneq$ and lower the osmotic pressure. However, stronger flow severs bonds and depletes the wake, removing the source of negative osmotic pressure and causing another reversal, this time to a positive flow-induced osmotic pressure that grows with Pe. In this strong-forcing limit, $Πneq$ approaches the high-Pe behavior of purely repulsive hard spheres, recovering the correct scaling of $Πneq∼Pe$ for freely draining systems.⁵³ 

Hydrodynamic interactions exaggerate this behavior: in Fig. 11(b), the flow-induced drop in osmotic pressure is more pronounced. Although it is already known that hydrodynamic stresses always act to suppress osmotic pressure,⁵⁴ the effect is stronger when combined with attractions. That is, an interplay between hydrodynamic coupling and attractions leads to more strongly negative contributions of both to osmotic pressure. Mechanistically, hydrodynamics enhance the bond durability, leading to a denser and more durable accumulation wake that comprises bonds stretched by flow. The stretched bonds produce negative osmotic pressure, and the relative motion of downstream wake particles coupled hydrodynamically to the probe produces a negative hydrodynamic pressure. Together, these further reduce $Πneq$ (the details of each osmotic pressure contribution are given in  Appendix G). In Fig. 11(b), the flow-induced osmotic pressure is dominated by hydrodynamic pressure $ΠneqH,ext$ in the intermediate-forcing regime of 1 ≲ Pe ≲ 100 for any sufficiently strong attractions, even when the initial, O(Pe²) osmotic pressure is positive.

Flow can become strong enough to detach the accumulation wake, but this threshold is higher when hydrodynamic interactions matter because lubrication couplings to downstream particles require strong flow to pull them away from the probe. When flow is strong enough to sever all downstream bonds, the high-Pe hard-sphere behavior for hydrodynamically interacting suspensions is recovered,⁵⁴ $Πneq∼Peδ$⁠. Overall, hydrodynamic interactions enhance flow-induced negative osmotic pressure in attractive suspensions. To understand whether this leads to phase instability under flow, we turn our attention next to the total osmotic pressure to determine if its total value can be made sufficiently negative by flow and attractions.

The pair-level osmotic pressure, $ΠP$⁠, reveals the tendency of a flowing suspension toward either expansion or contraction, and here, we look at conditions that can drive $ΠP$ low enough to produce phase instability. The Vliegenthart–Lekkerkerker (VL) critical point⁸ corresponds to a value of the reduced second virial coefficient $B2*=−1.5$⁠. This critical value of $B2*$ is associated with an equilibrium osmotic pressure, Π^eq/n_akT = 1 − 6ϕ_b. Values of $B2*$ more negative than this threshold value robustly predict equilibrium phase separation in colloids.⁷ Here, we seek an expanded definition of this criterion using the nonequilibrium osmotic pressure.

The particle-phase osmotic pressure $ΠP$⁠, above and beyond the ideal (single-particle) osmotic pressure, may be defined as the total osmotic pressure with the ideal (n_akT) contribution removed. Following Eqs. (24) and (26), it comprises both equilibrium and nonequilibrium contributions,

This particle-phase osmotic pressure is plotted in Figs. 12(a) and 12(b) as a function of flow strength Pe at varying attraction strengths for freely draining and hydrodynamically interacting suspensions, respectively. The dashed horizontal line corresponds to the value of the VL equilibrium critical point, which is valid only for Pe = 0, but is shown across all values of Pe to suggest that any points on the plot that fall below this line indicate conditions favorable for aggregation or phase separation in more concentrated conditions.

The left-most point of each curve on the plots recovers the equilibrium osmotic pressure {recalling its connection to $B2*$ [Eq. (23)]}, and at low Pe, the equilibrium contribution dominates the particle-phase osmotic pressure. As flow strength grows, $B2*$⁠, which is a constant, no longer adequately describes osmotic pressure. Instead, $Πneq$ begins to dominate at Pe ≳ O(1). A notable feature—and the major qualitative difference between the no-hydrodynamic (a) and strong-hydrodynamic (b) cases—is a dramatic reduction in osmotic pressure at moderate flow strength (1 < Pe ≲ 10) when hydrodynamics are strong [Fig. 12(b)]. Flow drives $ΠP$ significantly below its equilibrium value, even crossing below the VL critical point (the dotted line in Fig. 12), i.e., $ΠP/nakT<−6ϕb$⁠, showing that flow combined with hydrodynamics and attractions reduces suspension stability. This result, coupled with the fact that negative osmotic pressure drives nonequilibrium phase separation in more concentrated colloidal systems,^26,36 suggests a possible extension of the VL critical point to flowing suspensions—a flow-dependent “critical curve” that accounts for the combined effects of interparticle attractions, flow, and hydrodynamic interactions, and can predict flow-induced changes in stability, as described by the osmotic pressure.

The prospect of extending the VL critical point to flowing suspensions motivates the construction of a nonequilibrium phase map. Figure 13 shows two heat maps of the total particle-phase osmotic pressure, $ΠP$⁠, as a function of flow (horizontal axis) and attraction (vertical axis) strength for freely draining [Fig. 13(a)] and hydrodynamically interacting [Fig. 13(b)] suspensions. The red regions correspond to high, positive osmotic pressure; the blue regions correspond to strongly negative osmotic pressure; and the green regions correspond to zero osmotic pressure. An isobar (black dotted line) is drawn in both plots to demarcate the proposed nonequilibrium extension to the Vliegenthart–Lekkerkerker critical-point boundary,

Hydrodynamic suppression of osmotic pressure is apparent in the comparison between Figs. 13(a) and 13(b): in the bottom half of Fig. 13(b), the critical curve isobar is shifted rightward to 500% higher flow strength for strong attractions (bottom of each plot, $1−B2*∼100$⁠). In addition, hydrodynamics drive the critical curve upward, showing that changes in the strength of hydrodynamic interactions by, e.g., addition of salt can trigger nonequilibrium phase instability. Similar phase maps for systems at arbitrary strength of hydrodynamic interactions and other attractive potentials are given in  Appendix H. Overall, the phase maps enable prediction of flow-induced instability that could trigger phase separation in denser suspensions.

In this work, we have examined the viscosity and osmotic pressure in flowing, attractive colloidal suspensions, with a view toward connecting this behavior with suspension stability, a challenge prevalent in materials ranging from personal-care products to biofilms to industrial slurries. Specifically, we interrogated how microscopic interactions alter rheology enhance or hinder suspension stability. We began with the premise that the second virial coefficient traditionally used to describe the strength of colloidal attractions and the colloidal phase envelope predicted by the Vliegenthart–Lekkerkerker⁸ criterion may have a flow-induced analog for predicting phase stability in a flowing suspension: the osmotic pressure. We chose the active microrheology framework owing to its primacy in the pair interactions that underlie relative motion at the heart of suspension stability. We constructed a theoretical framework for describing the microstructural evolution of a hard-sphere suspension of hydrodynamically interacting colloids that exert attractive forces on one another and solved the resulting equations over a range of strengths of flow, attractive forces, and hydrodynamic interactions.

With probe motion as the driver of flow interacting with a nearby particle, we found that attractions condense or even reverse the microstructural distortion around the probe in the linear-response limit (Fig. 3). The rotating doublets underlying this reversal are quantitatively impacted by lubrication forces in this limit, recovering prior observations in the linear response of sheared suspensions.⁷⁹ When hydrodynamics are strong, attractions simply drive up the low-Pe Newtonian plateau (Fig. 7). However, as flow gets stronger, the flow-thinning seen in hard-sphere suspensions at Pe ≲ 1 is nearly absent in the presence of attractions because attractions smooth the microstructural asymmetry that underlies shear-thinning. Hydrodynamic interactions qualitatively change the impact of attractive forces as flow becomes weakly nonlinear, inducing low-Pe flow-thickening [Fig. 8(a)]. Without hydrodynamic interactions, this attraction-induced flow-thickening does not occur; it arises because hydrodynamics increase bond durability to the point of aggregation. This pair-level behavior signals potential phase instability, similar to the VL criterion, which also predicts phase transition based on pair interactions. The osmotic pressure is also affected: attractions suppress the weakly nonlinear osmotic pressure, and when hydrodynamic interactions also matter, this leads to negative flow-induced osmotic pressure [Fig. 9(a)]. A comparison to the VL critical value of the second virial coefficient, a proxy for equilibrium osmotic pressure, shows that the equilibrium contribution still dominates and the suspension stability is dictated by the VL criterion.

In the opposite limit of strong flow, we found that three distinct structural and dynamical regions form upstream from the probe (Fig. 6), which explains mechanistically the delay in the high-Pe flow-thickening: attractions drain the upstream particle density that, in purely repulsive suspensions, leads to flow-thickening. As flow gets even stronger, a competition between attractions and advection frustrates the attainment of structural symmetry expected in the pure-hydrodynamic limit. This competition also impacts structural asymmetry in freely draining suspensions, delaying the formation of a clean wake and thus attainment of the high-Pe plateau.

In between these two limits, in the intermediate-flow regime, 1 ≲ Pe ≲ 100, hydrodynamics qualitatively alter non-Newtonian rheology by exaggerating the flipped, distorted microstructure to the point that two-stage flow-thinning can emerge—attractive flow-thinning due to bond stretching, followed at stronger flow by hydrodynamic flow-thinning where flow breaks bonds and drains away particle density from the trailing wake (Fig. 8). The cooperative actions of attractions and hydrodynamic interactions also induce strongly negative osmotic pressure, which can drive the particle-phase pressure below values associated with the equilibrium VL point (Figs. 12 and 13). These conditions and thus this phenomenology are predicted to occur in biological complex fluids, for example. In protein solutions such as Ubiquitin (UBQ) and the IgG-binding domain of protein G (GB3), the second virial coefficients have been reported to be in the range⁹² $−1<B2*<0$⁠, values that fall within the range of attraction strengths where we expect negative flow-induced osmotic pressure when Pe < O(100). These flow conditions are realized in vitro during preparation of these solutions, as well as in vivo when cytoplasmic streaming and flow produced by motor proteins reach Pe ≈ 10 and Pe ≈ 5, respectively.⁹³ 

The idea of an extension of the Vliegenthart–Lekkerkerker⁸ criterion for the critical point to detect phase stability in flowing suspensions resulted in a non-equilibrium phase map, where we find that the particle-phase osmotic pressure, $ΠP$⁠, is qualitatively altered by the combination of hydrodynamic interactions, interparticle attractions, and flow—to the point that the onset of negative osmotic pressure that tends to condense the particle phase is driven to weaker attractions. That is, flow can induce instability in an otherwise stable suspension, especially if hydrodynamics are strong. We identified an isobar to demarcate a proposed nonequilibrium extension to the VL critical-point boundary (Fig. 13). Hydrodynamic suppression of osmotic pressure can shift the critical isobar to much higher flow strength and drive the critical curve upward, showing that changes in the strength of hydrodynamic interactions by, e.g., addition of salt can trigger nonequilibrium phase instability.

These results for microrheology have a straightforward connection to sedimentation and stability during gravitational settling. It is natural to ask whether these results for microrheology can be generalized to sheared suspensions. Prior work in microrheology has demonstrated its ability to recover many important qualitative behaviors observed in shear rheology, including shear-thinning and thickening,^46,61,62 flow-induced diffusion,^{63–65,94,95} and normal stress differences,^51,53,54,67 in excellent agreement with both dynamic simulation⁹⁶ and experimental studies.^{52,70–72,97,98} In a microrheology experiment, we predict that a probe particle of a stickiness $B2*<0$⁠, forced at a value of Pe within the regime of negative flow-induced osmotic pressure, would tend to aggregate particles as it moves through the suspension, effectively increasing its size and, by extension, its Pe. Our model also predicts that above a threshold Pe ∼ 100, flow-induced osmotic pressure reverses and begins to increase with flow strength. The macrorheological analog of this system is shear-induced aggregation;^99–101 by analogy with our results, a suspension of sticky particles sheared at intermediate Pe would be expected to form larger aggregates than observed at equilibrium, with potential phase separation occurring at higher concentrations.