The shear-curvature viscosity

For highly non-uniformly ﬂowing ﬂuids, there are contributions to the stress related to spatial variations of the shear rate, which are commonly referred to as non-local stresses. The standard expression for the shear stress, which states that the shear stress is proportional to the shear rate, is based on a formal expansion of the stress tensor with respect to spatial gradients in the ﬂow velocity up to leading order. Such a leading order expansion is not able to describe ﬂuids with very rapid spatial variations of the shear rate, like in micro-ﬂuidics devices and in shear-banding suspensions. Spatial derivatives of the shear rate then signiﬁcantly contribute to the stress. Such non-local stresses have so far been introduced on a phenomenological level. In particular, a formal gradient expansion of the stress tensor beyond the above mentioned leading order contribution leads to a phenomenological formulation of non-local stresses in terms of the so-called “shear-curvature viscosity”. We derive an expression for the shear-curvature viscosity for dilute suspensions of spherical colloids and propose an effective-medium approach to extend this result to concentrated suspensions. The validity of the effective-mediumpredictionisconﬁrmedbyBrowniandynamicssimulationsonhighlynon-uniformlyﬂowingﬂuids.


I. INTRODUCTION
Complex fluids can exhibit flow profiles in which there are unusually large spatial gradients in the local shear rate. For these highly non-uniformly flowing systems, there are contributions to the stress that arise from spatial variations of the shear rate, which are referred to as non-local stresses. Such large spatial gradients in the shear rate occur, for example, in micro-channel fluidics devices. Due to these large gradients, flow profiles in micro-channels cannot be described on the basis of the standard expression for the stress tensor. It is shown in Ref. 1 that non-local stresses are a necessary ingredient to describe the flow profiles of worm-like micellar systems in micro-channels with sufficiently large dimensions to describe the micellar solution as a continuum. Highly nonuniform flow profiles also occur when a fluid exhibits an instability that leads to gradient-banded flows. In the stationary state, two regions (the "bands") with spatially constant shear rates coexist. The shear rates within the bands are different for the two bands. The bands are connected by a sharp interface, where gradients in the shear rate are very large. There are two types of gradient-banding instabilities. A uniform flow profile is unstable when the stress decreases with increasing shear rate. In such cases, a transition occurs to a a) Electronic mail: w.j.briels@utwente.nl b) Electronic mail: j.k.g.dhont@fz-juelich.de stable gradient-banded flow profile. Gradient-banding can also occur due to shear-gradient induced mass transport. Migration of particles from regions of high shear rate to regions of low shear rate leads to concentration variations that in turn give rise to gradient-banded flow profiles. Two recent theoretical developments on shear-gradient induced diffusion of colloids and polymers can be found in Refs. 2 and 3, respectively. For both types of shear-banding scenarios, it is essential to include non-local stresses in the constitutive modeling of gradient-banding transitions, in order to account for the rapid spatial variations of the shear rate within the interface that connects the two bands. There are many different types of complex fluids which exhibit gradient-banding transitions. The most extensively studied systems that exhibit gradient banding are worm-like micellar systems (the early studies of these systems include Refs. [4][5][6][7][8][9]. Gradient banding also occurs in entangled polymers, [10][11][12][13][14] micellar cubic phases, 15,16 supra-molecular polymer solutions, 17 transient networks, 18 thermotropic liquid-crystalline polymers, 19 hexagonal phases of surfactant solutions, 20 hard-sphere suspensions, 2,21 soft glasses of charged colloidal rods, 22 metallic glasses, 23 and foams. 24 Gradient banding is sometimes the result of a shearinduced second phase, which has been observed in lamellar surfactant systems, [25][26][27] in a semi-flexible thermotropic liquid-crystalline polymer, 28 and in poly-crystalline colloids. [29][30][31][32][33] Review papers that address shear-banding phenomena include Refs. [34][35][36][37][38][39][40][41][42][43][44]. For all the above mentioned systems, non-local stresses are essential in the constitutive modeling of the interface between the shear-bands.
The standard expression for the deviatoric part of the stress tensor for isothermal incompressible fluids reads where p is the pressure,Î is the identity tensor, η is the viscosity, and E = (1/2)[∇v + (∇v) T ] (where "T " stands for transpose) is the symmetric part of the velocity gradient tensor ∇v of the fluid velocity v. This standard expression is obtained by formally expanding the stress tensor to first-order in gradients of the fluid velocity. There are two ways in which non-local stresses can be formulated. (i) A diffusive term can be added to the equation of motion for the stress tensor (see, for example, Refs. [45][46][47][48], where the corresponding diffusion coefficient is referred to as the stress-diffusion coefficient. Such a constitutive approach has been used to determine the stress diffusion coefficient for a micellar system from the kinetics of band formation. 49 This formulation of non-local stresses has also been applied to analyze the stability of the interface between gradient-bands, 50 where an undulation instability of the interface can give rise to vorticity banding. 51 (ii) Another possibility to include non-local stresses is to simply extend the formal expansion of the stress tensor with respect to spatial gradients to include the next higher-order spatial derivative of the flow velocity, as compared to the leading order expansion in Eq. (1). 48,52,53 For incompressible systems, this leads to where κ is referred to as the shear-curvature viscosity. 52 The minus sign in the non-local stress contribution renders κ positive. This is the constitutive equation that has been successfully used in Ref. 1 to describe the flow profiles of worm-like micellar systems in micro-channels. Note that the shear-gradient contribution to the stress tensor can be formulated within a general linear response approach for inhomogeneous systems as where the "viscous-stress kernel" S is zero for distances |r − r | larger than the range of correlations between fluid elements. Gradient expanding the velocity gradient tensor as identifies the viscosity and the shear-curvature viscosity in terms of moments of the viscous-stress kernel, The present theory describes non-local stresses in bulk, away from boundaries. Non-local stress/velocity-gradient relations, as formally given in Eq. (3), can also arise in confinement due to local modifications of the bulk viscosity. A detailed study of nanoscopically confined water by two sheared hydrated surfaces leading to such non-local stresses can be found in Ref. 54.
So far there are no (semi-)microscopic considerations to derive the constitutive relation in Eq. (2), which allow us to predict the magnitude of non-local stresses in bulk and, in particular, to predict or estimate the numerical value of the shear-curvature viscosity. It is the purpose of this paper to verify that the non-local stress contributions are of the form, as given in the constitutive equation (2), and to show that the shear-curvature viscosity κ is proportional to the shear viscosity η. Brownian dynamics (BD) simulations are performed on a relatively simply system of spherical particles. To derive an approximate theoretical expression for the shear-curvature viscosity, we invoke an effective-medium argument in which each particle in the system is considered as being immersed in a continuum that represents the rest of the system. In preparation to apply this effective-medium argument, we solve the corresponding hydrodynamic problem in Sec. II of a single sphere in a continuum fluid, where the velocity of the fluid prior to insertion of the sphere is highly non-uniform. This analysis confirms the proposed structure for the constitutive equation in Eq. (2) and leads to an explicit expression for the shearcurvature viscosity. The effective-medium approximation is further specified in Sec. III. The structure of the constitutive equation as well as the proportionality of the shear-curvature viscosity to the shear viscosity is verified by BD simulations of a relatively simple system of Brownian particles in Sec. IV.

II. NON-LOCAL STRESSES FOR DILUTE SUSPENSIONS: THE EINSTEIN ANALOG FOR THE SHEAR-CURVATURE VISCOSITY
Consider a sphere that is inserted in a fluid (hereafter referred to as the "solvent") that undergoes, prior to insertion of the sphere, significant spatial variations over a length scale comparable to the size of the sphere. On insertion of the sphere, stresses will be generated that not only depend on the local shear rate of the fluid at the position of the sphere before insertion but also on the local spatial variations of the local shear rate. The latter contributions to the stress are the non-local stresses.
In order to evaluate the stress tensor for a system that consists of a strongly inhomogeneously flowing solvent in which an assembly of spheres is embedded, we depart from an expression for the divergence of the stress tensor that is valid for inhomogeneous systems, as derived in Ref. 55 by two of the present authors, where the shear stress 2η 0 E (with η 0 as the shear viscosity of the solvent) and the pressure P ss are due to solvent-solvent molecule interactions (hence the subscript "ss"). The last term is the contribution of the colloidal particles to the stress, where N is the number of colloids under consideration, the integral ranges over the surface ∂V p of the colloids, f h p (r ) is the local force per unit area that the solvent exerts onto the surface elements of colloid p, and δ(·) is the Dirac delta function. The brackets · · · denote averaging with respect to the position coordinates (and the orientations in the case of non-spherical particles) of the N colloids. This expression generalizes Batchelor's expression 56,57 for the stress of homogeneously sheared systems to non-uniformly flowing systems. Batchelor's classic expression for the stress is obtained from Eq. (6) in the case of a constant shear rate and concentration. The interpretation of Eq. (6) for the stress tensor is rather straightforward: the integrals represent (minus) the force with which the surfaces of the colloids act onto the solvent, and the Dirac delta function counts only those colloids whose surfaces are located at position r where the divergence of the stress is evaluated. This is what is expected for the force per unit volume (which is by definition the divergence of the stress tensor) produced by the colloids. In this section, we consider an assembly of spheres embedded in a strongly inhomogeneously flowing solvent, where the spheres do not mutually interact. For such very dilute suspensions, each sphere can be considered as being embedded in an otherwise unbounded solvent since the presence of the remaining spheres does not affect the flow for the sphere under consideration.
Note that when the solvent velocity varies over distances comparable to the size of the colloids, their Brownian motion results in quite large temporal fluctuations of the local stress. The brackets · · · in Eq. (6) represent the thermal average over these fluctuations and are thus only applicable for macroscopic time-dependent flows which vary sufficiently slow, as compared to the time required for the colloids to diffuse over distances comparable to their size.
In order to obtain an expression for the shear-curvature viscosity, the stress tensor needs to be expanded up to thirdorder in the spatial gradient of the suspension flow velocity. Such a gradient expansion of the stress tensor with respect to spatial gradients is obtained by expanding the Dirac delta function in Eq. (6) according to where the " " denotes the contraction with respect to the polyadic products (r − r p ) n = (r − r p )(r − r p ) · · · (r − r p ) and ∇ n = ∇∇· · · ∇. Substitution into Eq. (6) leads to the following spatial-gradient expansion of the divergence of the stress tensor: where E = (1/2)[∇v + (∇v) T ] is the velocity gradient tensor of the suspension and where (with n ≥ 0) where the force moments T (n) p are defined as The calculation of the surface force density on the sphere and the resulting force moments for the inhomogeneous flow under consideration is a technical mathematical problem (the mathematical details are given in Appendix A). The resulting expressions for the force moments, up to contributions of fourth-order in spatial gradients of the flow velocity, are formulated in terms of the flow velocity u of the solvent prior to insertion of the sphere and the corresponding velocity gradient tensor E (0) = (1/2) ∇u + (∇u) T . As shown in Appendix A, it is found that (indices refer to the components of vectors and tensors) where F Br = −k B T ∇ ln ρ is the Brownian force (with k B as Boltzmann's constant, T as the temperature, and ρ as the colloid number density). The Brownian force and velocities are understood to be evaluated at the position r p of the inserted colloidal particle. Substitution of Eq. (11) into Eq. (9) leads to ∇ · Σ (0) (r) = −∇ · ρ k B TÎ , where the volume fraction of colloids is introduced, The leading gradient contributions to the divergence of the stress tensor are proportional to first-order gradients in the concentration and to fourth-order in the flow velocity, in accordance with Eq. (2). We will therefore keep only such leadingorder gradient contributions to the non-local stress so that mixed terms like ∇ϕ · ∇E (0) are neglected. Adding all terms in Eq. (12) for the various contributions to the body force originating from the colloids thus leads to As a last step, the velocity-gradient tensor E (0) corresponding to the solvent velocity prior to insertion of the sphere must be expressed in terms of the suspension's velocitygradient tensor E. The macroscopic flow velocity v of the suspension is the volume average of the solvent velocity u s after immersion of the sphere and the local velocity V p of the colloids, where, as before, ϕ is the volume fraction and where the brackets · · · denote averaging with respect to the positions of the non-interacting colloids. Since for dilute dispersions of noninteracting particles the stress is linear in the concentration, the explicit volume fraction dependence on the right-hand-side in Eq. (15) can be neglected: retaining the explicit volume fraction dependence leads to terms ∼ϕ 2 in the expression (14) for the stress tensor, which are terms that are beyond the validity of that expression. Hence, as far as the evaluation of the relation between E (0) and E is concerned, Eq. (15) reduces to v(r) = u s (r). Furthermore, there is a difference between u s (r) and u(r) only when there is a colloid in the vicinity of r, implying that the difference between u s (r) and u(r) is of the order of the volume fraction (remember that u s and u are the velocities after and prior to insertion of a sphere, respectively). Hence, for the evaluation of the difference between E (0) and E, we have v(r) = u(r), which implies that we can replace E (0) in Eq. (14) by the velocity-gradient tensor E of the suspension.
It is thus finally found that the stress tensor takes the form (2) that was phenomenologically suggested in Ref. 52, where the pressure P is equal to the viscosity η is given by and the shear-curvature viscosity is equal to Equation (17) reproduces the ideal gas law for the osmotic pressure of dilute suspensions, Eq. (18) reproduces Einstein's expression for the shear viscosity, and Eq. (19) is the corresponding analogous low-concentration expression for the shear-curvature viscosity. Hence, on insertion of a sphere in a solvent with a constant shear rate [see Fig. 1(a)], the additional stress caused by the presence of the sphere is equal to 2η 0 (5/2)ϕE [see Eqs. (16) and (18)]. On insertion of a sphere in a solvent with a curved flow profile [see Fig. 1(b)], the curvature of the flow leads to yet another stress contribution equal to −2η 0 (a 2 /10)ϕ∇ 2 E [see Eqs. (16) and (19)], which is the non-local contribution to the stress.

III. AN EFFECTIVE-MEDIUM APPROXIMATION FOR THE SHEAR-CURVATURE VISCOSITY
Within an effective-medium approach, the inhomogeneously flowing solvent in the considerations in Sec. II is replaced by the suspension of colloidal spheres which are identical to the sphere that is inserted. The suspension is thus considered as an effective medium that behaves like a solvent for the inserted sphere. The viscosity η 0 of the solvent in Eq. (19) is thus simply replaced by the shear viscosity η of the suspension at the given volume fraction ϕ and shear rateγ. According to Eq. (19), κ/N ∼ η 0 (with N being the number of colloidal spheres), that is, the contribution to the shear-curvature viscosity per colloidal sphere is linearly related to the shear viscosity. Since in concentrated suspensions each sphere on average contributes equally to the stress [corresponding to the particle-surface integrals in Eq. (6)], the effective-medium approach thus states that the additive contribution of each sphere to the shear-curvature viscosity is proportional to the suspension viscosity. This leads to the following effective-medium approximation for the shear-curvature viscosity: where the volume fraction [see Eq. (13)], the number concentration ρ, and the shear rate are understood to be evaluated at position r. This expression for the shear-curvature viscosity quantifies the non-local stress due to spatial variations of the shear rate and the concentration and will be compared to simulations in Sec. IV. Note that the above effectivemedium approximation predicts a shear-rate dependence of the shear-curvature viscosity that is similar to that of the shear viscosity.
It is important to realize that the radius a in the effectivemedium approach differs from the hard-core radius of the sphere. A sphere that is inserted in a suspension does not behave like a sphere with stick boundary conditions precisely at the geometrical boundary of that sphere. The volume fraction in the first equation in (20) is therefore not the volume fraction that corresponds to the hard-core radius of the sphere. In a comparison with simulations in Sec. IV, it is therefore necessary to make the distinction between a in Eq. (20) and the hard-core radius of the spheres. We shall hereafter refer to a as the hydrodynamic-insertion radius, which is the effective radius of a sphere that leads to an increase in the stress on insertion into the suspension when the suspension for that sphere behaves like a solvent with stick boundary conditions. The hydrodynamic-insertion radius will be larger, but of the same order of magnitude, as the hard-core radius, which will be discussed in more detail later.
The effective-medium result (20) for the shear-curvature viscosity will be validated against Brownian dynamics simulations in Sec. IV.

IV. BROWNIAN DYNAMICS SIMULATIONS
The effective-medium approach discussed above is an approximation, the accuracy of which will be tested against Brownian dynamics (BD-)simulations. In principle, one may perform simulations where the interactions between the colloids and their coupling to the stress lead to gradientbanding. [58][59][60] Since non-local stresses are only significant within the sharp interface that separates the two shear-bands and only a small fraction of the colloids reside within the interface, such simulations would suffer from bad statistics. We therefore chose to simulate a flow between parallel plates, which may be considered as a two-dimensional analog of the type of flow encountered in micro-fluidics devices. Essentially all the colloidal spheres now experience a highly non-linear flow. In order to probe the non-local stress with sufficient statistics for a quantitative comparison to the effective-medium prediction, a large number of particles must be simulated. Including full hydrodynamic interactions between the colloidal particles would therefore require unrealistically long computation times. The simulated fluid of Brownian spheres plays the role of the effective medium into which the sphere is immersed. On a continuum level, this Brownian fluid acts as a hydrodynamic medium. As a proof-of-principle, we therefore chose to simulate Brownian spheres that are subject to a spatially sinusoidally varying driving body force in the direction parallel to the confining plates in the absence of a solvent. The imposed body force is nothing but a means to induce a non-linear flow profile, which does not affect the measured shear-stress directly, but only indirectly through inter-particle interactions. The fundamental propagator for the BD-simulations thus reads where β = 1/k B T (with k B as Boltzmann's constant and T as the temperature), D 0 is the diffusion coefficient that sets the time scale of Brownian motion, and Θ i is a Gaussian variable with average zero and variance unity. Furthermore, F imp i is the above mentioned imposed force on particle number i in the x-direction with its gradient in the y-direction. The imposed force gives rise to an imposed velocity v with x i and y i as the xand y-component of the position coordinate r i of particle i andê x as the unit vector in the x-direction. The first term corresponds to a simple shear velocity with shear rateγ, and the second term is a superimposed sinusoidally varying velocity with wave number k and amplitudeĠ 0 /k. F p i is the force on particle i due to direct interactions with other particles, with V (r) being the inter-colloidal pair-interaction potential.
For the pair-interaction potential, we choose a hard-core potential corresponding to a sphere with radius R, augmented with a standard DLVO potential (details are given in Appendix B). The simple linear shear flow is maintained by Lees-Edwards boundary conditions, while the wave number k = 2πn/L with n = 1, 2, . . . and L is the height of the simulation box.
To calculate ensemble averaged quantities at a given position r, the simulation domain is divided into 200 slabs that are stacked along the gradient y-direction, spanning the flowvorticity plane. The ensemble averaged macroscopic velocity profiles are thus obtained from v(y) = 1 N(y) where the summation ranges over particles which reside within the slab with its position at y j and with a velocity v p (y j ) and where N(y) is the number of spheres within the same slab. The local macroscopic stress σ(y) is determined in the same way, where stresses arising from interactions with particles outside the slab are accounted for the averaging procedure, where inter-particle forces of particles whose line segment intersects with the slab under consideration are included in the stress computation, and account for contributions to the stress originating from spatial inhomogeneities. [61][62][63] The BD-simulations are also used to obtain the shear-rate dependence of the viscosity, for whichĠ 0 = 0. All simulations start from equilibrated states, which are achieved from simulations without flow over an extended period of time (more than 20 times R 2 /D 0 ). The stationary state under flow conditions is subsequently reached after an additional simulation under flow conditions over a time period of at least 100 × R 2 /D 0 for homogeneous shear flow and 1000 × R 2 /D 0 for inhomogeneous shear flow. Flow-and stress profiles are obtained once the stationary state is attained in the way described above. The time interval dt * = (R 2 /D 0 )dt is chosen as 5 × 10 −5 , which is verified to be sufficiently small. The number of BD-simulation time steps for an inhomogeneous flow is at least equal to 2 × 10 7 . The system consists of 16 000 spheres, dispersed within a rectangular box. The length of the box along the shear gradient direction is twice larger than that in the other two directions, in order to achieve a sufficiently finely divided set of values of the wave numbers k = 2nπ/L. The volume fraction corresponding to the hard-core radius R of the spheres is 0.30 so that the size of the cube along the gradient direction is about 100 times the particle radius. The simulation box size (100 R) is much larger than the distance over which particles are correlated (about 4-5 R). Furthermore, box sizes similar to ours have been shown to be sufficiently large to simulate shear-banded systems. [58][59][60] True bulk properties are therefore probed, assuring that the shearcurvature viscosity as obtained from our simulations is a true material function. Additional evidence will be presented in Sec. VI.

V. RESULTS AND DISCUSSION
In this section, we present the results of BD-simulations and compare them with the effective-medium approximation presented in Sec. III. All results will be presented in the dimensionless form. Lengths are non-dimensionalized with respect to the sphere radius R, times are non-dimensionalized with respect to R 2 /D 0 , velocities are non-dimensionalized with respect to D 0 /R, and forces are non-dimensionalized with respect to k B T /R: The expression for the shear-stress σ, that is, the xy-component of the stress tensor in Eq. (2), is nondimensionalized by the same characteristic variables mentioned above so that σ * = η − κ ∇ 2 γ × 6 π R 3 /k B T , and hence, where the dimensionless gradient operator ∇ * is equal to R∇ and the "relative viscosity" is equal to Note that for a Brownian particle in a solvent, η * = η/η 0 , with η 0 being the shear viscosity of the solvent. Furthermore, the local Peclet number is given bẏ while the dimensionless shear-curvature viscosity (per unit volume fraction) is defined as where is the hard-core volume fraction. Note that this volume fraction differs from the volume fraction ϕ in Secs. II-IV, which is based on the hydrodynamic-insertion radius [see Eq. (13) and the discussion at the end of Sec. III]. As a first step, the viscosity of a homogeneously sheared system (for whichĠ 0 = 0) is simulated as a function of the shear rate, for a given volume fraction of ϕ * = 0.3. The relative viscosity as a function of the Peclet numberγ * is given in Fig. 2. The black solid line is a fit of the blue simulation data points to the empirical Carreau-Yasuda viscosity equation, which has been shown before to accurately describe viscosity data, 64,65 where η q is the zero-shear viscosity, η ∞ is the high-shear viscosity, λ is a relaxation time, and m is a power-law index. A similar shear-thinning behavior and numerical values for the relative viscosity have been obtained, for example, in Refs. 66 and 67 for hard-sphere suspensions without hydrodynamic interactions with Brownian dynamics simulations and in Ref. 68 for Lennard-Jones fluids with Molecular dynamics simulations. The Carreau-Yasuda fit will be used later in order to analyze the numerical results for inhomogeneously sheared systems.
Since the imposed velocity in Eq. (22) for the inhomogeneously sheared systems is sinusoidal, the corresponding resulting particle velocities will be similarly sinusoidal, with an amplitude of approximatelyĠ 0 /k. The local shear rate is therefore equal toγ +Ġ 0 cos{k y}. In the simulation,Ġ 0 is a constant, in order to fix the values of local shear rates for each different wave number. The amplitude of the sinusoidal flow-velocity perturbation is chosen to be small as compared to the linear flow velocity corresponding to the shear rateγ. The reason for this is to clearly separate the local and non-local stresses. In case the amplitude of the sinusoidal flow velocity would be relatively large, there would be a significant sinusoidal variation of the stress due to the spatial variation of the shear viscosity itself on top of the non-local stress. The lower limit of the amplitudeĠ 0 to obtain accurate results is about D 0 /(10 R 2 ). The local shear-rate amplitude is not exactly equal toĠ 0 since direct inter-colloidal interactions change the flow profile. The local velocity is therefore written as v(y) =γ y +Ġ k sin{k y}, where the amplitudeĠ/k is determined from a fit of the simulated flow profile. The local shear rate is obtained by the differentiation of Eq. (31), γ(y) =γ +Ġ cos{ky}, (32) or in dimensionless variableṡ γ * (y * ) =γ * +Ġ * cos{k * y * }, where dimensionless quantities are introduced as before. An example of a velocity and shear-rate profile is shown in Fig. 3 for a Peclet numberγ * = 0.5 and a wave number of k = 4π/L. The velocity profile is plotted in Fig. 3(a), from which the amplitude V =Ġ * /k is determined by a fit to Eq. (31), which is indicated by the red solid curve. The corresponding shearrate profile as obtained from numerical differentiation of the simulation data is plotted in Fig. 3(b). Substitution of Eq. (33) for the inhomogeneous shear rate into Eq. (25) for the stress leads to where the relative shear viscosity is evaluated at the local shear rate given in Eq. (33). Examples of inhomogeneous shearstress profiles are given in Fig. 4(a) for an average shear rate ofγ * = 0.5 and for n = 1, 2, 3, 4 (the black, green, blue, and red data points and fits, respectively). The solid lines in Fig. 4(a) are the result of a global leastsquare fitting of all the stress profiles simultaneously. One fit  parameter is the dimensionless parameter, According to the effective-medium approximation (20), this dimensionless parameter should be equal to As before, a is the hydrodynamic insertion radius (see the discussion at the end of Sec. III), while R is the hard-core radius of the spheres. The effective-medium approximation thus predicts that Λ 0 is independent of the applied mean shear rateγ. Two more fitting parameters are introduced, which are necessary to correct for the errors involved in the independently obtained simulated viscosity data, for which the Carreau-Yasuda fit in Eq. (30) is used. We thus minimize the square ε(Λ 0 , C 1 , where the first sum ranges over all the different wave numbers [corresponding to all the stress profiles shown in Fig. 4(a)] and the second sum ranges over all positions where the stress is evaluated. Here, σ * ,sim (y * m |γ * , k * ) is the simulated stress, and with corresponds to the stress in Eq. (34), except for the two additional fitting parameters C 1,2 which account for the inaccuracy of the viscosity BD-simulation results. The fitting values of these two parameters should be close to unity. A fit for each separate wave vector (using the values for C 1,2 , as found in the global fits) gives numerical values for A(k * ) in Eq. (38), where the linear dependence on k * 2 in Eq. (39) is the prediction from the theory. Figure 4(b) shows that such a linear dependence is indeed reproduced by the simulations. Figure 4(c) shows the global fitting result for one of the stress profiles that is also given in Fig. 4(a) (with n = 4). The   FIG. 5. (a) Results of the least-square fits for Λ 0 , C 1 , and C 2 (the red stars, green circles, and blue triangles, respectively) for three values of the Peclet numberγ * = 0.2, 0.5, and 1.0. The horizontal red line is the weighted average of the three data points for Λ 0 . (b) The relative viscosity (in blue) as compared to the relative viscosity η * κ as obtained from the numerical values for Λ 0 = 3.8 and a/R = 2.07 ± 0.06 (in red) according to the effective-medium approximation in Eq. (20). inset shows the difference between the fitting results including the non-local stress and that without the contribution from the non-local stress. The comparison between these fits shows that non-local stresses are essential to describe the simulated stress profiles.
Results for the fitting parameters for three different values of the Peclet numberγ * are given in Fig. 5(a). As can be seen from this figure, the values of C 1,2 are quite close to unity, as they should, while the value of Λ 0 = 3.8 ± 0.5 is constant and independent ofγ * to within simulation errors, as predicted by the effective-medium approximation [see Eq. (36)]. According to Eq. (36), the numerical value of Λ 0 found in Fig. 5(a) corresponds to (a/R) 5 = 38 ± 5, and hence a/R = 2.07 ± 0.06. That is, the hydrodynamic-insertion radius of a sphere in the simulated system in the absence of hydrodynamic interactions is about two times larger than its hard-core radius. This is of the same order as the distance from a sphere over which the hydrodynamic behavior of a suspension sets in as found in simulations of Brownian hard-spheres without hydrodynamic interactions in 2D. 69 The relative viscosity η * κ that follows from the effective-medium prediction in Eq. (20) with the numerical value of Λ 0 = 3.8 and a/R = 2.07 is plotted in Fig. 5(b). This plot confirms the proportionality κ ∼ η, as predicted by the effective-medium expression for κ in Eq. (20), within a shear-rate range where the viscosity significantly shear-thins.

VI. SUMMARY AND CONCLUSIONS
For dilute suspensions, we derived a constitutive equation that includes non-local stresses, which reproduces the earlier phenomenologically proposed expression for the stress tensor 48,52,53 [see Eqs. (2), (16)- (19)]. This derivation is based on the evaluation of the additional stress that is generated by insertion of a sphere into an inhomogeneously flowing solvent, in combination with a general expression for the stress tensor that includes non-local stresses 55 [see Eq. (6)]. The analytical derivation leads to an explicit expression for the so-called shear-curvature viscosity that measures the nonlocal stress 52 [see Eq. (19)]. An effective-medium approach is proposed [see Eq. (20)], where the suspension itself is considered as "the solvent" into which a sphere is immersed. The effective-medium approach necessitates the introduction of the "hydrodynamic-insertion radius", which is the effective radius of a sphere that leads to an increase in the stress on insertion into the suspension when the suspension behaves like a solvent with stick boundary conditions. This radius is larger, but of the same order, as the hard-core radius of the sphere.
Brownian dynamics (BD-)simulations are performed on a highly non-uniformly flowing system of Brownian spheres with the neglect of hydrodynamic interactions. The simulations confirm the effective-medium prediction that the shearcurvature viscosity is proportional to the shear viscosity: the shear-rate dependence of the shear-curvature viscosity scales quite accurately with that of the shear viscosity. The hydrodynamic-insertion radius is found to be twice as large as the hard-core radius of the spheres. This is a quite reasonable result since the distance from a given sphere where a concentrated suspension (again with the neglect of hydrodynamic interactions) behaves as a macroscopic fluid is also of that order. 69 Note, however, that the hydrodynamic-insertion radius a is concentration dependent: it varies from about a ≈ 2R, as found in the present simulations at a high concentration, to a = R in dilute dispersions (with R as the hard-core radius). The effective-medium approach is not able to accurately predict the concentration dependence of the shear-curvature viscosity. Further analytical theory and/or simulations will be necessary to elucidate the concentration dependence of the hydrodynamic-insertion radius.
In simulations of the shear viscosity, through imposed wavelength dependent body forces, the box is considered to be sufficiently large when the inferred viscosity does not depend on the wavelength. Similarly, we find that the wave-vector dependence of the simulated stress agrees with the theoretical continuum prediction [see Eq. (34)], in which λ 0 does not depend on the wavelength [see Fig. 4(b)]. This confirms that the obtained shear-curvature viscosity is a true material property.
This work is limited to suspensions of spherical particles. It would be interesting to extend the approach described in this work to systems consisting of more complex macromolecules, like colloidal rods, polymers, and worm-like micelles. The development of an analytical theory, like for spheres in this work, may not be feasible. One approach could be to verify the relation between the shear-curvature viscosity and the shear viscosity as given in Eq. (20) by means of similar simulations as for spheres in this work and identify the length scale connected to the hydrodynamic-insertion radius for these more complex systems. The effective-medium relation between the shear-curvature viscosity and the shear viscosity in Eq. (20) may be more generally valid, where the challenge lies in the interpretation of the length scale corresponding to the hydrodynamic-insertion radius.
Similar to the experiments on worm-like micellar systems in Ref. 1, the experimental validation of the effectivemedium results in Eq. (20) for the shear-curvature viscosity could be established by probing velocity profiles in narrow micro-fluidics channels. Such experiments on suspensions of (spherical) colloids have not yet been performed. The effective medium approach for the non-local stress is also relevant for time-dependent flows, for example, oscillatory flows in narrow channels. For large-amplitude oscillatory flows, the non-linear response functions consist of additive local and non-local viscous contributions which are similarly proportional to each other as for the linear response functions.

APPENDIX A: EVALUATION OF FORCE MOMENTS
In this Appendix, the mathematical details are given for the calculation of the force moments in Eq. (10) for a single sphere immersed in an inhomogeneously flowing solvent.
Let u(r) denote the flow field that exists prior to insertion of a sphere. This incident flow field is gradient expanded up to third-order in spatial gradients around the instantaneous position r p of the sphere (which is taken to be at the origin), u(r) = u (0) (r) + u (1) (r) + u (2) (r) + u (3) (r) + u (4) (r) + · · · , (A1) where u (0) (r) = u(r p ), Hereafter, we will refer to u (n) as the nth-order flow field. Note that the prescribed velocity u of the solvent requires an external body force. The Stokes equations are still fulfilled, but with this body added to the divergence of the stress. Due to the linearity of the governing equations in the creeping flow regime, the induced forces are simply the sum of the surface forces induced by each of the separate flow fields in Eq. (A1). These contributions will be denoted as f (n) , The force densities f (n) (r) (withr as the outward unit normal vector on the spherical surface) are even functions ofr for even n and odd functions for odd n.
A little thought shows that the divergence of the stress tensor in Eq. (4) in the main text, up to fourth-order in spatial gradients, depends on the various moments of the nth-order force densities in Eq. (A3) as (2) (r), The flow field v(r) that exists after insertion of the sphere is written as v(r) = u(r) + u r (r), where u r (r) is the "reflected flow field," that is, the field that is due to the surface forces with which the sphere acts onto the fluid after insertion. The method of calculating reflected flow fields that we will employ here has been published in a monograph by one of the present authors (JKGD) 70 and is therefore not easily accessible. We will therefore first discuss this method.
The reflected flow field satisfies the incompressibility condition as well as the creeping flow equation, Without loss of generality, it is assumed that the sphere resides instantaneously with its midpoint at the origin. We seek a solution of these equations under the boundary conditions u r (r) = u 0 (r), for r ∈ ∂V 0 , = 0, for r → ∞, with u 0 as a known flow field (which will be related to the incident flow field u later) and where ∂V 0 is a spherical surface of radius a with its center at the origin (as indicated by the index "0"). The solution of this set of equations can be constructed as follows. 70 First expand the known flow field u 0 (r) around r = 0, where (∇ ) l is the lth-order polyadic product of ∇ and stands for the full contraction. The solution of the above set of equations is obtained by replacing the polyadic products r l by "hydrodynamic connectors" U (l+2) (r), These hydrodynamic connectors U (l+2) are tensors of rank l + 2. It is easily seen that this representation of the reflected flow field solves the above boundary-value problem, when The hydrodynamic connectors can be constructed from the m-rank tensors, The first few connectors are equal to The higher-order connectors vary like 1/r 2 or decay faster for large distances. In Ref. 70, higher-order connectors are given as well, where contributions ∼ ∇ 2 ∇ 2 u 0 are omitted. Note, however, that such terms are in the present case non-zero, as an external non-zero body force is generally exerted to sustain the incident flow. Up to the fifth-order connector, however, for which the fourth-order polyadic products of ∇ in Eq. (A9) do not appear, contributions ∼ ∇ 2 ∇ 2 u 0 are not present so that the expressions (A12) can nevertheless be employed for our purposes. The advantage of formulating the solution in terms of hydrodynamic connectors is that a general solution can now be constructed for arbitrary flow fields u 0 . That the expressions in Eq. (A12) solve the posed boundary value problem is easily verified using the identities, These relations will turn out to be quite useful in the evaluation of the pressure and induced force densities later on. The tensors H (m) that we will need, in terms of their components, are equal to H (3) αβγ = 1 r 4 3 r γ δ αβ +r α δ βγ +r β δ αγ − 15r αrβrγ , H (4) αβγp = 1 r 5 3 δ γp δ αβ + δ αp δ βγ + δ βp δ αγ + 105r prαrβrγ − 15 r prγ δ αβ +r prα δ βγ +r prβ δ αγ +r αrβ δ γp +r αrγ δ βp +r βrγ δ αp , H (5) αβγpq = 1 r 6 −15 δ αβ r q δ γp +r γ δ pq +r p δ γq + δ αγ r q δ βp +r β δ pq +r p δ βq + δ αp r q δ βγ +r γ δ βq +r β δ γq + δ βγ r α δ pq +r p δ αq + δ γp r β δ αq +r α δ βq + δ βp r γ δ αq +r α δ γq + 105 r βrγrp δ αq +r αrγrp δ βq +r αrβrp δ γq +r αrβrγ δ pq +r γrprq δ αβ +r αrprq δ βγ +r βrprq δ αγ +r αrβrq δ γp +r αrγrq δ βp +r βrγrq δ αp − 945r αrβrγrprq , wherer α are the components of the unit vector in the direction of r.
Once the reflected flow field is obtained via the above described method, the corresponding pressure is obtained from the creeping-flow equation The additional force density that originates from the incident field u(r) = u (1) in Eq. (A5) is found from (with p (1) 0 as the pressure resulting from the incident field) to be equal to The required force moments are thus equal to γ .

The second-order field
Contrary to the zeroth-and first-order surface force densities, the second-order force density is as yet not known and requires the full exploitation of the method of reflections in the form discussed before. For convenience of notation, the second-order contribution v (2) 0 to the incident flow field incident field is written as where Since the finite translation and rotational velocities of the sphere are fully accounted for in the calculation of the zeroth and first order contributions, in the calculation of the higher order moments, the sphere may be considered to have zero velocities. For stick boundary conditions, we therefore have that u r (r) = −(1/2) rr : F, for r ∈ ∂V 0 , which identifies the field u 0 in Eq. (A7), The gradient expansion in Eq. (A8) contains only a single term with l = 2 so that Eq. (A9) now reads where the contraction is with respect to three indices. Explicitly denoting the contractions gives (where summations over repeated indices is understood) Although we do not need the explicit expression for the reflected field, it is readily found from the above explicit expression for U (4) The pressure p r is again obtained from Eq. (A15). From the first and fourth expressions in Eq. (A13), it is readily found that It immediately follows that The explicit expression for the pressure is thus The stress contribution ∼∇ j u r ,i from the reflected field at the surface of the sphere, for which r = a, follows from Eq. (A43), where, as before,r p is the pth component of the unit vector r = r/r. From Eqs. (A16), (A46), and (A47), it is thus found that the force density f (2) r on the surface of the sphere originating from the reflected field is equal to The additional force density that originates from the incident field u(r) = u (2) in Eq. (A5) is found from (with p (2) 0 as the pressure resulting from the incident field) The additional force density that originates from the incident field u(r) = u (3) in Eq. (A5) is found from (with p (3) 0 as the pressure resulting from the incident field) while the body force B (3) = −η 0 ∇ 2 v (3) 0 +∇p (3) 0 that is necessary to maintain the third-order incident flow field is equal to The additional force density is thus equal to f 0,i = a 2 η 0 1 2r αrβrγ G βαγi The first-order moment of the total surface force density that is needed in Eq. (A4) is thus found to be equal to a 4 η 0 3 ∇ 2 ∇ p u i (r) + ∇ i u p + 5 ∇ 2 ∇ p u i (r) − ∇ i u p r=r p . (A66)

The fourth-order force moment
For the flow field u (4) , we need only the zeroth-order moment of the surface force density. In order to obtain that force moment, it is not necessary to calculate the reflected flow field and pressure as in Subsections 2-4 in this Appendix. Instead, this force moment can be obtained from the integral expression for the total flow field v that exists after immersion of the sphere in the flow field u (4) , v(r) = u (4) where T is the Oseen tensor. Note that the surface force density of corresponding to the reflected field occurs in the integral. For a sphere that is fixed, that is, not translating nor rotating, stick boundary conditions imply that u (4) (r) = ∮ ∂V 0 dS T(r − r ) · f (4) r (r ) , for r ∈ ∂V .
Integrating both sides over the surface of the sphere, using that where, as before, u is understood to be evaluated at the midpoint of the sphere and it is thus found that ∮ ∂V 0 dS f (4) r (r) = 6 π η 0 a a 4 120 ∇ 2 ∇ 2 u.
The contribution from the incident field u (4) can be found as follows. Since where it follows that the corresponding pressure is equal to while the body force is equal to The force density arising from the incident field is therefore equal to f (4) 0,i = 1 6 η 0 a 3 r αrβrγrq H αβγqi +r αrβrγrq H αβγiq −3r irβrγrq H ααβγq (A76) so that ∮ ∂V 0 dS f (4) 0 (r) = 0.
The total force is thus found to be equal to ∮ ∂V 0 dS f (4) (r) = 6 π η 0 a a 4 120 This force is non-zero since there is an external body force acting on the solvent to sustain the non-uniform flow.

The translational and rotational velocities
The velocities of the sphere can be obtained from force balance on the diffusive time scale, where the total hydrodynamic force F h balances the Brownian force F Br = −k B T ∇ ln ρ (with ρ as the number density of spheres) on the sphere, Since F h = ∮ ∂V 0 dS f (0) (r) + f (1) (r) + f (2) (r) + f (3) (r) + f (4) it follows from Eqs. (A23), (A51), and (A77) that V (p) = u(r) + a 2 6 ∇ 2 u(r) + a 4 120 ∇ 2 ∇ 2 u(r) This is Faxén's law for translational motion except for the last term between the square brackets, which is due to the fact that the flow field u is not body-force free (otherwise ∇ 2 ∇ 2 u and all other higher-order derivatives are zero). The above expression is valid for flow fields that are sustained by a non-zero body force up to order (a∇) 6 . Since for spherical particles (at infinite dilution) there is no Brownian torque, the hydrodynamic torque is zero, Since T h = ∮ ∂V 0 dS r × f (0) (r) + f (1) (r) + f (2) (r) + f (3) (r) + f (4) (r) and from Eqs. (A37) and (A66), we have ∮ ∂V 0 dS r × f (1) (r) = 8 π η 0 a 3 1 2 ∇ × u − Ω (p) , while the integrals over the even-superscripted forces are zero, it thus follows that Ω (p) = 1 2 ∇ × u(r) + a 2 12 ∇ × ∇ 2 u(r) This is the rotational Faxén's theorem (in the absence of an external torque), where the last term is due to the non-zero body force to sustain the field u (otherwise η 0 ∇ × ∇ 2 u and all higher-order derivatives are zero). This relation is valid up to order ∇( ∇ 2 ) 2 .