Equilibrium phase instability of colloids is robustly predicted by the Vliegenthart–Lekkerkerker (VL) critical value of the second virial efficient, but no such general criterion has been established for suspensions undergoing flow. A transition from positive to negative osmotic pressure is one mechanical hallmark of a change in phase stability in suspensions and provides a natural extension of the equilibrium osmotic pressure encoded in the second virial coefficient. Here, we propose to study the non-Newtonian rheology of an attractive colloidal suspension using the active microrheology framework as a model for focusing on the pair trajectories that underlie flow stability. We formulate and solve a Smoluchowski relation to understand the interplay between attractions, hydrodynamics, Brownian motion, and flow on particle microstructure in a semi-dilute suspension and utilize the results to study the viscosity and particle-phase osmotic pressure. We find that an interplay between attractions and hydrodynamics leads to dramatic changes in the nonequilibrium microstructure, which produces a two-stage flow-thinning of viscosity and leads to pronounced flow-induced negative osmotic pressure. We summarize these findings with an osmotic pressure heat map that predicts where hydrodynamic enhancement of attractive bonds encourages flow-induced aggregation or phase separation. We identify a critical isobar—a flow-induced critical pressure consistent with phase instability and a nonequilibrium extension of the VL criterion.

## I. INTRODUCTION

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,^{19} 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*_{2}, which describes the effect of pair-level interactions on the osmotic pressure in a virial expansion.^{39} Attractive forces decrease *B*_{2}, 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*_{2}, sets the onset of equilibrium phase separation. The threshold value of *B*_{2} is given by the Vliegenthart–Lekkerkerker (VL) criterion,^{8} which establishes the colloidal gas–liquid critical point at *B*_{2} = −6*v*_{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.^{42} The structural distortion induced by weak flow and the viscosity in the Newtonian plateau are quantitatively predictable through *B*_{2} irrespective of the detailed functional form of the interparticle attraction.^{42}

However, nonlinear rheology arising from stronger flows cannot be predicted by the second virial coefficient because *B*_{2} encodes only the equilibrium structure, while strong flow shifts the balance of microscopic forces and displaces the structure from equilibrium.^{43} 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.^{26} 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*_{2} can predict the onset of thermodynamic colloidal phase instability via the VL criterion,^{8} 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*_{2} 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^{22} 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,^{46} 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.

## II. MODEL SYSTEM

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.

### A. Modeling interparticle repulsion

The excluded-shell model is a simplified means of representing a no-flux thermodynamic surface surrounding the no-slip hydrodynamic surface of a particle^{83} (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*_{1} and a bath particle (particle 2) at *x*_{2}, giving a relative configuration ** r** =

*x*_{2}−

*x*_{1}(Fig. 1). The interparticle force exerted on a bath particle by the probe is then given by $F2P=\u2212\u2207rV(r)$, where

**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**

*r**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.

### B. Modeling interparticle attraction

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.^{84}

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} ∼ 2*kT*/*a*. We model an AO potential *V*^{attr} with attraction range 2*a* ≤ *r* ≤ 2*a*(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* = 2*a*)/(2*kT*/*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=\u2212\u2207rVattr$, and the force exerted on particle 1 by particle 2 is $F1attr=+\u2207rVattr$. Given Eq. (3), the attractive force $F1attr$ within the AO attraction range is

Here, $r\u0302$ 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.

## III. THEORETICAL FRAMEWORK

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.^{61} 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.

### A. Steady-state particle microstructure

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*_{2}(*x*_{1}, *x*_{2}) of finding the probe at a position *x*_{1} and any other bath particle at a position *x*_{2}. This probability function obeys a Smoluchowski equation, which in the pair limit reads

where the temporal evolution of particle distribution *∂P*_{2}/*∂t* is balanced by the sum of all particle fluxes, $\u2211\alpha =12j\alpha =\u2211\alpha =12U\alpha P2$, where *U*_{α} is the velocity vector of particle *α*. Here, *j*_{1} = *U*_{1}*P*_{2} denotes the probability flux of the probe and *j*_{2} = *U*_{2}*P*_{2} 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\alpha P$), and entropic (Brownian) forces [$F\alpha B\u2261\u2212kT\u2207\alpha lnP2(x1,x2)$], where particle motion and forces are coupled by the hydrodynamic mobility $M\alpha \beta 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\alpha \beta \u2261kTM\alpha \beta UF$ as required by the fluctuation–dissipation theorem.^{85} 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\alpha \beta A(r)$ and $y\alpha \beta 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*_{1}and a bath particle at

**=**

*r*

*x*_{2}−

*x*_{1}. 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*_{1}(** z**) denotes the probability of finding a probe at a given absolute position

**, and**

*z**P*

_{1/1}(

**|**

*r***) denotes the probability of finding a bath particle at relative position**

*z***, given a probe position**

*r***. The interparticle forces may be expressed as gradients of an interparticle potential, $F\alpha P=\u2212\u2207\alpha V$. This potential includes both attractive and repulsive contributions,**

*z**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,

**, to obtain**

*z*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\u2212x12a$ and $H(r/ah)=y11a\u2212y12a$ 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*_{11} − *M*_{12} = −*D*_{r}/2*kT*. With these scalings, the governing equation and the boundary conditions—rigid spheres permit no relative flux at contact, *j*_{1} − *j*_{2} = 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}/(2*kT*/*a*), emerges naturally from the scalings in Eq. (11). The unit vector $F\u0302ext\u2261Fext/Fext$ points along the applied force in $\u2212e\u0302z$.

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*.

### B. Average probe motion and microviscosity

Stokes’s drag law sets the velocity *U*^{0} 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,^{61} ⟨*U*_{1}⟩,

where the angle brackets around $U1$ denote an average of *U*_{1} 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^{62} and, subsequently, Mohanty and Zia,^{75}

where $\u27e8U1H\u27e9$, $\u27e8U1B\u0304\u27e9$, and $\u27e8U1P\u27e9$ signify the volume-averaged contributions of hydrodynamic, Brownian, and interparticle forces to the mean probe motion, respectively, averaged over all configurations. The overbar on $\u27e8U1B\u0304\u27e9$ 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(\kappa r+r)$, $y11A(\kappa 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.

### C. Nonequilibrium osmotic pressure

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,^{39}

where Π^{eq} denotes the equilibrium osmotic pressure and *n* is the particle number density. The sign of *B*_{2} 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\pi 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.

^{67}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), $\Sigma 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 $\Sigma H,ext$, Brownian drift-induced hydrodynamic stress $\Sigma B$, and dissipative interparticle stress $\Sigma P,dis$—is given by the probe number density times the volume-averaged stresslet: $\Sigma =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 $\Sigma $ 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.^{67}

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\alpha \beta P$, $X\alpha \beta G$, $Y\alpha \beta G$, and $Y\alpha \beta 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\alpha \beta A$, $y\alpha \beta A$, and $y\alpha \beta 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 $\Pi 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.

## IV. RESULTS

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.

### A. Microstructure

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,^{47} 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*=\u22120.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.

^{78}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.^{43} 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*^{−1})-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*=\u22121.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.

### B. Microviscosity

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^{42} [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,^{42} 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*_{1}(*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\u22122<2\Delta ;B2*\u226a0;\theta =\pi )<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}.

### C. Osmotic pressure

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 $\Pi neq$ in a flowing attractive suspension cannot be described simply by the second virial coefficient $B2*$ because, first, hydrodynamic interactions influence $\Pi 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*^{2}). Previous theoretical work established that for hard-sphere suspensions, $\Pi 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*^{2}) structural distortion, presented in Appendix B, from which we then obtained the effects of weak flow on osmotic pressure, $\Pi neq$. The evolution of this weak flow-induced osmotic pressure, $\Pi neq/Pe2$, with attraction strength is plotted in Fig. 9 for various strengths of hydrodynamics, showing that $\Pi neq$ evolves non-monotonically with attraction strength, $1\u2212B2*$. Weak-to-moderate attractions suppress nonequilibrium osmotic pressure [Fig. 9(a)], but this trend reverses near $B2*\u2248\u22120.5$. When attraction strength is increased such that $B2*<\u22121$, $\Pi neq$ becomes linear in the second virial coefficient, $\Pi neq\u223c\u2212B2*$, as shown in Fig. 9(b). This behavior is reminiscent of our previous results for microviscosity in the linear-response regime.^{42} However, we show in Appendix C that unlike the linear-response structure and viscosity, the strong-attraction scaling of $\Pi 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, $\Pi neq\u223cPe2\u226a|4B2*|$, but it reveals a tendency for flow to aggregate particles, which may induce phase separation in denser suspensions. Negative $\Pi neq$ appears at $B2*$ close to the Vliegenthart–Lekkerkerker critical-point criterion,^{8} 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, $\Pi 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\u2212B2*<3$. In the weakly nonlinear regime, negative osmotic pressure is only *O*(*Pe*^{2}) 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*\u21920$, attractions grow stronger and almost completely negate the contributions of repulsive forces to $\Pi neq$. Beyond this threshold attraction strength, attractive forces become the dominant contribution to osmotic pressure. Near $B2*=\u22120.5$, $\Pi neqrep$ and $\Pi neqattr$ cross over and change sign; this corresponds to a reversal in the *O*(*Pe*^{2}) distorted structure, as we show in Appendix B.

When hydrodynamic interactions are strong, they produce dissipative osmotic pressure contributions $\Pi neqP,dis$ that counteract elastic interparticle forces [Eq. (35)], weakening $\Pi neqattr$ and eliminating $\Pi neqrep$, which is replaced by an entropic contribution from the Brownian drift,^{54} $\Pi neqB$. In this strong-hydrodynamic limit, the weakened attractive and the Brownian osmotic pressure contributions behave qualitatively similar to the freely draining $\Pi neqattr$ and $\Pi neqrep$, respectively [Fig. 10(a) vs Fig. 10(b)]. However, now, a direct hydrodynamic osmotic pressure arises from probe forcing, $\Pi neqH,ext$, which is not counteracted by the other microscopic forces. These hydrodynamic forces increase $\Pi neq$ when particles accumulate upstream and decrease $\Pi neq$ when particles accumulate downstream. Hydrodynamic interactions enhance downstream particle retention, and therefore, they reverse the sign of $\Pi neqH,ext$, reducing the total flow-induced osmotic pressure $\Pi 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 $\Pi neq$ is explained by the separate contributions from attractions and entropic forces, $\Pi neqattr$ and $\Pi neqrep$ or $\Pi neqB$, which are nearly equal in magnitude and opposite in sign and eventually cross over and reverse sign because attractions flip the *O*(*Pe*^{2}) distorted structure. Strong hydrodynamic interactions add an external force-induced hydrodynamic contribution $\Pi neqH,ext$ that is not balanced by the other forces and becomes negative near the minimum $\Pi neq$, producing negative flow-induced osmotic pressure, $\Pi 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 $\Pi 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 $\Pi 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 $\Pi 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, $\Pi neq$ approaches the high-*Pe* behavior of purely repulsive hard spheres, recovering the correct scaling of $\Pi neq\u223cPe$ for freely draining systems.^{53}

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,^{54} 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 $\Pi 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 $\Pi neqH,ext$ in the intermediate-forcing regime of 1 ≲ *Pe* ≲ 100 for any sufficiently strong attractions, even when the initial, *O*(*Pe*^{2}) 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,^{54} $\Pi neq\u223cPe\delta $. 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, $\Pi P$, reveals the tendency of a flowing suspension toward either expansion or contraction, and here, we look at conditions that can drive $\Pi P$ low enough to produce phase instability. The Vliegenthart–Lekkerkerker (VL) critical point^{8} corresponds to a value of the reduced second virial coefficient $B2*=\u22121.5$. This critical value of $B2*$ is associated with an equilibrium osmotic pressure, Π^{eq}/*n*_{a}*kT* = 1 − 6*ϕ*_{b}. Values of $B2*$ more negative than this threshold value robustly predict equilibrium phase separation in colloids.^{7} Here, we seek an expanded definition of this criterion using the nonequilibrium osmotic pressure.

The particle-phase osmotic pressure $\Pi P$, above and beyond the ideal (single-particle) osmotic pressure, may be defined as the total osmotic pressure with the ideal (*n*_{a}*kT*) 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, $\Pi 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 $\Pi P$ significantly below its equilibrium value, even crossing below the VL critical point (the dotted line in Fig. 12), i.e., $\Pi P/nakT<\u22126\varphi 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, $\Pi 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\u2212B2*\u223c100$). 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.

## V. CONCLUSIONS

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^{8} 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.^{79} 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^{92} $\u22121<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.^{93}

The idea of an extension of the Vliegenthart–Lekkerkerker^{8} 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, $\Pi 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^{96} 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.

The work presented here focuses on semi-dilute suspensions, a simple model aimed at gleaning insight into the impact of attractions and flow on the viscosity and nonequilibrium osmotic pressure; by analogy with the VL criterion, which can predict phase behavior in dense systems based on a simple pair interaction and which has been empirically supported,^{40,41} our results may provide some insight into flow-induced phase instability of flowing, concentrated suspensions. A future study should include dynamic simulation in more concentrated systems to confirm the possibility of flow-induced phase separation from our semi-dilute theory. The effects of other attraction profiles, especially longer-range or asymmetric (patchy) attractions, are another possible area of future study that could further shed light on the effects of flow on phase transitions in biological cells, for example. Finally, our system considers particles that are forced by an external body force; the cooperative effect of hydrodynamics and attractions on structure and rheology in the active microrheology system raises questions of whether similar phenomena arise in active motion.

## ACKNOWLEDGMENTS

This research was funded, in part, by the Office of Naval Research, Director of Research Early Career Award (N000141812105). D.E.H. thanks Nicholas Hoh, Benjamin Dolata, Henry Chu, and Ritesh Mohanty for many fruitful conversations.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.

### APPENDIX A: HYDRODYNAMIC AND MOBILITY FUNCTIONS AND THE STRESS TENSOR

Values for the resistance^{90} and mobility^{47,86} functions introduced in Sec. III are obtained analytically from the literature.^{87,91,102} Our model is fully general for arbitrary size of probe relative to bath particles, but we restrict our analysis to same-size probe and bath particles. In this case, the scalar hydrodynamic functions depend only on the separation between particles, *r*, and the dimensionless repulsion range, *κ*, which sets the minimum approach distance (and therefore, the strength of hydrodynamic interactions). The hydrodynamic resistance function $R\alpha \beta SU$ couples the translational velocity of particle *β* to a stresslet on particle *α*, and the resistance function $R\alpha \beta S\Omega $ couples the rotational velocity of particle *β* to the stresslet on particle *α*. Following the notation of the work of Jeffrey *et al.*,^{91,102} these are given by

where $r\u0302i$ is the unit vector along the particles’ line of centers, *δ*_{ij} is the identity tensor, *ϵ*_{ijk} is the Levi-Cività tensor, and *X*_{αβ} and *Y*_{αβ} govern the motion of particles *α* and *β* parallel and perpendicular to their line of centers, respectively.

The hydrodynamic mobility functions $M\alpha \beta UF$ and $M\alpha \beta \Omega F$ couple a force on particle *β* to the induced translational and rotational velocity of particle *α*, respectively. These are given by

Here, *x*_{αβ} governs motion along the line of centers and *y*_{αβ} governs motion transverse to the line of centers of particles *α* and *β*. The scalar hydrodynamic functions depend on the center-to-center separation, *r*, and the dimensionless repulsion range, *κ*: in the Smoluchowski equation governing structure and in all the integral expressions for viscosity and nonequilibrium osmotic pressure given in Sec. III, these hydrodynamic functions are evaluated as functions of *r*(*κ* + 1).

### APPENDIX B: LINEAR-RESPONSE AND WEAKLY NONLINEAR SMOLUCHOWSKI EQUATION

In this section, we derive the weakly nonlinear Smoluchowski equation that we solved in the present work to obtain the *O*(*Pe*^{2}) structural distortion. We inserted the resulting microstructure into Eq. (35) to compute the effects of weak flow on nonequilibrium osmotic pressure. As defined in Eq. (13), the structural distortion *f* is the departure from the equilibrium structure *g*^{eq}, scaled by forcing strength, *Pe*,

The Smoluchowski equation governing structural distortion with arbitrary (spherically symmetric) potential *V*(*r*) is

where $Dr\u2261Gr\u0302r\u0302+H(I\u2212r\u0302r\u0302)$ is the dimensionless relative diffusivity tensor. There is no flux at contact and no long-range order,

In the weak-forcing limit, we expand the distortion in *Pe*,

where the functions *f*_{1}, *f*_{2}, and *h*_{2} depend on radial separation *r* only. Inserting Eq. (B4) into Eq. (B2) and gathering the leading-order terms give the Smoluchowski equation for linear-response distortion *f*_{1},

Here, the hydrodynamic function *W* is defined by the divergence of relative diffusivity,

The weakly nonlinear distortion consists of quadrupolar (*f*_{2}) and monopolar (*h*_{2}) components. Insertion of Eq. (B4) into Eq. (B2) and gathering the *O*(*Pe*) terms quadratic in forcing ($F\u0302extF\u0302ext:r\u0302r\u0302$) give the governing equation for *f*_{2},

and the monopolar distortion, which is obtained by gathering the *O*(*Pe*) terms with no angular dependence, is given by

For purely repulsive hard spheres, the interparticle potential is zero away from contact, and these equations match previous results for weakly nonlinear theory of purely repulsive hard-sphere suspensions.^{62,65} Solutions for these microstructures can be obtained numerically in Matlab utilizing a finite-difference scheme. The mobility functions *G*, *H*, and *W* are evaluated numerically using the twin-multipole expansions of Jeffrey and Onishi^{87} and the lubrication approximations of Kim and Karrila.^{48} These are plotted for freely draining systems in Fig. 14 and for strong hydrodynamic interactions in Fig. 15.

The Smoluchowski equations governing linear-response, *O*(*Pe*) distortion were laid out and solved in our previous work^{42} and are shown in Figs. 14(a) and 15(a). Solution of the Smoluchowski equations governing weakly nonlinear distortion reveals that when hydrodynamics are negligible, the transition point for the *O*(*Pe*^{2}) structure, where *f*_{2} and *h*_{2} reverse in sign, roughly coincides with the $B2*=\u22120.5$ threshold observed for the *O*(*Pe*) structure, *f*_{1}. This supports our observations from Fig. 10(a) in Sec. IV C.

We had previously observed that strong hydrodynamic interactions enhance attractive forces, shifting the reversal of low-*Pe* structure to weaker attraction strengths.^{42} This can be seen in the comparison between Figs. 14(a) and 15(a), where the latter transition occurs near $B2*=0.25$, i.e., repulsions are still stronger than attractions, in agreement with the contour plots shown in Fig. 3. The *O*(*Pe*^{2}) structures in Figs. 15(b) and 15(c) show that strong hydrodynamics also shift the transition point in the *O*(*Pe*^{2}) structure to weaker attraction strengths. However, the transitions in *f*_{1}, *f*_{2}, and *h*_{2} no longer coincide; most notably, the reversal in monopolar structure *h*_{2}(*r*) lags structural reversals in the *O*(*Pe*) dipole (*f*_{1}) and the *O*(*Pe*^{2}) quadrupole (*f*_{2}). As shown in Figs. 10 and 16, this lag is responsible for the different transitions observed for $\Pi neqH,ext$ and $\Pi neqattr$ in the low-*Pe* limit.

The total weakly nonlinear structural distortion at contact, *f*_{2}(2) + 3*h*_{2}(2), is plotted as a function of $B2*$ in Fig. 16, with the stagnation-point linear-response distortion *f*_{1}(2) for strong hydrodynamic interactions plotted on the opposing axis for comparison. When the weakly nonlinear stagnation-point distortion *f*_{2}(2) + 3*h*_{2}(2) reverses, weakly nonlinear flow dislodges bonded particles and reduces nearby particle density. In this strongly attractive regime, attractions dominate the osmotic pressure, so when weakly nonlinear flow reduces particle density, it, in turn, reduces the source of negative osmotic pressure, thus increasing $\Pi neq$, as shown in Figs. 9 and 10. The linear-response distortion *f*_{1}(*r*) contributes to $\Pi neq$ only when hydrodynamic interactions are strong, and because the transition in weakly nonlinear distortion lags that for *f*_{1}(*r*) (Fig. 16), hydrodynamic interactions enable negative flow-induced osmotic pressure, $\Pi neq<0$ (Figs. 9 and 10).

### APPENDIX C: IS NONLINEAR-RESPONSE OSMOTIC PRESSURE “MODEL AGNOSTIC”?

In Sec. IV, we investigated the weakly nonlinear osmotic pressure compared to the hard-sphere value for an attraction range Δ = 0.1. Here, we present results for longer-range attractions. The purpose of this extended parameter study is to understand whether the potential-agnostic response of the linear-response regime extends to the weakly nonlinear-response regime.

We use the expressions for the stress derived in Sec. III, inserting the weakly nonlinear microstructure from Appendix B. The *O*(*Pe*^{2}) osmotic pressure $\Pi neq$ is plotted in Fig. 17 as a function of attraction strength, $1\u2212B2*$, for various AO attraction ranges. For both weak and strong hydrodynamic interactions, variation in the attraction range drives pronounced differences in osmotic pressure that $B2*$ alone cannot resolve. Notably, the longest-range attractions result in the least suppression of osmotic pressure, and although the strong-attraction [Fig. 17(b)] scaling of $\Pi neq$ with $B2*$ remains linear, the scaling coefficient is no longer consistent, ranging from −2.0 to as high as −7.5 for strong hydrodynamic interactions. These results demonstrate that the *O*(*Pe*^{2}) structure and rheology are much more dependent on the specific attraction profile, especially when hydrodynamic interactions are strong.

### APPENDIX D: FINITE-DIFFERENCE SCHEME

Computation of volume-averaged quantities such as the microviscosity and the nonequilibrium osmotic pressure at arbitrary forcing and attraction strengths requires numerical solution of Eq. (12) for the microstructure, *g*(** r**). To this end, we employ a centered second-order finite-difference scheme to compute the microstructural distortion,

*f*(

**), which is defined relative to the total microstructure**

*r**g*as

*g*(

**) ≡**

*r**g*

^{eq}[1 +

*Pef*(

**)]; this aids numerical stability during computation by removing the steep gradients present in the equilibrium structure**

*r**g*

^{eq}for strongly attractive suspensions. In axisymmetric spherical coordinates, the Smoluchowski equation governing structural distortion is given by

where the hydrodynamic functions *G*(*s*), *H*(*s*), and *W*(*s*) are evaluated in terms of *s* ≡ (*κ* + 1)*r*. The boundary conditions are

As forcing grows strong, steep gradients in *f*(** r**) arise along the radial coordinate near probe contact, which requires special treatment of the gridpoints in the boundary layer. To faithfully represent the physics in this region, we map the radial coordinate to a finite domain, following the work of Hoh and Zia.

^{51,64}To do so, we define a new radial coordinate,

*T*,

We additionally define a parameter, *ω*, for exponential mapping of *T*,

where *T*_{∞} is the largest stretched radial distance calculated and Δ*y* is the increment for the finite radial coordinate, *y*,

In this work, we solve the Smoluchowski equation, thus transformed, for *f*(*y*, *θ*), utilizing the LaPack direct banded solver in Matlab. We use 1200 logarithmically spaced grid points for the radial coordinate *r* and 1200 uniformly spaced grid points for the polar coordinate *θ*.

### APPENDIX E: ADDITIONAL MICROSTRUCTURE PLOTS

By solving Eq. (12) numerically via the finite-difference method described in Appendix D, we obtained the microstructural distortion *f*(** r**) at arbitrary strengths of hydrodynamic interactions, attractions, and external forcing. Here, we present plots of our solutions to Eq. (12) for a broad range of

*Pe*and

*A*. Figure 18 contains contour plots of

*f*(

**) at varying strengths of the Asakura–Oosawa potential,**

*r**A*, and Péclet numbers,

*Pe*, for an attraction range of Δ = 0.1 and a relative repulsion range of

*κ*= 10

^{−4}, i.e., strong hydrodynamic interactions. Going from left to right, the forcing strength increases, and structural asymmetry grows as a result. Going down Fig. 18, attraction strength increases, and the structural reversal caused by attractions becomes apparent. A similar figure for the microstructure in the absence of hydrodynamic interactions is given in our previous work.

^{43}In comparison with the no-hydrodynamic case, strong hydrodynamic interactions facilitate structural reversal, i.e., downstream accumulation, whether in a dipole or an accumulation wake, at weaker attraction strengths. For strongly attractive suspensions under strong forcing, hydrodynamic interactions dramatically enhance the durability of the downstream accumulation wake, which persists at significantly higher

*Pe*compared to the freely draining case

^{43}before splitting due to advective bond breakage.

### APPENDIX F: TRAJECTORY ANALYSIS

Our study examines probe motion through a microstructural and statistical-mechanics perspective. However, analyzing relative trajectories between the probe and bath particles provides further insight into the role of attractions on probe motion and viscosity. In the steady, Stokes-flow regime, with zero Brownian motion, particle velocities may be written as linear combinations of external and interparticle forces, where the mobility tensors $M\alpha \beta UF$ couple external forces $Fext$ and interparticle forces *F*^{P} to the translational velocity,

The motion of a bath particle relative to a probe is thus given as

To determine complete particle trajectories, we integrated the equations over time utilizing a Runge–Kutta method. Trajectories for the purely repulsive hard-sphere case and relative attraction strengths *F*^{P}/*F*^{ext} = 0.1, 0.3, and 0.5 are plotted in Figs. 19 and 20.

When hydrodynamic interactions are negligible, bath particles passing at an initial vertical offset of $y\u03020\u22652+2\Delta $, beyond the range of either repulsions or attractions, are completely unaffected by the probe. When the particles are purely repulsive, as in Fig. 19(a), bath particles with an initial vertical offset less than the particle diameter, $y\u03020\u22642$, collide with the probe before shifting downstream along its surface due to advection. Upon reaching the downstream face ($x\u0302\u22640$), bath particles detach immediately at vertical or *θ* = *π*/2. The addition of attractive forces changes the nature of these two-particle encounters. Within the range of attractions, $0\u2264r\u0302\u22122<2\Delta $, the bath particle is pulled toward the probe center. When attractions are much weaker than advection, as in Fig. 19(b), the downstream wake is narrowed slightly and the duration of a two-particle encounter increases by a small amount. Increasing attraction strength, as shown in Fig. 19(c), narrows the wake significantly; bath particles passing almost anywhere inside the attraction range are advected to the probe surface and the downstream trajectories merge. A strong enough attraction strength prevents a bath particle from ever detaching, as seen in Fig. 19(d).

In the opposite limit of strong hydrodynamics, entrainment between the probe and a bath particle affects relative motion at much longer distances. Figure 20 plots the relative trajectories of particles with a very low repulsion range, *κ* ≪ 1. Relative trajectories for purely repulsive particles are shown in Fig. 20(a); in comparison with the purely repulsive trajectory map shown in Fig. 19(a), fore-aft symmetry is preserved, and hydrodynamic interactions cause closely approaching bath particles to separate farther downstream, i.e., *θ* < *π*/2. The microstructural analog, as seen in Fig. 5, is a narrowing of the trailing wake. When both hydrodynamic interactions and attractions are present, the downstream wake is narrowed even more dramatically, as shown in Figs. 20(b) and 20(c), illustrating the tandem action of lubrication forces and attractions in delaying pair separation. In Figs. 19(d) and 20(d), attractive forces are strong enough to keep bath particles permanently attached despite external forcing.

Although the plots of relative trajectories help illustrate the hydrodynamic enhancement of attractive bonds and its impact on microstructure, the most dramatic effects are shown by plotting encounter times, as shown in Fig. 21. Here, we show the dimensionless, advectively scaled encounter time, *t*_{pass}/*τ*_{adv}, as a function of the initial vertical offset, *y*_{0}, for varying strengths of hydrodynamics, both with and without attractions—here, set at *F*^{attr}/*F*^{ext} = 0.3 to enable advective separation. When bath particles pass the probe outside its attraction range, there is no observed difference in encounter time (i.e., the solid and dashed lines coincide). As bath particles approach within the probe’s attraction range, however, the pair forms a bond, and the bath particle is advected to the downstream face of the probe, attractions delay detachment, increasing the observed encounter time; this is seen in Fig. 21 as the sharp upturn in the curves where *F*^{attr}/*F*^{ext} = 0.3 (solid lines with symbols). In the complete absence of Brownian motion, attractions only slightly increase the encounter time observed for negligible (*κ* → ∞). An intermediate (*κ* = 0.1) strength of hydrodynamics, where the dimensionless repulsion range is reduced to *κ* = 0.1, results in a significant increase in encounter times as the bath particle approaches closely enough to experience attractions, as evidenced by the larger upturn in the green curve (square symbols); hydrodynamics are strong enough to reinforce attractive bonds, but repulsions still shield the particles from experiencing strong lubrication interactions. However, when repulsion range is reduced to *κ* = 10^{−5}, an even more dramatic increase in encounter time with attractions is observed: as a bath particle moves within the attraction range, it becomes bonded strongly enough to resist separation. Without any attractions, a bath particle can form a permanently joined doublet due to lubrication interactions if *y*_{0} < 0.5, but attractive forces increase this “filter range” to *y*_{0} = 1.0. For all bath particles passing within this vertical range of −1 < *y*_{0} < 1, attractions combine with advection upstream to rapidly move bath particles to the probe surface; the pair forms a permanently joined doublet that advection cannot separate, even when attractions or hydrodynamic interactions alone would not suffice to keep the particles bonded.

It is important to note that exact trajectories can only be determined in the non-colloidal limit of *Pe*^{−1} ≡ 0, where Brownian motion is negligible and particle trajectories are fully deterministic. Even relatively weak Brownian motion can help to detach bath particles downstream; however, strong hydrodynamic interactions reduce the ability of Brownian motion to detach particles, causing the dramatically increased bond durability seen in the microstructure at finite *Pe* (Fig. 5).

### APPENDIX G: MECHANISTIC ORIGIN OF NEGATIVE FLOW-INDUCED OSMOTIC PRESSURE

In Sec. IV, we demonstrated the qualitative influence of hydrodynamic interactions on flow-induced osmotic pressure. Here, we present the mechanistic origin of this behavior by examining the various contributions to nonequilibrium osmotic pressure—attractions, repulsions, Brownian drift, and hydrodynamic drag [Eq. (35)]—as they evolve with flow. In this section, we describe attraction strengths both in terms of dimensionless AO force, *A*, and in terms of the decrease in the reduced second virial coefficient from the hard-sphere value, $B2attr\u22611\u2212B2*$.

The attractive and repulsive contributions to the osmotic pressure in freely draining (*κ* → ∞) are plotted in Figs. 22(a) and 22(b), respectively. Most notably, the attractive contribution to osmotic pressure plateaus as flow gets very strong, but, in contrast, the repulsive contribution grows without bound. For *A* > 17 ($B2attr>1.5$), the attractive osmotic pressure is dominant when *Pe* < 1, but strong flow [*Pe* ≳ *O*(100)] eventually breaks attractive bonds, draining away particle density within the attraction well and removing the source of negative osmotic pressure. Meanwhile, the constant growth of repulsive flow-induced osmotic pressure with flow strength eventually drowns out the influence of attractions. The value of *Pe* at which this occurs sets the location of the critical isobar in Fig. 13.

Mechanistically, the negative flow-induced osmotic pressure that occurs in suspensions with strong hydrodynamics is hydrodynamic in origin: the external force-induced osmotic pressure $\Pi neqH,ext$, which arises due to entrainment between no-slip surfaces [cf. Eqs. (31) and (35)]. This is deduced from Fig. 23(a), which shows the hydrodynamic force-induced osmotic pressure, $\Pi neqH,ext$, scaled both advectively and attractively, i.e.,

as it evolves with forcing strength *Pe*, at various attraction strengths. As attractions grow strong, prominent troughs arise in $\Pi neqH,ext$ in Fig. 23(a), mirroring the troughs we observed as flow-reduced osmotic pressure $\Pi neq<0$ in Fig. 11. The structural origin of this negative hydrodynamic osmotic pressure $\Pi neqH,ext<0$ is a reversal in the distorted structure, as previously shown in Sec. IV A, which persists at intermediate-forcing strengths of 1 < *Pe* < 10. The hydrodynamic contribution $\Pi neqH,ext$ is positive when particles accumulate upstream and negative when particles accumulate downstream [cf. Eq. (35)], as evidenced by the direct link between structural reversal and negative $\Pi neqH,ext$ shown in Figs. 10(b) and 16. Overall, the prominent flow-induced reduction in osmotic pressure is the direct result of external force-induced hydrodynamic pressure $\Pi neqH,ext$ and the indirect result of attraction-induced structural reversal, *f*(** r**).

For hard-sphere suspensions with strong hydrodynamic interactions, it is known that although $\Pi neqH,ext$ dominates $\Pi neq$ under strong forcing, the osmotic pressure contribution from Brownian drift continues to grow with *Pe*. In Fig. 23(b), we turn to the other osmotic pressure contributions and show for strongly attractive systems that, as external forcing grows strong, the sum of the interparticle (both equilibrium plus nonequilibrium) and Brownian contributions to osmotic pressure vanishes: $(\Pi eq/nakT\u22121)+\Pi neqattr+\Pi neqrep+\Pi neqB\u21920$ as *Pe* → ∞, and the osmotic pressure becomes almost purely hydrodynamic. This is shown in Fig. 23(b), where the black dotted line shows the evolution of the Brownian drift term $\Pi neqB$ (plus an insignificant repulsive contribution, $\Pi neqrep, el+\Pi neqrep, dis\u22480$) with *Pe*. The addition of attractive forces suppresses this growth in $\Pi neqB$ as forcing grows strong; when attractions are weak, the attractive osmotic pressure $\Pi neqattr$ contributes negatively to $\Pi neq$, counteracting the growth in $\Pi neqB$. However, as attractions grow strong, they become the dominant contribution to $B2*$ and qualitatively different behavior arises: the sum of entropic and interparticle contributions approaches a constant value at high *Pe*. When $B2*\u226a\u22121$, this limiting value approaches the *negative* of the equilibrium osmotic pressure, i.e.,

Overall, strong flows break attractive bonds, removing the source of negative interparticle osmotic pressure and causing the attractive osmotic pressure to reach a non-negative plateau, but for intermediate-flow strengths, hydrodynamic enhancement of attractive bonds leads to negative hydrodynamic osmotic pressure, which causes the flow-induced reduction in osmotic pressure that shifts the critical-point isobar upward in Fig. 13.

### APPENDIX H: OSMOTIC PRESSURE MAPS

The strength of hydrodynamic interactions can be varied from weak to moderate to strong, and here, we select several finite values of the repulsion range *κ* that sets this strength and plot its influence on the osmotic pressure map, Fig. 24. The primary result is that once the repulsion range is smaller than of order of the particles’ size, the critical isobar changes qualitatively. These maps were generated from the particle-phase osmotic pressure, $\Pi P\u2261\Pi tot/nakT\u22121$, following Eq. (36) from Sec. IV C. The most dramatic contrast is between the limiting cases of no hydrodynamics and strong hydrodynamics [Figs. 24(a) and 24(d)], as described in Sec. IV C. Figures 24(b) and 24(c) show intermediate repulsion ranges of *κ* = 1 and *κ* = 0.1, respectively. When *κ* = 1, the evolution of nonequilibrium osmotic pressure with attractions and forcing strength is nearly identical to the no-hydrodynamic case in Fig. 24(a). Upon shortening the repulsion range to *κ* = 0.1, however [Fig. 24(c)], the negative flow-induced osmotic pressure becomes visually apparent; an isotherm drawn along the Vliegenthart–Lekkerkerker^{8} boundary (the critical curve) curves upward near *Pe* = 10. The most dramatic effects of hydrodynamic interactions arise as the repulsion range is reduced below *κ* ≤ 0.1, such as the flow-induced osmotic pressure minimum—the upward shift of the critical-point Π boundary shown in Fig. 13(b)—and the increase in forcing strength *Pe* required to raise osmotic pressure above zero in strongly attractive suspensions—the rightward expansion of the blue, negative-pressure region in Fig. 24(d).