Inhomogeneous steady shear dynamics of a three-body colloidal gel former

Gelation in soft matter is a complex and important phenomenon that has many practical applications ranging from the use in household materials to advanced technological processes.^1–6 A common property of gels is their ability to sustain weak external stresses due to the formation of persistent long-range network structures. From a microscopic point of view, it is the nontrivial correlation of particles arising from their internal interactions that gives gels their characteristic mechanical response.^7–9 

However, the route to the generation of the network topology can be diverse.¹⁰ One common path to the gelation of colloids involves the crossing of a liquid–gas spinodal, e.g., by a sudden quench in temperature and the subsequent dynamical arrest of heterogeneous dense regions. This arrested spinodal decomposition is non-equilibrium in nature, and it bears similarity to the glass transition although it is driven by interparticle attraction rather than by repulsion.^11–13 On the other hand, an equilibrium route to gel formation lies open by careful choice of the interparticle interactions in order to prevent macroscopic liquid–gas phase separation and to favor instead the local arrangement of particles into interconnected clusters or chains. In this spirit, a multitude of interaction potentials have been investigated, which incorporate, for example, limitation of particle connectivity,^14–17 competing short-range attraction and long-range repulsion,^18,19 and anisotropy^20–22 via, e.g., “patchy” interaction sites.^23–25 The liquid–gas spinodal can sometimes be pushed to very low temperatures and densities, which enables large parts of the phase diagram to be governed by the percolation into dilute networks, as in so called “empty liquids.”^26–28 

A further class of particle models for colloidal gels that the present work focuses on is based on the inclusion of three-body interactions to the interparticle interaction potential, which consists otherwise only of isotropic pair-interactions.^29–36 It has been shown that an appropriate choice of the three-body term reproduces the distinctive network topology³⁷ as well as the characteristic non-linear response to homogeneous shear, including strain hardening and yielding.³³ Especially under external load, the dynamics of such a gel can be intricate, e.g., exhibiting cooperative restructuring of particle bonds³¹ and shear banding.³³ 

While many studies have considered the response of gels to a linear shear profile up to their breaking point, not much is known about their viscous flow behavior in inhomogeneous external force fields. However, it can be expected that the intrinsic features of a gel former, such as the tendency of particles to percolate, have substantial ramifications in such out-of-equilibrium scenarios. Specifically, one is tempted to assume that some of the genuine non-equilibrium effects^38,39 already reported for simple fluids (such as shear migration) might even be amplified by additional three-body interactions.

In this work, we show that gels modeled via a modified Stillinger–Weber⁴⁰ potential with a preferred three-body angle of 180° as proposed by Saw, Ellegaard, Kob, and Sastry (SEKS)^29,30 are indeed highly susceptible to these non-equilibrium effects when sheared by a sinusoidal external force profile. For this, we numerically investigate the behavior of the SEKS model with adaptive Brownian dynamics⁴¹ (adaptive BD), which is a stable and efficient method for the simulation of many-body systems governed by the overdamped Langevin equations of motion. We find that the properties of the emerging stationary state vary strongly with temperature and with the amplitude of the external force profile. Different behaviors occur in the shape of both the density and the internal force profiles as compared to simple fluids. In particular, we show that the superadiabatic (i.e., genuine out-of-equilibrium) contribution to the internal force is substantial in magnitude and that it is responsible for the structural and viscous behavior of the stationary shear flow. This is discussed from a microscopic point of view as well as in a coarse-grained fashion, where we use power functional theory^42,43 (PFT) to develop a quantitative model for the superadiabatic viscous force. Besides representing a generic situation, the sinusoidal shear flow profile could be seen as a toy model for a mesoscopic convection roll. Convection typically occurs in sedimentation as upward streams alternate with downward streams.⁴⁴ 

This work is structured as follows: In Sec. II A, the modified Stillinger–Weber potential as well as details for its efficient computation is given. The adaptive BD method and its advantages for our non-equilibrium simulations are laid out in Sec. II B. In Sec. II C, the protocol for the simulation of the stationary flow state is described. In Secs. III A and III B, we show one-body profiles of the density as well as the parallel and perpendicular components of the internal force for a range of simulation parameters and discuss their behavior and interplay. An analogous interpretation on the level of internal stresses is given in  Appendix A. In  Appendix B, we showcase results for different values of the three-body angle of the Stillinger–Weber potential, and in  Appendix C, the unusual non-equilibrium response of the gel is contrasted with numerical results for the simple Lennard-Jones fluid. In Sec. III C, the description of superadiabatic forces with PFT is illustrated and the results are compared with those from simulation. We conclude in Sec. IV and give an outlook to the investigation of further dynamical phenomena observed in our simulations and to a more extensive analysis with PFT.

The Stillinger–Weber potential⁴⁰ has originally been used for the simulation of solid and liquid silicon, and it has since been optimized and adapted to other particle types.^45,46 The interparticle interactions consist of a two-body potential u₂(r) that models both isotropic attraction and repulsion depending on the distance r between two particles as well as a three-body contribution u₃(r, r′, Θ). This three-body term imposes an energetically favorable angle Θ for three particles where a central particle is separated by the pairwise distances r and r′ to two other particles. The directionality of internal interactions is therefore only realized via u₃. Crucially, there is no need to explicitly incorporate orientational degrees of freedom, which is an advantage both in simulations as well as in a theoretical treatment.

In total, the internal energy potential possesses the form

with

for a certain particle configuration r^N = {r⁽ⁱ⁾, …, r^(N)} of the many-body system with N particles.

The parameters p, q, A, B, a, γ, λ, and Θ₀ can be tuned to alter the shape of the potential. A choice for these quantities, which is used in the present work and varies in some aspects from the one used originally by Stillinger and Weber,⁴⁰ is given in Table I. In particular, following previous studies of Saw et al.,^29,30 we tune Θ₀ to obtain a gel former, which is described in more detail in the following. The formulation in Eqs. (2) and (3) refrains from using absolute units and only involves intrinsic energy (ϵ) and length (σ) scales. In an overdamped system with friction coefficient ζ, all physical quantities can therefore be expressed in a reduced form.

We note that the parameter a sets the cutoff distance since both u₂(r) and u₃(r, r′, Θ) as well as their gradients vanish smoothly for r → aσ and r′ → aσ. The potential is therefore inherently short-ranged (cf. the small value of a in Table I). This is a favorable property for the treatment in computer simulations since it enables the use of neighbor-tracking algorithms to avoid superfluous evaluations for particles beyond the cutoff distance, which substantially reduces the computational cost in large systems.

The parameter Θ₀ in the three-body term u₃(r, r′, Θ) sets the preferred angle of a certain particle triplet (note that u₃ vanishes for Θ = Θ₀ and that it is otherwise strictly positive for particles within the cutoff distance). Most commonly, as discussed in the following, tetrahedral configurations are desired for which one chooses cos Θ₀ = −1/3. The strength of the three-body interaction term is adjusted via λ, which is often referred to as the tetrahedrality⁴⁷ for the above choice of Θ₀.

A further computational optimization is employed, which makes use of the concrete structure of u₃(r, r′, Θ) as given in Eq. (3). Via a rewriting of the three-body sum and the introduction of accumulation variables, an evaluation of the total energy and of all particle forces is possible by only iterating twice over all interacting particle pairs. In contrast, a naive implementation would require an iteration over particle triplets. For details of this exact reformulation, which leads to a substantial speedup in our simulations,⁴⁸ consult Ref. 30.

In summary, the versatility and computational efficacy of the Stillinger–Weber potential make it applicable to a wide range of problems. An important example, which conveys its use as an effective interaction potential for more complex particle types, is the monatomic water model of Molinero and Moore.⁴⁷ It has been shown by these authors that thermodynamic and structural properties of water (e.g., for the study of interfacial phenomena⁴⁹) can be captured accurately by this model via an appropriate choice of the absolute values of ϵ and σ as well as the tetrahedrality λ. By comparison with the melting temperature of water, they determined an optimal value of λ = 23.15, which lies between the respective tetrahedralities of silicon and carbon and which is adopted in our simulations.

While Eq. (3) has initially been conceptualized as a model for tetrahedrally coordinated particles, it is entirely conceivable to alter the preferred three-body angle Θ₀. A variation of Θ₀ has significant consequences for the spatial correlations of the fluid since the formation of droplets might become energetically unfavorable and the self-assembly into interconnected chains that form open networks is enforced. Therefore, the careful choice of the values of Θ₀ and λ is a means to reduce the effective valency and to suppress the liquid–gas phase transition, making the Stillinger–Weber potential (1) a suitable model for colloidal gels. In the following, we set Θ₀ = 180° although other values of Θ₀ have been shown to support gelation as well, e.g., as reported in Refs. 29 and 30, where a detailed investigation of the phase diagram and percolation behavior was carried out for various choices of λ and Θ₀. (In  Appendix B, illustrative results are presented for lower values of Θ₀, which shows that its precise value has little impact on the sheared steady state as long as network formation can occur.) It is worth noting that gelation has also been investigated for other choices of two- and three-body interaction terms u₂ and u₃ apart from those given in Eqs. (2) and (3); see, e.g., Refs. 31–36. For instance, to yield stronger angular rigidity, the cosine difference in Eq. (3) can been exponentiated.^31–33 

An important property of gels is their mechanical response to externally imposed strain. As particle bonds within the network are capable of sustaining substantial forces and torques without breaking, a gel exhibits elastic behavior before stiffening⁵⁰ as well as yielding at intermediate and large shear strain due to bending and breaking of bonds, respectively.³³ Numerically, these results can be obtained, e.g., by performing a linear deformation of the simulation box and measuring the stress tensor. When using non-equilibrium molecular dynamics, adequate thermostatting is required,^51,52 which is not straightforward if spatially inhomogeneous deformations are considered. This is even more problematic if cause and effect are reversed, and an external force profile is applied, which generates a macroscopic net flow that is, hence, not known a priori. We circumvent these issues by considering overdamped dynamics, where thermostatting is intrinsic.

Furthermore, an advanced numerical integration scheme known as adaptive BD⁴¹ is applied, which improves upon conventional BD simulations as described in the following. We consider the overdamped Langevin equations

i = 1, …, N, as the relevant equations of motion to obtain particle trajectories r^N(t) in our system consisting of N identical particles. Here, f⁽ⁱ⁾(r^N(t)) is the total force acting on particle i, which can be split into external and internal contributions, $fext(i)(r(i)(t))$ and $fint(i)(rN(t))=−∇iU(rN(t))$⁠, respectively. The friction coefficient ζ is the same for each (identical) particle, and the over dot denotes a time derivative. The vectors R⁽ⁱ⁾(t), i = 1, …, N, are Gaussian distributed and must therefore satisfy ⟨R⁽ⁱ⁾(t)⟩ = 0 and ⟨R⁽ⁱ⁾(t)R^(j)(t′)⟩ = Iδ_ijδ(t − t′). Here, the angular brackets denote an average over realizations of the random process, I is the 3 × 3-unit-matrix, δ_ij is the Kronecker delta, and δ(⋅) is the Dirac delta function. Thermostatting irrespective of applied external forces and any (possibly inhomogeneous) net flow is inherent in overdamped Brownian dynamics as the temperature-dependent prefactor of the Gaussian random vectors in Eq. (4) determines the average magnitude of the random displacements.

Because Eq. (4) is a set of coupled stochastic differential equations, its numerical treatment requires particular care. Specifically, the use of the Euler–Maruyama method,⁵³ which is usually employed in conventional BD simulations, has serious drawbacks regarding both its stability and accuracy. This is primarily due to using a constant timestep interval Δt, which may lead to faulty particle displacements and erroneous force evaluations when particle collisions are not resolved with the required precision (i.e., with a small enough Δt).

Within adaptive BD,⁴¹ the automatic choice of an appropriate timestep length Δt_k is ensured in each iteration k → k + 1 by the evaluation of an embedded Heun–Euler integrator,

which yields two estimates $r̄k+1N$ and $rk+1N$ for the new particle positions at time t_k + Δt_k. If large discrepancies of $r̄k+1N$ and $rk+1N$ are detected, the timestep Δt_k is reduced and the step k → k + 1 is retried. In such a case of a rejected trial step, one must carefully choose appropriate discrete random increments $Rk(i)$ to retain the Gaussian nature of the target random process R⁽ⁱ⁾(t). For this, adaptive BD utilizes Rejection Sampling with Memory (RSwM),⁵⁴ which is an efficient algorithm to counteract the rejection of previously drawn random increments. RSwM, hence, guarantees the correct generation of a specified random process. With the numerical treatment of Eq. (4) via adaptive BD, stable and accurate long-time simulations of overdamped many-body systems are possible in equilibrium but also under extreme non-equilibrium conditions, such as when driving the system with a large external force f_ext(r). A detailed description of adaptive BD is given in Ref. 41.

Hydrodynamic interactions that are mediated by the implicit solvent are neglected in the equations of motion (4). From a computational standpoint, the performance of the simulations turns out to be crucial to obtain accurate results for the quantities of interest, as described in Sec. II C. The numerical treatment of hydrodynamic interactions, which requires both the evaluation of long-ranged forces as well as the generation of appropriately correlated random displacements (e.g., via a Cholesky decomposition of the diffusion matrix⁵⁵), would have a significant impact on the required computational effort. Additionally, within adaptive BD, the Heun–Euler pair (5) and (6) is conceived to handle only additive noise in the underlying stochastic differential equation. Instead of the Heun method (6), one would have to resort to an integration scheme with a sufficient strong order of convergence for general noise terms.⁵³ Moreover, from a physical point of view, we argue that the omission of hydrodynamic interactions simplifies the analysis of the results in Sec. III as all observations are ensured to stem solely from the properties of the Stillinger–Weber particle model. In particular, we show that its ability to form networks is the crucial mechanism that causes the reported out-of-equilibrium response. As hydrodynamic interactions tend to support anisotropic coagulation and transient network states,^56–58 we expect no significant qualitative change in the reported observations.

In this work, we model an inhomogeneously sheared system by imposing an external force profile f_ext(z) that is parallel to the x axis and modulated in the z-direction. The x-component of the force field is sinusoidal with amplitude K such that it complies with the periodic boundary conditions of the cubic simulation box with side length L, i.e.,

With this choice, Lees–Edwards boundary conditions⁵⁹ as used in simulations with linear shear profiles⁶⁰ are not required because f_ext(r) and its derivatives are continuous at the periodic boundaries. The application of the time-independent but spatially inhomogeneous shear force (7) is not to be confused with time-dependent oscillatory shear.⁶¹ The considered external force constitutes the lowest-order Fourier mode within the simulation box, and it can, hence, be taken as a generic model for experimentally relevant scenarios, as occur, e.g., in convection⁴⁴ or when inducing inhomogeneous forces with a laser tweezer.

The simulation procedure is as follows: We set L = 30σ and initialize N = 1000 particles on a regular lattice, which yields a mean number density of ρ_b ≈ 0.037σ⁻³. This configuration is randomized for a short time (10⁴ steps) at a high temperature of k_BT = 10ϵ using the adaptive BD method before instantaneously reducing the temperature to the desired value and imposing the sinusoidal external force profile (7). At this point, no particle bonds have formed yet and a flow in the x-direction sets in immediately. From here, the actual production run begins and the respective observables are sampled, which is described in more detail in the following. Due to the nature of the external force profile (7), the system retains translational invariance in the x–y-plane and forms a flow channel in the upper and lower half of the simulation box. Since a transient from the randomized particle distribution into this stationary flow occurs initially, we partition the sampling of the production run into consecutive sections of 10⁶ steps. Thus, the sections where a stationary flow has not been reached yet can be discarded, and the remaining ones are averaged over. During individual runs, asymmetric channel populations that persist for a long time are observed. Rather than performing longer simulation runs to yield better time-averages, we average over ∼50 distinct realizations of a simulation until symmetric profiles are obtained. The typical simulation time of the stationary flow in each individual run is then of the order of 10⁵τ with the Brownian timescale τ = σ²ζ/ϵ.

Using this protocol, a range of external modulation amplitudes K and temperatures T is investigated. For each set of parameters, we obtain the density profile ρ(r) as well as the force density profile F(r) from sampling of the density operator

and force density operator

respectively. Thus, $ρ(r)=⟨ρ̂(r)⟩$ and $F(r)=⟨F̂(r)⟩$⁠, where angular brackets denote an average over configurations of the stationary flow state obtained according to the above simulation procedure. Specifically, for the force density profile, we focus on its internal contribution

to better reveal how the stationary state is stabilized by the internal interaction (1). The internal force density profiles are then normalized by the density to acquire the internal force profile

A sufficiently large number of samples is necessary to yield accurate results for the internal force profile as its convergence is slower than that of the density profile.⁶² To investigate possible finite-size effects, which might occur in gels specifically due to their long-range effective correlations, we have conducted additional simulations where the side length L of the box has been doubled while extending the external potential (7) in the z-direction by an additional shear period. No significant impact was found on the behavior of the sheared system compared to the results shown in Sec. III for the original choice of L.

A sketch of the system and of the flow velocity profile resulting from the applied shear force (7) is depicted in Fig. 1, where we also show characteristic snapshots of the quiescent and of the sheared gel. Additionally, the spatial variations of the one-body profiles are illustrated, and we indicate locally by arrows the directions of the one-body force contributions. Actual simulation results are presented and analyzed in the following, and we highlight the labels of the one-body profiles in subsequent figures according to the colors used in Fig. 1.

For certain state points and values of the amplitude of the external force, large variations in the one-body profiles of density ρ(r) = ρ(z) and internal force f_int(r) = f_int(z) are observed while the system retains translational symmetry in the x- and y-directions. To investigate the onset and origin of these inhomogeneities, we first vary the temperature T and maintain a large constant amplitude K = 5ϵ/σ of the external force profile. The results are shown in Fig. 2.

The one-body profiles remain almost featureless for k_BT = 0.3ϵ. At k_BT = 0.2ϵ, an inhomogeneous structure begins to appear in the internal force profiles f_int,x(z) and f_int,z(z). This becomes more clearly visible as variations of the density profile ρ(z) from its bulk value for k_BT = 0.15ϵ. For k_BT = 0.1ϵ, remarkable modulations occur in all three quantities with spatial density variations of the order of the mean bulk density ρ_b itself.

The emergence of structural features in the one-body profiles when decreasing temperature is rapid and continuous. The spatial modulations are significantly stronger than those observed in (non-percolated) simple fluids,³⁸ and we illustrate this in  Appendix C via a comparison to results for the dilute Lennard-Jones fluid. Therefore, this effect can be linked to the percolation transition in equilibrium, which sets in at similar thermodynamic state points for the considered particle model.⁶³ We support this reasoning by an investigation of the cluster size distribution C(n), which gives the probability of finding a random particle in a cluster of size n. As is standard, we define the agglomeration of particles into clusters to be transitive with two particles belonging to the same cluster if their distance is below the cutoff distance aσ of the interparticle potential. In Fig. 3, C(n) is shown for varying temperatures in a system sheared according to Eq. (7) with K = 5ϵ/σ. One recognizes that the mean cluster size grows with decreasing temperature and that clusters span up to half of the system for k_BT = 0.1ϵ.

Additionally, to better reveal the internal structure of the clusters, we monitor the probabilities of the coordination numbers P_n, i.e., the proportion of particles having n neighboring particles within the cutoff distance aσ. The behavior of the coordination numbers n = 0, 1, 2, 3 is shown in Fig. 4 as a function of inverse temperature. It is apparent that for low temperatures, the structure of the network is dominated by particle chains. Branching still occurs, which interconnects the chains within the flow channels. This shows that even in strongly sheared systems, microscopic correlations are dominated by the three-body contribution to the internal interaction potential (1).

While the formation of finite-size clusters can be understood as a relic of the equilibrium percolation transition, its effect on the concrete structure of ρ(z) and f_int(z) turns out to be substantial and can only be explained if genuine non-equilibrium dynamics are considered. In the following, the response of the three-body gel over a range of external force amplitudes K at constant (low) temperature k_BT = 0.1ϵ is investigated, whereby the system is driven further away from equilibrium with increasing K. The corresponding one-body profiles are shown in Fig. 5.

To rationalize the results, it is instructive to work on the level of forces and to consider the one-body force balance⁴³ 

The above relation is exact for arbitrary many-body Hamiltonians, which can be shown, e.g., via an integrating-out of the Smoluchowski equation or in equilibrium, where v(r) = 0, by an application of Noether’s theorem.⁶⁴ The external force f_ext(r) is imposed in our system via Eq. (7), and ρ(r) as well as f_int(r) is accessible from their microscopic definitions given in Eqs. (8)–(11). The term −k_BT∇  ln ρ(r) = f_id(r) on the right-hand side of Eq. (12) is the force arising from ideal-gas diffusion. Furthermore, the time dependence has been dropped as we consider a stationary state where v(r) = J(r)/ρ(r) is the time-independent one-body velocity, which can be obtained from the one-body current J(r). The current obeys the continuity equation ∂ρ(r, t)/∂t = −∇ · J(r, t), and it is therefore divergence-free in the present case since the density profile is stationary. We recall that J(r) is an average of the microscopic operator $Ĵ(r)=∑iv(i)δ(r−r(i))$⁠, and it is thus directly accessible from simulation.⁶⁵ 

We proceed similar to Ref. 39 and distinguish between adiabatic and superadiabatic contributions to the internal force profile f_int(r) = f_ad(r) + f_sup(r). The adiabatic force f_ad(r) is defined to be that of a reference equilibrium system, which is constructed to have the same density profile as the original non-equilibrium state. Only the superadiabatic part f_sup(r) consists of purely out-of-equilibrium forces, which, hence, determine both (inhomogeneous) structure and flow of a driven colloidal suspension. By specializing to our planar geometry, we now analyze the force profiles in Fig. 5 to determine adiabatic and superadiabatic contributions.

Parallel to the flow, the density remains homogeneous due to translational symmetry such that the x-component of its gradient vanishes. This also implies a vanishing adiabatic force f_ad,x(z) = 0. The respective component of the internal force therefore only consists of the superadiabatic contribution, i.e., f_int,x(z) = f_sup,x(z). Hence, the x-component of the force balance Eq. (12) simplifies to

which clarifies that f_sup,x(z) plays the role of a viscous force and that it is readily available via the simulation results for f_int,x(z) shown in Fig. 5.

In the z-direction, the density varies inhomogeneously, and the internal force therefore consists of both adiabatic and superadiabatic contributions. To distill f_sup,z(z) from the data of f_int,z(z) in Fig. 5 would require the construction and simulation of an appropriately chosen equilibrium system,^65,66 which will be considered in the future work. Nevertheless, due to the absence of driving and flow in the z-direction, i.e., f_ext,z(z) = 0 and v_z(z) = 0, the force balance along the z axis reduces to

Therefore, the non-equilibrium force component f_sup,z(z) is necessary to stabilize the density gradient, and it can thus be referred to as a structural superadiabatic force. Equation (14) also reveals that the internal force density is straightforwardly related to the derivative of the density profile due to F_int,z(z) = f_int,z(z)ρ(z) = k_BT∂ρ(z)/∂z, which can be utilized as a means to “force sample”^67–69 the density profile with a reduced variance. Additionally, an analogous description of viscous and structural effects on the level of internal stresses is given in  Appendix A.

In simple fluids, where the constituent particles only interact via an isotropic pair-potential, non-equilibrium viscous and structural forces have been reported to occur both in an analogous sinusoidal shear profile³⁸ and in more complex two-dimensional flows.³⁹ However, the emerging features of density and force profiles—while being measurable and conceptually important—are rather frugal especially in the quasi-one-dimensional case (cf.  Appendix C for results of the sheared Lennard-Jones fluid). The relative variation in density is comparatively small even for moderate external force, and particles consistently accumulate in regions of low shear rates, i.e., at the center of the flow channels. The superadiabatic forces possess a sinusoidal shape such that the structural force drives particles to the center of the channels. The viscous force is Stokes-like in a broad range of shear amplitudes, and it is always opposed to the flow direction. In the following, these observations are compared to the markedly different one-body profiles of the sheared three-body gel illustrated in Fig. 1. The results of the variation of K are shown in Fig. 5.

For small values of the external force amplitude K, the density is sinusoidal in shape but the amplitude is inverted as compared to the simple fluid scenario such that particles accumulate in regions of large velocity gradient. When K is increased, the density maxima shrink while the depletion at the center of the channels remains pronounced. For large values of K, we observe that particles now tend to flee the regions of high velocity gradient. However, the migration is not simply directed toward the center of the flow channels where the local shear rate vanishes, as would be the case in simple fluids. Instead, the density profile develops a double-peak and retains a depletion zone right at the location of maximum flow velocity. This behavior is reflected in the form of the internal force f_int,z(z), which progresses from a sinusoidal profile for small K to a rapidly varying quantity for large K, thereby promoting and maintaining the observed double-peak structure of ρ(z) within the flow channels. The purely superadiabatic viscous force f_int,x(z) counteracts partially the flow for low to intermediate K similar to the behavior found in simple fluids. For large external force amplitudes, however, f_int,x(z) locally acts in the same direction as the flow velocity at the sides of the channels, which is anomalous phenomenology for a viscous force.

The striking signal in both structural and viscous forces can be explained as a consequence of the three-body interaction (1). As the system is weakly sheared, particles can still percolate into a large network for the chosen temperature. With increasing K, bonds are first broken in regions of maximum external force such that particles become mobile and evade these regions—thus, a density depletion zone develops. At even larger K, the formation of an extensive network cannot be maintained and bonds break and dynamically rejoin across the whole system. However, driven by the three-body term in Eq. (1), particles still tend to develop finite-size chains, which then align parallel to the flow direction. The mobility of the individual chains enables the migration to regions of low velocity gradient, and the density profile, hence, inverts. Within the flow channels, the chains organize into two lanes that are slightly offset from the center and thus lead to a double-peak structure in ρ(z). This is because their alignment parallel to the flow is driven by inhomogeneous shear rate, and it can therefore only occur if ∂v_x(z)/∂z ≠ 0. Yet, at the extrema of the external force profile, the gradient of the resulting flow vanishes and particle bonds are not aligned. This explains the spatial offset of the chain formation to regions of finite velocity gradient; cf. Fig. 5. The arrangement of particles into aligned chains also clarifies the anomalous behavior of the viscous force f_int,x(z) that is encountered in this case and that can, hence, be understood as a dynamical “drag-along.” In summary, the inclusion of three-body terms in the interaction potential greatly affects the response of colloidal suspensions to inhomogeneous shear and results in collective effects, which influence and amplify structural and viscous forces.

We next turn to a theoretical description of the simulation results with PFT^42,43 and give a brief summary of its core concepts in the following. PFT is based on an exact variational principle that reproduces the time-dependent force balance equation

Thereby, the nontrivial contributions f_ad(r, t) and f_sup(r, t), which together constitute the internal force profile f_int(r, t), are made accessible via universal generating functionals of the density profile ρ(r, t) and current profile J(r, t). Together with the external and diffusive forces (right-hand side), they are balanced by the friction of the overdamped system (left-hand side).

More precisely, the adiabatic force f_ad(r, t) incorporates the functional derivative of the intrinsic excess Helmholtz free energy F_exc[ρ],

where brackets denote functional dependencies. If one uses Eq. (16) in Eq. (15) and neglects f_sup(r, t), classical dynamical density functional theory (DDFT)⁷⁰ is recovered as an uncontrolled approximation. In the sheared three-body gel, as was shown in Secs. III A and B, the dynamics are governed by genuine out-of-equilibrium effects. Being a purely adiabatic theory by construction, DDFT is strictly unable to reproduce or describe the observed behavior in our system.⁷¹ 

Instead, in order to go beyond an adiabatic description, superadiabatic forces f_sup(r, t) have to be taken into account. Within PFT, this is made possible by functional differentiation of the superadiabatic excess power functional P_exc[ρ, J],

The force balance Eq. (15) can then be written as

and it involves both adiabatic and superadiabatic interparticle forces as systematically generated via the respective functionals.

Therefore, if P_exc[ρ, J] and F_exc[ρ] are known, PFT enables the dynamical description of a system subjected to an arbitrary external force profile f_ext(r) via Eq. (18) and the continuity equation. This reformulation, which reduces the many-body problem to a variational principle on one-body quantities, is exact, in principle. Crucially, both P_exc[ρ, J] as well as F_exc[ρ] are intrinsic functionals that depend only on internal interactions and further intrinsic properties of the system (e.g., temperature and density), but not on the externally applied force profile f_ext(r, t). In practice, for a certain interparticle interaction potential, approximations for P_exc[ρ, J] and F_exc[ρ] must be found, which poses a nontrivial problem. For P_exc[ρ, J], the functional dependence will, in general, be non-local both in space and in time (i.e., non-Markovian) as the history of ρ(r, t) and J(r, t) has to be considered to obtain an accurate dynamical theory for time-dependent problems.

In the following, we use the framework of PFT to develop a model that is capable of reproducing the found anomalous behavior of the viscous force profile in the sheared three-body gel. We focus on the viscous part because it is directly accessible in simulation (we recall that f_int,x(z) is purely superadiabatic), which, hence, simplifies the following considerations. Recall also that a stationary state is considered, which implies ∇ · v(r) = 0 due to the continuity equation and the chosen geometry. Since ρ(r) is time-independent, we perform a change of variables to formulate P_exc[ρ, v] as a functional of the velocity profile and use δ/δJ(r) = ρ(r)⁻¹δ/δv(r). To yield an approximate explicit expression for P_exc[ρ, v], a semi-local velocity gradient expansion⁷² is assumed. Due to being in a stationary state, we can further specialize to a Markovian model such that

with a suitable integrand ϕ(ρ(r), ∇v(r)).

As in Refs. 38 and 39, we perform the general expansion up to second order in ∇v(r) to obtain an expression for the integrand ϕ(ρ(r), ∇v(r)). Assuming a local dependence in space and imposing rotational invariance, this expression can be reduced to

where η is the coefficient of the superadiabatic viscous response. This coefficient depends on intrinsic properties, such as interparticle potential, density, and temperature, but it crucially is independent of the imposed external force profile. Furthermore, the value of η only alters the magnitude of the superadiabatic force profile resulting from Eq. (20) with its shape being fully determined by the forms of ρ(r) and v(r).

The functional minimization (17) of Eq. (19) with the model integrand (20) results in the superadiabatic force profile

where the right-hand side can be evaluated for the sheared gel. For this, we approximate v(r) ≈ f_ext(r)/ζ since the magnitude of the viscous force is small compared to ζv(r) and take the density profile from the simulations as the input. The x-component of Eq. (21) then yields the viscous force, which can be compared to the actual simulation data f_int,x(z) as shown in Fig. 5. To obtain a quantitative comparison, the value of the transport coefficient η is fitted to match the magnitude of the simulation results universally for the considered shear amplitudes. The superadiabatic force profiles for the viscous force within this PFT description are shown in Fig. 6, where a value of η = 5 has been used for the viscous coefficient for all considered values of K.

It is apparent that f_sup,x(z) displays a double-peak within the flow channels and therefore differs from the simulation results, where only a single peak is observed. However, this inaccuracy is not surprising since the model functional (19) and (20) is obtained merely by an expansion in gradients of the velocity profile. The density enters the functional only locally, and the model is thus expected to fail in regions where higher derivatives (e.g., the curvature) of ρ(r) are significant, such as in the center of the flow channels. To achieve better results in these regions, the integrand (20) of P_exc[ρ, v] could be augmented by an expansion in ρ(r), which will be considered in future work.

In between the flow channels, the viscous force profiles obtained from Eqs. (19) and (20) match the simulation results across the range of investigated shear amplitudes K. Particularly, the anomalous change in the sign of the viscous force, which we attribute to a dynamical drag-along of particles, is captured by the PFT model as well, and it shows the same K-dependent behavior as in the simulation. The successful reproduction of this phenomenon exemplifies that even simple model functionals for the excess power are capable of resolving nontrivial superadiabatic effects and that PFT is a concise framework for their systematic investigation.

In this work, we have studied the behavior of a colloidal gel modeled by the Stillinger–Weber potential (1), where the three-body interaction (3) has been modified similar to Refs. 29 and 30. The gel is subjected to a sinusoidal external shear profile. For the numerical investigation, we have utilized adaptive BD,⁴¹ which facilitates to carry out efficient and stable long-time simulation runs to accurately obtain the density and internal force profile in the stationary flow state. Markedly different behaviors have been encountered depending on the chosen temperature T and the amplitude K of the external force profile.

The simulations over a range of temperatures revealed that the effect of the equilibrium percolation transition—which leads to the formation of an extended and dilute network under quiescent bulk conditions—transfers to situations far from equilibrium. Thus, while a system-spanning network is not formed for sufficiently strong inhomogeneous shear, the local arrangement of particles into finite-size chains is still viable, which we have shown via the cluster size distribution C(n). An investigation of the probabilities of the coordination numbers P_n revealed that the clusters are dominated by chains, which interconnect via branching. This clustering effect is crucial to describe the emergence of structure in the density profile ρ(z) and in the parallel and perpendicular component of the internal force f_int(z) with respect to the flow direction. Note that the global temperature acts as a control parameter for the network formation in our system. In depletion-induced gels, a similar effect could be achieved by a variation of the concentration of the depletion agent to tailor the effective attraction between colloids.^11–13 

For an in-depth analysis, we have further split the internal force into adiabatic and superadiabatic contributions with the latter being the driving mechanism for genuine out-of-equilibrium effects. Due to the chosen planar geometry, the parallel component of the internal force could be associated directly with a superadiabatic viscous force. The perpendicular component consists of both adiabatic and superadiabatic contributions instead, where the latter is needed to stabilize the emerging density inhomogeneity.

When comparing the found results of the three-body gel with known observations of colloids consisting of simpler particle types,^38,39 we found anomalous behavior for both viscous and structural effects. This could be attributed to be a direct consequence of the internal three-body contributions. The emerging density modulation is much larger in magnitude and shows a richer phenomenology than in simple fluids, as we have illustrated via a comparison to the Lennard-Jones fluid in  Appendix C. In particular, the accumulation of particles can occur both in regions of high and low velocity gradient depending on the applied external force. For large amplitudes of the latter, the formation of particle chains occurs within a double-lane near the center of the flow channels. The superadiabatic viscous force, which generally opposes the flow direction in simple fluids, has been shown here to flip its usual counteracting direction for large K in some regions of the channels. We deduced this “drag-along” to be another consequence of the formation of particle chains. Therefore, in both components of the internal force profile, collective effects are involved, which substantially amplify the non-equilibrium response of the system. As we have shown, colloidal gels are very susceptible to out-of-equilibrium phenomena, and they can, hence, be taken as a prototypical model for future study.

By utilizing PFT, a possible route to a coarse-grained description of the found results was given. This was exemplified for the viscous force profile, where we have shown that a simple excess power functional suffices to reproduce the simulation results and capture the anomalous drag-along in the three-body gel. In the future work, more sophisticated model functionals will be investigated in order to alleviate some deficiencies of this simple description. Building upon the found results, a similar analysis of the structural force profile will be considered. This requires, however, the construction of an equilibrium reference state to perform the splitting of the respective internal force component into adiabatic and superadiabatic contributions.

In the conducted simulations, it was observed that asymmetric channel populations that persist over long time scales occur especially for intermediate values of the shear amplitude. Hence, another objective for future work is a study of their statistics and stability, possibly being indicative of a dynamical phase transition as reported already in dense colloidal suspensions of simpler particles that exhibit flow-induced ordering or layering phenomena.^73,74 Further interesting research could incorporate a variation of other parameters of the Stillinger–Weber potential besides Θ₀ to study their impact on the response of the driven system. This is especially important from a practical perspective as the tuning of microscopic interactions to yield desired material properties is a central concept of material science, which has also been applied to colloidal gels under shear.⁷⁵ For a quantitative prediction, hydrodynamic interactions might become relevant, and it would be useful to augment adaptive BD in this regard, possibly accompanied by efficient evaluation schemes of then correlated random increments.^76,77 Additionally, going beyond the steady state and investigating time-dependent situations, such as transients in a switching protocol of the external force,⁷⁸ could reveal the nature of non-equilibrium memory effects. This is especially interesting from the view point of PFT as memory kernels can be directly incorporated in the theory, such that time-dependent phenomena may provide further assistance in the development of accurate functionals.

See the supplementary material for an animation of the steady shear flow of the three-body gel at temperature k_BT = 0.1ϵ and shear amplitude K = 5ϵ/σ.

We thank Tobias Eckert and Matthias Fuchs for useful comments. This work was supported by the German Research Foundation (DFG) via Project No. 436306241.

The authors have no conflicts to disclose.

Florian Sammüller: Formal analysis (lead); Methodology (equal); Software (lead); Writing – original draft (lead). Daniel de las Heras: Conceptualization (equal); Methodology (equal); Writing – review & editing (equal). Matthias Schmidt: Conceptualization (equal); Methodology (equal); Project administration (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Instead of working on the level of the force balance Eq. (12) directly, one can consider a similar decomposition of the stress tensor as a generator of the respective force profiles. For this, we use the definition

of the total stress tensor σ(r, t). To identify its internal contribution, we multiply Eq. (12) by the density profile ρ(r, t), which yields the force density balance

Insertion of Eq. (A2) into Eq. (A1) and an analogous splitting then gives rise to the definition

for the internal stress tensor σ_int(r, t).

In the considered stationary state, the time dependence can be dropped. Obtaining σ_int(r) from the sampled force density profile F_int(r) requires an integration of its spatial components according to Eq. (A3). Pressure-like contributions (corresponding to integration constants) are not accessible from the force density profiles alone and require further suitable measurements in simulation.^79,80 (The standard Irving–Kirkwood⁸¹ treatment is only valid for pair-potentials.) We omit such constants in the following and only consider relative inhomogeneities of the internal stress. Additionally, a non-unique⁸² divergence-free part of σ_int(r) remains undetermined from the integration of Eq. (A3) and is set to zero.

We specialize to the planar geometry of our system, which enables a straightforward integration to obtain two relevant components of σ_int(z) via

Analogous to Fig. 5, where the x- and z-component of the internal force is depicted, we show results for the components σ_int,zz(z) and σ_int,zx(z) of the internal stress tensor as obtained by Eqs. (A4) and (A5) in Fig. 7. Here, the integration constants were chosen such that σ_int,zz(z) vanishes at the boundaries of the box and σ_int,zx(z) is anti-symmetric under motion reversal (v(r) → −v(r)).

It is observed that σ_int,zz(z) reproduces the shape of the density profile, which is consistent with the considerations in the main text; cf. Eq. (14). For σ_int,zx(z), a sinusoidal shape is obtained at low shear amplitudes. When increasing K, the zx-component of the internal stress tensor develops a secondary structure. This is indicative of the non-linear response of a colloidal gel to applied shear, which manifests itself for inhomogeneous shear in an anomalous behavior of the viscous contribution.

In Fig. 8, we show illustrative results of the sheared three-body gel for different values of the preferred three-body angle Θ₀. For lower values of Θ₀, it is increasingly difficult to obtain symmetric profiles. We choose Θ₀ = 150° as the lowest value to keep away from the liquid–gas binodal and to prevent the formation of droplets within the flow channels, which hinder an accurate sampling. It is apparent from the results that the choice of Θ₀ = 180° in the main text is not artificial and that similar behavior can be achieved also for lower values of Θ₀ as long as gelation is enforced. When decreasing Θ₀, one even observes larger local forces at the sides of the flow channels [cf. f_int,z(z) and f_int,x(z) in Fig. 8] as the desorption of particle strands is enhanced due to the increased ability of branching. We refer to Refs. 29 and 30 for an investigation of the equilibrium behavior of the three-body gel for different values of the three-body angle Θ₀ and the three-body interaction strength λ.

For comparison, we show the behavior of the truncated Lennard-Jones fluid under an analogous shear protocol as for the three-body gel. The Lennard-Jones interaction potential only consists of the radially isotropic pairwise contribution

with the cutoff distance r_c = 2.5σ, and it can, hence, be taken as an example of a simple fluid or colloidal suspension.

In Fig. 9, the density profile as well as the parallel and perpendicular contribution of the internal force profile is shown for a temperature of k_BT = 1.5ϵ and for various (large) shear amplitudes K. All other system parameters are adopted from the simulations of the sheared gel, which yields the same low mean density of ρ_b ≈ 0.037σ⁻³. One recognizes that the superadiabatic response of the Lennard-Jones fluid differs starkly from that of the three-body gel; cf. 5. The density inhomogeneity of the simple liquid is orders of magnitude smaller and possesses a sinusoidal shape that does not change qualitatively for different shear amplitudes. Note that despite driving the Lennard-Jones system with much stronger external forces, the onset of notable superadiabatic effects occurs only for sufficiently large inhomogeneous shear as opposed to the three-body gel, where a substantial density inhomogeneity develops also for low values of K. In particular, no inversion of the extrema in the density profile ρ(z) is observed, as was the case for the three-body gel when transitioning from low to high shear. The internal force components reflect this situation with both f_int,z(z) and f_int,x(z) being much smaller and showing less features than in the three-body gel. Especially for f_int,x(z), no anomalous drag-along is observed as the superadiabatic viscous force in the Lennard-Jones fluid always counteracts the flow.