# Heterogeneous shear in hard sphere glasses

###### Abstract

There is growing evidence that the flow of driven amorphous solids is not homogeneous, even if the macroscopic stress is constant across the system. Via event driven molecular dynamics simulations of a hard sphere glass, we provide the first direct evidence for a correlation between the fluctuations of the local volume-fraction and the fluctuations of the local shear rate. Higher shear rates do preferentially occur at regions of lower density and vice versa. The temporal behavior of fluctuations is governed by a characteristic time scale, which, when measured in units of strain, is independent of shear rate in the investigated range. Interestingly, the correlation volume is also roughly constant for the same range of shear rates. A possible connection between these two observations is discussed.

## Introduction.—

Heterogeneous flow and shear banding are central to the rheology of complex fluids and are widely observed in many industrial and natural materials, such as foams, emulsions, pastes, or even rocks olmsted_rheol_2008 (); schall_review2010 (). Despite its ubiquitous appearance, many aspects of this phenomenon are still not well-understood. In the simplest case, shear banding can be captured by a nonmonotonic dependence of the shear stress, , on the shear rate, FieldingOlmstedCates (); Fielding2009 (). For certain complex fluids, like colloidal gels, shear banding can be associated to the competition between a structural phase transition and a shear Moller (). However, in other systems, such as dense hard sphere (HS) colloidal suspensions Besseling_SCC () and granular materials Losert2000 (), flow heterogeneity is often observed without such accompanying structural changes. The rheological response of these systems is essentially determined by the competition between an inherent slow dynamics and the acceleration caused by the external drive Sollich1997 (); Berthier2002 (). This may lead to a spatially and temporally heterogeneous flow if the system is close to the yielding threshold VarnikPRL_03 (); VarnikJCP_static_yields (); Varnik2008 ().

Recently, it has been proposed that flow localization in dense hard sphere suspensions can be alternatively rationalized in terms of shear-concentration coupling (SCC) Besseling_SCC (), a well-known feedback mechanism for flow instability in complex fluids SchmittPRE95_banding (). In this model, regions of high (low) shear rate are associated with high (low) diffusivity, giving rise to a shear-induced flux from high to low shear-rate regions. This leads to an increase (decrease) of density and thus viscosity in low (high) shear-rate regions, thereby enhancing shear-rate fluctuations further. This, in turn, enhances the migration of particles from high to low shear-rate regions, unless balanced by diffusive counter flux. As pointed out recently furukawa (), enhancement of density fluctuations and the associated formation of voids can also be essential for the occurrence of fracture.

However, while the experimental data in Besseling_SCC () are interpreted consistently within the proposed macroscopic picture, no test of the basic underlying assumptions, such as the presence of a correlation between shear-rate and concentration fluctuations or the growth of emergent velocity fluctuations, has been provided so far for colloidal hard sphere glasses. This is not surprising, as the relevant density fluctuations that trigger the initial instability are quite small and hardly accessible to experiments. Furthermore, the available experimental time window for the observation of velocity fluctuations is limited to a few hundred percent strain, thus making a temporal analysis rather difficult.

Here, we study these and related issues via event-driven molecular dynamics (MD) simulations of a polydisperse hard sphere system. It is explicitly shown that fluctuations of the local volume fraction, , are correlated to the fluctuations of the local shear rate . More precisely, a decrease of local density is accompanied by an increase of the local shear rate and vice versa. However, while the relative amplitude of shear-rate fluctuations increases upon decreasing shear rate, the magnitude of the fluctuations of volume fraction remain constant. The temporal behavior of the fluctuations is characterized by a dominant time scale which is identical for both observables and . When measured in units of strain, this time scale is of the order of a few hundred percent strain and fairly independent of the imposed shear rate in the studied range. Interestingly, the correlation volume, as determined from the maximum of the four-point susceptibility of density fluctuations, , is also roughly constant in the investigated range ( is the particle number and is the incoherent scattering function at wave vector ). This suggests a possible link between an inhomogeneous shear and a dynamically heterogeneous environment with a given spatial correlation.

## Simulation setup.—

We perform event driven MD simulations of a polydisperse (11%) hard sphere system in 3D rapaport_book (). The volume fraction is varied from below to above the glass transition point, which, for the present polydisperse system, is located at a volume fraction of Pussey2009 (). The quiescent properties of this system have been studied extensively in Williams (). The temperature is fixed at via velocity rescaling. In all the simulations reported here, particles ^{1}^{1}1We have explicitly tested that, in the investigated parameter range, the obtained results are independent of the particle number. are placed in a random configuration between two walls. The walls are made of the same kind of particles as the bulk, but have infinite mass.
A constant shear is applied by moving the walls with velocities in the direction (planar Couette flow). The shear gradient points in the direction. The (overall) shear rate is defined as , where is the distance between the walls. Local quantities, such as velocity and density profiles, are computed as an average over particles within distinct layers of finite width parallel to the walls. Note that, for simplicity of notation, we will drop the subscript “av.” Local quantities will be identified by their arguments, e.g., and .

## Theoretical description.—

Central to the theory of shear-concentration coupling are a non-Newtonian constitutive relation for the shear stress and the notion of a nonequilibrium particle pressure SchmittPRE95_banding (). Both quantities are coupled through the Navier-Stokes equations. In Figs. 1(a) and 1(b), we show the dependence of the shear stress on shear rate and volume fraction in the yield stress regime, as obtained from our simulations. In the glassy phase (), our data can be well described by a Herschel-Bulkley expression

(1) |

where , , and are constants to be determined by a fit. is the reduced volume fraction ( corresponds to random close packing). The first term in Eq. (1) represents the dynamic yield stress Varnik2006 (), while the second term accounts for the effect of shear. The dynamic yield stress can be naturally associated with shear heterogeneity, if the rigid regions are understood to be locally below the yielding threshold, while the liquidlike regions are above CoussotPRL02_shearbanding (); VarnikPRL_03 (); schall_review2010 (); Besseling_SCC (). However, it does not explain the mechanism of the initial flow instability.

The particle pressure obtained from our simulations for different densities and shear rates is shown in Figs. 1(c) and 1(d). Noting that the pressure data show similar behavior to the shear stress, it appears natural to write (, , and are fit parameters)

(2) |

The first term constitutes the pressure contribution responsible for ordinary concentration diffusion Brady1995 (). The shear-rate dependence of the particle pressure is a manifestation of shear-induced particle migration, also known as “dilatancy” leighton_acrivos (); Deboeuf (); Haxton2011 ().

Within SCC, the feedback responsible for heterogeneous flow occurs if an initial density excess leads, via , to a locally lower shear rate, which in turn drives, via , further particle migration towards the low shear region. The flow instability develops as soon as these fluxes cannot be overcome anymore by the concentration and viscous momentum diffusion, described by and , respectively. A detailed calculation Besseling_SCC (); SchmittPRE95_banding () shows that this happens if

(3) |

In the interesting lower shear-rate regime, reduces to

(4) |

As seen from Figs. 1(a) and 1(c), for the glassy phase in the limit of low shear rates, shear stress as well as pressure show quite a similar dependence on shear rate. Indeed, trying various fit procedures revealed that leads to consistent fit results for all the simulated data. This suggests that is practically independent on . Using this information in , the critical shear rate for the onset of instability is obtained:

(5) |

## Heterogeneous shear.—

Figure 2(a) shows the velocity profile at successive strain intervals, illustrating the occurrence of flow heterogeneity in the studied polydisperse HS model in the glassy phase (). Note the qualitative similarity to experiments of colloidal hard spheres Besseling_SCC () and to simulations of a binary Lenard-Jones glass VarnikPRL_03 (). Here, we go a step further and perform a detailed survey of temporal fluctuations of both local shear rate, , and volume fraction, [Fig. 2(b)]. The data suggest that a positive is often accompanied by a negative and vice versa. Moreover, even though the temporal analysis is performed at considerably higher (due to limited computer time), the time scale of fluctuations in Figs. 2(a) and 2(b) seems to be roughly the same in units of strain.

We compute the statistical average of the instantaneous correlations between fluctuations of the local shear rate and volume fraction, . A typical result on is shown in Fig. 3(a) for an average shear rate of . We first note that most data points are distributed around a constant value indicated by a dashed line in the plot. This value is small but definitely nonzero (see the error bars, representing a standard deviation determined from independent runs). It is worth mentioning that the computational effort to obtain the present statistical accuracy is quite significant (40 independent runs, each with a duration of 60 days on a 3GHz CPU). The large values of at are a consequence of wall effects. Indeed, the volume fraction close to the walls is slightly (by about ) but systematically larger than the average volume fraction, leading to a corresponding decrease of wall shear rate. A detailed study of this issue will be presented elsewhere. For the reminder of this Letter, we will be concerned with bulk properties only.

To test the stability phase diagram, Eq. Eq. (5), via our simulations, we identify a heterogeneous flow by requiring max. As can be seen from Fig. 3(b), simulation and theoretical predictions are in reasonable agreement. Thus, SCC seems to describe at least the onset of instability quite well for the present system.

## Time scale of fluctuations.—

To elucidate the nature of the instability, we study layer-resolved fluctuations, and , and determine, via a Fourier analysis, the peak frequency and the associated characteristic time scale of fluctuations at for shear rates ranging from homogeneous to inhomogeneous flow regimes [dashed line in Fig. 3(b)]. We first note that this time scale is independent of if the distance to the walls is above a few particles in diameter (data not shown). More importantly, it is practically identical for both the observables and [Fig. 4(a)]. Moreover, when measured in the units of strain, it is independent of applied shear in the investigated range. Using the same set of data, we also determine the relative amplitude of fluctuations as a function of shear rate [Fig. 4(a)]. In marked contrast to the behavior of the characteristic time scale, this quantity depends on the specific observable. In particular, increases significantly with decreasing imposed shear rate, while is perfectly constant in the interesting limit of low shear rates, i.e., when entering the inhomogeneous flow regime.

The above observations have a number of consequences for a possible interpretation of flow heterogeneity. First, the presence of a correlation between fluctuations of local shear rate and volume fraction [Fig.3(a)] and the analysis of the stability diagram [Fig. 3(b)] are, at least in the first sight, in favor of SCC theory. However, the basic assumption of SCC theory is that fluctuations of shear rate give rise to fluctuations of volume fraction and vice versa. This implies that, as the instability is approached, both and should grow to some extent. This is at odds with the data shown in Fig. 4(a). To see this, we first note that, in contrast to the supercooled state, where a slight change of volume fraction can lead to considerable variations of viscosity as the glass transition is approached, the viscosity of a glass changes only little with volume fraction. This can easily be inferred from a comparison of shear stress for volume fractions below and above in Fig. 1(a). It is thus not plausible that tiny (in the sense of hardly detectable) changes of volume fraction can account for the observed increase of shear-rate fluctuation amplitude.

In this light, it is tempting to search for a link between the observed temporal behavior of and and dynamic heterogeneity. Without any claim for rigor, we start from the intuitive idea that some regions in the material are more mobile (i.e., appear more fluidlike) than others and therefore can support larger shear. These regions are not static. Rather, they continuously form and dissociate. The time scale of the fluctuations in the flow velocity is thus expected to reflect a time scale inherent to this microscopic dynamics. This inherent time scale should in turn be related to the characteristic size of these dynamically correlated regions.
If this idea is consistent, the constant time scale ^{2}^{2}2Here we assume that, in the present case, the strain is the relevant unit of time (note that thermal fluctuations are negligible). observed in our simulations would then imply a constant correlation volume. This is indeed borne out in Fig. 4(b), where we plot , whose maximum is a measure of the correlation volume Martens2011 (), as a function of shear rate. We remark that dynamic correlations are not visible in the fluctuations of one-particle density so that, based on this idea, no conclusion can be made as to how local density fluctuations should behave. Therefore, the different behaviors of the fluctuation amplitudes of shear rate and density shown in Fig. 4a are not a priori in conflict with this picture.

In conclusion, we find that, although the theory of shear-concentration coupling describes some aspects of flow heterogeneity in hard sphere glasses, serious inconsistencies remain, such as a significant increase of the relative amplitude of shear-rate fluctuations upon approaching the instability, as opposed to constant fluctuations of volume fraction. A more stringent test of the SCC theory would be to numerically solve the underlying Navier-Stokes equations using the experimental (or simulated) and as input benzi_glasses (). An important question to be answered would be whether the original SCC theory predicts stable shear bands or a fluctuating behavior as observed in our simulations. Furthermore, keeping in mind that solvent effects on the system response are negligible in the limit of high and low , it would also be interesting to work out the consequences of shear-concentration coupling by solving the compressible Navier-Stokes equation (complemented by a non-Newtonian constitutive law) only for the colloids. When seeking for alternative routes, on the other hand, it must be kept in mind that, in its present form, the proposed SCC theory for a HS glass is a hydrodynamic approach with a fully local constitutive law. Dynamic correlations, however, give rise to nonlocal effects. These could possibly be taken into account (e.g., by means of a correlation length of some measure of “fluidity” Goyon2008 ()) by introducing generalized hydrodynamic equations.

###### Acknowledgements.

S.M. is financially supported by the Max-Planck Society. We thank R. Besseling and A. Donev for useful discussions. Boundary-driven simulations are performed using the code by D. C. Rapaport rapaport_book (). We are also very grateful to M. Bannerman for adding polydispersity to DynamO DynamO (), which allowed us the study of homogeneous flow.## References

- (1) P. D. Olmsted, Rheol. Acta 47, 283, (2008).
- (2) P. Schall and M. van Hecke, Ann. Rev. Fluid Mech. 42, 67, (2010).
- (3) S. M. Fielding and P. D. Olmsted, Eur. Phys. J. E 11, 65 (2003); M. E. Cates and S. M. Fielding, Adv. Phys. 55, 799 (2006).
- (4) S. M. Fielding, M. E. Cates and P. Sollich, Soft Matter 5, 2378 (2009).
- (5) P. C. F. Moller, S. Rodts, M. A. J. Michels and D. Bonn, Phys. Rev. E 77, 041507 (2008).
- (6) R. Besseling, L. Isa, P. Ballesta, G. Petekidis, M. E. Cates and W. C. K. Poon, Phys. Rev. Lett. 105, 268301 (2010).
- (7) W. Losert, L. Bocquet, T. C. Lubensky and J. P. Gollub, Phys. Rev. Lett. 85, 1428 (2000).
- (8) P. Sollich, F. Lequeux, P. Hebraud and M. E. Cates, Phys. Rev. Lett. 78, 2020 (1997).
- (9) L. Berthier, J.-L. Barrat and J. Kurchan, Phys. Rev. E 61, 5464 (2000).
- (10) F. Varnik, L. Bocquet, J.-L. Barrat and L. Berthier, Phys. Rev. Lett. 90, 095702 (2003).
- (11) F. Varnik, L. Bocquet and J.-L. Barrat, J. Chem. Phys. 120, 2788 (2004).
- (12) F. Varnik and D. Raabe, Phys. Rev. E 77, 011504 (2008).
- (13) V. Schmitt, C. M. Marques and F. Lequeux, Phys. Rev. E 52, 4009 (1995).
- (14) A. Furukawa and H. Tanaka, Nature 443, 434 (2006); Nat. Mater. 8, 601 (2009).
- (15) D. C. Rapaport, The Art of Molecular Dynamics Simulation (Cambridge University Press, 1995).
- (16) P. N. Pusey, E. Zaccarelli, C. Valeriani, E. Sanz, W. C. K. Poon and M. E. Cates , Philos. Trans. R. Soc. London, Ser. A 367, 4993 (2009).
- (17) S. R. Williams, I. K. Snook and W. van Megen, Phy. Rev. E 64, 021506 (2001).
- (18) F. Varnik and O. Henrich, Phys. Rev. B 73, 174209 (2006).
- (19) P. Coussot, J. S. Raynaud, F. Bertrand, P. Moucheront, J. P. Guilbaud, H. T. Huynh, S. Jarny and D. Lesueur, Phys. Rev. Lett. 88 218301 (2002).
- (20) J. F. Brady and M. Vicic, J. Rheol. 39, 545 (1995).
- (21) A. Deboeuf, G. Gauthier, J. Martin, Y. Yurkovetsky and J. F. Morris, Phys. Rev. Lett. 102, 108301 (2009).
- (22) D. Leighton and A. Acrivos, J. Fluid. Mech. 177, 109 (1987).
- (23) T.K. Haxton, M. Schmiedeberg and A. Liu, Phys. Rev. E 83 031503 (2011).
- (24) K. Martens, L. Bocquet and J.-L. Barrat, Phys. Rev. Lett. 106 156001 (2011).
- (25) R. Benzi, M. Bernaschi, M. Sbragaglia and S. Succi, EPL 91, 14003 (2010); Phys. Rev. Lett. 106, 164501 (2011).
- (26) J. Goyon, A. Colin, G. Ovarlez, A. Ajdari, L. Bocquet, Nature 454, 84 (2008).
- (27) M. Bannerman, R. Sargant and L. Lue, J. Comp. Chem. 32, 3329 (2011).