Connection between the anisotropic structure and nonlinear rheology of sheared colloidal suspensions investigated by Brownian dynamics simulations

Using Brownian dynamics simulations, we investigate the connection between the shear-induced microstructural distortion and nonlinear rheology of charged colloidal suspensions subject to steady shear. We demonstrate that their rate-dependent flow behavior is a consequence of localized elastic response, which we define as transient elastic zone (TEZ), generated by particle interaction. The body of colloids under shear behaves like an elastic solid in short distances but like a fluid at long distances. The short-lived, localized elastic region, i.e. transient elastic zone, plays a crucial role in determining the observed rheological behaviors. Our findings shed new light on understanding the nature of nonlinear rheology of soft matters with strong interactions.


Introduction
Colloidal particles suspended in a liquid under flow are used for numerous and varied kinds of applications. For example, collodial suspensions are found in paints, food, and cosmetics. Improved understanding of their flow behavior could potentially be of great benefit to these industries [1,2].
The nature of the relation between the microstructure of colloids and the macroscopic rheological properties of colloidal suspensions, such as viscosity, has been a subject of intense research in soft condensed matter physics [3][4][5][6].
It is of particular interest to understand the microscopic mechanisms of their non-Newtonian behavior, in which the viscosity of the suspension depends strongly upon on flow rate [4][5][6][7][8][9]. For dense suspensions of charged colloids or hard spheres, subtle changes in microstructure due to shear result in significant decreases to steady shear viscosity with shear rate. That is to say, shear thinning occurs without discontinuous jump [4][5][6][7]. This continuous shear thinning behavior is in contrast to some cases where shear may lead to sudden and dramatic changes in the long-range order of a complex fluid. Despite extensive studies, the underlying physics behind the phenomenon of shear thinning are still under debate [4,[7][8][9][10][11][12].
In recent decades, computer simulations employing various model systems such as hard spheres [11][12][13] or charged particles [14] have been performed to understand the exact microscopic mechanism responsible for the observed macroscopic nonlinear flow behavior. However, no conclusive signature has yet been identified.
For example, although it is known that discontinuous shear thinning of colloidal suspensions is associated with formation of particle layers [2,3], whether such ordering is necessary for continuous shear thinning remains unresolved [15]. If the two-dimensional ordering is not the main driving mechanism of this effect, what then is the possible physical origin? To answer this question, the approach from a structural point of view is of great significance, and it becomes necessary to construct a structure-property relationship in colloidal suspension in a shear flow.
An important clue comes from a recent simulation of a sheared metallic glass, which has revealed a close connection between local connectivities in atomic structure network and steady shear viscosity by a detailed analysis of structural anisotropy [16]. Also similar results have been found in recent flow-SANS experiments on suspension of charged colloids [8], thus suggesting a universal nature in shear thinning for a variety of strongly interacting disordered materials. Therefore it is important to elaborate the concept developed in references [8,16].
In colloidal suspension it has been previously considered that the continuous shear thinning is mainly attributed to the distortion of sheared microstructure away from the quiescent state [17]. As a result the contribution of Brownian stress due to thermal fluctuations decreases with increasing shear rate, leading to shear thinning. This picture provides a qualitative explanation for shear thinning of moderately concentrated suspension, but it does not apply to such systems as those with an elastic response or in which a signature dominantly starts to emerge, including deeply supercooled liquids and glasses or highly-charged suspensions. To understand strongly interacting systems, an explanation beyond a simple interplay between shear and thermal forces is needed. Therefore we suggest an alternate mechanism that may describe the microscopic origin of shear thinning in a charged colloidal suspension in terms of the anisotropic structure of colloids.
In this paper we study the correlation between the mobility of colloidal particles and structural anisotropy during the shear-thinning regime. This is accomplished using Brownian Dynamics (BD) simulation of a binary charged colloidal suspension in a simple steady-state shear flow, and we show that shear thinning at high shear stress is associated with decreases in the mesoscopic zone where elasticity develops.
This paper is organized as follows: in section 2 we briefly discuss the model used in this paper, the BD technique and the calculation of rheological and structural properties. The simulation results of the structure under shear are described in section 3, followed by the summary in section 4.

Details of BD simulations
The interaction between binary colloidal particles is described by a pairwise repulsive hard-core Yukawa potential, which has been used for a typical model of charged colloids [18]. The potential between particles i and j, ( ) u r , ij can be expressed as where σ i is the diameter of particle i, z the screening length, and K ij the interaction strength between particles i and j. Here r is the inter-particle distance and b = k T 1 .

B
We define the shear velocity to lie in the x direction, the shear gradient to lie in the y direction, and the vorticity in the z direction. In order to avoid crystallization we chose a Kob-Andersen [19] mixture consisting of 80% A particles and 20% B particles with s s = 1.5,

A B
K AA =K, K BB =0.5 K and K AB =0.75 K. The parameters of potential K, z, and f were determined a priori from SANS measurement of a charged colloidal suspension [1]. These parameters are K=9.69, z=4.86, and f=0.42. We performed Brownian dynamics simulations of this binary mixture of 4000 particles. The algorithm used in this work was proposed by Ermak and McCammon [20], without hydrodynamic interactions and with an additional shear flow in the x direction. The resulting equation for the motion of colloid particles can be described by is the short-time self-diffusion of the particle. Here the positional vector of particle i at time t is = ( ) ( ( ) ( ) ( )) r t x t y t z t , , is extracted from Gaussian distributions with zero mean and variance given by The simulations were performed in a cubic box under periodic boundary cons. Shearing was incorporated into the boundary conditions using the 'sliding brick' boundary conditions of Lees and Edwards [21]. The initial configurations were generated by the following procedure: first, particles were placed randomly in the simulation box at the desired density, and then the overlap between the particles was reduced or eliminated. Once the initial configuration was constructed, several thousand cycles were performed to achieve thermal equilibrium, followed by at least two million cycles where the data was collected. We found that the results obtained were almost the same for each cycle, and that there is no systematic effect of system size in the calculation of these properties.
Throughout this paper we use the following reduced units for length * 0 which is the time required for the particle to diffuse over a distance of σ A The equation (2) in reduced units, for particle A reads Note that equation (3) is governed by a single parameter Pe only, and the Pe measures the strength of shear force relative to thermal force. Similarly the equation (4) for the motion of particle B in reduced units reads During the simulation, we calculated the potential part in macroscopic shear stress as [14] åå This term is the contribution to the stress from the direct interactions between the colloidal particles. The kinetic part in shear stress was ignored because the effect is negligibly small in this study. At a finite shear rate, the steady shear viscosity η is given by The potential contribution to the infinitely-high frequency shear modulus, G ∞ in the formalism of Zwanzig and Mountain [22] is calculated as,

Results and discussions
In this section we present the simulation results for applying shear to a charged colloidal suspension. We begin by presenting the flow curve and anisotropic structure of particles in terms of radial distribution function. We then follow with the results for non-affine displacement parameters, where we observe a connection between the viscosity and the structure at a particle level. Figure 1(a) shows the dependence of shear stress on different Péclet numbers using equation (5). At high shear rates curves merge and can be fit using the Herschel-Buckley model [23].

Flow curve
where τ 0 is the dynamic yield stress below which there is no flow within the accessible range of shear rates, α is the consistency (which is the same value as η evaluated at g =  1), and n is the flow index which defines the extent of non-Newtonian behavior. Our simulations show n=0.4, which indicates shear thinning [1]. Figure 1(b) shows the steady shear viscosity versus Peclet number using equation (6). Compared with the shear stress, we observed that the viscosity at low Peclet number (Pe<0.5) remains constant. This value corresponds to the viscosity at equilibrium. As we increase the Peclet number (Pe >0.5), we notice that viscosity starts to decrease, consistent with what has been reported in several studies [4,[7][8][9][10][11][12]. At a high shear rate, we can use the Herschel-Buckley model to fit the viscosity, similarly to shear stress.

Pair distribution function (PDF)
When a system is in thermal equilibrium without external forces, the symmetry of the system will be isotropic. The isotropic pair distribution function (PDF), g(r), has been the commonly used tool for characterization of the structure of these systems. However if we apply an external force, such as shear stress, to the system, the isotropic symmetry breaks down and structural anisotropy appears, depending upon the type of deformation [16,24]. Such an anisotropic structure can be described by the anisotropic PDF [24], defined by expanding the three dimensional PDF in terms of the spherical harmonics function, where l and m are integers and represent deformation modes, and the coefficient ( ) g r l m defines the anisotropic PDF. It has been recognized [16] that the flow-velocity gradient (x-y) plane is the most relevant to shear viscosity in the simulated geometry. From symmetry arguments, the relevant anisotropic PDFs are those for the second order spherical harmonics of l=2, as these are simply the lowest that satisfy the real space parity conditions and the rotational symmetries of applied shear. Examples include the anisotropic components, To capture an overall picture of the anisotropic structure, we show the behavior of these PDFs. Note that the components -( ) g r 2 2 and ( ) g r 2 2 are related to the shear stress and the first normal stress difference, respectively, as decrease. This observation indicates that structures at high shear rates resemble those of liquids at high temperature, and that shearing raises the effective temperature of the system. There has been much interest in defining the effective temperature under this non-equilibrium condition [16]. Figure 3 shows -( ) g r 2 2 at different Pe. In contrast to ( ) g r , 0 0 we observe a systematic increase in the amplitude of -( ) g r 2 2 even at small values for Pe. The increase in the amplitude directly reflects an increase in shear stress as shown in figure 1(a). These behaviors are consistent with the results obtained from simulation of metallic glasses subject to steady shear [16].
The influence of particle interaction on the macroscopic flow curve may be further elucidated by investigating the dependence of -( ) g r 2 2 on shear rate in the context of rheological behavior. To characterize the anisotropic structure of colloidal suspensions under shear, we used the concepts described in the reference [16]. For a solid subjected to homogenous deformation with a simple shear strain of a dense soft-sphere fluid under shear [25]. This finding suggests that, in the liquid, the strong many-body interaction acts as a free energy barrier for an applied shear to perturb the equilibrium structure. As demonstrated in figure 3, the disagreement between -( ) g r closely resembles equation (13), indicating that, for low shear rates, the system is essentially elastically deformed. This discovery demonstrates that the inter-particle interaction acts like a free energy barrier to resist the applied strain, and the charged colloidal suspension responds elastically against the imposed steady shear. Our results show that in low steady shear field, colloid particles are dynamically correlated over a certain distance, as they are in solids. Therefore we show that the concept of elasticity can be applied not only to elastic solids, but also quite directly to the local response of colloidal suspensions in steady flow at low shear rate as well. It is important to point out that the localized elasticity has also been observed in simulations on metallic liquid [16], but it has not been reported in prior computational and theoretical studies of nonlinear rheology of charge colloids [26,27].
Normal stress difference is an important factor in addition to as shear stress to monitor rheological behavior. Recently, normal stress has been considered to be the key to understanding marked shear thickening in a non-Brownian suspension. Here we show how ( ) g r 2 2 develops during shear thinning. In figure 4, we present anisotropic ( ) g r 2 2 at different values of Pe. For Pe<0.25, we observe little change in the amplitude of ( ) g r . Here parameter γ 0 is calculated from the data of shear stress t xy and shear modulus G ∞ .
Pe>0.25, we notice an increment in the amplitude of ( ) g r . Compared to -( ) g r , 2 2 the peak height of ( ) g r 2 2 becomes even smaller, and the shape near the first peak is found to be asymmetric. This behavior may reflect a higher order structural anisotropy because external force is imposed into the x-y shear plane. is very small, while for Pe>0.5 we notice a pronounced increase in the amplitude of ( ) g r . We found that around the first peak in ( ) g r , has an almost positive value, indicating that associated stress difference is also positive.

Structural analysis of anisotropic PDF
Here let us embark on the structural analysis of -( ) g r 2 2 as already shown in figure 3. For Pe>0.25, the differences in the amplitudes between -( ) g r   heterogeneous response which causes the system to deviate gradually from equation (13). In order to describe these heterogeneities at high shear rate, a generalized equation has been proposed which takes into account the heterogeneous behavior [16]: Equation (14) retains the general expression of equation (13), but the strain is allowed to depend on distance r. In this way, we introduce the shear-induced heterogeneous strain γ(r) into the anisotropic PDF, which make it possible to analyze the structure beyond the continuum approximation.
The profile of γ(r) defines a range over which colloidal particles, on time average, are considered to deform homegenously. We define this zone within which particles are elasticatically deformed as the transient elastic zone (TEZ) [8]. In figure 6, we compare -( ) g r 2 2 and equation (14) using γ(r), which was, extracted simply by the  . The dependence of shear strain γ on r as a function of Pe for a sheared charged colloidal suspension using equation (14). σ A is the particle diameter. We see that the agreement between -( ) g r 2 2 and equation (14) is improved. Figure 7 shows the obtained γ(r) for different values of Pe and represents how the size of TEZ changes with Pe. It clearly illustrates that elastic deformation persists only up to a limited spatial range. Upon increasing the shear rate, the size of TEZ is seen to decrease progressively. A connection between viscosity and the shearinduced structural change at particle level emerges from the results of figure 7. On physical ground, it is reasonable to regard the shear-thinning phenomenon as the consequence of the diminishing size of TEZ. Within the lifetime of TEZ, the colloid particles are highly correlated. With a larger TEZ, the particles can resist the flow more effectively, resulting in higher viscosity. We therefore identify the importance of the influence of the local elastic coherency, manifested by TEZ, upon the rheological behavior of interacting colloidal suspensions.

Investigating the shear-induced structural heterogeneity
We tracked the Brownian motion of an individual colloid under applied shear. The deviation from affine deformation along a mean flow is measured by mapping the motion of nearest-neighbor particles on an affine transformation zone. These zones are identified by quantity D , min 2 as originally proposed in [28]. The D min 2 represents the local deviation from the affine zones, or the non-affine zones, during the interval .This quantity has been reported to be an excellent metric of plasticity that detects local irreversible shear transformations [29,30].
The results for D min 2 reveal an interesting link between shear viscosity and microscopic transport processes at the particle level. Figure 8 shows the 2D map of D min 2 in the x -y plane for Pe=0.16 and 3.3. Here, the time interval δt was set to be the Maxwell relaxation time, h ¥ G . The blue area represents small magnitude of D , min 2 while red area represents larger magnitude. Surprisingly, at the particle level, local distortions in structure under steady shear are highly heterogeneous in space. It turns out that as Pe is increased from 0.16 and 3.3, the number of red zones clearly increases. We will show later that the highly distorted zones can lead to disagreement between -( ) g r 2 2 and the derivative of ( ) g r  observed at the higher shear rate as shown in figure 6. To complement, in figures 9 and 10, we calculated the nonaffine displacement for the planes 1-3 and 2-3, respectively. In both Péclet numbers, we observe small nonaffine motion, maintaining a homogenous flow as we increase the shear rate from 0.16 and 3.3. This implies that heterogeneities observed in the plane 1-2 are related with the -( ) g r 2 2 and shear stress t .
xy Figure 11 gives the comparison between the projection of the -( ) g r 2 2 on the flow (1)-velocity gradient (2) plane of sheared colloidal solutions and the derivative of their respective isotropic angular-averaged ( ) g r . The shear rate is equivalent to P e =3.3. In blue regions, where particles are characterized with small non-affine motion, these two quantities are essentially in quantitative agreement with a properly scaled γ 0 (bottom left panel). However, discernible difference is seen for particles with large non-affine components (bottom right panel) even when properly scaled with γ 0 . The information of figures 8 and 11 suggests that in the blue regions the majority of particles deforms elastically. On the other hand, in the red regions, increasing amount of deformation is dominated by the plastic flow, which involves topological rearrangement of neighboring particles. One can therefore envision that the non-equilibrium structure of sheared liquids is generally  (dashed lines). The difference between them is more significant for the particles characterized with higher non-affine parameter. characterized by both elastically and plastically deformed zones at a microscopic level. This finding suggests that flows in sheared colloidal suspensions are highly heterogeneous.

Reciprocal space analysis
In many experimental techniques, such as neutron scattering, the information is obtained in reciprocal (momentum) Q space. Therefore, it is important to have a procedure to analyze the previous results in this space.
To investigate the crossover between the strong shear rate dependence of the structure factor along the gradient direction, and the weak shear rate dependence along the vorticity direction, it is useful to visualize the structure factor in some planes in the wave-vector space. Evidence of the anisotropy may be obtained by investigating the static structure factor S(Q) along specific directions in Q space: In figure 12, we present the S(Q) with the contribution of the form factor P(Q) in three planes: the flowvelocity gradient plane (1-2), the gradient-vorticity (1-3) plane, and the velocity-vorticy plane (2-3). Plane 1-2 (figures 12(a)-(c)) shows that the increase of the structure factor along the velocity direction is quite localized in the flow-velocity gradient plane. Moreover, we note that at the highest shear rate, the whole ring-like pattern becomes strongly deformed. For plane 1-3 (figures 12(d), (e)), we also observe an elliptical shape when we increase the shear rate. However, they have different angular dependence. In addition we show the plane 2-3 (figures 12(h), (i)). The 1-3 and 2-3 planes indicate particles are preferentially aligned along the vorticity plane during shear with increasing shear rate. Moreover, in figure 12(c), it is observed that the two additional spots are developed along the y-axis when P e =3.3. As demonstrated in figure A1 in the appendix, we show a picture of the particle trajectory at different planes at P e =3.3, and no discernible formation of layers along the flow direction is observed. We therefore conclude that the origin of the two spots is the difference in the density fluctuation along and perpendicular to the flow direction.
In the SANS experiment, the outputs are 2D patterns similar to those we plot in figure 6. In order to analyze such data, we adopt the spherical harmonics approach to describe the orientational dependence of the measured spatial distribution/correlation function. For the static structure factor, S(Q) is expanded by spherical harmonics as, One way is the direct Fourier transformation of ( ) g r 0 0 and -( ) g r .

2
The second one is to extract the information from the 2D patterns as in figure 6, which is explained in detail in the appendix.
In figure 13, we plot To correct this difference, it is necessary to add the contribution of ( ) In this way, the results obtained from the two different methods are the same.
The other component we need to obtain is -( ) S Q ,   where j 2 (Qr) is the spherical Bessel function with = ℓ 2. Although the exact mathematical form of γ(r) remains to be explored, based on the previous results of figure 7, it may be approximately described as: where γ 0 defines the maximum strain of the affinely deformed zone. The overall size of the affinely deformed zone is collectively determined by p, which is the position of the maximum of the affinely deformed zone and δ, which gives its width from p. Equations (20) and (21) provide the critical link between the macroscopic rheological responses due to microscopic structural evolution triggered by shear. The optimized γ(r), given by the solid curves in figure 15, is seen to be quantitatively consistent with that obtained from trajectory analysis. This agreement demonstrates the validity of this approach to determine the shear-induced flow heterogeneity. is due a factor of 15/8. Figure 15. The dependence of shear strain γ on r for several Péclet numbers for a sheared charged colloidal suspension. Solid lines represent a model strain calculated using equation (21). σ A is the particle diameter. Figure 16 gives the instantaneous shear strain g t = ¥ ( ) G , xy which is compared with the level of strain at the nearest neighbor γ NN as a function of Pe. When Pe is low, they are very similar, but this does not hold true when Pe is large. This indicates that the level of strain at the nearest neighbor, determined by anisotropic PDF, indeed represents the macroscopic instantaneous strain.

Summary
In summary, Brownian dynamics (BD) simulation on a concentrated colloidal suspension under shear indicates a non-uniform flow in the nonlinear shear thinning regime. An increasingly heterogeneous deformation is identified from the stress-dependent anisotropic structure factor. The observed structural behavior appears to be the manifestation of the transient elastic zone (TEZ) at mesoscopic length scale. A strong correlation between the average size of TEZ and decrease in the steady state shear viscosity is discovered. The connection of our observation to the origin of viscosity is also addressed. Now the phenomenological description of TEZ has been clarified based on our BD observation, perhaps the more difficult problem of the influence of hydrodynamic interaction on the viscoelastic properties of flowing colloids can be addressed by incorporating the mobility tensor into our BD or via a more explicit approach of dissipative particle dynamics. This implementation could extend our methodology to provide insights into nonlinear rheological responses that occur in a wide variety of highly interacting soft matter systems.
Clearly we can see that in x-z plane the ( ) S Q