Transient stress evolution in repulsion and attraction dominated glasses
aa r X i v : . [ c ond - m a t . s o f t ] J un Transient stress evolution in repulsion and attraction dominated glasses
Christian P. Amann a) and Matthias Fuchs b) Fachbereich Physik, Universität Konstanz, 78457 Konstanz, Germany (Dated: February 18, 2018)
We present results from microscopic mode coupling theory generalized to colloidal dispersions under shear inan integration-through-transients formalism. Stress-strain curves in start-up shear, flow curves, and normalstresses are calculated with the equilibrium static structure factor as only input. Hard spheres close totheir glass transition are considered, as are hard spheres with a short-ranged square-well attraction at theirattraction dominated glass transition. The consequences of steric packing and physical bond formation onthe linear elastic response, the stress release during yielding, and the steady plastic flow are discussed andcompared to experimental data from concentrated model dispersions.PACS numbers: 82.70.Dd, 83.60.Df, 64.70.PfKeywords: Colloids, Nonlinear rheology, Glass transition
I. INTRODUCTION
Colloidal dispersions have been established as model systems in materials science. They behave like fluids at highdilution, and form condensed phases if particle interactions dominate over entropic disorder. Equilibrium statisticalmechanics explains equilibrium phases and their near-equilibrium properties based on direct particle interactions.While purely repulsive interactions lead to fluid and crystalline solids, attractions can also give rise to liquids. Ascolloidal dispersions offer the unique possibility to tailor the depth and range of the attractions, the conditionswhen liquid phases become stable could be established in studies e.g. by Lekkerkerker et al. (1992). Yet colloidaldispersions also form various metastable solid states like fractal networks, particle gels, and glasses, which can notbe described with purely equilibrium statistical mechanics. They are important in industrial processes and products,and in biological systems, which generally are far from thermal equilibrium. The specific mechanical properties of theresulting soft solids are often tailored adjusting the competition between attractive and repulsive interactions. Varioussolid states, with very different elastic properties, have been prepared; recent studies include works by Lu et al. (2008),Gibaud et al. (2011), Kim et al. (2013), and others.Mixtures of colloids with non-adsorbing polymers constitute one of the simplest systems, where effective interactionsamong the (larger) colloidal constituents can be adjusted in a controlled way. Two different glass states, viz. amorphoussolids formed by the freezing-in of the cooperative structural dynamics, could be observed conclusively by Pham et al. (2002) and Eckert and Bartsch (2002), as had been predicted by Fabbian et al. (1999) and Bergenholtz and Fuchs(1999) using mode coupling theory. Controlling the range of the attraction proved crucial for the formation of’attraction driven glasses’, where particles form physical bonds to their neighbors. Attraction driven glasses requireshort ranged attractions so that particles become tightly localized, as can be clearly seen in computer simulations(Puertas et al. , 2002). In ’repulsion driven glasses’, the cooperative behavior of the neighbors forming the cagelocalizes particles less tightly. The localization length there corresponds to Lindemann’s length (Lindemann, 1910)which is roughly a tenth of the average particle separation (Hansen and McDonald, 2009). Lindemann had foundthat atomic displacements increase up to this value when (crystalline) solids are heated to their melting temperature.Mode coupling theory as developed by Götze (2009) and others established that Lindemann’s length also characterizesthe frozen-in structure of repulsion driven glasses, as was verified experimentally by Pusey and van Megen (1987) andvan Megen and Underwood (1993) in colloidal hard sphere dispersions. As attractions change the structure of theseglasses only if they are of short range — mode coupling theory predicts that qualitative changes require attractionranges (somewhat) shorter than Lindemann’s length — and as attractions in molecular systems act across longerranges, molecular glass transitions fall into the class of repulsion driven transitions. This has made hard spherecolloids a model system for studying the glass transition in molecular, metallic, and simple supercooled liquids, whichwas recently reviewed by Hunter and Weeks (2012). Attraction driven glasses require attraction strengths of the orderof only a few k B T (thermal energies) and can form at particle concentrations lower than repulsive ones (Dawson et al. ,2000). The aspect that increasing the strength of attraction at first destabilizes and melts the (repulsion driven) glass, a) Electronic mail: [email protected] b) Electronic mail: [email protected] provides insight into the mechanism of vitrification and its theoretical description. For any strength, short rangedattractions increase the equilibrium probability of close-particle contacts, viz. the contact value of the pair correlationfunction. Yet, this increased stickiness reduces the medium-ranged order in the disordered fluid as measured in theprincipal peak of the equilibrium structure factor. Mode coupling theory predicts that this decrease destabilizesthe repulsion driven cages so that a hard sphere glass melts upon turning on a narrow attraction. This causes anappreciable shift of the glass transition to higher concentrations. In a recent investigation, Willenbacher et al. (2011)achieved a shift of the transition packing fraction by more than 10 %. For fixed concentration, short ranged attractionsthen cause a glass transition when at increasing attraction strength large wavevector contributions in the equilibriumstructure factor contribute strongly. Inbetween, a reentrant fluid region lies where equilibrium is approached at longtimes. The continuation of this states diagram to lower concentrations, where colloidal gels are observed (Lu et al. ,2008), remains an active area of research. E.g. computer simulations investigate it by following the temporal evolutionof the average localization length of individual particles as function of concentration (Zaccarelli and Poon, 2009).Metastable glassy systems exhibit a complex interplay between their structural relaxation and flow causing astrongly non-Newtonian response. Moreover, the structural and mechanical changes during or after processing usingflow or other external fields also are crucial for achieving desired material properties. Thixotropy and ageing areimportant, and are moving into the reach of first principles theories only recently (Ballauff et al. , 2013). While anumber of studies exists of the mechanical response of dispersions approaching the repulsion driven glass transition— a recent one (Siebenbürger et al. , 2009) where the transition density could be approached very closely is reviewedby Siebenbürger et al. (2012) — there are far fewer studies of the nonlinear rheology in the complete region coveringrepulsion and attraction driven glass transitions. Besides the already mentioned work by Willenbacher et al. (2011)combining rheology and dynamic light scattering, especially the seminal investigations by Pham et al. (2006) andPham et al. (2008) provided insights into the nonlinear mechanical behavior under different rheological protocols andvarying repulsive and attractive effects by varying concentration and attraction strength. As also seen in computersimulations by Puertas et al. (2007), the linear shear moduli are far larger at the attraction driven glass transitionthan at the repulsion driven one. In mode coupling theory this arises from the dominance of large wavevectorsin the approximated Green-Kubo relation for the linear moduli, which predicts an increase with the square of theinverse of the relative attraction range in the sticky limit (Bergenholtz et al. , 2000). Another intriguing observationby Pham and coworkers concerns the yielding of glasses under applied strain. A typical strain of around 10% to 20 %characterizes the shear-induced yielding of repulsion driven glasses when experiencing step strain and start-up flow,while attraction driven glasses yield in a two-step process. The former observation, in agreement with light scatteringstudies under large amplitude oscillatory shearing by Petekidis et al. (2002), nicely ties to the picture of the cage effectand its characteristic length following Lindemann’s criterion; yet this connection has not been established theoreticallyup to now. The second observation, already visible in oscillatory shear experiments by Gadala-Maria and Acrivos(1980) and supported by a detailed investigation under start-up shear flow by Koumakis and Petekidis (2011), cannot be interpreted so easily by the cage picture because the two characteristic strain values are around 10% and 100%. The second value is far larger than the local cage picture would imply, and the attraction range, which was around5% in the experimental system, appears not to characterize the stress-strain relations.Aim of the present contribution is to determine the stress-strain relations close to repulsion and attraction drivenglass transitions from mode coupling theory (MCT) as generalized to sheared colloidal dispersions in the integration-through-transients (ITT) framework. We will consider the quintessential repulsive glass transition, viz. the one in ahard sphere fluid, and a typical attraction driven one for a narrow square-well pair potential. Specifically, we willconsider start-up shear flow with fixed shear rate and determine the transient shear stress as function of accumulatedstrain. In the generalization of the microscopic mode coupling theory developed by Brader et al. (2012) based on thework by Fuchs and Cates (2002), the structural relaxation under arbitrary, homogeneous, and incompressible flows isdeduced from the equilibrium structure factor so that caging and bond formation, as in the quiescent situation, canbe discussed.Our work bears similarity to the study by Henrich et al. (2009) where hard disks were considered in two dimensions.Here we present the first calculations within microscopic MCT-ITT for hard spheres in three dimensions and addition-ally consider attractive glasses in the second part. Our work also bears some similarity to the one by Miyazaki et al. (2004), who, however, concentrated on fluid states under shear and on time-dependent fluctuations around the steadystate. We focus here on the transient dynamics of yielding glass states and their stationary, time-independent prop-erties which are not accessible to the theory by Miyazaki et al. (2004). Our work also bears similarity to the studyby Priya and Voigtmann (2014) in the present volume of the Journal of Rheology, who also consider the non-linearrheology of repulsion and attraction dominated glass transitions. They use a simplified MCT-ITT, where shear de-formations are isotropically averaged, which enables them to study wider variations in the attraction ranges andstrengths than possible in our solution of MCT-ITT without additional approximations.In Sect. II the pertinent equations of mode coupling theory are summarized. Section III gives an overview of thestudied systems and the glass states diagram, while Sects. IV and V describe the results for hard spheres without andwith square-well attraction, respectively. Section VI concludes with a comparison of the findings with experimentaldata.
II. NONLINEAR RHEOLOGY WITH MODE COUPLING THEORY
The integration-through-transients (ITT) formalism which generalizes mode coupling theory (MCT) to drivensystems provides a method to calculate the complete time evolution of a concentrated dispersion under homogeneousstrain deformation. This yields more information than e.g. just obtaining the steady state properties as the evolutionfrom elastic to plastic response can be observed. The MCT-ITT approach by Brader et al. (2012) is presently restrictedto incompressible flows and neglects hydrodynamic interactions. We will consider start-up shear flows, which is asimple time-dependent deformation protocol where the former condition is obeyed, and high particle concentrations,where the dominance of structural correlations motivates our neglect of hydrodynamics. Solvent effects are presumedto renormalize the hydrodynamic radii and short time diffusion coefficients.Microscopic starting point is the Smoluchowski equation for interacting Brownian particles in a given shear flow. Theparticles’ time evolution follows from affine motion with the flow and random motion causing non-affine displacements,both combined in the Smoluchowski operator:
Ω = N X i =1 ∂∂ r i · (cid:20) ∂∂ r i − F i − ˙ γ y i ˆ x (cid:21) . (1)The force on particle i derives from a potential, where the chosen pair interaction will enter in later sections. Dimen-sionless units are used, where length, energy and time are measured in units of particle diameter d , thermal energy k B T , and d /D , respectively. The effect of shear relative to Brownian motion is measured by the bare Péclet numberPe = ˙ γd /D , which in these units agrees with the shear rate. Note that the Weissenberg number Wi = ˙ γτ , with τ an intrinsic α -relaxation time scale of a fluid, is also called (dressed) Péclet number Pe. Viscoelastic response can beobserved at Pe ≪ and Wi & , while the response for Wi ≪ is the one of a Newtonian fluid.An equation of motion for a transient density correlator Φ q ( t ) encodes rapid local motion without structuraldecay, and elasticity and plasticity owing to structural rearrangements. The transient density correlator Φ q ( t ) = h ρ ∗ q ρ q ( t ) ( t ) i /N S q , is the correlation function built with density fluctuation, ρ q = P Nj =1 exp { i q · r j } . Their timeevolution is given by the adjoint of the Smoluchowski operator from Eq. (1), ρ q ( t ) = e Ω † t ρ q . In ITT, the average canbe performed over the equilibrium Gibbs-Boltzmann ensemble, because the system is assumed to be in equilibriuminitially. Thus the normalization, giving Φ q (0) = 1 , is done with the equilibrium structure factor S q . The shear-advected wavevector q ( t ) = ( q x , q y − ˙ γtq x , q z ) T appearing in the definition accounts for the affine particle motion withthe flow, and gives Φ q ( t ) ≡ in the absence of non-affine motion. Random motion, affected by the shear flow, causes Φ q ( t ) to decay. In MCT-ITT this is given by an equation of motion containing a retarded friction kernel which arisesfrom the competition of particle caging and shear advection of fluctuations ˙Φ q ( t ) + Γ q ( t ) (cid:26) Φ q ( t ) + Z t dt ′ m q ( t, t ′ ) ˙Φ q ( t ′ ) (cid:27) = 0 . (2)The initial decay rate contains Taylor dispersion as Γ q ( t ) = q ( t ) /S q ( t ) . The generalized friction kernel m q ( t, t ′ ) isan autocorrelation function of fluctuating stresses. Based on the insights of quiescent MCT as described by Götze(2009), it is approximated by a quadratic polynomial in the density correlators m q ( t, t ′ ) = Z d k (2 π ) nS q ( t ) S k ( t ′ ) S p ( t ′ ) q ( t ) q ( t ′ ) V qkp ( t ) V qkp ( t ′ ) Φ k ( t ′ ) ( t − t ′ ) Φ p ( t ′ ) ( t − t ′ ) (3)with abbreviation p = q − k , and n the particle density. The vertex function is given by V qkp ( t ) = q ( t ) · (cid:0) k ( t ) c k ( t ) + p ( t ) c p ( t ) (cid:1) (4)where c k is the Ornstein-Zernicke direct correlation function, c k = (1 − /S k ) /n . These equations of motion werederived in detail by Fuchs and Cates (2009) using a Zwanzig-Mori projection-operator formalism together with modecoupling approximations. The equilibrium structure factor, S k , encodes the particle interactions and introducesthe experimental control parameters like density and temperature. Again generalizing quiescent MCT to flow, thepotential part of the stress σ αβ ( ˙ γ ) = h σ αβ i ( ˙ γ ) /V is approximated assuming that stress relaxations can be computedfrom integrating the transient density correlations σ xy ( t ) = ˙ γ Z t dt ′ Z d k π ) " k x k y ( − t ′ ) k y kk ( − t ′ ) S ′ k S ′ k ( − t ′ ) S k Φ k ( − t ′ ) ( t ′ ) . (5)Equation (5) is the central equation which we will evaluate for different systems. It can also be used to derive aconstitutive equation for the shear stress of the form σ xy ( t ) = ˙ γ R t dt ′ g xy ( t ′ , ˙ γ ) with a generalized shear modulus g xy ( t, ˙ γ ) = 12 Z d k (2 π ) " k x k y k y ( t ) k ( t ) k S ′ k ( t ) S ′ k S k ( t ) Φ k ( t ) . (6)Here a shift of the integration variable was performed collecting the wavevector advection in a time-dependent vertex(the term in the square bracket) which is weighted with the square of the density correlator at wavevector k . Withoutwavevector advection, g xy ( t, ˙ γ = 0) recovers the quiescent MCT expression for the stress autocorrelation function(Götze, 2009); see a review by Fuchs (2010) for more discussions. For an ideal elastic solid, the modulus g xy wouldbe constant and stress and accumulated strain ˙ γt would be proportional. If g xy ( t ) does not depend on ˙ γ and decayson an intrinsic time scale τ , the finite time integral over g xy ( t ) is the long time viscosity η xy of a Newtonian fluid,stress and shear rate are then proportional. This recovers Maxwell’s model of linear response. Viscoelastic mediaexhibit a non-linear behavior in ˙ γ , because of a ˙ γ functionality of g xy ( t, [ ˙ γ ]) . MCT can obviously provide a microscopicdescription of such viscoelasticity. Equation (6) will be studied in the following sections specifying different systemsby implementing different equilibrium structure factors. III. OVERVIEW OF CONSIDERED REPULSIVE AND ATTRACTIVE SYSTEMS AND THEIR QUIESCENT GLASS STATES
In MCT-ITT, the structural relaxation is determined by the direct correlation function c q (equivalently the equilib-rium structure factor S q ) which describes the effective interaction between density fluctuations. In density functionaltheory it gives the quadratic term in the interaction part of the free energy functional, in MCT it enters the interac-tion vertices when fluctuating stresses are connected to density fluctuations. As discussed in detail by Dawson et al. (2000), there are two major mechanisms for vitrification in the self-consistency equations of MCT for pair potentialsconsisting of an excluded volume core and a short-ranged attraction. They can be discerned from the wavevectorrange where the major contributions arise in the memory kernel of Eq. (3).The normal situation is that the principal peak in S q dominates, and that the glass transition is crossed whenit becomes high. This glass transition originates in the local ordering caused by the excluded volume constraintsin dense fluids, the transition leads to a ’repulsion dominated glass’ (RDG), and the width of the glass form factor f q in reciprocal spaces is set by the localization length following Lindemann’s criterion. The inset of Fig. 1 showsthe glass form factor f cq at the transition of a fluid of monodisperse hard spheres (HS), which is the most simpleexample for this transition (Götze, 2009). Solving Eqs. (2) to (4) numerically at vanishing shear rate and usingthe Percus-Yevick approximation for the structure factor, the transition lies at packing fraction φ c = 0 . ,which is somewhat below the value φ c = 0 . established in experiment for slightly polydisperse hard sphere colloids(Pusey and van Megen (1987), Hunter and Weeks (2012)). The present numerical result for φ c and especially f cq (e.g. the small wiggles close to the second and third peaks) differ from more precise calculations reviewed by Götze(2009) because the chosen discretization of the wavevector integrals in Eq. (3) is optimized to handle the anisotropicdistortions under shear, including with attractions. Presently we cannot chose a finer discretization, because a singlestress vs strain curve (from Fig.4) takes 90 hours running time on a modern CPU.A second mechanism causes arrest of the structural relaxation when a short-ranged attraction is strong enough ina colloidal dispersion. A square-well potential (SWP) of width ∆ and attractive depth u is chosen here to exemplifythis. It acts outside the excluded volume core. This study extends the one by Dawson et al. (2000) by includingshear flow. The potential depth will be made dimensionless using thermal energy k B T , i.e. U = u /k B T . Relativeattraction ranges δ = ∆ /d smaller than Lindemann’s ratio, δ < . are required which cause strong contributions inthe memory kernel of Eq. (3) at large wavevectors. The ’attraction dominated glass’ (ADG) transition takes placewhen wavevectors of the order k ∼ / ∆ dominate. The resulting glass form factors f q extend to large wavevectors andtheir q -dependent width corresponds to a localization length of the order of the attraction range ∆ . This is shown inthe inset of Fig. 1, where the form factors f cq are shown at the attraction driven glass transition at U c = 4 . and δ = 0 . at the same packing fraction φ c as the HS transition; the different widths at the HS and ADG transitionsare apparent.The complete glass transition lines (states diagram) in Fig. 1 for two attraction ranges verify that the first effectwhen turning on a short-ranged attraction is a destabilization of the repulsion dominated glass. The glass transitionline moves to higher packing fractions initially with increasing attraction strength U . The effect is stronger for smallerrelative attraction range δ . A reentrant fluid region emerges which extends up to a maximal packing fraction, aroundwhere the two glass transition lines merge or intersect. The scenario depends on the precise attraction range. For δ = 0 . both transition lines merge in an A singularity, while for shorter ranges both glass transition lines intersectat a crossing point, and the ADG ends at an A singularity (Dawson et al. , 2000). U f d = 0.035 d = 0.0465 q [1/ d ] f q c U = 0 f = f c U = U c f = f c Figure 1.
Main panel:
States diagram showing the critical glass-transition interaction strength U c = u /k B T and packingfraction φ c for hard spheres with square-well attraction. The solid lines are taken from Dawson et al. (2000); see the legendfor the relative attraction range of the SWP δ = ∆ /d , with ∆ attraction range and d particle diameter. A red star marks thehigher order glass transition point A at U = U A . The symbols mark ADG calculations, with q max = 79 . ( black squares andgreen diamonds ) or q max = 119 . ( red circles ). Calculations (marked by symbols × ) along two paths crossing a glass transitionwill be discussed latter: Crossing the HS transition increasing φ for U = 0 , and crossing the ADG increasing U at φ c . The inset shows (isotropic) non-ergodicity parameters f cq at the two transitions, viz. at the critical packing fraction φ c for U = 0 (hard-sphere case) and U = U c (attraction driven transition). In the present study, the transient stress evolution close to two typical transitions of both kinds shall be explored.Numerical calculations will be performed for two paths, one varying the concentration at fixed (vanishing) attractionstrength, and the other one at fixed packing fraction for increasing well depth. Equilibrium structure factors, whichare the only input to the theory, will be taken from the work by Dawson et al. (2000). While the choice of U = 0 is quite natural and leads to the pure HS system, qualitatively similar results would hold for all RDG transitionsat fixed U (much) smaller than the attraction strength of the A singularity. The second transition to the ADGis chosen for the same packing fraction as the HS transition in order to eliminate density dependent differences.According to MCT qualitatively similar results would hold for all φ below the density of the A singularity. In orderto achieve a clear time scale separation between the slow structural relaxation of interest and the fast local dynamics,small relative separations from the glass transitions shall be considered. They will be expressed using the relativeseparations ε = ( φ − φ c ) /φ c and ε U = ( U − U c ) /U c , respectively. IV. HARD-SPHERE TRANSITION
After the presentation of the MCT equations and a short overview of the quiescent glass states diagram, the yieldingof hard sphere glasses under applied shear strain shall be discussed first. It gives the example relevant to concentrateddispersions without especially short ranged attractions, and to molecular and metallic glass-forming liquids.Equation (5) is a constitutive equation for the shear stress undergoing simple shear after start-up at t = 0 . Thevertex (the term in the square bracket) depends on time only via the accumulated strain, γ = ˙ γt . The squaredcorrelators depend on time and accumulated strain independently in general and thus the stress-strain relationsobtain different forms depending on the bare Péclet number Pe and the distance to the glass transition ε . Thenumerical calculation proceeds by first solving the self-consistency equations (2) to (4) for the density correlators Φ q ( t ) at given shear rate ˙ γ , and second by integrating Eq. (6) to obtain the generalized shear modulus. If desired,flow-curves and stress strain relations can then be obtained.Figures 2 and 3 provide a good qualitative illustration of the transient stress regime, as well as already the quanti- Pe = 10 -1 -2 -3 -4 -5 -6 PSfrag replacements a) ˙ γt σ x y (cid:2) k B T / d (cid:3) G c ∞ = 10 -1 -2 -3 -4 -5 -6 PSfrag replacements b) ˙ γt σ x y (cid:2) k B T / d (cid:3) G c ∞ = 10 -1 -2 -3 -4 -5 -6 PSfrag replacements c) ˙ γt σ x y (cid:2) k B T / d (cid:3) G ∞ = 10 -4 e = 10 -2 -3 -3 -2·10 -3 -10 -2 PSfrag replacements d) ˙ γt σ x y / σ s t x y Figure 2. Transient shear stress σ xy as function of accumulated strain ˙ γt for several separation parameters ε = ( φ − φ c ) /φ c andbare Péclet numbers Pe calculated for a hard sphere system (see legends for Pe ). Panel a) is in the fluid, ε = − − , b) atthe transition, ε = 0 + , and c) in the glass, ε = 10 − . The elastic moduli G ∞ and G c ∞ , calculated from the plateaus of Fig. 3,are shown as gray, dash-dotted lines . Panel d) shows stress curves rescaled by the steady state value varying ε at constantPe = 10 − . tative data. Figure 2 shows the shear stress vs strain. Figure 3 shows the generalized shear modulus. Three regimesof viscoelasticity can be identified from the curves: ( i ) First the stress grows linearly with strain proportional to an elastic shear modulus G ∞ . This expresses Hooke’slaw, σ xy = G ∞ γ . Figure 3 shows that this corresponds to a plateau in g xy ( t, ˙ γ ) for intermediate times. It is reachedafter some shear-rate independent short-time relaxation, and holds up to the final decay time. This plateau has thevalue G ∞ . The index ∞ can be understood to refer to an infinite intrinsic relaxation time, which characterizes anelastic solid. In fluid states, the intrinsic final structural or α -relaxation causes a decay of g xy ( t, ˙ γ = 0) , which softensthe stress-strain relations. See Fig. 2 a), where σ xy ( t ) becomes smaller at low shear rates, i.e. it becomes plasticbecause ˙ γτ < . There is a loss of memory in the system, causing deviations from the linear elastic and leading toa viscous response. At times long compared to the intrinsic α -relaxation time τ , a Newtonian viscosity is observed, σ xy ( t ≫ τ ) → η xy ˙ γ . Above the glass transition, MCT predicts that the glass structure is persistent and that theplateau does not decay (ideal glass), which is an idealization which is not observed experimentally. Nevertheless,the time-window where the shear modulus is nearly time-independent can be made arbitrarily large by supercoolingfurther. In MCT, the frozen-in glass structure is the reason for solid elasticity. Increasing the shear rate, the shear-distorted structure of the hard spheres inside their yielding structural cages stores additional stress. This increase of σ xy ( γ ) at small strains is an ’anelastic’ effect in the β -process of MCT, which describes the instability of the cagestrapping the particles (Voigtmann et al. , 2012); this can be seen best in panels b) and c) of Fig. 2. ( ii ) Shear induced decay melts the glass, viz. the G ∞ -plateau decays after a decay time proportional to the inverseof the shear rate; see Fig. 3. Consequently the time integral in Eq. (5), viz. the area under g xy ( t ) vs time, does notincrease anymore. A steady state value of σ xy ( t ) is reached in Fig. 2 at strains of order one, denoted as flow curvevalue σ stxy ( ˙ γ ) ; see Fig. 4 for the flow curve. The approach to a steady flow curve requires plastic effects, which causethe final decay of g xy ( t ) . For increasing Pe , the plateau decay starts earlier and the short-time relaxation becomesmore important. ( iii ) Between the steady-state regime of the flow curve and the elastic regime occurs a transient plastic regime. FromFig. 3 can be verified that g xy ( t ) takes negative values prior to its final decay to zero. This negative area under the g xy ( t ) curve adds a negative portion to σ xy ( t ) so that the stress decreases onto the steady-state plateau. The emergingbump in Fig. 2 is called stress overshoot. The stress is maximal for the peak strain value γ ∗ , i.e. Eq. (5) identifies γ ∗ as the zero of g xy ( γ = ˙ γt ) . If, in the fluid phase, the structural relaxation time τ is smaller than the shear induced one,the stress overshoot vanishes; Fig. 2a) illustrates this. This agrees with a result from linear response theory, viz. thatthe equilibrium shear modulus g xy ( t, ε < is completely monotone (Götze, 2009), and is observed in experiments oncolloidal dispersions by Koumakis et al. (2012a), Koumakis et al. (2012b), and Amann et al. (2013). -3 -2 -1 G ∞ c G ∞ ( ε = -3 ) c Pe = 10 -1 -2 -3 -4 -5 -6 G ∞ c G ∞ ( ε = -3 ) ε = -10 -3 -3 t (cid:2) d /D (cid:3) g x y ( t , ˙ γ ) (cid:2) k B T / d (cid:3) ˙ γt g x y ( t , ˙ γ ) (cid:2) k B T / d (cid:3) Figure 3. Generalized shear modulus g xy ( t, ˙ γ ) , Eq. (6), as function of time (main panel) and accumulated strain ( inset ; setsfor every second Pe left out for clarity). The legend provides color coded the strain rates and line-style coded the separationparameters ε ; a letter c labels the critical g xy ( t, for ε = 0 + . In the main panel, fluid curves at Pe = 0 and − overlap,in the inset glass curves at Pe = 10 − and − overlap. Elastic shear moduli G ∞ ( ε ) can be read off from quiescent curves(Pe = 0 ), with G ∞ (0 + ) = G c ∞ = 18 . k B T /d and G ∞ (10 − ) = 21 . k B T /d . The inset shows subtle differences in thetypical strain γ ∗ , where g xy = 0 , when shear-driven and internal relaxation in Φ q ( t, ˙ γ ) interfere. This becomes most clear forPe = 10 − in the fluid phase, where the internal relaxation dominates and the undershoot disappears. Figure 4 shows the flow curve values σ stxy ( ˙ γ ) = σ xy ( t → ∞ , ˙ γ ) of the steady-state regime of the stress-strain curvesfrom Fig. 2 and also the long-time shear viscosity η xy = σ stxy / ˙ γ . The difference between fluid phase and glass phasein the context of MCT can be seen. In the glass phase and for Pe → , the area under g xy ( t ) becomes proportionalto / ˙ γ . In consequence, a constant dynamic yield stress σ + xy can be read off directly from the flow curve for vanishingshear rate, σ + xy = σ stxy ( ˙ γ → , ε ≥ . This stress is necessary to keep the glass yielding and flowing at infinitesimallysmall shear rates. Because this yield stress is non-zero, the Newtonian viscosity diverges in the glass. Increasing theshear rate, processes on intermediate time scales (including the β -process in MCT studied by Götze (2009)) cause thestress to become larger. The β -process is slowest close to the glass transition, and therefore the flow curve at ε = 0 varies sensitively with shear rate already at very small Pe . For a schematic model of MCT-ITT, Hajnal and Fuchs(2009) deduced a Herschel-Bulkley law at the transition, which unfortunately cannot be tested in our calculationsbecause of the coarse discretization. Deeper in the glass, the flow curve stays rather constant, because the β -processhas become faster, and σ + xy can be observed for a wider window in Pe . It must be distinguished from a static yieldstress, which a static load would have to overcome to fluidize a glass. One could within this context try to identifythe peak-stress observed during the overshoot as static yield stress, but the corresponding shear-protocol is shear-rateand not stress controlled. The static yield stress could be (and is) different if one increases the stress in a controlledmanner until the glass starts to flow. The issue of creep does then arise, which can be modeled with a stress controlledMCT as shown by Siebenbürger et al. (2012).If in the fluid phase ( ε < ), the shear-rate independent structural decay characterized by the quiescent α -time τ takes place much earlier than the shear-induced one at time / ˙ γ , the integral in Eq. (5) becomes shear-rate independentand leads to the Newtonian viscosity η xy . The linear relation between stress and shear rate defines a Newtonian fluid, -3 -2 -1 -6 -5 -4 -3 -2 -1 s s t xy [ k B T / d ] Pe h ¥g ˙ e = -10 -2 -10 -3 -3 -2 slope 110 -5 -3 -1 Pe h xy [ k B T / ( D d ) ] e eff = -0.050.0050.020.03 10 -5 -3 -1 Pe h xy [ k B T / ( D d ) ] e eff = -0.050.0050.020.03 Figure 4. The main panel shows flow curves σ xy ( t → ∞ ) = σ st xy vs bare Péclet number Pe for separation parameters ε =( φ − φ c ) /φ c as given in the legend. The colored symbols denote numerically obtained values for six shear rates Pe = 10 {− ... ; − } which are connected with straight lines as guides to the eye. For ε = 0 + , higher Pe were calculated up to the failure of thenumerics at Pe ≥ ; the Pe = 10 point was computed without friction kernel in Eq. (2) (i.e. just Taylor dispersion).The open black symbols with effective separation parameters ε eff given in the legend are experimental data obtained byCrassous et al. (2008); the separation parameters ε eff were deduced from loss G ′′ ( ω ) and elastic G ′ ( ω ) spectra. The measuredeffective packing fractions were φ eff = 0 . , . , . and 0.622. A gray dash-dotted line illustrates a slope of 1 and thus theNewtonian regime. The high shear viscosity η ˙ γ ∞ indicates the hydrodynamic contribution ( η ˙ γ ∞ = 1 . k B T / ( D d ) as measuredby Crassous et al. (2008)) neglected in the MCT-ITT calculations. The comparison is discussed in Sect. VI. The inset showsthe stationary viscosity η xy = σ st xy / ˙ γ . and more generally the Newtonian viscosity is defined in the limit Wi = ˙ γτ → . A first Newtonian plateau can beidentified for fluid states in Fig. 4 inset. A second Newtonian plateau would arise if the initial decay Γ q ( t ) in Eq. (2)gave a rapid shear-rate independent decay of the transient correlators at high shearing. This would also make g xy ( t ) independent from ˙ γ for high shear rates, viz. Pe . The present MCT-ITT cannot address this, as it describesthe physics of structural arrest at long times and uses for short-times the quiescent S q and Brownian motion with D as input without further considering how shear might affect them prior to structural arrest. Another problemarises, because the actual numerical iteration algorithm becomes unreliable for Pe > . The flow curve becomesnon-monotonic for stronger shear rates, and can even turn negative for some parameter values.How a stress overshoot emerges in microscopic MCT has been discussed already for an simplified model with isotropicshear-distortions by Zausch et al. (2008) and for schematic MCT by Amann et al. (2013). It provides insights into thephysical mechanisms involved in the yielding of glass. These will be discussed in more detail in context with Fig. 12comparing effects from repulsion and attraction. Here, the pertinent results for hard spheres shall be summarized.The peak in the transient stress indicates a characteristic strain value γ ∗ for the yielding process of a glass because itonly arises for Wi ≫ . Another quantity, the relative peak amplitude, σ pkxy /σ stxy − , characterizes the stress built-upduring the linear response of the glass, which is released during the later stage of the yielding process.Numeric evaluations for the present hard sphere system, lead to the dependence of the peak strain γ ∗ , viz. the zeroof g xy ( t, ˙ γ ) , on shear rate and packing fraction, which is shown in Fig. 5. One verifies that in the glass ( ε > ) and forsmall Pe , the critical yield strain γ ∗ is independent of shear rate. The reason is that intrinsic time scales then playno role in Eq. (5), because the transient density correlators decay with accumulated strain, as does the vertex in anycase. If the transient density correlators, which encode the structural relaxation, vary with time (and wavevector)even for negligible accumulated strain, then the value of γ ∗ changes due to the dependence of the integral on thewhole k -range. When the linear response regime of the fluid dispersion is approached, the stress overshoot vanishes,so that the zero of g xy ( t ) still can be used to define γ ∗ , but it loses its role as position of a noticeable overshoot in thestress-strain curve. The inset in Fig. 5 shows the vanishing of the relative overshoot height for negative separationparameters, ε < . Only if the the Weissenberg number Wi = ˙ γτ exceeds unity, a well developed stress overshootallows to observe γ ∗ in the stress-strain curve directly. For the fluid state at ε = − − , fitting a Kohlrausch lawto the quiescent shear modulus g xy ( t, gives τ = 8063 d /D (and stretching exponent β = 0 . ), so that Wi = 1 holds at bare Péclet number Pe = 1 . · − . This corresponds well to the inset in Fig. 5, where one estimates arelative amplitude 0.238 (lin. interpolated) of the stress overshoot at this Péclet number. Decreasing Pe further, theovershoot vanishes quickly. Except for this rapid variation when the quiescent fluid is approached, the position γ ∗ ofthe stress overshoot increases smoothly with packing fraction for all shear rates. This can be seen from the main panelof Fig. 5, and arises from the slight deepening of the undershoot in the shear modulus g xy ( t, ˙ γ ) , which can be noticedin the inset of Fig. 3. The γ ∗ also is an increasing function of shear rate for all separations to the glass transition. -6 -5 -4 -3 -2 -1 e = 10 -2 e = 10 -3 e = 0 e = -10 -3 e = -10 -2 = 10 -1 Pe = 10 -2 Pe = 10 -3 Pe = 10 -4 Pe = 10 -5 Pe = 10 -6 PSfrag replacements Pe γ ∗ ε · σ p k x y / σ s t x y − Figure 5. The main panel shows peak-strain values γ ∗ as function of bare Péclet number Pe . The grouping in separationparameters ε is given by the outside legend. The γ ∗ s are read off from the zeros of g xy ( t, ˙ γ ) , Eq. (6). Symbols are connectedwith straight lines as guide for the eye. The inset shows the relative overshoot height σ pk xy /σ st xy − of the MCT stress overshootsas function of separation parameters ε , grouped in bare Péclet numbers Pe ; see inset legend. Both, peak-positions andpeak-heights, increase with shear rate and with packing fraction. The overshoot’s shape and height remains invariant, as long as internal relaxations play no role and parametersremain close to the glass transition, where τ is large. The inset of Fig. 5 shows that the relative overshoot height σ pkxy /σ stxy − , with σ pkxy the maximal stress value, increases with packing fraction and with Pe outside the asymptoticregime. The former effect is due to the growth of the vertices and correlators in Eq. (6) with packing fraction, thelatter due to the growing role of the β -decay with increasing shear rate. Note that the increase with density differsfrom the findings by Koumakis et al. (2012a), where the opposite was measured. A decrease of the relative overshootheight with packing fraction was attributed to approaching random close packing (RCP). Besides speculations byAmann et al. (2013) if ageing played a role, this difference remains unresolved. It may indicate that MCT-ITT doesnot describe the regime close to RCP, but around the glass transition, and that MCT’s predictions fail at much higheror lower densities than at φ c .To illustrate the tensorial nature of the stress tensor σ ( t ) and of the 3d MCT-ITT approach, Fig. 6 shows acalculation of the first and second normal-stress differences N ( t ) = σ xx ( t ) − σ yy ( t ) and N ( t ) = σ yy ( t ) − σ zz ( t ) .Already Brader et al. (2009) showed tensorial, numerical evaluations of σ ( t ) via a schematic MCT. However, inschematic MCT, N = 0 under shear by construction. More recently, Farage et al. (2013) calculated for small shearrates the leading quadratic order of both normal stresses for stationary flows in fluid states. Figure 6 shows resultsfor the nonlinear regime of yielding glasses complementing their study.Within MCT-ITT under simple shear, all components of the stress tensor can be calculated analogously toEq. 4, by simply changing the components of some k vectors. Because of the neglect of an isotropic contribu-tion (Fuchs and Cates, 2009), however, the shear-dependent pressure can not be calculated presently. Figure 6 showsthe stress on a plane-element perpendicular to the flow direction, viz. σ xx ( t ) , and the normal-stress differences asdefined above, which are of distinguished rheological interest. The calculation has been done at the glass transition( ε = 0 + ) and for Pe = 10 − − − , i.e. for a genuine, critical glass behavior with small Pe and large Wi. Thefirst observation is that they all exhibit a transient regime, with stress overshoots, which look qualitatively like thoseof the shear stress, cf. Fig. 2. The overshoots however occur at strains larger than 0.4, which is larger than all γ ∗ determined from the shear stress. At all times t > , it holds σ xx > σ zz > σ yy > , which renders N > and N < . Comparing more quantitatively, Brader et al. (2009) showed that the Lodge-Meissner relationship holds inthe elastic regime. It reads N ( t ) /σ xy ( t ) = γ , and states that the slope of the 1st normal-stress difference is quadratic0 -3 -2 -1 s xx : Pe = 10 -2 -4 -6 N = s xx - s yy : Pe = 10 -6 - N = s zz - s yy : Pe = 10 -6 PSfrag replacements ˙ γt σ xx , N , − N (cid:2) k B T / d (cid:3) -2 -1 N / s xy : Pe = 10 -2 -4 -6 - N / s xy : Pe = 10 -2 g ˙ t g ˙ t PSfrag replacements ˙ γt ± N / σ x y Figure 6.
Left panel: transient normal stresses σ xx ( solid lines ), first normal-stress difference N = σ xx − σ yy ( dot-dashed lines ),and second normal-stress difference N = σ yy − σ zz ( dashed lines , plotted with negative sign) as functions of accumulated strain ˙ γt at the critical packing fraction ( ε = 0 + ) and for the Pe s given in the legend (Pe s are color coded for all line-styles ). Oneverifies that all σ ii > and σ zz > σ yy . Stress overshoots are apparent; they lie at 25% higher peak strains than in the shearstress.The right panel shows the normal-stress differences divided by shear stress, N /σ xy and − N /σ xy , to illustrate the propor-tionality of this ratio to the accumulated strain in the elastic region. Grey dashed and gray dotted lines are linear fits withprefactors 1 and 0.28 (see legend). in accumulated strain in the elastic regime with prefactor G ∞ , the shear modulus. The full 3d numerics recovers thisrelation and gives a corresponding one for the 2nd normal-stress difference: N ( t ) /σ xy ( t ) ≈ − . γ ; both relationsare shown in the right panel of Fig. 6. Outside the elastic regime, the simple linear relation between normal stressesand shear stress breaks down, as is clear from the different stress overshoot positions. V. ATTRACTION DRIVEN GLASS
Turning on a square-well potential (SWP) among the hard spheres, the attraction depth U relative to the thermalenergy becomes an additional control variable. Increasing U the states diagram depends on the relative attractionrange δ . We fix δ = 0 . , far smaller than Lindemann’s ratio, in order to explore the consequences of physical bond-formation in concentrated dispersions. At this δ , the two different types of glass transitions, repulsion dominated glass(RDG) and attraction dominated glass (ADG), merge in an higher order singularity of MCT, denoted A bifurcationpoint (Götze and Sperl, 2002). This occurs around U A ≈ . ; see Fig. 1. For U below this value, the major effectof the attraction is to destabilize the RDG and to induce a (reentrant) fluid phase where equilibrium is reachedat long times (Pham et al. (2002) and Eckert and Bartsch (2002)). For attractions close to and slightly above U A (and concentrations below it), the ADG transition takes place where the short-ranged attraction causes cooperativebonding among caging particles. Even though the attraction strength is of the order of a few k B T only, the ’physicalbonding’ stabilizes a second non-ergodic state with quite different mechanical properties than the RDG.Figure 7 shows the (quiescent) shear moduli G c ∞ , viz. the elastic coefficients under volume conserving deformations,along the glass transition lines as functions of attraction strength U ; note that φ varies non-monotonically along thered curve in Fig. 1 for these calculations. For U = 0 the modulus of HS indicated in Figs. 2 and 3 is recovered. Atfirst, the attraction little changes the elastic shear modulus relative to the RDG value until around U A , G c ∞ increasesrapidly. For even shorter attraction range, e.g. δ = 0 . , the modulus would jump discontinuously at the crossing ofglass transition lines in Fig. 1 (Dawson et al. , 2000). The mechanism causing the increase is the tighter localization ofparticles in the ADG than in the RDG. It remains an entropic effect like in the RDG, as G c ∞ ∝ k B T , but asymptoticallyfor small ranges MCT predicts an increase scaling like: G c , ADG ∞ ∝ (1 /δ ) G c , RDG ∞ . Note that for increasing U along theglass transition line (outside the range of Fig. 7) the elastic modulus decreases non-monotonically, which was shownby Bergenholtz et al. (2003). This happens because the particle concentration strongly decreases along the transitionline (and with it the elastic modulus) for increasing interaction strength.Turning on shear destroys the elasticity of the amorphous solids and causes plastic deformations also in the presenceof attractions. We choose the packing fraction φ c given by the RDG transition of hard spheres (HS) in order to study1 PSfrag replacements U G c ∞ ( φ , U ) (cid:2) k B T / d (cid:3) U G ∞ ( φ c , U ) (cid:2) k B T / d (cid:3) Figure 7. Shear moduli of hard spheres with square-well attraction of relative range δ = 0 . as function of attractionstrength U = u /k B T . Squares mark calculations along the glass transition line and correspond to the black squares in Fig. 1;the packing fraction φ c varies non-monotonically here. The inset shows shear moduli varying over a wider range. Symbols × mark G ∞ calculated at fixed packing fraction φ c and increasing U deeper into the attraction driven glass. the interplay of attractions and shear. Figure 8 gives an overview of the stress-strain curves for various U . TheHS curve is included and for the solid states the elastic law with the independently calculated linear shear moduli G ∞ ( φ c , U ) from Fig. 7. Intermediate shear rates are chosen, in order to see strong effects but to remain well belowthe limit of applicability of MCT-ITT; recall that for HS, instabilities emerged for Pe > . Starting the discussionof panel a) in Fig. 8 in the elastic regime and at U = 0 , the linear region in the HS stress-strain curve lies somewhatabove the linear elastic asymptote G c , RDG ∞ γ because the bare Péclet number is not asymptotically small; compareFig. 2. For attraction strengths in the reentrant fluid region, the stress-strain relations clearly exhibit fluid behavior.At U = U c / the glassy structure on intermediate time scales relaxes so fast that the stress-strain curve becomesmonotonous as holds in linear viscoelastic response, where g xy ( t ) is independent of the shear rate. We estimateWi = ˙ γτ ≈ . there, while it is infinite for U = 0 and U = U c . For the critical attraction strength U c of the ADG,the stress-strain curve first exhibits a linear elastic regime, then an overshoot and finally a steady state. The linearregime is again well described by the asymptote G c , ADG ∞ γ , yet, the linear shear modulus at the ADG transition exceedsthe HS one by roughly a factor . . The linear elastic response holds over a far smaller strain range than for HS,and the ensuing stress peak is far broader at the ADG than at the RDG transition. The physical bonds apparentlystart to get broken already at strains comparable to the attraction width so that plastic rearrangements occur andthe linear elastic limit is left early. The position γ ∗ of the stress peak, however, is unexpectedly somewhat larger atthe ADG than at the RDG transition. Yet bonds — possibly by rotating and stretching — still manage to bear stressup to strains even somewhat larger than characteristic for the RDG. Thus the stress overshoot is broad. Deep in theADG, at twice the attraction strength than at the transition, the stress strain curve exhibits the same regimes. Yet,the elastic modulus is larger and the position of the stress overshoot has shifted to noticeably smaller values. Also, therelative height of the overshoot has increased very strongly; see panel b) in Fig. 8 and Fig. 11 below. Consequently,the steady state stress σ st xy ( γ → ∞ ) is higher because of the steeper linear elastic increase, but lower than this effectwould imply because of the large stress release after the maximal stress σ pk xy = σ xy ( γ ∗ ) . This is brought out clearlyin panel d) of Fig. 8, where the stress divided by the shear modulus is shown for a representative set of attractionstrengths. All curves overlap in the linear regime by construction. States affected by the short ranged attractionleave the linear regime at very small strains, while the hard sphere curve follows the linear elastic response far longer.The large fraction of stress released during yielding in states deep in the ADG is apparent, and the broadness of thestress overshoots with attractions. The transient stress curves at the high U also exhibit a broad maximum, whichapproaches the long-time limit only around strains of order unity. On top of this slow transient, the calculations in2 = 10 -3 e U = 0-0.5-0.867-1 PSfrag replacements a) ˙ γt σ x y (cid:2) k B T / d (cid:3) ADG: G c ∞ HS: G c ∞ e U = 1 e U = 0.4 e U = 0.2 PSfrag replacements b) ˙ γt σ x y (cid:2) k B T / d (cid:3) G ∞ = 10 -1 -2 -3 -4 -5 -6 PSfrag replacements c) ˙ γt σ x y (cid:2) k B T / d (cid:3) G c ∞ Pe = 10 -6 : e U = -1 e U = 0 e U = 0.4 e U = 1 PSfrag replacements d) ˙ γt σ x y / G ∞ G ∞ Figure 8. Overview of the transient shear stress σ xy vs accumulated strain ˙ γt curves at fixed density when varying the attractionstrength given by the relative separation to the critical attraction strength, ε U = ( U − U c ) /U c ; the HS φ c is chosen. In panela), ε U = − , − . , − . , and 0 correspond to the HS case, two values in the reentrant fluid region, and at the attraction-dominated transition; the shear-rate is fixed at Pe = 10 − . In panel b), ε U = 0 . , . , and 1 (line style coded) enter deeperinto the ADG. Here, calculations for Pe between − and − overlap (master functions shown). Note, the oscillations uponreaching the steady state at the strongest U are due to the limitations of the rough q -grid discretization. Panel c) varies theshear-rate as given in the legend at the ADG-transition, ε U = 0 + . Panel d) shows transient stresses divided by the correspondingshear moduli, σ xy ( t ) /G ∞ for some glass states (labeled by attraction strength separation ε U ). All curves overlap in the linearregime by construction as the shear rate is low, Pe = 10 − . In all panels, the elastic moduli G ∞ and G c ∞ , calculated from theplateaus of Figs. 3 and 9, or from Fig. 7, are shown as gray, dashed, dotted, or dash-dotted lines . Fig. 8 exhibit noticeable oscillations which we consider artifacts of the coarse discretization of Eq. (6).After this overview of the yielding and plastic deformation when turning on a short ranged attraction, the differenteffects shall be explored in more detail. At first the region close to the critical attraction strength U c shall be exploredvarying ε U slightly, and second the behavior for large U is investigated. Finally, the steady state flow curves arepresented. Recall that ε U = ( U − U c ) /U c , i.e. ε U = − in the HS case, ε U = 0 at the ADG transition line, and ε U = O (1) in the ADG glass region.Panel b) in Fig. 8 and Fig. 9 present a detailed look at the transient stress evolution close to the ADG transitionwhere shear rate, non-dimensionalized as bare Péclet number, and relative separation ε U are varied. For attractionstrengths close to the ADG, the linear elastic regime generally is followed by a stress maximum. Here, the linearelasticity only holds for rather small strains. Starting at strains below 5%, the stress grows sub-linearly with strain.As the final approach to the steady state asymptote requires strains of order unity, in general a rather broad stressovershoot can be seen in the ADG. Comparing the relative magnitude of the stress overshoot and its strain-position3 -3 -2 -1 G ∞ c G ∞ ( ε U = -2 ) c Pe = 10 -1 -2 -3 -4 -5 -6 ε U = -10 -2 -2 t (cid:2) d /D (cid:3) g x y ( t , ˙ γ ) (cid:2) k B T / d (cid:3) ˙ γt g x y ( t , ˙ γ ) (cid:2) k B T / d (cid:3) Figure 9. Generalized shear modulus g xy ( t, ˙ γ ) of hard spheres with square-well attraction close to the ADG transition, asfunction of time (main panel) and accumulated strain ( inset ; curves for every second Pe left out for clarity). The legendprovides color coded the strain rates and line-style coded the relative attractions ε U = ( U − U c ) /U c ; all curves at the HS φ c . Aletter c labels the critical g xy ( t, for ε U = 0 + . Elastic shear moduli G ∞ ( ε U ) can be read off from quiescent curves (Pe = 0 ),with G ∞ ( ε = 0 + ) = G c ∞ = 89 . k B T /d and G ∞ ( ε = 10 − ) = 122 k B T /d . The inset shows subtle differences in the peakposition γ ∗ , where g xy ( t, ˙ γ ) = 0 , which are caused by ˙ γ independent β and α decays; glass curves for Pe = 10 − and − overlap. For Pe = 10 − in the fluid, the undershoot (almost) disappears. γ ∗ , similar results are observed as in the HS case.To support this comparison, Fig. 10 shows stress-peak strains γ ∗ and relative overshoot magnitudes σ pk xy /σ st xy − for different ε U and Pe around the ADG and RDG transition. They, correspond to the zeros and negative areasin Fig. 9, respectively. A noteworthy and at first unexpected effect is the larger characteristic strain value at theADG than at the RDG. For shear rates giving Pe ≤ . , the γ ∗ ’s at the transition are somewhat larger for a glasswhere particles feel strong bonds to their neighbors than for a glass where repulsive interactions dominate. Closeto the ADG, the characteristic strain γ ∗ does not correlate with the localization length, which is far smaller at theADG than at the RDG as shown in the inset of Fig. 1. The cause of this subtle effect, which also depends on shearrate, is given by the extreme stretching of the quiescent α -process at this ADG transition. The physical cause wasanalyzed in detail by Götze and Sperl (2002): The scenario of two different glass transition lines causes very broadrelaxation curves, culminating in logarithmic relaxation close to higher order glass transitions. At the present choiceof density, attraction range and strength, the quiescent shear modulus happens to exhibit logarithmic decay for morethan five decades; see Fig. 9. This causes a shear-distortion of the shear modulus around its zero which pushes γ ∗ to larger values. The inset of Fig. 9 shows the shallow negative region caused by the broad α -process for small shearrates. Thus γ ∗ increases for low Pe . For larger shear rates the relaxation curves become steeper with larger U ,and thus γ ∗ decreases entering the ADG. The extremely slow intrinsic relaxation also explains the broadness of thestress-overshoot because shear affects the internal relaxation in a wide time window. This can also be deduced fromthe inset of Fig. 10, which shows that the stress-overshoot exists for a broad window of shear rates at the ADG. Againthis holds because the intrinsic α -process is characterized by a broad distribution of relaxation times. The broadercrossover between a characteristically harder elastic regime and the steady state, reached after a larger plastic stressrelease, thus is the qualitative difference of the stress-strain curves at an attraction driven compared to a repulsiondriven glass transition.Figure 11 continues the investigation of Fig. 10, i.e. peak strain γ ∗ and relative overshoot magnitudes as functionof ε U and Pe , but deep inside the ADG phase. There, the transient stress evolution changes strongly and in acharacteristic way. First, the elastic constants G ∞ increase dramatically. The inset of Fig. 7 shows that at anattraction strength twice as large as the critical value of the ADG transition, the shear constant has increased bya factor around one hundred relative to the HS one. This can also be seen directly from the linear elastic regimein the stress-strain curves in Fig. 8 b). Non-linearities set in at strain values comparable to the attraction range,as was observed at the ADG transition already. Thus, deep inside the ADG phase where the overshoot peak isdominated completely by the elastic energy stored in the short range particle bonds (which is discussed in more detailtogether with Fig. 12 below), a bond reordering at smaller strain values shifts the peak position to much smaller γ ∗ for4 -6 -5 -4 -3 -2 -1 e U = 10 -2 -2 -0.5-0.818-0.867-100.10.20.30.40.50.6 e U -1 -0.8 -0.6 -0.2 0 Pe = 10 -1 Pe = 10 -2 Pe = 10 -3 Pe = 10 -4 Pe = 10 -5 Pe = 10 -6 PSfrag replacements Pe γ ∗ ε θ σ p k x y / σ s t x y − Figure 10. The main panel shows peak-strain values γ ∗ as function of bare Péclet number Pe going from the hard sphere glasstransition ( ε U = − ) through the reentrant fluid state( − < ε U < ) to the attraction dominated glass transition ( ε U = 0 ); thepacking fraction is fixed at φ c . The upper outside legend gives the relative separations in attraction strength ε U = ( U − U c ) /U c .The γ ∗ s are read off from the zeros of g xy ( t, ˙ γ ) , Eq. (6). Symbols are connected with straight lines. The inset shows the relativeovershoot height σ pk xy /σ st xy − of the stress overshoots as function of the relative separation in attraction strength ε U and forseveral Pe (lower outside legend). The variations mainly result from the shift of the broad α -process through the range set bythe shear rate. with inverse temperature (similar to the HS case). increasing U . The stress-overshoot retains its width in strain values, as the steady flow curve values is approached ataccumulated strains of order unity (as we find for all considered states; recall that the wiggles in Fig. 8 at high U arenumerical artifacts). This implies that the steady state stress then arises from a competition of shearing and repulsiondominated caging which is not very strongly affected by attractions. In consequence, a very prominent feature ofthe stress-overshoot deep in the attraction driven glass is its relative magnitude. The lower panel of Fig. 11 showsthat the the relative magnitude σ pk xy /σ st xy − increases by around a factor six for the present path into the ADG. Alarge portion of the stress stored elastically for intermediate times is released during the late stage of the transient.The flow-curve deep in the ADG is raised by roughly a decade for the present parameters relative to the HS flowcurve. This increase is smaller by a decade than the attraction-driven increase of the elastic constant G ∞ . Becausethe internal β -relaxation at the states deep in the glass domain is quite rapid, no shear-rate dependence is observedin the stress-strain curves. This holds because the generalized shear modulus, which in general depends on time andstrain independently, has become a function of accumulated strain only, g xy ( t, ˙ γ ) = g xy ( t ˙ γ ) . Thus the stress-strainrelations for different shear rates collapse onto a common master curve, which, for different U , are shown in Fig. 8panel b).The structural rearrangements involved in the transient stress release around γ ∗ can be recognized from the wavevec-tor dependent contributions in the generalized shear modulus. Figure 12 shows the wavevector dependent vertex,viz. the square bracket in Eq. (6), for some relevant attraction strengths U and densities. The strain values are (closeto) the position γ ∗ of the stress overshoot. Two specific directions are chosen along the extensional ( ϕ k = 45 o ) andcompressional ( ϕ k = 135 o ) axis, where the distortions of the structure are maximal under shear (Koumakis et al. ,2012a). The correlators Φ q ( t ) in Eq. (2) depend on time and accumulated strain independently. The former due tothe intrinsic relaxation, the latter due to wavevector advection. In Eq. (6), the squared correlators are a weight forthe purely strain/advection dependent vertex, and a stress overshoot emerges depending on the different weighting ofwavevector contributions. If the product of structure factor deviations S ′ k S ′ k ( t ) becomes negative for enough k values,the k -space integral and thus the shear modulus become negative. At the HS transition we recover the result whichZausch et al. (2008) obtained by isotropic averaging that negative contributions arise close to the principal peak of thestatic structure factor; for HS it lies close to kd = 7 . The dominance of the principal peak in S ( k ) at the RDG verifiesthat the RDG originates in the local steric packing always present in dense fluids. The origin of the stress-overshootbeing negative contributions from the local order peak in S ( k ) , also holds deep in the RDG, viz. for HS at ε = 0 . ,where the vertex has grown with density and varies more rapidly with wavevector, and where the correspondingdensity correlators are more glass like, viz. have higher plateau amplitudes. Figure 12 includes the vertex for ε = 0 . to exemplify this; other results at this ε are not shown as they can be extrapolated based on the data presented in5 e U Pe = 10 -1 Pe = 10 -2 Pe = 10 -6 e U Pe = 10 -1 Pe = 10 -2 Pe = 10 -3 Pe = 10 -4 Pe = 10 -5 Pe = 10 -6 PSfrag replacements γ ∗ σ p k x y / σ s t x y − Figure 11. Peak-strain values γ ∗ (upper panel) as function of the relative separation in attraction strength ε U = ( U − U c ) /U c spanning from the hard sphere glass transition, ε U = − , to the attraction dominated glass transition, ε U = 0 , until deep inthe ADG, ε U = 1 . The lower panel shows the relative overshoot height σ pk xy /σ st xy − of the stress overshoots as function of ε U .Data for several shear rates as shown as labeled with Pe ; they overlap in the explored range. Sect. IV: E.g. the characteristic strain has changed little relative to the HS transition at ε = 0 , and takes the value γ ∗ ( ε = 0 .
1) = 0 . at Pe = 10 − . Somewhat unexpectedly, the same wavevector range as at the RDG transitiondominates the stress integral at the stress maximum of the ADG transition. This is shown in the panel at ε U = 0 of Fig. 12. One notices that larger wavevector contributions have grown, but the dominant contributions remainsclose to the peak in S q . The large wavevector contributions in MCT-stress kernels capture the formation of physicalbonds which result from the increased stickiness of the attractive square-well potential. It increases the equilibriumstructure factor at large k . These high- k modes are responsible for the early breakdown of the linear elastic regime inADG states; see Fig. 8. Nevertheless, the dominant contribution at the ADG transition remains connected with thelocal ordering of neighbor shells. Apparently, this is a universal characteristic of the yielding process of a glass at thetransition in MCT-ITT. This holds even though the elastic modulus is characteristically larger at the ADG than atthe RDG as shown in Fig. 7. The situation changes deep in the ADG at U = 2 U c (i.e. ε U = 1 ), where Fig. 12 indicatesthat negative contributions in the stress relaxation arise dominantly at large wavevectors. Then the characteristicstrain γ ∗ becomes smaller, see Fig. 11, and the fraction of released stress grows strongly. This holds because thecontributions at large wavevectors, viz. local rearrangements of the physical bonds, are rather rapid. When bond-formation and breakage dominate the stresses deep in the ADG, the advected wavevector changes quickly with time,and the vertices of MCT, which are positive only in the quiescent state, become negative rapidly. The presence ofcontributions from the main peak in S q apparently cause that the bonded glass can rearrange (quasi-) elastically untilthe packing-dominated cages yield. Then the stress is released, which was stored elastically in the physical bonds.As the bonds become distorted starting from very small strains, a characteristically broad stress-overshoot results.The vertex at U = 2 U c in Fig. 12 is not decayed to zero at the upper cut-off in k of our integration. This indicatesthat the results at ε U = 1 are not completely converged. Yet, we expect only quantitative corrections because of the6evolution of results from smaller ε U . -20-15-10-5 0 5 10 15 0 10 20 30 40 50 60 70ADG: ε U = 0 ↑ ϕ k = 45 °ϕ k = 135 ° ε = 0 ↑ ↓ -49.3 ϕ k = 45 °ϕ k = 135 ° HS: ε = 0.1 ↑ ↓ -199 ϕ k = 45 °ϕ k = 135 ° ε U = 1 ϕ k = 45 °ϕ k = 135 ° k [1 /d ] v x y ( ϑ k = ◦ , g x y = ) [ k B T ] Figure 12. Wavevector k dependent vertex of the generalized shear modulus in Eq. (6); viz. the contents of the square bracketthere. The directions of extensional ( ϕ k = 45 o , red lines) and compressional ( ϕ k = 135 o , green lines) strain are chosen and theaccumulated strain is (close) to γ ∗ . Four different glasses are considered as labeled and discussed in the text. Figure 13 shows the flow curves for a few representative states spanning from the hard sphere to the attractiondriven glass transition and beyond, deep in the ADG. The packing fraction is kept fixed at the critical value of theHS transition, so that the ε U = − curve corresponds to the critical HS curve at ε = 0 of Fig. 4, where flow curves ofhard spheres were shown. The reentrant fluid region lies at small U , which correspond to negative separations ε U < .The states from ε U = − . to ε U = − . exhibit Newtonian viscosities which decrease with increasing U . Raising U up to close to the value of the ADG transition, the Newtonian viscosity increases again strongly. At e.g. ε U = − . the Newtonian regime lies outside the window of Fig. 13 at lower bare Péclet numbers. Thus, crossing the reentrantliquid region, the Newtonian viscosity varies non-monotonically as observed experimentally by Willenbacher et al. (2011). Entering the ADG, a yield stress σ + xy arises as holds universally at MCT-ITT glass transitions. Comparingthe values of σ + xy at the HS and at the ADG transition at the same packing fraction φ c , one notices an increasecaused by the attractions. The yield stress increases by roughly the same factor as does the shear modulus G ∞ atboth transitions; compare Fig. 7. Entering into the ADG, the steady stresses increase yet again. Because the localdynamics of caging and bonding has become quite fast according to MCT-ITT deep in the ADG, the flow-curvebecomes shear-rate independent. At ε U > . the numerical curves indicate no dependence on Pe in the windowof Fig. 13. It is noteworthy that the increase of the steady stress in the bonded glass relative to the hard sphereglass is far smaller than the increase of the corresponding elastic constant. The inset of Fig. 7 indicates that G ∞ hashardened by around two orders when going from U = 0 to U ≈ , while the yield stress increases only from around σ + xy ( U = 0) ≈ k B T /d to σ + xy ( U = 2 U c ) ≈ . k B T /d . The reason behind the comparatively weak increase in thesteady stress lies in the different stress recovery after imposing flow in the ADG and the RDG. A large fraction of thestress which is elastically stored in the physical bonds is released when the ADG fluidizes for strains of the order of γ ∗ ; see the overview of stress-strain curves in Fig. 8 and the detailed analysis of the magnitude of the stress-overshootin Fig. 10.7Converting the flow curves into viscosities, η xy = σ xy / ˙ γ , straightens out the curves over a wide range along theordinate. The subtle sigmoidal shape of the flow curves of MCT-ITT close to a glass transition thus get ironed-out. The inset of Fig. 13 shows the corresponding viscosities which exhibit a Newtonian plateau in fluid states forsmall Weissenberg numbers and then cross over to shear-thinning with asymptotic exponent -1. Restrictions in thenumerical code prevent us also for the ADG to address the question of the existence of a second Newtonian plateauat high shear rates. The MCT-ITT calculations continue to decrease for the numerically accessible range of ˙ γ . -3 -2 -1 -6 -5 -4 -3 -2 -1 ε U = 110 -2 -2 -0.5-0.818-0.867-110 -5 -3 -1 ADG: φ = 0.5998HS: φ = 0.5949 Pe σ s t x y (cid:2) k B T / d (cid:3) Pe η x y [ k B T / ( D d ) ] Figure 13. The main panel shows flow curves σ xy ( t → ∞ ) = σ st xy vs bare Péclet number Pe for relative attractions ε U =( U − U c ) /U c as given in the legend. The curves span from the hard sphere glass transition, ε U = − , to the attractiondominated glass transition, ε U = 0 , until deep in the ADG, ε U = 1 . The colored Symbols were calculated for six shear ratesPe = 10 {− ... ; − } and are connected with straight lines as guides for the eye. The open black symbols with packing fractions φ given in the legend are experimental data obtained by Pham et al. (2008). The comparison is discussed in Sect. VI. The inset shows the stationary viscosity η xy = σ st xy / ˙ γ . Figure 14 shows the first and second normal-stress differences N = σ xx − σ yy , N = σ yy − σ zz . Choices in thenumerical algorithm aimed at capturing strong flows give errors in the coefficients when they become too small influid states. Hence we cannot compare with the calculations by Farage et al. (2013) who considered the prefactorof the quadratic scaling at low Péclet numbers. As in the case of the shear stress, the steady state normal stressesdeep in the ADG are above the ones of the HS at the same packing fraction. Yet, during the transient, again a largeamount of the stress built-up during the linear elastic response is released upon yielding. The build-up of normalstresses during the deformation of physical bonds in the linear regime again obeys the Lodge-Meissner relationship astested in Fig. 6 for hard spheres; the tests for the ADG states are not shown. VI. CONCLUSIONS AND COMPARISON WITH EXPERIMENTAL DATA
We presented the first fully quantitative solutions of the MCT-ITT equations for the nonlinear response of sheardriven colloidal dispersions. The results on the transient shear stress response, the steady flow curves, and normalstresses at colloidal glass transitions exhibit the central qualitative features which have previously been discussedusing schematic MCT-ITT models. They recover and quantify phenomena like the initial elastic response, which islinear in the accumulated strain, the yielding of glasses characterized by a dynamic yield stress, and the transientstress overshoot, which defines a strain γ ∗ characterizing the yielding process. Solutions of the MCT-ITT equations byHenrich et al. (2009), Krüger et al. (2011), and Amann et al. (2013) considering hard disks in two dimensions showedthat these phenomena are universal at MCT-ITT glass transitions under shear and also arise in two-dimensionalsystems confined to a plane. These authors also discussed transient density correlators, tagged particle motion, anddistorted structures, which, for three dimensions, can only be presented in future. The present three-dimensionalsolutions, however, enable comparison with experimental data, which are available for dispersions of colloidal hardspheres and of colloids mixed with non-adsorbing polymers which induce a short-ranged attraction of Asakura-Oosawaform.8 -7 -6 -5 -4 -3 -2 -1 -6 -5 -4 -3 -2 -1 e U = 110 -2 -2 -0.5-0.818-0.867-110 -5 -3 -1 PSfrag replacements Pe N (cid:2) k B T / d (cid:3) Pe Ψ (cid:2) k B T d / D (cid:3) -6 -5 -4 -3 -2 -1 -2 -1 e U = 110 -2 -2 -0.5-0.818-0.867-110 -5 -3 -1 PSfrag replacements Pe − N (cid:2) k B T / d (cid:3) Pe − Ψ (cid:2) k B T d / D (cid:3) Figure 14. The left main panel shows stationary first normal-stress differences N ( t → ∞ ) = σ st xx − σ st yy vs bare Pécletnumber Pe for six shear rates Pe = 10 {− ... ; − } and for relative attractions ε U = ( U − U c ) /U c as given in the legend. The Symbols are connected with straight lines as guides to the eye. The left inset shows the first long-time normal-stress coefficient Ψ = ( σ st xx − σ st yy ) / Pe . The right main panel shows stationary second normal-stress differences − N ( t → ∞ ) = σ st zz − σ st yy , andits inset the corresponding long-time normal-stress coefficient − Ψ = ( σ st zz − σ st yy ) / Pe . A quite stringent comparison of the theoretical results for hard spheres in Sect. IV can be made with the experimentsby Crassous et al. (2008) who studied core-shell particles consisting of a polystyrene core and a thermosensitive,crosslinked PNIPAM shell. The microgel dispersions with size polydispersity of 9.3% can be mapped onto the phasediagram of monodisperse hard spheres using the freezing density φ F = 0 . , and exhibit the strongest tendency tocrystallize around φ eff = 0 . . The rheology at higher packing fractions closer to vitrification appears little affectedby crystallization. The effective packing fraction can be adjusted by changing the weight percentage of particles orby changing the effective size R H (viz. the hydrodynamic radius taken to be d/ here) by temperature, which makesit possible to approach the glass transition packing fraction of hard spheres φ c = 0 . quite closely. Moreover, thefrequency-dependent linear response moduli G ′ ( ω ) and G ′′ ( ω ) were already quantitatively analyzed using MCT so thatthe mapping of the theory onto the rheology is known. Crassous et al. (2008) found that roughly the same value of thecritical packing fraction as established by (Pusey and van Megen, 1987) studying PMMA hard sphere colloids usingdynamic light scattering rationalizes the microgel rheology. The separation parameters obtained by Crassous et al. (2008) and listed in the caption of Fig. 4 refer to this experimental hard sphere glass transition density. Clearly,the theoretical results, which are calculated (not fitted) for comparable separation parameters from the theoretical φ c show quite comparable flow curves. The need to use the separation from the critical packing fraction instead ofthe actual value of the packing fraction in order to compare MCT and experiment is well established (Götze, 2009).It arises from the approximate nature of MCT which misses the precise value of φ c , while capturing the sensitivedependence of the viscoelasticity on the separation to the glass transition. The comparison in Fig. 4 can in principlebe done without adjustable parameter, because D sets the time scale in Eq. (2) and can be calculated from the solventviscosity following Stokes, Einstein and Sutherland (Einstein, 1905; Sutherland, 1905). Also the stresses are calculateddirectly. Yet, hydrodynamic interactions are neglected by MCT-ITT. They affect the short time diffusion coefficient,which at high concentrations is a more relevant scale than the Stokes-Einstein-Sutherland diffusion coefficient atinfinite dilution. Additionally, hydrodynamic interactions add a contribution to the viscosity at high shear rates.Crassous et al. (2008) found that the linear moduli agree best when assuming that the hydrodynamic interactionsslow down the short time diffusion to . D . (The viscosity η ˙ γ ∞ observed at high shear rates is indicated in Fig. 4to describe the second origin of hydrodynamic deviations.) Also they observed that MCT underestimates stressesby 40% (a rescaling factor c Gy = 1 . was used). The data in Fig. 4 were rescaled by the given ratio of the diffusioncoefficients, but by a different stress-rescaling factor: c σy = 0 . . Quite satisfactorily MCT deviates by less than 50%from either experiment. The aspect that theory underestimates the linear elastic stress but overestimates the steadystate stress of the yielding glass can be traced back to the error of MCT-ITT in determining the characteristic strainvalue γ ∗ . While transient stress-strain curves are not available for the microgel dispersions by Crassous et al. (2008),stress-overshoots were measured in more polydisperse microgel samples. Amann et al. (2013) measured γ ex ∗ ≈ . while our MCT-ITT calculation gives γ mct ∗ ≈ . . Apparently, MCT-ITT underestimates the speed-up by shearof stress fluctuations. While experiments close to the glass transition in hard sphere colloids, including the largeamplitude oscillatory measurements by Petekidis et al. (2002), find characteristic strain values around 10%, MCT-ITT overestimates it by a factor around three. Consequently the steady stresses are somewhat overestimated, eventhough the linear elastic regime is somewhat underestimated. Reassuringly, the deviations by MCT-ITT in threedimensions are appreciably smaller than the deviations in two dimensional hard disk systems, where Henrich et al. c σy = 0 . . Considering the densitydependence of the stress-overshoot peak-strain γ ∗ , one notices a similar increase than found by Petekidis et al. (2002)for the strain where irreversible rearrangements first appear in states close to the glass transition. Intriguingly, theirlight scattering measurements of this characteristic strain reveal a maximum at intermediate densities followed bya decrease when approaching random close packing. The dependence of γ ∗ on packing fraction thus is richer thanthe dependence of the localization length on φ . The latter is expected to be monotonically decreasing with higherpacking fraction as particles get localized more tightly. While both length scales thus characterize the cage effectin repulsion dominated glass transitions, and take comparable values right at the hard sphere glass transition, theirprecise relation is non-linear and not straightforward.Turning on the short-ranged attractions by adding polymer to the colloids, first the states diagram can be tested.It consists of two different glass states, a reentrant fluid region, and possibly glass-to-glass transitions and higherorder glass singularities. Pham et al. (2002) and Eckert and Bartsch (2002) found qualitative agreement concerningthe transition lines, with the ones from theory shifted, but tracking the experimental ones. Willenbacher et al. (2011)extended these studies by pushing the reentrant fluid region to higher packing fractions. The enhanced elastic stiffnessof the glasses with physical bonds was convincingly seen by Pham et al. (2006) and Pham et al. (2008). While stressvs strain curves after shear start-up are not available, Pham et al. (2008) present and discuss as equivalent stressesafter step strain deformations. The observed characteristic strain γ ∗ ≈ . for hard spheres corresponds well to theabove discussion. For glasses with a polymer-induced attraction of roughly 6% range, the stress-strain relations showtwo characteristic differences to the ones of hard spheres. First, nonlinear deviations to the linear elastic response setin at rather small strain values. This, considering the differences in attraction potential, agrees well with our findingsfrom MCT-ITT. The origin of the nonlinearities lies in the high-wavevector contributions of the memory kernels whichcause a broadly stretched quiescent structural relaxation. They also are very susceptible to wavevector advection,the mechanism by which shear affects the structural relaxation in MCT-ITT. Thus small strains suffice to soften theelastic response. The second experimental finding is a second maximum of the stress at strain values beyond one,which is not observed by MCT-ITT. The stress maximum at large strains is higher than the maximum at strainscomparable to the repulsive case, distorting it to a shoulder. While MCT-ITT appears to capture the phenomena atthe first characteristic strain γ ∗ which remains close to the hard sphere value, the second stress maximum is missed.Presumably it arises from structural correlations in the bonded glassy state which reach beyond the cage-effect lengthscale. Going to the final steady state, Fig. 13 contains the experimental flow curves obtained by Pham et al. (2008)for a repulsion and an attraction dominated glass state. The experimental data are mapped onto the theoreticalcalculations by estimating the diffusion coefficient D = k B T / (3 πη solv d ) by using η solv = 1 mPa s with the particlesize d = c σy = 3 . gives the best agreement for the hard sphere data. Afterthis mapping of the hard sphere data, the increase of the yield stress deep in an attraction dominated glass statecan be addressed in Fig. 13. While the differences in packing fraction, attraction range and strength prevent a moredetailed comparison, the hardening of the flow curves at (roughly) fixed packing fraction upon increasing the attractionstrength agrees qualitatively with the MCT-ITT calculation. A more detailed comparison appears justified, whichwould require improved equilibrium structural input.Based on the encouraging comparisons of the nonlinear stress strain relations from microscopic MCT-ITT withexperiments on model colloidal systems it appears worthwhile to consider shear distorted structure and transientdensity correlations in order to gain a deeper understanding of the mechanisms of plastic deformation and yielding ofcolloidal glasses. Work along these lines is underway. VII. ACKNOWLEDGEMENTS
We thank Th. Voigtmann for valuable discussions, M. Siebenbürger for help in interpreting the experimental data,and the Deutsche Forschungsgemeinschaft (DFG) for financial support in the initiative FOR 1394, project P3.
Appendix A: Numerical implementation
In this appendix, the numerical implementation of MCT-ITT in three dimensions is summarized. The standardchallenge to solve MCT-equations over around ten decades in time is made more difficult by the requirement tocompute · d -dimensional wavevector integrals for the memory kernels. Desirable requirements to the numerics are torecover an isotropic quiescent solution (Franosch et al. , 1997), while choosing a sufficiently close q -grid discretizationin Fourier space, in order to minimize discretization errors.0The numerical evaluation of Eq. (2) depends mainly on the discretization of the friction kernel m q ( t, t ′ ) in Eq. (3).Under simple start-up shear, the dependence on two times t and t ′ < t can be simplified on the dependence on thetime interval τ = t − t ′ , yielding m q ( t ′ ) ( τ ) , which is explained in detail by Fuchs and Cates (2009). This yields a socalled single-time MCT, which is far simpler to solve than two-times MCT considered by Voigtmann et al. (2012). Thetemporal evaluation of Eq. (2) follows the standard scheme of previous MCT-ITT calculations (also schematic MCT)and is described in detail e.g. by Voigtmann et al. (2012) and Amann (2013). The temporal discretization is performedon a time grid consisting of blocks, which each consist of N t = 64 linearly spaced time instances. A straightforwarddiscretization backwards in time (see Brader-Voigtmann algorithm (Amann, 2013)) and Φ q ( t ) is iteratively solved,depending on all time instances t ′ t . The first time block is initialized by Φ q = exp( − Γ q t ) . As second input, thePercus Yevick (PY) structure factor for hard spheres and square well potentials is used (Hansen and McDonald, 2009;Baxter, 1970; Dawson et al. , 2000). Subsequent time blocks are generated by doubling the time step and using thearithmetic average of Φ q and m q ( t ) of two instances of the preceding time block to initialize the first N t / instancesof the new time block. Via this decimation procedure, a large time range can be covered.The discretization of wavevectors in q space (and the k integral in the friction kernel, Eq. (3)) has been performedusing standard spherical coordinates, q = ( q x q y q z ) = q (cos ϕ q sin ϑ q sin ϕ q sin ϑ q cos ϑ q ) , discretizing modulus q , az-imuthal angle ϕ q , and inclination angle ϑ q . This helps to keep the isotropy in quiescent calculations and simplifiesidentifying spherical symmetries in the computed observables. Under shear however, the q grid becomes completelyanisotropic in three dimensions. The V qkp ( t ) is · d dimensional. The computation of k ( t ) respects shear advectionvia coordinate transformations and depends on the according q ( t = 0) , which has been chosen as k z axis; see Amann(2013) for more details. Table I shows the generic parameter choice used for the computation of the results of thiswork as best trade-off between computation time and precision. HS q ADG q ∆ q π/ ∆ ϕ q π/ ∆ ϑ q HS φ c [0.2;39.8] [0.2;79.8] 0.4 24 24 0.515712(1)Table I. Generic parameters of the q grid used for this work. φ c is rounded up in the seventh digit. With
OpenMP parallelization on 32 CPUs á . GHz and GB RAM, the calculation of one stress-strain curvein the ADG takes up to 90 hours. Only the repeated calculation of m q ( t ′ ) ( τ ) can be parallelized effectively to gaincomputation speed. Hence, a compromise in precision and computation time must be accepted. A method of savingcomputation time in the iteration of Φ q ( t ) , which needs a repeated calculation of m q ( t ′ ) ( τ ) , is to store the vertex for theyoungest time instance t , which consumes a relevant fraction of the available RAM and limits the grid discretization.Thus, RAM consumption and computation time increase approximately with the square of the grid-point density.1 REFERENCES
Amann, C. P.,
Time dependent flows in arrested states , Ph.D. thesis, Universität Konstanz (2013).Amann, C. P., F. Weysser, M. Fuchs, M. Siebenbürger, M. Ballauff and M. Krüger, “Overshoots in stress strain curves: Colloid experimentsand schematic mode coupling theory,” J. Rheol. , 149–175 (2013).Ballauff, M., J. M. Brader, S. U. Egelhaaf, M. Fuchs, J. Horbach, N. Koumakis, M. Krüger, M. Laurati, K. J. Mutch, G. Petekidis,M. Siebenbürger, Th. Voigtmann and J. Zausch, “Residual stresses in glasses,” Phys. Rev. Lett. , 215701 (2013).Baxter, R. J., “Ornstein-Zernike Relation and Percus-Yevick Approximation for Fluid Mixtures,” J. Chem. Phys. , 4559 (1970).Bergenholtz, J. and M. Fuchs, “Nonergodicity transitions in colloidal suspensions with attractive interactions,” Phys. Rev. E , 5706(1999).Bergenholtz, J., M. Fuchs and Th. Voigtmann, “Colloidal gelation and non-ergodicity transitions,” J. Phys.: Condens. Matter , 6575(2000).Bergenholtz, J., W. C. K. Poon and M. Fuchs, “Gelation in model colloid-polymer mixtures,” Langmuir , 4493–4503 (2003).Brader, J. M., M. E. Cates and M. Fuchs, “First-principles constitutive equation for suspension rheology,” Phys. Rev. E , 021403 (2012).Brader, J. M., Th. Voigtmann, M. Fuchs, R. G. Larson and M. E. Cates, “Glass rheology: From mode-coupling theory to a dynamicalyield criterion,” Proc. Natl. Acad. Sci. U.S.A. , 15186–15191 (2009).Crassous, J. J., M. Siebenbuerger, M. Ballauff, M. Drechsler, D. Hajnal, O. Henrich and M. Fuchs, “Shear stresses of colloidal dispersionsat the glass transition in equilibrium and in flow,” J. Chem. Phys. , 204902 (2008).Dawson, K., G. Foffi, M. Fuchs, W. Götze, F. Sciortino, M. Sperl, P. Tartaglia, Th. Voigtmann and E. Zaccarelli, “Higher-order glass-transition singularities in colloidal systems with attractive interactions,” Phys. Rev. E , 011401 (2000).Eckert, T. and E. Bartsch, “Re-entrant glass transition in a colloid-polymer mixture with depletion attractions,” Phys. Rev Lett. ,125701 (2002).Einstein, A., “Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendiertenTeilchen,” Ann. Phys. , 549 (1905).Fabbian, L., W. Götze, F. Sciortino, P. Tartaglia and F. Thiery, “Ideal glass-glass transitions and logarithmic decay of correlations in asimple system,” Phys. E , 1347 (1999).Farage, T. F. F., J. Reinhardt and J. M. Brader, “Normal-stress coefficients and rod climbing in colloidal dispersions,” Phys. Rev. E ,042303 (2013).Franosch, T., M. Fuchs, W. Götze, M. R. Mayr and A. P. Singh, “Asymptotic laws and preasymptotic correction formulas for the relaxationnear glass-transition singularities,” Phys. Rev. E , 7153–7176 (1997).Fuchs, M., “Nonlinear Rheological Properties of Dense Colloidal Dispersions Close to a Glass Transition Under Steady Shear,” Adv. Polym.Sci. , 55–115 (2010).Fuchs, M. and M. E. Cates, “Theory of nonlinear rheology and yielding of dense colloidal suspensions,” Phys. Rev. Lett. , 248304 (2002).Fuchs, M. and M. E. Cates, “A mode coupling theory for brownian particles in homogeneous steady shear flow,” J. Rheol. , 957–1000(2009).Gadala-Maria, F. and A. Acrivos, “Shear-Induced Structure in a Concentrated Suspension of Solid Spheres,” J. Rheol. , 799–814 (1980).Gibaud, T., F. Cardinaux, J. Bergenholtz, A. Stradner and P. Schurtenberger, “Phase separation and dynamical arrest for particlesinteracting with mixed potentials–the case of globular proteins revisited,” Soft Matter , 857 (2011).Götze, W., Complex Dynamics of Glass-Forming Liquids - A Mode-Coupling Theory , Oxford University Press (2009).Götze, W. and M. Sperl, “Logarithmic relaxation in glass-forming systems,” Phys. Rev. E , 011405 (2002).Hajnal, D. and M. Fuchs, “Flow curves of colloidal dispersions close to the glass transition Asymptotic scaling laws in a schematic modelof mode coupling theory,” Eur. Phys. J. E , 125–138 (2009).Hansen, J. and I. McDonald, Theory of Simple Liquids , Academic Press, London, rd edn. (2009).Henrich, O., F. Weysser, M. E. Cates and M. Fuchs, “Hard discs under steady shear: comparison of brownian dynamics simulations andmode coupling theory,” Phil. Trans. R. Soc. A , 5033–5050 (2009).Hunter, G. L. and E. R. Weeks, “The physics of the colloidal glass transition,” Rep. Prog. Phys , 066501 (2012).Kim, J. M., J. Fang, A. P. R. Eberle, R. Castaneda-Priego and N. J. Wagner, “Gel transition in adhesive hard-sphere colloidal dispersions:The role of gravitational effects,” PRL , 208302 (2013).Koumakis, N., M. Laurati, S. U. Egelhaaf, J. F. Brady and G. Petekidis, “Yielding of hard-sphere glasses during start-up shear,” Phys.Rev. Lett. , 098303 (2012a).Koumakis, N., A. Pamvouxoglou, A. S. Poulos and G. Petekidis, “Direct comparison of the rheology of model hard and soft particleglasses,” Soft Matter , 4271–4284 (2012b).Koumakis, N. and G. Petekidis, “Two step yielding in attractive colloids: transition from gels to attractive glasses,” Soft Matter ,2456–2470 (2011).Krüger, M., F. Weysser and M. Fuchs, “Tagged-particle motion in glassy systems under shear: Comparison of mode coupling theory andBrownian dynamics simulations,” Eur. Phys. J. E , 1–22 (2011).Lekkerkerker, H., W. Poon, P. Pusey, A. Stroobants and P. Warren, “Phase behaviour of colloid + polymer mixtures,” Europhys. Lett. , 559 (1992).Lindemann, F. A., “Über die Berechnung molekularer Eigenfrequenzen,” Physik. Z. , 609–612 (1910).Lu, P. J., E. Zaccarelli, F. Ciulla, A. B. Schofield, F. Sciortino and D. A. Weitz, “Gelation of particles with short-range attraction,” Nature , 499 (2008).Miyazaki, K., D. R. Reichman and R. Yamamoto, “Supercooled liquids under shear: Theory and simulation,” Phys. Rev. E , 011501(2004).Petekidis, G., A. Moussaid and P. Pusey, “Rearrangements in hard-sphere glasses under oscillatory shear strain,” Phys. Rev. E , 051402(2002).Pham, K. N., G. Petekidis, D. Vlassopoulos, S. U. Egelhaaf, W. C. K. Poon and P. N. Pusey, “Yielding behavior of repulsion- andattraction-dominated colloidal glasses,” J. Rheol. , 649–676 (2008).Pham, K. N., G. Petekidis, D. Vlassopoulos, S. U. Egelhaaf, P. N. Pusey and W. C. K. Poon, “Yielding of colloidal glasses,” Europhys. Lett. , 624–630 (2006).Pham, K. N., A. M. Puertas, J. Bergenholtz, S. U. Egelhaaf, A. Moussaid, P. N. Pusey, A. B. Schofield, M. E. Cates, M. Fuchs and W. C. K. Poon, “Multiple glassy states in a simple model system,” SCIENCE , 104–106 (2002).Priya, M. and Th. Voigtmann “Nonlinear rheology of dense colloidal systems with short-ranged attraction: A mode-coupling theoryanalysis,” J. Rheol. same issue same volume as actual, please insert (2014).Puertas, A. M., M. Fuchs and M. E. Cates, “Comparative simulation studyof colloidal gels and glasses,” PRL , 098301 (2002).Puertas, A. M., C. D. Michele, F. Sciortino, P. Tartaglia and E. Zaccarelli, “Viscoelasticity and stokes-einstein relation in repulsive andattractive colloidal glasses,” , 144906 (2007).Pusey, P. N. and W. van Megen, “Observation of a glass transition in suspensions of spherical colloidal particles,” Phys. Rev. Lett. ,2083 (1987).Siebenbürger, M., M. Fuchs and M. Ballauff, “Viscoelasticity and shear flow of concentrated, noncrystallizing colloidal suspensions: Com-parison with mode-coupling theory,” Soft Matter , 4014–4024 (2012).Siebenbürger, M., M. Fuchs, H. Winter and M. Ballauff, “Viscoelasticity and shear flow of concentrated, noncrystallizing colloidal suspen-sions: Comparison with mode-coupling theory,” J. Rheol. , 707–726 (2009).Sutherland, W., “Dynamical theory of diffusion for non-electrolytes and the molecular mass of albumin,” Phil. Mag. , 781–785 (1905).van Megen, W. and S. M. Underwood, “Glass transition in colloidal hard spheres: Mode-coupling theory analysis,” Phys. Rev. Lett. ,2766 (1993).Voigtmann, Th., J. M. Brader, M. Fuchs and M. E. Cates, “Schematic mode coupling theory of glass rheology: single and double stepstrains,” Soft Matter , 4244–4254 (2012).Willenbacher, N., J. S. Vesaratchanon, O. Thorwarth and E. Bartsch, “An alternative route to highly concentrated, freely flowing colloidaldispersions,” Soft Matter , 5777–5788 (2011).Zaccarelli, E. and W. C. K. Poon, “Colloidal glasses and gels: The interplay of bonding and caging,” Proc. Natl. Acad. Sci. U.S.A. ,15203 (2009).Zausch, J., J. Horbach, M. Laurati, S. U. Egelhaaf, J. M. Brader, Th. Voigtmann and M. Fuchs, “From equilibrium to steady state: thetransient dynamics of colloidal liquids under shear,” J. Phys.: Condens. Matter20