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

We investigate the stationary flow of a colloidal gel under an inhomogeneous external shear force using adaptive Brownian dynamics simulations. The interparticle forces are derived from the Stillinger-Weber potential, where the three-body term is tuned to enable network formation and gelation in equilibrium. When subjected to the shear force field, the system develops remarkable modulations in the one-body density profile. Depending on the shear magnitude, particles accumulate either in quiescent regions or in the vicinity of maximum net flow, and we deduce this strong non-equilibrium response to be characteristic of the gel state. Studying the components of the internal force parallel and perpendicular to the flow direction reveals that the emerging flow and structure of the stationary state are driven by significant viscous and structural superadiabatic forces. Thereby, the magnitude and nature of the observed non-equilibrium phenomena differs from the corresponding behavior of simple fluids. We demonstrate that a simple power functional theory reproduces accurately the viscous force profile, giving a rationale of the complex dynamical behavior of the system.


I. INTRODUCTION
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][2][3][4][5][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 give gels their characteristic mechanical response [7][8][9] .
However, the route to the generation of the network topology can be diverse 10 .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 nonequilibrium in nature and it bears similarity to the glass transition, although it is driven by interparticle attraction rather than by repulsion [11][12][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][15][16][17] , competing short-range attraction and long-range repulsion 18,19 , and anisotropy [20][21][22] via, e.g., "patchy" interaction sites [23][24][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][27][28] .
a) Electronic mail: Matthias.Schmidt@uni-bayreuth.deA 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][30][31][32][33][34][35][36] .It has been shown that an appropriate choice of the three-body term reproduces the distinctive network topology 37 as well as the characteristic non-linear response to homogeneous shear, including strain hardening and yielding 33 .Especially under external load, the dynamics of such a gel can be intricate, e.g.exhibiting cooperative restructuring of particle bonds 31 and shear banding 33 .
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 40 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 41 (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 behavior occurs in the shape of both the density and the internal force profiles as compared to simple fluids.In particu-lar, we show that the superadiabatic (i.e.genuine out-ofequilibrium) 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 44 .
This work is structured as follows.In Sec.II A, the modified Stillinger-Weber potential as well as details for its efficient computation are 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 nonequilibrium 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.

A. Particle model
The Stillinger-Weber potential 40 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 2 (r) that models both isotropic attraction and repulsion depending on the distance r between two particles, as well as a three-body contribution u 3 (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 3 .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 for a certain particle configuration r N = {r (i) , . . ., r (N ) } of the many-body system with N particles.The parameters p, q, A, B, a, γ, λ and Θ 0 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 40 , is given in Table I.In particular, following previous works of Saw et al. 29,30 , we tune Θ 0 to obtain a gel former, which is described in more detail below.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 2 (r) and u 3 (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 propery 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 Θ 0 in the three-body term u 3 (r, r , Θ) sets the preferred angle of a certain particle triplet (note that u 3 vanishes for Θ = Θ 0 and that it is otherwise strictly positive for particles within the cutoff distance).Most commonly, as discussed below, tetrahedral configurations are desired, for which one chooses cos Θ 0 = −1/3.The strength of the three-body interaction term is adjusted via λ, which is often referred to as the tetrahedrality 47 for the above choice of Θ 0 .
A further computational optimization is employed, which makes use of the concrete structure of u 3 (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 48 , consult Ref. 30.
In summary, the versatility and computational efficacy of the Stillinger-Weber potential make it applica-ble 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 47 .It has been shown by these authors that thermodynamic and structural properties of water (e.g. for the study of interfacial phenomena 49 ) 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 Θ 0 .A variation of Θ 0 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 Θ 0 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 Θ 0 = 180 • , although other values of Θ 0 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 Θ 0 .(In Appendix B, illustrative results are presented for lower values of Θ 0 , 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 2 and u 3 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][32][33] .TABLE I.In eqs.( 2) and (3), we adopt the parameters p, q, A, B, a and γ of the original Stillinger-Weber 40 potential and choose the three-body strength λ as determined in Ref. 47.In accordance with Saw et al. 29,30 , a preferred three-body angle of Θ0 = 180 • then leads to the percolation of inter-connected chains, enabling colloidal gelation in equilibrium.

B. Adaptive Brownian dynamics
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 50 as well as yielding at intermediate and large shear strain due to bending and breaking of bonds respectively 33 .Numerically, these results can be obtained, e.g., by performing a linear deformation of the simulation box and measuring the stress tensor.When using nonequilibrium 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 41 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 (i) (r N (t)) is the total force acting on particle i, which can be split into external and internal contributions, , respectively.The friction coefficient ζ is the same for each (identical) particle and the over dot denotes a time derivative.The vectors R (i) (t), i = 1, . . ., N , are Gaussian distributed and must therefore satisfy 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 53 , 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 41 , 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 r(i) k+1 = r which yields two estimates rN k+1 and r N k+1 for the new particle positions at time t k + ∆t k .If large discrepancies of rN k+1 and r N k+1 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 R (i) k to retain the Gaussian nature of the target random process R (i) (t).For this, adaptive BD utilizes Rejection Sampling with Memory (RSwM) 54 , 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 nonequilibrium 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 the next section.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 55 ), 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 53 .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][57][58] , we expect no significant qualitative change in the reported observations.

C. Simulation protocol
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 59 as used in simulations with linear shear profiles 60 are not required, because f ext (r) as well as 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 timedependent oscillatory shear 61 .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 44 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σ −3 .This configuration is randomized for a short time (10 4 steps) at a high temperature of k B T = 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 below.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 respectively.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 6 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 approximately 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 in the order of 10 5 τ with the Brownian timescale τ = σ 2 ζ/ .
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 FIG. 1. Characteristic snapshots of an equilibrium gel (left) and a sheared gel in steady-state (right) a .The particles are colored according to the cluster to which they belong.The simulation box is a cube of side length L and the temperature is set to kBT = 0.1 .For the sheared state, an external force amplitude of K = 5 /σ is chosen.The forces acting on the sheared gel are schematically represented: A sinusoidal external force pointing along the x-direction drives the flow of particles, forming two flow channels with the velocity profile (gray) closely following the external force profile (7).In steady state, a superadiabatic viscous internal force (cyan) emerges that locally either opposes or supports the flow.A strong density modulation (orange) develops along the z-direction.The ideal gas diffusive force (olive) that tends to homogenize the density profile is balanced by an internal force (pink) along the z-direction, which incorporates both adiabatic and superadiabatic structural components.
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 62 .
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 the next section 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 onebody 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.

A. Variation of temperature
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 The density profile ρ(z) (a) as well as the component fint,z(z) (b) and fint,x(z) (c) of the internal force (7) is shown.A constant shear amplitude of K = 5 /σ is maintained and the temperature is varied with values of kBT / = 0.1, 0.15, 0.2, 0.3 (indicated by ticks on the color scale).While fint,x(z) acts parallel to the flow direction, fint,z(z) constitutes a force perpendicular to the flow that leads to the observed density inhomogeneity.This is illustrated by arrows, which accentuate in particular the alternating direction in both the parallel and the perpendicular internal force component for low temperature.The onset of structural inhomogeneities in the one-body profiles is continuous and occurs rapidly for decreasing T when the equilibrium percolation transition is encountered.
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 B T = 0.3 .At k B T = 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 vari- For low temperatures, individual particles (n = 0) as well as particle pairs (n = 1) are desorbed into the network, as both P0 and P1 decrease.The network structure is dominated by chains (P2 is large) and branching is still viable, as can be deduced from the moderate value of P3.
ations of the density profile ρ(z) from its bulk value for k B T = 0.15 .For k B T = 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 38 , 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 63 .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 temperature 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 B T = 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).

B. Variation of external force amplitude
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 B T = 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 43 The above relation is exact for arbitrary many-body Hamiltonians, which can be shown e.g.via an integratingout of the Smoluchowski equation or in equilibrium, where v(r) = 0, by an application of Noether's theorem 64 .The external force f ext (r) is imposed in our system via eq.( 7), and ρ(r) as well as f int (r) are accessible from their microscopic definitions given in eqs.( 8) to (11).The term −k B T ∇ ln ρ(r) = f id (r) on the right-hand side of eq. ( 12) is the force arising from ideal-gas diffusion.Further, 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) = i v (i) δ(r − r (i) ) and it is thus directly accessible from simulation 65 .
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-ofequilibrium 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 xcomponent 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 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.Eq. ( 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 B T ∂ρ(z)/∂z, which can be utilized as a means to "force sample" [67][68][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 38 and in more complex two-dimensional flows 39 .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 rate, 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 towards 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 doublepeak 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 threebody 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.

C. Power functional theory
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 (15) 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) 70 is recovered as an uncontrolled approximation.In the sheared three-body gel, as was shown in the previous section, the dynamics are govered 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 71 .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, 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) −1 δ/δv(r).To yield an approximate explicit expression for P exc [ρ, v], a semilocal velocity gradient expansion 72 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.Further, 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 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 of sign in 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.

IV. CONCLUSION AND OUTLOOK
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 41 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 behavior has 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 transitionwhich leads to the formation of an extended and dilute network under quiescent bulk conditions -transfers to situations far from equilibrium.Thus, while a systemspanning 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][12][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 threebody 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 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 which 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 Θ 0 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 75 .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 78 , 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.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.
FIG. 8.The density profile ρ(z) (a) as well as the perpendicular (b) and parallel (c) component of the internal force profile fint(z) are shown for the sheared three-body gel with modified preferred three-body angles of θ0 = 175 • , 170 • , 160 • , 150 • (indicated by ticks on the color scale).We set a temperature of kBT = 0.1 and a shear amplitude of K = 5 /σ.As network formation also occurs for the above values of θ0 and as it is the driving mechanism for the strong superadibatic response, one can observe similar behavior as for the choice of θ0 = 180 • in the main text.Below a value of θ0 = 150 • , an accurate sampling of the steady state was hindered by the formation of droplets in the flow channels.

Appendix A: Internal stress tensor
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 ∇ • σ(r, t) = ζJ(r, t) (A1) of the total stress tensor σ(r, t).To identify its internal contribution, we multiply eq. ( 12) by the density profile for the internal stress tensor σ int (r, t).
In the considered stationary state, the time dependence can be dropped.To obtain σ int (r) from the sampled force density profile F int (r) requires an integration of its spatial components according to eq. (A3).

FIG. 3 .FIG. 4 .
FIG. 3.The cluster size distribution C(n) is shown for different values of the temperature (indicated by ticks on the color scale) at constant shear amplitude K = 5 /σ.Particles tend to form chains when driven by the shear force, cf.fig.1, and the mean size of the chains grows when temperature is decreased.Additionally, for kBT = 0.1 , the occurence of large clusters that span across a flow channel and include up to half of the particles in the system is observed.

FIG. 5 .
FIG. 5.Similar to fig.2, the one-body profiles of velocity vx(z) (a), density ρ(z) (b), perpendicular internal force fint,z(z) (c) and parallel internal force fint,x(z) (d) are depicted.The temperature is now fixed to kBT = 0.1 and the amplitude of the external force is varied with values of Kσ/ = 0.1, 0.5, 1, 1.5, 2, 3, 4, 5, which are indicated by ticks on the color scale (K = 0 corresponds to the bulk state and is not shown).Arrows indicate the local direction of the internal force components for low (blue) and high (yellow) shear amplitude.

FIG. 6 .
FIG.6.The superadiabatic viscous force fsup,x(z) is shown as obtained from the model (19) and(20) for Pexc[ρ, v] with η = 5.For this, the expression (21) is evaluated with the density profiles from the simulations in Sec.III B. The velocity profiles are approximated analytically by v(r) ≈ fext(r)/ζ.When comparing the shown PFT results (b) with the adaptive BD simulation data (a) for fint,x(z), good agreement is found.Solely in the center of the flow channels, the simple model for Pexc[ρ, v] leads to deficiencies due to the complex behavior of ρ(z) in these regions.

FIG. 9 .
FIG.9.The steady state behavior of a sheared low-density Lennard-Jones fluid is shown for a temperature of kBT = 1.5 and for various shear amplitudes Kσ/ = 5, 10, 20, 50, 100 (indicated by ticks on the color scale).The superadiabatic response of this representative simple fluid is much weaker than in the sheared three-body gel.The migration of particles always occurs towards the center of the flow channels, where the velocity gradient vanishes, as can be deduced from the density profile ρ(z) (a) and the perpendicular internal force profile fint,z(z) (b).The viscous superadiabatic force fint,x(z) (c) counteracts the flow direction and unlike in the three-body gel, no drag-along is observed.