Particles suspended in a Newtonian fluid raise the viscosity and also generally give rise to a shearrate dependent rheology. In particular, pronounced shear thickening may be observed at large solid volume fractions. In a recent article [R. Seto et al., Phys. Rev. Lett. 111, 218301 (2013)], we have considered the minimum set of components to reproduce the experimentally observed shear thickening behavior, including discontinuous shear thickening. We have found frictional contact forces to be essential and were able to reproduce the experimental behavior by a simulation including this physical ingredient along with viscous lubrication. In the present article, we thoroughly investigate the effect of friction and express it in the framework of the jamming transition. The viscosity divergence at the jamming transition has been a well known phenomenon in suspension rheology, as reflected in many empirical laws for the viscosity. Friction can affect this divergence, and in particular the jamming packing fraction is reduced if particles are frictional. Within the physical description proposed here, shear thickening is a direct consequence of this effect: As the shear rate increases, friction is increasingly incorporated as more contacts form, leading to a transition from a mostly frictionless to a mostly frictional rheology. This result is significant because it shifts the emphasis from lubrication hydrodynamics and detailed microscopic interactions to geometry and steric constraints close to the jamming transition.
I. INTRODUCTION
A. Shear thickening
Suspensions (solid particles immersed in a fluid) exhibit a wide range of rheological behaviors, including shear thinning, shear thickening, and finite normal stress differences. Shear thickening [Barnes (1989); Mewis and Wagner (2011); Brown and Jaeger (2014)], where the viscosity increases with shear rate even if the suspending fluid is Newtonian, is a particularly intriguing phenomenon. In the extreme case of discontinuous shear thickening [DST; for early descriptions see the work of Williamson (1930); Williamson and Hecker (1931); Freundlich and Roder (1938)], which is observed at high volume fractions of solid material, the viscosity can increase by several orders of magnitude at a critical shear rate. Note that we do not include here the irreversible shear thickening occurring due to particle aggregation {as seen for example in fumed silica [Crawford et al. (2013)]} under the name DST.
The variety of systems showing a DST suggests that this may be a universal behavior that is obtained with a minimum set of physical ingredients. Even though most of the data available are for Brownian suspensions (that is, submicrometer particles) [Metzner and Whitlock (1958); Hoffman (1972); Bender and Wagner (1995, 1996); Frith et al. (1996); Fagan and Zukoski (1997); Boersma et al. (1990); D'Haene et al. (1993); O'Brien and Mackay (2000); Maranzano and Wagner (2001b,a, 2002)], thermal motion does not seem necessary to observe DST, as experiments with micrometer scale [Boersma et al. (1990); Lootens et al. (2004); Lootens et al. (2005); Larsen et al. (2010)] or larger particles [Freundlich and Roder (1938); Bertrand et al. (2002); Brown and Jaeger (2009, 2012); Fall et al. (2010, 2012)] show. Inertia has also been associated with the existence of shear thickening in suspensions [Lemaître et al. (2009); Fall et al. (2010); Trulsson et al. (2012); Fernandez et al. (2013); Kawasaki et al. (eprint)], based on the initial ideas of Bagnold (1954), but the Bagnoldian scaling of a viscosity proportional to the shear rate is clearly milder than the abrupt shear thickening observed in experiments. More importantly, the thickening due to inertia arises at Stokes number $ S t \u2261 \rho a 2 \gamma \u0307 / \eta 0 $ (with $ \gamma \u0307 $ being the shear rate, η_{0} being the viscosity of the fluid phase, and ρ and a are the mass density and size of the solid particles, respectively) of order $ S t \u2248 1 $ in the simulations [Trulsson et al. (2012); Fernandez et al. (2013); Kawasaki et al. (eprint)], which contrasts with the values of at most $ S t \u2248 10 \u2212 3 $ found in experiments {see for example [Maranzano and Wagner (2001b); Fall et al. (2010); Brown and Jaeger (2012)]}. Thus inertia is not necessary for DST in the conditions probed by rheometric flows. (Flows involving significantly larger particles, say in the mm range, or significantly larger shear rates $ \gamma \u0307 \u226b 1 \u2009 s \u2212 1 $, can probe the regime where inertia matters).
Indeed, the apparently simple experimental system of nearly rigid nonBrownian neutrally buoyant particles immersed in a Newtonian fluid exhibits DST in the Stokes regime [Brown and Jaeger (2009, 2012)], indicating that the phenomenon stems from a rather restricted set of simple ingredients. The puzzle has very few pieces. In the Stokes regime, however, nonBrownian neutrally buoyant hard spheres in a Newtonian fluid will create a suspension whose rheology is independent of shear rate, as there is only one force scale, the hydrodynamic one (we will further explain this fact in Sec. II B 1). So there cannot be too few pieces.
B. Fluid mechanics and granular physics perspectives
The flow of suspensions has historically been studied from a fluid mechanics perspective, and the emphasis has been on a description based on hydrodynamic interactions. Suspensions are usually described as hard particles immersed in a Newtonian fluid, interacting through hydrodynamics (including Brownian forces) and sometimes an additional soft repulsive potential (approximating an electric doublelayer or mimicking polymer coating, for instance). They are typically studied in the Stokes flow regime, at vanishing particle Reynolds number $ R e \u2261 \rho 0 a 2 \gamma \u0307 / \eta 0 $ (with ρ_{0} being the mass density of the fluid phase) and Stokes number. The key point is that the hard cores of the particles are treated as boundary conditions for the Stokes equations that describe the fluid phase; the particles never directly generate forces through contacts. This treatment is selfconsistently justified by the fact that the Stokes flow between two rigid surfaces leads to a lubrication force whose resistance coefficient diverges at contact, effectively preventing two particles from colliding. Within this framework, shear thickening is explained by the creation at large shear rates of locally denser clusters of particles (hydroclusters), which are highly dissipative due to the singular lubrication flows between the particles [Brady and Bossis (1985); Bender and Wagner (1996); Melrose and Ball (2004a); Wagner and Brady (2009)]. While this purely fluid mechanical point of view is able to describe the flow and rheology of moderately concentrated suspensions, it seems unable to explain the abrupt shear thickening observed in highly concentrated suspensions. Simulations by Stokesian Dynamics give a weak logarithmic shear thickening [Brady and Bossis (1988); Bossis and Brady (1989); Phung et al. (1996); Foss and Brady (2000); Melrose and Ball (2004b)], thereby raising the issue of the validity of this approach at high volume fractions. Moreover, DST is also observed in athermal suspensions, where the shear rate dependence introduced by the Brownian force scale disappears: In this regime, a purely hydrodynamic perspective would predict a shearrate independent rheology.
In the past few years, new ideas have emerged from the granular rheology perspective. Boyer et al. (2011) developed an analogy between the rheology of suspensions and the rheology of granular materials. In dry granular flow, the rheology depends on the ratio between inertia and particle pressure [Deboeuf et al. (2009)]: If the inertia dominates, particles bounce around and the stresses are essentially due to momentum exchange based on collisions (a “granular gas” regime), whereas if the pressure dominates, particles are forced to stay at contact and the stresses are dominated by contact force chains (a “granular packing” regime). This perspective can be incorporated in a dimensionless inertial number $ I \u2261 \gamma \u0307 d \rho / \Pi $ [da Cruz et al. (2005)], where d is the diameter of the solid particles, ρ their density, and Π is the particle pressure. Similarly, for suspensions in the Stokes regime, there is a viscous number $ I v \u2261 \eta 0 \gamma \u0307 / \Pi $, where η_{0} is the viscosity of the suspending fluid, which compares viscous dissipation to the particle pressure [Boyer et al. (2011)]. The stresses can be dominated either by viscous dissipation if $ I v \u226b 1 $ (when particles are far apart) or by a contact network if $ I v \u226a 1 $. This results in constitutive laws for the viscosity and other stress components and for the volume fraction that are unique functions of I_{v}. One can anticipate that the regimes $ I \u226a 1 $ and $ I v \u226a 1 $, both of which are dominated by the approach to the jamming transition, will share some similarities. Indeed, they do share a powerlaw scaling (typical of the jamming transition) between particle pressure and distance to the jamming point $ \Pi \u223c ( \varphi J \u2212 \varphi ) \u2212 \lambda $ [Forterre and Pouliquen (2008); Boyer et al. (2011)]. The scaling is the same for the other stress components, including the important case of shear stress $ \sigma \u223c ( \varphi J \u2212 \varphi ) \u2212 \lambda $. These power laws are now argued to be a generic behavior close to jamming [Lerner et al. (2012); Andreotti et al. (2012)], though the value of the exponent λ might be system dependent (it might depend on shape or friction for instance). With these ideas, the emphasis has shifted from hydrodynamics and detailed microscopic interactions to geometry and steric constraints close to the jamming transition.
While those new ideas have proven to be successful in explaining some of the rheology at high volume fractions [Boyer et al. (2011); Trulsson et al. (2012); Lerner et al. (2012)], they cannot account for nonNewtonian behavior in the Stokes regime. The case of DST is particularly puzzling: Whereas intuition suggests that DST is a manifestation of jamming [Cates et al. (1998); Bertrand et al. (2002); Lootens et al. (2003); Hébraud and Lootens (2005); Brown and Jaeger (2009)], DST is not captured by the I_{v}based rheology, which predicts no shear rate dependence in the Stokes regime [Lerner et al. (2012); Trulsson et al. (2012)], and reduces to a $\varphi $dependent rheology [Boyer et al. (2011)]. Simply stated, this rheology predicts the correct volume fraction dependence of the viscosity at fixed shear rate but completely fails at predicting the shear rate dependence at fixed volume fraction. Again, it is worth noticing that this absence of shear rate dependence hints at a missing force (or time) scale in the description of dense suspensions.
C. This work
In this article, which extends and complements a previous publication [Seto et al. (2013)], we show how jamming is related to shear thickening and how to escape the apparent contradiction described above. The starting point is to note that the jamming transition, while it exhibits some generic features that are independent of the microscopic interactions like the scalings described above, retains a clear signature of the microscopic details in the volume fraction at which it occurs. For instance, it is well known that jamming can occur anywhere between $ \varphi J \mu = \u221e \u2248 0.55 $ and $ \varphi J \mu = 0 \u2248 0.64 $ depending on the friction coefficient μ between the grains [Silbert et al. (2002); Jerkins et al. (2008); Song et al. (2008); Silbert (2010)]. This means that the relation $ \varphi ( I v ) $ actually depends on μ, which in turn implies that the viscosity $ \eta ( \varphi ) $ depends on μ. Therefore, if there exists a mechanism such that μ (or an effective friction coefficient) varies with the shear rate $ \gamma \u0307 $, one can obtain a nonNewtonian rheology within the framework of Boyer et al. (2011). This mechanism has to be associated with an extra force scale besides the hydrodynamic one in order to provide the shearrate dependence.
We introduced such a mechanism in the recent article [Seto et al. (2013)], and here we thoroughly explore its consequences through numerical simulations of dense (i.e., highly concentrated) nonBrownian suspensions under simple shear flow. We note that two other recent papers have independently introduced similar but somewhat different mechanisms to achieve a ratedependent friction coefficient $ \mu ( \gamma \u0307 ) $ in simple shear flow simulations [Fernandez et al. (2013); Heussinger (2013)]. Both of them prove that interparticle friction is essential to shear thickening of dense suspensions. In contrast to those works, the present work however shows that inertia is inessential to this phenomenon, as we obtain it in an overdamped simulation (see Sec. II). This is important because results coming from simulations including inertia give shear thickening at shear rates that are several orders of magnitude larger than the experimental values (see previous Sec. I A). We show here that our inertialess model gives a quantitative description of shear thickening (see Sec. V).
The present work shows the importance of frictional contacts in dense suspension rheology. Of course, particles must form contacts during flow in order for the rheology to depend on the friction coefficient, and here we come up against one consequence of the fluid mechanics perspective. In fact, it is known that particles do come in contact under flow, even in the dilute limit, perhaps due to surface roughness [Davis (1992); Zhao and Davis (2002); Lootens et al. (2005); Blanc et al. (2011)]. One simply has to accept that those contacts, though not essential to the rheology of dilute or semidilute suspensions (as the successes of Stokesian Dynamics demonstrate), become important at large volume fractions. But contacts are a necessity if one thinks of the large volumefraction limit of jamming $ ( I v \u2192 0 ) $: There is no more room in this limit for lubrication films between particles, and contacts must proliferate and dominate the rheology. Besides enforcing the geometric constraint of no overlapping particles, those contacts bring a new restriction to reorganization at the microscopic level: They carry significant tangential forces due to friction. It is worth noting that even though the closerange hydrodynamics (i.e., lubrication forces) also generate some tangential forces, they are of a fundamentally different nature. For a relative motion between two nearby particles, the normal lubrication force diverges as the inverse of the interparticle gap, whereas the tangential lubrication force diverges only logarithmically with vanishing gap, which means that the effective friction coefficient vanishes in the limit of zero gap for smooth bodies. Thus lubrication, as in its literal meaning, provides very little resistance to relative tangential motion compared to relative normal motion. This stands in contrast to the frictional contact forces, for which the friction coefficient is finite; i.e., tangential and normal contact forces are of comparable scale. It follows that the constraint introduced by contact friction is much more efficient at increasing the viscosity than is lubrication.
Describing the physics of contact and motion when the gap between two particles in a suspension is smaller than, say, 1 nm, is not an easy task. Moreover, one cannot expect to find generic behavior, as interactions will vary depending on the nature of the solid particles and of the suspending fluid. However, this level of detail may not be essential to understanding the physics that emerges at a macroscopic level. At high volume fractions, it is reasonable to expect that the singularity associated with the jamming transition will dominate the rheology of suspensions. Our approach in this work is thus to study the rheology of minimal model systems that include hydrodynamics, a repulsive interaction, and frictional contacts between particles. These interactions, while being realistic, will be simplified to their essence.
II. MODELS AND METHODS
In this section, we will describe our numerical model of a dense nonBrownian suspension under shear flow [Denn and Morris (2014)]. We intend to have a minimal model containing as few physical elements as possible while exhibiting a DST. Even if the model is intended to be as generic as possible, a few parameter choices are necessary. For those choices, we try to set values that correspond to a realistic interpretation of a suspension of hard (say, silica) spheres with radius in the 1–10 μm range, density matched (or such that the sedimentation occurs over times much larger than the inverse shear rate) and charge stabilized. The only parameter that takes an unrealistic value is the particles' stiffness, which is smaller in the simulations than in a typical experimental system. A stiffness comparable to experimental values for silica (or even Polymethyl methacrylate (PMMA)), would require the use of very small time steps in the simulation in order to resolve the dynamics of the contacts which correspond to an overlap in the nm range. We however circumvent this problem by tuning the particle stiffness with applied shear rate such that we get rid of particle deformation effects in the rheology (see Sec. II A 3 and Appendix B for details).
A. Merging hydrodynamic interactions and granular contact models
1. Equations of motion
In a suspension, the flow of the fluid is described by the Navier–Stokes equation, while the motion of the particles is given by Newton's equations of motion
where $m$ is the mass/momentofinertia matrix, $U$ and $\Omega $ are the translational and rotational velocity vectors, and $ F \alpha $ and $ T \alpha $ are the force and torque vectors, respectively. The righthand side consists of two different types of forces and torques: Some depend only on the configurations $ r N $ of the particles (e.g., forces derived from a potential), but some also depend on the velocities $U$ and $\Omega $, including inelastic contact forces (e.g., contacts with dashpot terms) and fluidparticle interactions (hydrodynamic interactions).
We will study the rheology of this suspension under an imposed simple shear flow field $ U \u221e ( r ) $ expressed using the vorticity $ \Omega \u221e $ and the rateofstrain tensor $ E \u221e $ as
At a shear rate $ \gamma \u0307 $, a simple shear flow corresponds to the following nonzero elements: $ \Omega 3 \u221e = \gamma \u0307 / 2 $ and $ E 12 \u221e = E 21 \u221e = \gamma \u0307 / 2 $.
In many experimental flows, including the ones for which shear thickening is typically observed, the particle Reynolds number and the Stokes number are very small due to the small size of the suspended particles. This means that the equation of motion for the particles can be studied in its overdamped limit, which is simply a quasistatic force balance equation. As will be discussed in Secs. II A 2, II A 3, and II B below, in this regime the velocity dependence coming from the hydrodynamic interactions and the contact forces has a linear form, and the force balance equation can be written as a linear algebraic relation having the form
Therefore, our simulation of a suspension under a simple shear flow in the Stokes regime requires solving (3) for the velocities, and obtaining the positions $ r N $ of the particles at any time t through time integration of these velocities.
2. Hydrodynamic interactions
In principle, hydrodynamic interactions can be determined by solving the Stokes equations, but this is extremely expensive from a computational view point. For dense suspensions, however, the dominant hydrodynamic interactions come from the fluid flow in the narrow gaps between nearby solid particles [Ball and Melrose (1997)], because the resistance to relative motion becomes singular when the gaps are vanishingly small. They give rise to pairwise shortrange lubrication forces (as opposed to the manybody nature of the full, longrange hydrodynamic interactions). As a consequence, the linear relations between velocities and hydrodynamic forces simply contain a contribution from the Stokes drag and a contribution from the lubrication
where $ R Stokes $ is a diagonal matrix giving Stokes drag forces and torques, and $ R Lub $ and $ R \u2032 Lub $ are sparse matrices [Ball and Melrose (1997)].
Consistent with our choice of keeping only physically relevant neardistance interactions, we use only the leading terms for the resistance matrices. Physically, the terms included in our model correspond to the squeeze, shear, and pump modes of Ball and Melrose (1997) (we do not consider their twist mode, as it is not associated with a divergence of the resistance at contact and is thus subdominant). Defining the nondimensional gap h^{(i,j)} between particles i and j having radii a_{i} and a_{j} as $ h ( i , j ) \u2261 2 ( r ( i , j ) \u2212 a i \u2212 a j ) / ( a i + a j ) $ and the centertocenter unit vector $ n i j \u2261 ( r ( j ) \u2212 r ( i ) ) / r ( i , j ) $ with $ r ( i , j ) \u2261  r ( j ) \u2212 r ( i )  $, the modes associated with relative displacements, respectively, along and tangential to $ n i j $ have singular behaviors with leading terms diverging as 1/h^{(i,j)} and $ log ( 1 / h ( i , j ) ) $ [Jeffrey and Onishi (1984); Jeffrey (1992)]. The consequence of the singularity in the squeeze mode (i.e., along $ n i j $) is that there should not be any contact between particles [Ball and Melrose (1995)]. However, as mentioned in the introduction, this statement is only true for perfectly idealized situations and not realistic even for model experimental systems consisting of spherical particles due to, for example, a finite surface roughness [Davis (1992); Zhao and Davis (2002); Lootens et al. (2005); Blanc et al. (2011)] or a breakdown of the continuity assumption for the fluid at the molecular mean free path scale [Ho and Tai (1998)]. In order to mimic this reality, we regularize the singularities arising in the lubrication resistances by inserting a small cutoff length scale δ [Trulsson et al. (2012)], which can be thought of as the length scale of the particle surface roughness; the leading terms we use for normal and tangential displacements then behave as $ 1 / ( h ( i , j ) + \delta ) $ and $ log ( 1 / ( h ( i , j ) + \delta ) ) $. In the simulations, we set δ = 10^{−3}. This value is typical for nonBrownian particles [Davis (1992); Blanc et al. (2011)]. The detailed expressions for the resistance matrices $ R Lub $ and $ R \u2032 Lub $ are given in Appendix A.
3. Contact model
There are several models that one can use to describe frictional contacts between particles. We use a stick/slide friction model employing springs and dashpots that is commonly used in granular physics [Cundall and Strack (1979); Luding (2008)]. The normal and tangential components of the force and the torque for particles having radii a_{i} and a_{j} are obtained as
and fulfill Coulomb's friction law $  F C , \u2009 tan \u2009 ( i , j )  \u2264 \mu  F C , nor ( i , j )  $. In the above expressions, k_{n} and k_{t} are the normal and tangential spring constants, respectively, and γ_{n} is the damping constant. The normal and tangential velocities are $ U n ( i , j ) \u2261 n i j n i j \xb7 ( U ( j ) \u2212 U ( i ) ) $ and $ U t ( i , j ) \u2261 ( I \u2212 n i j n i j ) \xb7 [ U ( j ) \u2212 U ( i ) \u2212 ( a i \Omega ( i ) + a j \Omega ( j ) ) \xd7 n i j ] $. Finally, the quantity $ \xi ( i , j ) $ is the tangential spring stretch. This contact model could be made more general, but at the price of numerical difficulties; see Appendix B.
The computation of the tangential spring stretch $ \xi ( i , j ) $, described in the following, requires some care, as we have to impose Coulomb's law. We apply an algorithm described in Luding (2008). At the time t_{0} at which the contact (i,j) is created, we set an unstretched tangential spring $ \xi ( i , j ) ( t 0 ) = 0 $. At any further time step t in the simulation, the tangential stretch $ \xi ( i , j ) ( t ) $ is incremented according to the value of a “test” force $ F C , \u2009 tan \u2009 \u2032 ( i , j ) ( t + d t ) = k t \xi \u2032 ( i , j ) ( t + d t ) $ with $ \xi \u2032 ( i , j ) ( t + d t ) = \xi ( i , j ) ( t ) + U t ( i , j ) ( t ) d t $. If $  F C , \u2009 tan \u2009 \u2032 ( i , j ) ( t + d t )  \u2264 \mu  F C , nor ( i , j ) ( t + d t )  $, the contact is in a static friction state and we update the spring stretch and force as
However, if $  F C , \u2009 tan \u2009 \u2032 ( i , j ) ( t + d t )  > \mu  F C , nor ( i , j ) ( t + d t )  $, the contact is in a sliding state and the spring and force are updated as
where the direction $ t i j $ is the same as the one for the test force, i.e., $ t i j \u2261 F C , \u2009 tan \u2009 \u2032 ( i , j ) ( t + d t ) /  F C , \u2009 tan \u2009 \u2032 ( i , j ) ( t + d t )  $. In this case, Coulomb's law is not violated, as $  F C , \u2009 tan \u2009 ( i , j ) ( t + d t )  = \mu  F C , nor ( i , j ) ( t + d t )  $.
B. Other interactions, full models
1. Shear rate dependence: Why Stokes hydrodynamics plus hard spheres is not enough
The shearrate dependence of a suspension requires the existence of a time scale distinct from the inverse shear rate. A well known case is the one of Brownian suspensions, where the inverse shear rate competes with the Brownian diffusion time to give a shearrate dependent rheology. This competition is adequately measured by a nondimensionalized shear rate, the Péclet number, which is the ratio of the two time scales: $ P e \u2261 6 \pi \eta 0 a 3 \gamma \u0307 / k B T $. This quantity can also be thought of as a ratio of two force scales: The typical hydrodynamic force and the typical Brownian force. Thus, a shearrate dependent rheology generically requires at least one other force scale besides the hydrodynamic one.
Such a time scale does not exist for nonBrownian suspensions of hard spheres in the Stokes regime. The reason is straightforward: The hard sphere force has no typical value, and a hardsphere contact between two particles can withstand any applied load. Therefore, there is nothing with which to compare the hydrodynamic force. This means that the solutions of the force/torque balance equations (3) (i.e., the particle trajectories) are independent of the shear rate, leading to a shearrate independent rheology. Notice that this conclusion holds whether the hard spheres are frictionless or frictional: Coulomb's law does not introduce a force scale. There must be another force scale in the system to provide the shearrate dependence.
An important point in our model is that, even if we try to mimic a hard sphere suspension, our contact forces actually are elastic forces, thus bringing a natural force scale k_{n} a. We preserve the shearrate independent hard sphere behavior by selecting the particle stiffness k_{n} such that the nondimensional number based on the stiffness $ 6 \pi \eta 0 a \gamma \u0307 / k n $ remains much smaller than unity; i.e., we stay in the asymptotic regime of nearly hard particles. We describe the techniques we use to achieve this in our simulation in Appendix B.
2. A minimal model: Criticalload friction
The first model we consider is an intentionally simplistic one, arguably the simplest model with an additional force scale besides the hydrodynamic force. Simplicity is the main motivation for this model, as it enables the clear physical discussion one can expect from a minimal model.
In this model (which we call critical load model, CLM), we introduce an extra force scale in the friction law itself, namely a threshold normal load F_{CL} to activate friction. When the normal load is smaller than this threshold, particles interact as frictionless hard spheres. When the load is beyond the threshold, friction is activated. Overall, the friction law reads
We may consider the CLM as the short Debye length limit of the electrostatic repulsion model (ERM) presented below (Sec. II B 3).
3. Electrostatic repulsion model
An electrostatic repulsion is a simple and plausible interaction. Many experimental systems have such a force, as it stabilizes the suspension. Thus, we will consider an ERM including this ingredient in the simplest form of a repulsive doublelayer electrostatic force $ F R ( i , j ) $ [Israelachvili (2011)]. If the Debye length κ^{−1} is small compared to the radius of the particles, the approximate form
can be used. In the simulations, we set κ^{−1 }= 0.05a, meaning a Debye length in the 0.1 μm range for particles of a few μm; see Mewis and Wagner (2011).
C. Stress tensor and bulk rheology
The mechanical stress applied to the suspension arises from the several interactions included in the model: Particles disturb the flow, creating hydrodynamic stresses, and they develop force chains via contacts (and/or electrostatic repulsion for the ERM).
The hydrodynamic stresslets acting on the particles are given by [Batchelor (1970); Brady and Bossis (1988)]
where $ S H \u2261 ( S H ( 1 ) , \u2026 , S H ( n ) ) $ and the matrices $ R Lub S $ and $ R Lub \u2032 S $ contain leading terms of the lubrication resistances [Jeffrey and Onishi (1984); Jeffrey (1992)] (see Appendix A) in a manner consistent with the hydrodynamic forces considered in Sec. II A 2.
The stress due to the contact or repulsive force between particles i and j is simply written as
respectively.
The bulk stress, in which the isotropic part of the fluid pressure is omitted, is the sum of the above contributions
where V is the volume of the simulation box (the electrostatic stress $ S R $ appears only for the ERM). Shear stress σ, normal stress differences N_{1} and N_{2}, and particle pressure Π [Yurkovetsky and Morris (2008)] are defined as $ \sigma \u2261 \Sigma 12 , \u2009 N 1 \u2261 \Sigma 11 \u2212 \Sigma 22 , N 2 \u2261 \Sigma 22 \u2212 \Sigma 33 $, and $ \Pi \u2261 \u2212 ( \Sigma 11 + \Sigma 22 + \Sigma 33 ) / 3 $, respectively. The relative viscosity is given by $ \eta r \u2261 \sigma / \eta 0 \gamma \u0307 $.
D. Additional points in the simulation model

Bidispersity: We consider a suspension of bidisperse frictional hard spheres immersed in a Newtonian fluid with viscosity η_{0}. The bidispersity is introduced to hinder the formation of the ordered phase observed for dense monodisperse suspensions under shear. An effective choice for the size ratio is a_{2}/a_{1} = 1.4 (where a_{1} = a), with the two populations occupying equal volumes; i.e., $ \varphi 1 = \varphi 2 $. Other effective size ratios are possible, as is polydispersity, with the same qualitative effects. The choice of size ratio 1:1.4 is however the most common one in the granular matter literature, where bidispersity is used to avoid crystallization [see for example O'Hern et al. (2003)]. Indeed, with a smaller size ratio 1:1.2, slow and weak strain thinning due to ordering is seen in the low viscosity state.

Dimensionless shear rate: We discussed in Sec. II B 1 that another force scale besides the hydrodynamic one is required to yield shearrate dependence. In CLM, the threshold value gives the force scale $ F \u2217 = F CL $, and in ERM, the force at contact gives $ F \u2217 = F ER $. Therefore, the shear rate dependence is given by the ratio $ \gamma \u0307 / \gamma \u0307 0 $, with $ \gamma \u0307 0 \u2261 F \u2217 / 6 \pi \eta 0 a 2 $.

Periodic boundary conditions and fixedvolume simulation: Rheology is a bulk property. If solid walls are used for the boundaries of the system, we need a large system size to reduce the influence of walls. We can avoid the use of solid walls in a simulation by using periodic boundary conditions. For the straincontrolled simple shear, we use the Lees–Edwards [Lees and Edwards (1972)] boundary condition.
Particle migration never develops under these periodic boundary conditions. Although some experimental observations suggest that global migration may be a cause of shear thickening [Fall et al. (2010)], our simulation is free of this effect.
There is no way for the system to dilate with these periodic boundary conditions. Such a fixed volume condition is expected in most rheology measurements. For shear thickening fluids, however, the open edges may have some influence [Brown and Jaeger (2012)]. In granular physics, systems are often sheared under a given normal stress, hence the volume is not fixed [Boyer et al. (2011)].
III. RESULTS
This section presents the results from the simulation of the two models, the “minimal” CLM and the ERM. In Secs. III A–III E, whenever possible, the two models are treated at the same time in the text, and we will show data plots for the CLM. However, when the two models give slightly different results, the results of the CLM will be described first, as it allows a simple and clear understanding of the underlying physics, and the more realistic model ERM will be described afterward. The differences are essentially in the low shearrate limit for the ERM (existence of a shearthinning regime), due to the repulsive potential.
Although friction is the key factor in this work, few experimental estimates are available for the friction coefficient. We select a friction coefficient μ = 1 [which is comparable to the measurements of Fernandez et al. (2013) for ∼10 μm quartz particles coated by polymer brushes] for most of the simulations, except when the dependence on μ is specifically investigated.
In the data plots shown in this section, the error bars represent the standard deviation, which we define for an observable A(t) as $ \u27e8 A 2 \u27e9 \u2212 \u27e8 A \u27e9 2 $, where $ \u27e8 \xb7 \u27e9 $ is the time (strain) average $ T \u2212 1 \u222b 0 T \xb7 \u2009 d t $ over the simulated units of strain. (We perform simulations over 50 strain units, hence T = 50 except when the time averages are performed over subsets of the whole simulation in the case of intermittent data around DST.)
Rheological data are plotted versus shear rates and shear stresses nondimensionalized by the characteristic shear rate $ \gamma \u0307 0 = F \u2217 / 6 \pi \eta 0 a 2 $ and stress $ \eta 0 \gamma \u0307 0 $, respectively, where η_{0} is the viscosity of the suspending fluid.
A. Frictionless and frictional rheologies
In the CLM, due to the threshold force, the friction between grains is absent at low shear rates and activated at high shear rates. Because of this, we expect that the low shearrate limit for concentrated suspensions will have a rheology $ \eta ( \varphi , \gamma \u0307 \u2192 0 ) $ typical of a system close to the jamming transition for frictionless particles, while the high shearrate limit shows a rheology $ \eta ( \varphi , \gamma \u0307 \u2192 \u221e ) $ typical of a system close to the jamming transition for particles with a friction coefficient μ = 1.
These two limiting viscosities are shown in Fig. 1, where we also show the high shearrate behavior for the infinite friction case (μ = ∞) for reference. Each diverges at a different volume fraction, thus friction shifts the jamming point [Otsuki and Hayakawa (2011)]. We fit our data with powerlaw divergences $ \eta \u221d C ( 1 \u2212 \varphi / \varphi J ) \u2212 \lambda $, with parameters $ ( \varphi J , \lambda , C ) $ as detailed in the caption of Fig. 1.
For the ERM, the situation is the same at high shear rates but differs at low shear rates. While friction is not felt for $ \gamma \u0307 \u2192 0 $, as particles do not contact, the finite range of the repulsive potential creates a shearthinning behavior from which we could not obtain the low shearrate limit $ \eta ( \varphi , \gamma \u0307 \u2192 0 ) $. The shear thinning comes from the fact that the system behaves essentially as soft particles at low shear stresses, with an apparent size that includes the hard sphere and a part of the surrounding soft repulsive potential. Indeed, a simple force balance argument at the scale of a particle says that if a particle is subject to a driving shear force σa^{2}, the minimum gap h_{min} between this particle and its neighbors will be such that $ F R ( h min ) \u223c \sigma a 2 $ [with the limitation that it must respect the geometrical constraint, i.e., $ h min < ( 1 \u2212 \varphi / \varphi J 0 ) 1 / 3 $]. Then, at this shear stress the minimum centertocenter distance between two particles is $ 2 a ( 1 + h min ( \sigma ) / 2 ) $, which means that the particles have an apparent radius larger than the one of their hard core. The jamming transition of these effectively bigger frictionless particles is at an apparent packing fraction (that is, based on the apparent diameter) $ \varphi \u2032 = \varphi J \u2032 0 \u2248 0.66 $, which corresponds to an actual packing fraction $\varphi $ significantly lower than $ \varphi J \u2032 0 $. Because of this low volume fraction jamming transition at $ \gamma \u0307 = 0 $, the viscosity is larger than at finite (but small) $ \gamma \u0307 $, leading to the shear thinning we observe. Such thinning is actually known to occur due to repulsive forces in charge stabilized suspensions [Krieger (1972); Maranzano and Wagner (2001a)].
B. Shear thickening, continuous, and discontinuous
We can switch from one rheology to the other by varying the shear rate. Physically, the transmitted stress increases as the shear rate increases, which triggers the formation of frictional contacts between particles. Thus, by increasing the shear rate, the viscosity interpolates between the frictionless and frictional rheology curves, which means we can observe shear thickening. All this should be a natural consequence of the existence of two distinct rheologies at $ \gamma \u0307 = 0 $ and $ \gamma \u0307 = \u221e $. What we cannot anticipate a priori is the way in which the system switches from the low viscosity state to the high viscosity one: Do we observe a continuous shear thickening (CST) or a DST?
The shear rate dependence of the viscosity, shown in Fig. 2 for the CLM, demonstrates the existence of both CST and DST in our system, depending on the volume fraction. As in experiments, when $ \varphi < \varphi c $ the shear thickening is continuous, getting steeper and steeper as we approach $ \varphi c $, at which point it becomes discontinuous and keeps this behavior for $ \varphi > \varphi c $ and up to $ \varphi J 0 $. Note that there appears to be a real discontinuity in our data for these volume fractions: The time series of the viscosity show an intermittent behavior switching between two states, which we address in Sec. III C. In Fig. 2, the intermittent data are split between low and high viscosity states: Two points that correspond to a separate time average for each of the two states appear at the same shear rate.
When plotted against stress in Fig. 2, the viscosity curves show another interesting feature, namely that the onset of shear thickening occurs at a stress $ \sigma ST / ( \eta 0 \gamma \u0307 0 ) \u2248 0.3 $ that is roughly independent of the volume fraction, as is observed in many experiments [Frith et al. (1996); Bender and Wagner (1996); Maranzano and Wagner (2001b,a); Lootens et al. (2005); Fall et al. (2010); Larsen et al. (2010); Brown and Jaeger (2012, 2014)]. Similarly, the transition toward the shearrate independent plateau at high viscosity also occurs at a stress independent of $\varphi $. The stress range over which thickening occurs remains constant from a mild shear thickening to a marked DST, as expected from a simple balance argument between the driving stress $ \sigma = \eta \gamma \u0307 $ and the stress required to create a frictional contact $ \u223c F CL / a 2 $ (for the CLM) with the area a^{2} meaning that we require almost every particle to have a contact. It implies the thickening should roughly take place around $ \sigma \u2248 F CL / a 2 = 6 \pi \eta 0 \gamma \u0307 0 $, which gives $ \sigma / ( \eta 0 \gamma \u0307 0 ) \u2248 6 \pi $. This value falls at the upper end of the thickening regime as shown in Fig. 2, which is consistent with our assumption that at this stress scale every particle should have contacts with its neighbors, hence the system is in the frictional rheology state. The inverse argument applied to the onset stress σ_{ST} finds that the area A_{ST} over which one unit of critical load applies is σ_{ST} = F_{CL}/A_{ST}, with $ A ST \u2248 ( 7 a ) 2 $. That means that at the onset of shear thickening, the contact chains must be sparsely distributed.
Friction is an essential ingredient for the shear thickening to occur, as is illustrated in Fig. 3, which shows the effect of a reduction of the friction coefficient on the $ \eta r ( \gamma \u0307 ) $ curve in the CST regime. If the friction is removed (μ = 0), the effect is completely suppressed (which is expected, as the threshold force scale F_{CL} disappears from the equations of motion for μ = 0, and the simulation is then rigorously the same for all $ \gamma \u0307 $). If the friction is only reduced (μ = 0.2), the shear thickening is milder than with μ = 1 friction at the same volume fraction. In that case, the CST however becomes more pronounced and turns to DST when one increases the volume fraction, just as for the μ = 1 case.
C. Discontinuity and hysteresis
The existence of two distinct flowing states suggests that one might expect to observe a hysteresis associated with DST. It is thus remarkable that, aside from suspensions that are hysteretic at the microscopic level (with attractive interactions, for example, since once particles are aggregated under flow they stay aggregated upon cessation of the flow), there are very few experimental reports of hysteresis in systems showing DST [one of the best examples is provided by Bender and Wagner (1996)].
What is commonly observed, however, is a switching behavior, where the time series of the stress shows that the system shares its time between the two states, erratically switching in what looks like activated events [Boersma et al. (1991); D'Haene et al. (1993); Bender and Wagner (1996); Lootens et al. (2003); Lootens et al. (2004); Lootens et al. (2005); Hébraud and Lootens (2005)]. We observe this switching behavior in both models that we studied. In Fig. 5, we show time series of the stress and the corresponding histograms around the DST. These data show a typical “coexistence” behavior: Significantly below and above shear thickening, the system, respectively, stays in the low and high viscosity state, but close to DST, the two states coexist in the same time series, the proportion of time spent in the high viscosity state increasing with $ \gamma \u0307 $. While it is tempting to use the analogy of DST with an equilibrium firstorder transition to conclude that the switchings are finitesize “activated” events, we need to be careful in doing so, because intermittency seems to survive for unusually large numbers of particles. The intermittency is still present in simulations with N = 1000 particles, and the switching frequency does not seem to decrease relative to simulations with N = 512 particles, whereas in a firstorder thermodynamic transition, one would expect the switching rate to be drastically reduced by a doubling of the system size. In fact, this effect might persist for much larger particle numbers, as is suggested by the experimentally observed intermittency [Boersma et al. (1991); D'Haene et al. (1993); Bender and Wagner (1996); Lootens et al. (2003); Lootens et al. (2004); Lootens et al. (2005); Hébraud and Lootens (2005)] in systems where N is considerably larger. The simulations of Heussinger (2013) also show intermittency for systems of several thousands of particles. Moreover, should this phenomenon disappear in the thermodynamic limit, the fact that this limit is not observed until there are very large numbers of particles might be valuable information about the nature of this outofequilibrium phase transition, just as strong finitesize effects have a profound physical significance near a critical point.
In order to observe hysteresis, one must have a system where the relaxation time τ_{relax} (the length of the transient) is much smaller than the activation time τ_{act} (the typical time between two switches). One then does a proper measurement with a shear rate sweep on a time scale τ_{sw} such that $ \tau relax \u226a \tau sw \u226a \tau act $: The first inequality ensures that the system is always in a steady state during the sweep and that the measurement is done with sufficient time averaging, while the second inequality ensures that we are not averaging data over two distinct states. In our simulations, as can be seen by looking at the time series in Fig. 5, a clear separation of time scales is not achieved (we have $ \tau relax \u2272 \tau act $, which does not permit finding a proper τ_{sw}) and hysteresis cannot be observed unambiguously. Many experimental data actually suffer from the same limitations, and this might be the reason for the very small number of hysteretic flow curves actually reported. One exception is the noted work of Bender and Wagner (1996).
D. Normal stresses
One common rheological characteristic of dense suspensions is the appearance of finite normal stress differences. Our measurements of the two normal stress differences N_{1} and N_{2} are shown in Figs. 7 and 6 for the CLM. The ERM behaves very similarly and is not shown.
We obtain for N_{2} a behavior consistent with most experimental data available: It is negative, comparable to the shear stress (and much larger than N_{1} that we discuss in the next paragraph), and its behavior is reminiscent of that of the shear stress. This comes from the fact that most of the stress is transmitted by forces in the shear plane (flowgradient), not along the vorticity direction. Indeed the ratio N_{2}∕σ is rather insensitive to the stress magnitude, increasing at most by a factor of two across shear thickening while the stresses increase by almost two orders of magnitude. What is remarkable is the independence of N_{2}∕σ on the volume fraction: All volume fractions investigated here collapse on a master curve for both models.
There has been debate in the community regarding the value (and even the algebraic sign) of N_{1}. The first notable feature of our data concerning N_{1} is that it is dominated by the fluctuations. In one time series, a cursory look indicates that its value fluctuates around zero. Long time averaging reveals more structure, however, as shown in Fig. 7 for the CLM. Prior to shear thickening, the average value of N_{1} is nearly zero (slightly negative for CLM, slightly positive for ERM). The shear thickening transition is marked by two different behaviors, depending on the volume fraction: At the lowest volume fractions studied, N_{1} decreases across shear thickening, while at larger $\varphi $ there is a clear upturn toward positive values at the shear thickening transition. For the CLM, the behavior can be systematized by plotting N_{1} as a function of the stress, as in Fig. 7: N_{1}∕σ is an increasing function of σ, negative at low stresses, with a qualitative behavior roughly independent of $\varphi $. The same plot for the ERM is similar but less systematic. In any case, even above shear thickening, the ratio N_{1}∕σ is always small, never exceeding 0.1.
Overall, the behavior we observe for N_{1} is in reasonable agreement with most experiments [Lootens et al. (2005); Lee et al. (2006); Larsen et al. (2010); Couturier et al. (2011); Dbouk et al. (2013)]. Zarraga et al. (2000); Singh and Nott (2003); Dai et al. (2013) report a negative value for N_{1}. Lee et al. (2006) agree with our trend at low shear rates but observe a subsequent change toward negative N_{1} at very high shear rates, while Dbouk et al. (2013) find a positive but significantly larger N_{1}. In any case, the amplitude of the fluctuations seems to be the relevant physical information about N_{1}, as the average, positive, or negative, is buried in the fluctuations in the time series, such that at any time N_{1} can be either positive or negative, even at large $\varphi $ and $ \gamma \u0307 $.
E. Microstructure
NonNewtonian behavior of Brownian suspensions is often discussed with shear induced microstructure; the particle configuration is nearly at equilibrium at low Péclet number, while an anisotropic microstructure is induced by shear at high Péclet number [Morris (2009)]. Ideas associated with the orderdisorder transition scenario [Hoffman (1972, 1974, 1998)] also predicted a clear structural change with shear thickening. Experimental observations, however, revealed that the low viscosity state was not always ordered and that the shear thickening can occur with only a subtle signature in the microstructure [Bender and Wagner (1996); Watanabe et al. (1997); Watanabe et al. (1998); Newstein et al. (1999); Maranzano and Wagner (2002)]. One suggestion comes from the hydroclustering scenario, which assumes the existence of hydrodynamically created extended density fluctuations. Some experimental data are indeed in qualitative agreement with this scenario [Maranzano and Wagner (2002); Cheng et al. (2011)].
In the scenario we suggest in this article, the main structural modifications across the shear thickening transition should be sought in the contacts between particles. At low shear rates, frictional contacts are avoided because the particle pressure is too small to overcome the threshold (for the CLM) or the repulsion between particles (for the ERM). Those changes should thus be detected by measures specifically sensitive to the contacts. The pair correlation function should show few dramatic modifications across shear thickening. This is what we show in this section.
We define the pair correlation function
The system being bidisperse, we can also define partial pair correlation functions restricted to, e.g., pairs of smallsmall particles $ g SS ( r ) \u2261 ( V / N S 2 ) \u2211 i , j \u2208 S S \delta ( r \u2212 r i j ) $ where $ S S $ is the subset of N_{S} small particles. In the same way, we define the structure factor (which is related to the pair correlation function via a Fourier Transform) as
In Figs. 8 and 9, we show the pair correlation function restricted to the shear plane (velocity/gradient) at $ \varphi = 0.57 $ for four values of the shear rate: $ \gamma \u0307 / \gamma \u0307 0 = 0.005 $, well below DST; $ \gamma \u0307 / \gamma \u0307 0 = 0.015 $, just below DST; $ \gamma \u0307 / \gamma \u0307 0 = 0.02 $, just above DST and $ \gamma \u0307 / \gamma \u0307 0 = 0.05 $, well above DST. As expected, the evolution of those plots seems weak around DST. The main feature is a loss of contrast, unveiling a more isotropic structure above shear thickening.
The structure factor shown in Fig. 10 reveals another interesting feature: There is a clear anisotropy in the shear plane in the low viscosity state, with a peak in the gradient direction that is strongly reduced (but not absent, see the right of Fig. 10) in the high viscosity state, above DST. This anisotropy is absent in the flowvorticity plane; see Fig. 11. Looking back at the pair correlation functions in Fig. 8, we see that the peak of $ S all ( k ) $ along the gradient direction is the signature of a slight ordering of the low viscosity phase: There are small but visible downstream ripples aligned with the flow direction at small shear rate.
The slight layering we observe is however quite far from the string order usually associated with the orderdisorder scenario. In supplementary material, we show movies of the system projected on the gradientvorticity plane (i.e., looking down the stream direction) for shear rates just below and just above shear thickening, where we reduced the size of the particles to a third of their actual radius to allow a better visualization. Those movies show warped layers along the flowvorticity plane in some parts of the system, while in other parts ordered layers are completely absent. In any case, there is no further string ordering within the layers. This is consistent with our structure factor data in the gradientvorticity plane (not shown), which exhibits only some structure in the gradient direction and no sixfold symmetry patterns typically observed when string ordering takes place [Laun et al. (1992); Kulkarni and Morris (2009)].
IV. DISCUSSION
A. Contact network
In our simulation, the shear thickening transition is due to the appearance of a growing number of frictional contacts as the shear rate (or more precisely the shear stress) increases. The contact network that is built across shear thickening is shown in Fig. 12 for the CLM for the two cases of CST and DST, respectively. (Two movies showing $ \gamma \u0307 / \gamma \u0307 0 = 0.015 $ and 0.018 of DST are available in supplementary material. The time in these movies is scaled with $ 1 / \gamma \u0307 $.)
At low shear rate, in the low viscosity state, frictional contacts only seldomly appear, being concentrated in small force chains along the compressional axis. At high shear rate, frictional contacts are the norm, creating a frictional contact network that is very close to being jammed; i.e., having only a few (collective) degrees of freedom left to reorganize under the applied stress.
Between those two extreme situations, a whole continuum of gradually denser frictional contact networks is seen across the CST transition in Fig. 12 (top). But in the case of DST in Fig. 12 (bottom), the contact network discontinuously changes from sparse to dense at the transition, never showing configurations of intermediate densities. Another interesting point is the occurrence of intermittency in the contact network, which is the immediate structural origin of the intermittent stress behavior shown in Fig. 5: Close to DST, the contact network suddenly switches from mostly frictionless to mostly frictional, showing a large sensitivity to fluctuations. When this happens, the viscosity immediately follows by switching to high/low viscosity when the contact network, respectively, densifies/loosens.
We can make the link between frictional contacts and shear stress explicit by looking at the fraction of frictional contacts f [Wyart and Cates (2014)]. This quantity is unambiguously defined in the CLM model as the ratio of the number of contacts $ N C f $ that are in the frictional state (those for which the normal force exceeds the critical load F_{CL}) divided by the total number of contacts N_{C}. This ratio has a direct relation to the shear stress, as the function f(σ) plotted in Fig. 13 shows. This single relation demonstrates that, more than structural aspects, the shear thickening is above all a manifestation of the proliferation of frictional contacts in the system. We note one remarkable aspect of the relation f(σ): It is independent of the volume fraction, at least in the range of $\varphi $ that we have studied (which is rather large, owing to the fact that it should be understood as a range of distance to jamming $ \varphi \u2212 \varphi J \mu $ rather than a bare range of $\varphi $). This peculiar aspect will be studied in a separate article.
B. Phase diagram, relation to jamming
While the jamming transition is the critical phenomenon underlying DST, it should be noted that the jamming transition and DST are distinct, and in particular $ \varphi c \u2260 \varphi J \mu $. According to a recent theoretical argument by Wyart and Cates (2014), a scenario based on two diverging rheologies like the one we present in this work implies under reasonable assumptions that $ \varphi c < \varphi J \mu $; i.e., DST could then occur between two flowable (unjammed) states for volume fractions $ \varphi c < \varphi < \varphi J \mu $.
In our simulations, we estimate $ \varphi J \mu = 1 \u2248 0.58 $ (see Fig. 1), while we observe DST for $ \varphi = 0.56 $ (see Fig. 2), which indeed implies $ \varphi c < \varphi J \mu = 1 $. For a volume fraction $ \varphi > \varphi J \mu = 1 $ only the low viscosity state would be visible under shear, because the high viscosity frictional system is jammed at this volume fraction. This sets a maximum shear rate for the shear of the system at this volume fraction.
Physically, the two qualitatively different behaviors CST and DST, and the fact that $ \varphi c < \varphi J \mu = 1 $ can be explained in our model as stemming from the level of contrast that exists between frictionless and frictional rheologies, as suggested by Wyart and Cates (2014) [part of this explanation, when the contrast is diverging, was also suggested by Nakanishi et al. (2012)]. Recalling that the fraction of frictional contacts f(σ) essentially depends only on the applied shear stress (see the previous Sec. IV A), and noting that the viscosity interpolates between the two rheologies according to the number of frictional contacts, we see that we can write a direct relation η(σ) between the viscosity and the applied stress. If the difference between the two rheologies is small—i.e., $ \eta ( \sigma \u2192 \u221e ) \u2212 \eta ( \sigma \u2192 0 ) $ is small—the curve η(σ) is a sufficiently mildly increasing function such that $ \gamma \u0307 = \sigma / \eta ( \sigma ) $ is also an increasing function of σ; this corresponds to a singlevalued curve $ \sigma ( \gamma \u0307 ) $, and thus to CST. But if the difference $ \eta ( \sigma \u2192 \u221e ) \u2212 \eta ( \sigma \u2192 0 ) $ becomes too large, η(σ) increases faster than σ in some interval, which means that $ \gamma \u0307 ( \sigma ) $ is a decreasing function in this same interval. This corresponds to a multivalued curve $ \sigma ( \gamma \u0307 ) $, which is unstable, and shows up experimentally as a discontinuous curve, i.e., DST.
An increasing difference $ \eta ( \sigma \u2192 \u221e ) \u2212 \eta ( \sigma \u2192 0 ) $ is naturally provided in our model by the fact that the frictional rheology diverges at a jamming volume fraction $ \varphi J \mu $ that is smaller than the frictionless rheology jamming point $ \varphi J 0 $. Then, at low volume fraction, the difference is small, but there must be a point $ \varphi c $ below $ \varphi J \mu $ where the difference becomes large enough to observe a DST.
These ideas can be summed up in a phase diagram. The rheology is essentially frictionless in the lower part of the diagram Fig. 14, as the stress is too small to activate friction between grains. This rheology diverges at the frictionless jamming point $ \varphi J 0 $. The rheology is frictional in the upper part of the diagram, as friction is activated under the applied stress. This rheology diverges at the frictional jamming point $ \varphi J \mu $ and thus shows a larger viscosity than the frictionless rheology. Thus two rheologies coexist on the shear ratevolume fraction plane. Those rheologies are separated by a shear thickening that is continuous for $ \varphi < \varphi c $ and discontinuous for $ \varphi c < \varphi < \varphi J \mu $.
The discontinuity is actually related to the coexistence of the two rheologies in the triangle delimited by a dashed line. In this region, we observe one consequence of the coexistence in the intermittency of the flow: The system switches from one state to the other through activation events, as is shown by the time series of the stress (see Fig. 5). For $ \varphi J \mu < \varphi < \varphi J 0 $, DST is strictly speaking no longer observed, as the high viscosity flowable state does not exist any more. In this region, DST is actually replaced by shear jamming [Bi et al. (2011)]: If one applies a shear stress larger than σ_{ST}, the system goes to a solid state and cannot flow. This implies the existence of a forbidden region, in gray, where no flow is possible for hard spheres. In the lower shearrate domain above $ \varphi J \mu $, the low viscosity state is stable, but the activated events responsible for the intermittency at lower $\varphi $ in small systems might lead to complete jamming in a finite time. Lastly, the dot in the upper left corner of the DST region is a critical point separating CST from DST. It shares some features with a critical point of a second order phase transition, for instance, a diverging susceptibility $ d \sigma / d \gamma \u0307 $.
Of course, as we numerically work with a shearrate controlled scheme, the system always flows at any shear rate and any volume fraction, even in the forbidden region of the phase diagram. But in the latter region flow is only achieved by creating large overlaps between the particles (i.e., compressing them), hence violating the criteria we set to mimic hard sphere suspensions. At high shear rate and for $ \varphi > \varphi J \mu $, a real hard sphere system would not flow, but it does in the simulation because of our inability to enforce the hard sphere condition in a high (possibly infinite) stress state. The spheres we simulate in that case cannot be considered hard any more, and their stiffness sets a cutoff to the stress scales under shear. This situation is somehow similar to the one observed in experiments, where it is observed that the measured stress in the high viscosity state is set by the weakest stress scale in the sample (which might be a large stress scale for the rheometer usage range), whether it is the particles' stiffness or the surface tension on the edge of the rheometer [Brown and Jaeger (2012)]. Therefore, in an experiment, even if the idealized system would be jammed under the investigated conditions, the real system still flows, perhaps thanks to dilation, whereas in our simulation, under the same conditions, it still flows thanks to the deformability of particles.
The difference between $ \varphi c $ and $ \varphi J \mu $ may have been observed in some stress drop experiments [O'Brien and Mackay (2000); Larsen et al. (2010)]. For $ \varphi c < \varphi < \varphi J \mu $, the phase diagram of Fig. 14 predicts that the stress in the thickened state would relax quickly upon flow cessation, as the contact network is unjammed. By contrast, for $ \varphi J \mu < \varphi < \varphi J 0 $ the stress would not relax entirely, as part of it would be stored elastically in the jammed contact network. This scenario is actually consistent with the data from [O'Brien and Mackay (2000); Larsen et al. (2010)].
V. CONCLUSION
We have introduced a frictionalviscous model of dense suspensions under shear that is built on the framework used by previous standard models of suspensions: Simulations are performed in the zero particleinertia and zero fluidinertia limits ( $ S t \u2192 0 $ and $ R e \u2192 0 $) and include relevant hydrodynamic interactions and a shortrange repulsive potential. The critical innovation is that, besides this purely hydrodynamic basis, it includes frictional contacts between solid particles. This model can be used to simulate sheared suspensions over a range of volume fractions, from mildly dense, where the shortrange hydrodynamic interactions dominate, to the jamming point, where contact forces dominate.
The effect of adding frictional contacts is striking at large volume fractions: Tangential forces due to friction restrict microscopic rearrangements in the system, resulting in a large shear stress. The number of contacts created during the flow is the result of a competition between applied shear forces and the shortrange repulsive forces. Therefore, alongside the contactless (hence frictionless) rheology at low shear rates a new contact dominated frictional rheology appears at high shear rates. Shear thickening is then simply the transition from the contactless rheology to the contact dominated rheology with increasing shear stress (and shear rate). The strongest evidence that this stressinduced friction scenario is correct is the ability of our model to reproduce shear thickening, including its discontinuous form, and its volumefraction dependence. Therefore, we conclude that friction is a key element for shear thickening.
The shear thickening mechanism we describe is qualitatively independent of our parameter choices. Quantitatively, we can summarize the effects of the various parameters as follows. The friction coefficient μ only affects the volume fraction range over which the phenomenon happens: If μ increases, the critical volume fraction $ \varphi c $ decreases. Any friction coefficient μ > 0 leads to a DST, with $ \varphi c ( \mu ) \u2192 \varphi J 0 $ in the limit $ \mu \u2192 0 $. The lubrication cutoff δ mainly sets the dissipation level in the systems and thus sets an overall prefactor in the entire curve η(σ): The viscosity is larger when δ decreases. For instance, we performed simulations with a cutoff $ \delta = 10 \u2212 2 $ and compared the results to the ones shown in the present work obtained with $ \delta = 10 \u2212 3 $: The curve $ \eta \delta = 10 \u2212 2 ( \sigma ) $ can be superimposed to $ \eta \delta = 10 \u2212 3 ( \sigma ) $ by a mere multiplicative factor $ \eta \delta = 10 \u2212 2 ( \sigma ) \u2248 0.8 \eta \delta = 10 \u2212 3 ( \sigma ) $. This has the effect of pushing the shear thickening to lower shear rates, as the onset stress σ_{ST} is independent of δ. The particle stiffness constants k_{n} and k_{t} are chosen such that the simulations are run in the asymptotic regime where the dimensionless numbers $ 6 \pi \eta 0 a \gamma \u0307 / k n , t $ are very small. In this asymptotic regime, the rheology is independent of these stiffness constants. The Debye length κ^{–1} governs the shearthinning behavior at low shear rates. Our CLM can be thought as having a vanishing Debye length. Comparing with the ERM with κ^{– 1 }= 0.05a (with a the radius of the small particles), we observe that besides the shearthinning regime, a finite κ^{–1} has virtually no effect on the shear thickening and the high viscosity state.
We now turn to a quantitative comparison with experimental data. The main tuning parameter in our description is the force scale F_{CL} for CLM (or F_{ER} for the ERM). This parameter sets the onset stress at which shear thickening occurs. We can estimate this force for a system of 1 μm silica spheres in water. We find an electrostatic repulsion at contact of order $ F ER \u2248 0.2 \u2009 nN $ [as evaluated in Behrens and Grier (2001) with a charge density $ \u2248 5 \xd7 10 \u2212 4 \u2009 C \u2009 m \u2212 2 $ and a Debye length $ \kappa \u2212 1 \u2248 0.1 \u2009 \mu m $]. This gives for our model the stress scale $ \eta 0 \gamma \u0307 0 = ( 6 \pi a 2 ) \u2212 1 F ER $, which together with the nondimensionalized onset stress $ \sigma ST / ( \eta 0 \gamma \u0307 0 ) \u2248 0.3 $ obtained the ERM model (see Fig. 4) gives a predicted onset stress of $ \sigma ST \u2248 8 \u2009 Pa $. For the very same system, Lootens et al. (2005) experimentally find an onset stress of $ \sigma ST \u2248 5 \u2009 Pa $, with which our simulation is thus in good agreement. Thus, our mechanism seems to give a quantitative description of shear thickening in dense suspensions, even if it would require comparison to more experimental data to be definitely conclusive. (Obtaining a reliable value of the microscopic force scale is rarely feasible in the published data on DST, and this restricts the number of comparisons we can make at the moment.)
Although our model assumes nonBrownian suspensions, we may expect that the same mechanism explains shear thickening of Brownian colloidal suspensions. Indeed, the Brownian forces qualitatively play a role similar to the repulsive potential that we use in our ERM: At low shear stress, the Brownian forces are large enough to randomly open gaps between contacting particles, while at high shear stress this becomes much less likely, as relative trajectories of pairs of particles become smooth at high Peclet numbers [Nazockdast and Morris (2013)]. The fact that a random noise can relax friction between grains is a known phenomenon in granular matter [Gao et al. (2009); Karim and Corwin (2014)]. Friction would then develop in a Brownian suspension in a manner very similar to the present work.
ACKNOWLEDGMENTS
The authors would like to thank Mike Cates and Matthieu Wyart for the many discussions concerning our respective work and Philippe Claudin for discussions about the μ(I_{v}) rheology. The authors' code makes use of the CHOLMOD library by Tim Davis (https://www.cise.ufl.edu/research/sparse/cholmod/) for direct Cholesky factorization of the sparse resistance matrix. This research was supported in part by a grant of computer time from the City University of New York High Performance Computing Center under NSF Grant Nos. CNS0855217, CNS0958379, and ACI1126113. J.F.M. was supported in part by NSF PREM (DMR 0934206).
APPENDIX A: HYDRODYNAMIC RESISTANCES
The hydrodynamic interactions arising from an imposed background flow and relative motions between nearby particles are given by
where a diagonal matrix $ R Stokes $ comes from the (onebody) Stokes drag and sparse matrices $ R Lub $ and $ R \u2032 Lub $ come from the (twobody) lubrication.
Using the basic units $ L 0 \u2261 a $ for lengths, $ U 0 \u2261 L 0 \gamma \u0307 $ for velocities, and $ F 0 \u2261 6 \pi \eta 0 L 0 U 0 $ for forces, the elements of $ R Stokes $ give the Stokes drag through
while $ R Lub $ and $ R \u2032 Lub $ consist of offdiagonal blocks giving the lubrication forces and torques for a pair (i, j) through
Following the notation of Jeffrey and Onishi (1984) and Jeffrey (1992), the matrices $ R Lub ( i , j ) $ and $ R Lub \u2032 ( i , j ) $ are obtained from the particle positions as (we give only the upper triangular part of the symmetric $ R Lub ( i , j ) $)
In these expressions, we have introduced the normal projection operator $ P n i j \u2261 n i j n i j $, the tangential projection operator $ P n i j \u2032 \u2261 I \u2212 n i j n i j $, and the “cross product” operator $ P n i j r $ defined as $ P n i j r q \u2261 n i j \xd7 q $. We also used the operators $ Q n i j , \u2009 Q n i j \u2032 $, and $ Q n i j r $ defined for an arbitrary matrix $M$ as
The scalar coefficients X and Y have an explicit dependence on the nondimensional interparticle gap $ h ( i , j ) \u2261 2 ( r ( i , j ) \u2212 a i \u2212 a j ) / ( a i + a j ) $, and we use only the terms of leading order
With $ \lambda \u2261 a j / a i $, the coefficients g^{X} and g^{Y} appearing in $ R Lub ( i , j ) $ are written
Similarly, the terms appearing in $ R Lub \u2032 ( i , j ) $ are
APPENDIX B: CONTACT MODEL
1. Tuning the spring constants
Our objective is to study the rheology of hardsphere suspensions. We have no way to mimic rigorously hard spheres by using a contact model with springs. To provide the best approximation to hard spheres, the elastic constants appearing in the contact model must satisfy a simple constraint: They should be large enough to generate as little geometric overlap as possible between the particles. (An equivalent way to state this is through the requirement that the contact based nondimensionalized shear rate $ 6 \pi \eta 0 a \gamma \u0307 / k n \u226a 1 $ introduced in Sec. II B 1 is sufficiently small.)
We therefore set a criterion for the overlap: The largest overlap between any two particles during the simulation should not be larger than h_{max} (5% of the particle radius in this simulation). As the overlap $  h  $ depends on the shear stress, an estimate being $ \u27e8  h  \u27e9 \u223c ( 1 / k n ) \sigma ( \varphi , \gamma \u0307 ) a 2 $, the spring constant has to be tuned for every $\varphi $ and $ \gamma \u0307 $ so that $ max (  h  ) < h max $. We introduce a similar criterion for the tangential spring stretch mimicking static friction: We pick $ k t ( \varphi , \gamma \u0307 ) $ so that $\xi $ is smaller than 5% of the particle radius. We detail below how we fulfill these two conditions and retain hardsphere like behavior in our simulation.
For a given $\varphi $, the largest stress is obtained in the high shearrate limit. Hence, we start by determining high shearrate spring constants $ ( k n \u2217 ( \varphi ) , k t \u2217 ( \varphi ) ) $ by running preliminary simulations at $ \gamma \u0307 \u2192 \u221e $, where these values are regularly updated with a certain interval until the criteria on the overlap and the tangential spring stretch are fulfilled.
Second, a trivial shearrate dependence comes in if the contact model introduces a force scale in addition to the hydrodynamic one. Avoiding this requires picking the shear rate dependence of $ k n ( \varphi , \gamma \u0307 ) $ and $ k n ( \varphi , \gamma \u0307 ) $ by scaling these parameters with the shear rate; i.e., $ k n ( \varphi , \gamma \u0307 ) = \gamma \u0307 k n \u2217 ( \varphi ) $ and $ k t ( \varphi , \gamma \u0307 ) = \gamma \u0307 k t \u2217 ( \varphi ) $. With this scaling, there is no competition between hydrodynamic and contact forces, as they are proportional, so an additional explicit force scale (as the one in Sec. II B) is required to have a shearrate dependence in our simulation. Also note that as the high shearrate limit is always the largest viscosity at a given $\varphi $, with this choice of scaling $ k n ( \varphi , \gamma \u0307 ) $ and $ k t ( \varphi , \gamma \u0307 ) $ always fulfill the criteria we set for the overlap and tangential spring stretch.
2. Combining the overdamped dynamics and the contact model
In this appendix, we present a slightly more general contact model than the one presented in the main text. While the contact model may not deserve a lengthy description in itself, immersing an orthodox contact model in overdamped dynamics leads to nontrivial difficulties, which justifies our use of a simpler version of the model.
A more general stick/slide friction model using springs and dashpots is the following [Cundall and Strack (1979); Luding (2008)]:
These forces must fulfill Coulomb's law of friction $  F C , \u2009 tan \u2009 ( i , j )  \u2264 \mu  F C , nor ( i , j )  $. In the above expression, k_{n} and k_{t} are, respectively, the normal and tangential spring constants, and γ_{n} and $ \gamma t $ are the damping constants. The normal and tangential relative velocities between two particles i and j are $ U n ( i , j ) \u2261 P n i j ( U ( j ) \u2212 U ( i ) ) $ and $ U t ( i , j ) \u2261 P n i j \u2032 [ U ( j ) \u2212 U ( i ) \u2212 ( a i \Omega ( i ) + a j \Omega ( j ) ) \xd7 n i j ] $. Finally, the quantity $ \xi ( i , j ) $ is the tangential spring stretch.
The computation of the tangential spring stretch $ \xi ( i , j ) $, described in the following, requires some care as we have to impose Coulomb's law, which is made difficult by the overdamped dynamics. At the time t_{0} at which the contact (i, j) is created, we set an unstretched tangential spring $ \xi ( i , j ) ( t 0 ) = 0 $. At any further time step t in the simulation, the tangential stretch $ \xi ( i , j ) ( t ) $ is incremented according to the value of a test force, $ F C , \u2009 tan \u2009 \u2032 ( i , j ) ( t + d t ) = k t \xi \u2032 ( i , j ) ( t + d t ) + \gamma t U t ( i , j ) ( t + d t ) $ with the tentative update of stretch $ \xi \u2032 ( i , j ) ( t + d t ) = \xi ( i , j ) ( t ) + U t ( i , j ) ( t ) d t $. If $  F C , \u2009 tan \u2009 \u2032 ( i , j ) ( t + d t )  \u2264 \mu  F C , nor ( i , j ) ( t + d t )  $, the contact is in a static friction state and we update the spring stretch as $ \xi ( i , j ) ( t + d t ) = \xi \u2032 ( i , j ) ( t + d t ) $ and the corresponding tangential force is $ F C , \u2009 tan \u2009 ( i , j ) ( t + d t ) = F C , \u2009 tan \u2009 \u2032 ( i , j ) ( t + d t ) $. However, if $  F C , \u2009 tan \u2009 \u2032 ( i , j ) ( t + d t )  > \mu  F C , nor ( i , j ) ( t + d t )  $, the contact is in sliding state and the spring is updated as $ \xi ( i , j ) ( t + d t ) = ( 1 / k t ) ( \mu  F C , nor ( i , j ) ( t + d t )  t i j \u2212 \gamma t U t ( i , j ) ( t ) ) $, where the direction is the same as the test force, i.e., $ t i j \u2261 F C , \u2009 tan \u2009 \u2032 ( i , j ) ( t + d t ) /  F C , \u2009 tan \u2009 \u2032 ( i , j ) ( t + d t )  $. Assuming that the velocities are continuous in time, this ensures that the contact force will at most weakly violate the Coulomb friction law for the next time step, as the force effectively used in the equation of motion will be $ F C , \u2009 tan \u2009 ( i , j ) ( t + d t ) = k t \xi ( i , j ) ( t + d t ) + \gamma t U t ( i , j ) ( t + d t ) = \mu  F C , nor ( i , j ) ( t + d t )  t i j + \gamma t U . t ( i , j ) d t $, so that $  F C , \u2009 tan \u2009 ( i , j ) ( t + d t )  = \mu  F C , nor ( i , j ) ( t + d t )  + O ( d t ) $.
Of course, if velocities happen to be discontinuous, the Coulomb's law might be violated significantly, and this is the source of one problem when merging this contact model with the hydrodynamic interaction model. Indeed, when contacts are created and destroyed during the flow (which happens very frequently), the overall resistance matrix changes discontinuously, as dashpots are switched on and off. But, as other positiondependent forces are continuous in time, solving the force balance equation will lead to a discontinuity in velocities when a contact forms or breaks. This, in turn, leads to large violations of Coulomb's law, and hence to numerical instabilities where contacts keep switching between static and dynamic cases at every time step.
Another problem occurs when a contact forms (i.e., $ \xi ( i , j ) ( t ) = 0 $): As by definition the normal load on the new contact $  F C , nor ( i , j ) ( t + d t )  $ vanishes (or is very small in the simulation owing to the finite time step), the finite test force $ F C , \u2009 tan \u2009 \u2032 ( i , j ) ( t + d t ) \u2243 \gamma t U t ( i , j ) ( t ) $ will lead to an immediate rescaling of the tangential spring as $ \xi ( i , j ) ( t + d t ) \u2243 \u2212 ( \gamma t / k t ) U t ( i , j ) ( t ) $; i.e., an entirely new finite force appears right at contact time.
These two problems actually also exist in simulations of dry granular matter [Walton (1993b)], but they are more acute with an overdamped dynamics because of the direct relation between forces and velocities. In this case, any discontinuity in forces or velocities rapidly causes numerical instabilities. A solution is to use a continuously varying damping $ \gamma t ( h i j ) $ such that $ \gamma t ( h i j = 0 ) = 0 $, which ensures the continuity of forces and velocities [Walton (1993a)] at the cost of increased complexity of the model. The solution that we prefer using here is to eliminate of the tangential dashpot by setting $ \gamma t = 0 $. In any case, this dashpot has no physical significance, and it is only used in granular matter simulations as an efficient numerical stabilizer. In our already overdamped context, as long as the hydrodynamic resistance associated with tangential motion is not too small, we do not need such an extra resistive stabilizer, and we do not face any major difficulty by simply dropping this dashpot.
We can quantify this assertion by looking at the relaxation times associated with a contact. Both normal and tangential contacts have relaxation times. For the normal part, it is the one of a spring and dashpot system, $ \tau n = \gamma n / k n $. For the tangential part, the damping is provided by the hydrodynamic resistance, which we here simply denote $ \gamma t H $, so that $ \tau t = \gamma t H / k t $. On the one hand, in order to have a stable numerical scheme, one should choose these relaxation times large enough compared to the time step $ \tau n , \u2009 \tau t \u226b d t $. On the other hand, contacts of hard spheres should react instantaneously to any external load, so physics imposes relaxation times smaller than other typical times $ \tau n , \u2009 \tau t \u226a 1 / \gamma \u0307 $. For the normal part of the contact, we are free to choose the damping γ_{n} to achieve this. For the tangential part however, we check that those scale separations are verified.