Shear banding, discontinuous shear thickening, and rheological phase transitions in athermally sheared frictionless disks

TL;DR

This study uses numerical simulations to explore rheological phase transitions, shear banding, and discontinuous shear thickening in athermal, frictionless disk systems under shear, revealing a first-order transition and the impact of rotational coupling.

Contribution

It uncovers a first-order rheological phase transition in athermal disks and highlights the role of rotational-translational coupling in shear banding and thickening behaviors.

Findings

Identifies a sharp Bagnoldian to Newtonian transition in the phase diagram.

Shows the system is always Newtonian at jamming.

Demonstrates the importance of rotational coupling at small Q.

Abstract

We report on numerical simulations of simple models of athermal, bidisperse, soft-core, massive disks in two dimensions, as a function of packing fraction $\phi$, inelasticity of collisions as measured by a parameter $Q$, and applied uniform shear strain rate $\dot\gamma$. Our particles have contact interactions consisting of normally directed elastic repulsion and viscous dissipation, as well as tangentially directed viscous dissipation, but no inter-particle Coulombic friction. Mapping the phase diagram in the $(\phi,Q)$ plane for small $\dot\gamma$, we find a sharp first-order rheological phase transition from a region with Bagnoldian rheology to a region with Newtonian rheology, and show that the system is always Newtonian at jamming. We consider the rotational motion of particles and demonstrate the crucial importance that the coupling between rotational and translational degrees of freedom has on the phase structure at small $Q$ (strongly inelastic collisions). At small $Q$ we show that, upon increasing $\dot\gamma$, the sharp Bagnoldian-to-Newtonian transition becomes a coexistence region of finite width in the $(\phi,\dot\gamma)$ plane, with coexisting Bagnoldian and Newtonian shear bands. Crossing this coexistence region by increasing $\dot\gamma$ at fixed $\phi$, we find that discontinuous shear thickening can result if $\dot\gamma$ is varied too rapidly for the system to relax to the shear-banded steady state corresponding to the instantaneous value of $\dot\gamma$.

Figures

Equations

$${\cal V}_{ij}(r_{ij})=\left\{\begin{array}[]{cc}\frac{1}{\alpha}k_{e}\left(1-r_{ij}/d_{ij}\right)^{\alpha},&r_{ij}<d_{ij}\\ 0,&r_{ij}\geq d_{ij}.\end{array}\right.$$

$${\mathbf{f}}_{ij}^{\mathrm{el}}=-\frac{d{\mathcal{V}}_{ij}(r_{ij})}{d{\mathbf{r}}_i}=\frac{k_e}{d_{ij}}{\left(1-\frac{r_{ij}}{d_{ij}}\right)}^{\alpha -1}{\hat{\mathbf{r}}}_{ij},$$

$$\mathrm{model}\mathrm{CD}:{\mathbf{f}}_{ij}^{\mathrm{dis}}=-k_d({\mathbf{v}}_i-{\mathbf{v}}_j),$$

$$m_i\frac{d^2{\mathbf{r}}_i}{d{t}^2}=\sum _j\left[{\mathbf{f}}_{ij}^{\mathrm{el}}+{\mathbf{f}}_{ij}^{\mathrm{dis}}\right],$$

$${\mathrm{model}\mathrm{CD}}_{\mathrm{rot}}:{\mathbf{f}}_{ij}^{\mathrm{dis}}=-k_d({\mathbf{v}}_i^{C_j}-{\mathbf{v}}_j^{C_i}),$$

$${\mathbf{v}}_i^{C_j}={\mathbf{v}}_i+\frac{d{\theta }_i}{dt}\hat{\mathbf{z}}\times {\mathbf{s}}_{ji},{\mathbf{s}}_{ji}\equiv {\mathbf{r}}_{ji}\left(\frac{d_i}{d_i+d_j}\right).$$

$$I_i\frac{d^2{\theta }_i}{d{t}^2}=\hat{\mathbf{z}}\cdot \sum _j{\mathbf{s}}_{ji}\times {\mathbf{f}}_{ij}^{\mathrm{dis}},$$

$${\tau }_e\equiv \sqrt{m{d}_s^2/k_e},{\tau }_d\equiv m/k_d.$$

$$Q\equiv {\tau }_d/{\tau }_e=\sqrt{m{k}_e/(k_d{d}_s{)}^2}$$

$${\tau }_0\equiv {\tau }_e^2/{\tau }_d=d_s^2{k}_d/k_e,$$

$$\varphi =\frac{\pi N}{2L^2}\left[{\left(\frac{d_s}{2}\right)}^2+{\left(\frac{d_b}{2}\right)}^2\right].$$

$${\mathbf{p}}^{\mathrm{el}}=\frac{1}{L^2}\sum _{i<j}{\mathbf{f}}_{ij}^{\mathrm{el}}\otimes {\mathbf{r}}_{ij},$$

$${\mathbf{p}}^{\mathrm{dis}}=\frac{1}{L^2}\sum _{i<j}{\mathbf{f}}_{ij}^{\mathrm{dis}}\otimes {\mathbf{r}}_{ij},$$

$${\mathbf{p}}^{\mathrm{kin}}=\frac{1}{L^2}\sum _i{m}_i\delta {\mathbf{v}}_i\otimes \delta {\mathbf{v}}_i,$$

$$\mathbf{p}={\mathbf{p}}^{\mathrm{el}}+{\mathbf{p}}^{\mathrm{dis}}+{\mathbf{p}}^{\mathrm{kin}}.$$

$$p=\frac{1}{2}\left[⟨p_{xx}⟩+⟨p_{yy}⟩\right],\sigma =-⟨p_{xy}⟩,$$

$$\mathbf{P}\equiv ({\tau }_e^2/m)\mathbf{p},$$

$$P\equiv ({\tau }_e^2/m)p,\Sigma \equiv ({\tau }_e^2/m)\sigma ,$$

$$B_p\equiv \frac{P}{(\dot{\gamma }{\tau }_e{)}^2}=\frac{p}{m{\dot{\gamma }}^2},B_{\sigma }\equiv \frac{\Sigma }{(\dot{\gamma }{\tau }_e{)}^2}=\frac{\sigma }{m{\dot{\gamma }}^2}.$$

$${\eta }_p\equiv \frac{P}{\dot{\gamma }{\tau }_0}=\frac{1}{k_d}\frac{p}{\dot{\gamma }},{\eta }_{\sigma }\equiv \frac{\Sigma }{\dot{\gamma }{\tau }_0}=\frac{1}{k_d}\frac{\sigma }{\dot{\gamma }}.$$

$${\tilde{\mathbf{v}}}_{i1}$$

$${\tilde{\mathbf{r}}}_{i1}$$

$${\mathbf{v}}_{i1}$$

$${\mathbf{r}}_{i1}$$

$$Z=2{\tau }_{\mathrm{dur}}{\nu }_{\mathrm{coll}}.$$

$${\mathbf{r}}_{ij}\equiv {\mathbf{r}}_i-{\mathbf{r}}_j,{\mathbf{v}}_{ij}\equiv {\mathbf{v}}_i-{\mathbf{v}}_j,$$

$${\hat{\mathbf{v}}}_{ij}\cdot {\hat{\mathbf{r}}}_{ij}=\cos \vartheta ,-180^{\circ }<\vartheta \le 180^{\circ }.$$

$$\mathrm{var}[{\dot{\theta }}_i/\dot{\gamma }]\equiv ⟨({\dot{\theta }}_i/\dot{\gamma }{)}^2⟩-{⟨{\dot{\theta }}_i/\dot{\gamma }⟩}^2,$$

$$⟨({\dot{\theta }}_i/\dot{\gamma }{)}^2⟩\equiv ⟨\frac{1}{N}\sum _i({\dot{\theta }}_i/\dot{\gamma }{)}^2⟩.$$

$${\mathbf{V}}_{ij}^C={\mathbf{V}}_{ij}+{\mathbf{V}}_{ij}^R,$$

$${\mathbf{V}}_{ij}^R\equiv \left[\frac{{\dot{\theta }}_i}{\dot{\gamma }}\frac{d_i}{d_s}+\frac{{\dot{\theta }}_j}{\dot{\gamma }}\frac{d_j}{d_s}\right]\hat{\mathbf{z}}\times {\hat{\mathbf{r}}}_{ji}$$

$$V_{ij}^R\equiv {\mathbf{V}}_{ij}^R\cdot {\hat{\mathbf{t}}}_{ji}$$

Full text

I Introduction

Simple models of spherical particles, interacting via soft- or hard-core repulsive contact interactions, have been used to model a wide variety of granular and soft-matter materials, such as dry granular particles, foams, emulsions, non-Brownian suspensions, and colloids. The rheology of such systems under applied shear strain has been a topic of much active investigation. For athermal systems ($T=0$), such models generally give a rheology with the following features: (i) At low density and low strain rate the system is in a flowing liquid-like state where the pressure $p$ and shear stress $\sigma$ scale with the shear strain rate $\dot{\gamma}$ according to either a Newtonian ($p,\sigma\sim\dot{\gamma}$) or Bagnoldian ($p,\sigma\sim\dot{\gamma}^{2}$) rheology, depending on the form of energy dissipation in the system VOT_rheology . Newtonian rheology is often associated with overdamped motion due to the Stokes drag on particles that are embedded in a host fluid OT_2007 . Bagnoldian rheology is often associated with dry hard granular particles, assuming that the frequency of collisions and the momentum change per collision are both proportional to the strain rate $\dot{\gamma}$ LLC . However Newtonian rheology may still be observed in models which neglect interactions with a host fluid VOT_rheology ; VOT_RDCD , and Bagnoldian rheology may still be observed even in the presence of a host fluid, as the strain rate increases Bagnold ; Fall . (ii) As the particle packing fraction $\phi$ increases, a shear driven jamming transition occurs at a well defined $\phi_{J}$, where the macroscopic transport coefficient of the liquid diverges OT_2007 ; VOT_RDCD ; Barrat ; OT_RD ; Lerner ; DeGiuli ; Ikeda ; VOT_CDn ; VOT_CDnQ . (iii) For soft-core particles, above $\phi_{J}$ the system develops a finite yield stress $\sigma_{Y}(\phi)$. For $\sigma<\sigma_{Y}$, the system is in a static jammed solid state, while for $\sigma>\sigma_{Y}$ the system is in a state of plastic flow, generally described by a Herschel-Bulkley rheology, $\sigma=\sigma_{Y}+C\dot{\gamma}^{b}$ OT_HB . When particles are frictionless, the jamming transition is continuous, with $\sigma_{Y}$ vanishing continuously as $\phi\to\phi_{J}$ from above. In this case, critical scaling is found to well describe behavior near the jamming transition OT_2007 ; Barrat ; OT_RD ; Lerner ; VOT_RDCD ; DeGiuli ; Ikeda ; VOT_CDn . When particles have inter-granular friction, the jamming transition may become discontinuous, with $\sigma_{Y}$ dropping discontinuously to zero at $\phi_{J}$; the rheology may now show hysteretic and re-entrant behavior near $\phi_{J}$ OH_friction ; Bi_Nature ; Heussinger ; Ciamarra ; Grob ; Grob2 ; Schwarz ; Jaeger_Nature .

However in many materials, the rheology is observed to be more complex. One such departure from the above scenario is the phenomenon of “shear banding,” in which the shear strain rate $\dot{\gamma}$ becomes spatially inhomogeneous, sometimes breaking into well-defined distinguishable bands with different values of $\dot{\gamma}$. Shear banding has been found in a wide range of materials, including complex fluids CatesFielding , non-Brownian suspensions Coussot ; Ovarlez0 , pastes, gels, and emulsions Ovarlez , foams Kabla , granular matter Jaeger_RMP ; Mueth , soils and rocks Rudnicki , metallic glasses Spaepen ; Li , and colloidal glasses Besseling ; Chikkadi . In some cases, shear banding is believed to result from inhomogeneous stress fields that vary spatially from values above to values below the yield stress of the material Schall . In other cases, shear banding is believed to result from the coupling of the system to a boundary wall Latzel ; Wang . However, in yet other cases, shear banding appears to exist as a bulk phenomenon, and in situations of spatially uniform stress Goddard ; Moller . In such cases, shear banding is believed to result from some inherent instability in the shear flow that leads to a non-monotonic stress-strain rate curve; shear bands then result as coexisting phases at a discontinuous transition, as in analogy to the van der Waals picture for the liquid-gas transition. In some of these cases, one of the bands is seemingly stationary, while the flowing band exists only above some critical strain rate $\dot{\gamma}_{c}$ Coussot ; Ovarlez0 ; Moller ; Dennin . In other cases, however, bands form which have different finite values of $\dot{\gamma}$ CatesFielding ; Besseling ; Chikkadi . Various theoretical models have been proposed Fielding ; Dhont ; Olmsted for shear banding, and specific microscopic mechanisms have been suggested for specific materials CatesFielding ; Besseling ; Dogic ; Bonn1 ; Bonn2 ; Bonn3 . For more complete reviews see Refs. Ovarlez and Schall . It remains of interest to see how general the phenomenon of shear banding might be, and whether relatively simple mechanisms can produce shear banding in simple homogeneous systems of repulsive spheres.

Another phenomenon that has received considerable recent attention is discontinuous shear thickening (DST). In dense non-Brownian suspensions, the canonical example being cornstarch in water, increasing the shear strain rate can lead to a sudden and dramatic increase of shear stress Brown . There seems to be a growing consensus that inter-particle frictional interactions are key to understanding DST. Several recent works have sought to explain DST as a transition from a state where hydrodynamic lubrication forces prevent particles from prolonged contact, to a state where such lubrication forces are overcome, particle contacts increase, and elastic Coulomb frictional forces at those contacts become important Fernandez ; Seto ; Mari ; Cohen . Other works OH_friction ; Heussinger ; Grob ; Grob2 ; Mari ; WyartCates ; Zippelius ; Hayakawa1 have described DST as a precursor of the discontinuous jump in stress associated with the shear-driven jamming of frictional particles OH_friction ; Bi_Nature . While both DST and shear banding have been associated with discontinuous transitions, it remains unclear if there is a direct connection between these two phenomena.

In this work we report on numerical simulations of simple models of athermal, bidisperse, soft-core disks in two dimensions. We consider massive particles and energy dissipation occurs via particle collisions only. While we take our particles to have no inter-granular elastic frictional force, we do include a tangential viscous dissipative force. This tangential dissipation is key to the novel behaviors we find. In one model, denoted as model “CD”, this dissipative force is taken proportional to the center of mass velocity difference between two particles in contact. In a second model, denoted as “CDrot”, it is taken proportional to the local velocity difference between the two particles at their point of contact. Model CDrot thus introduces a non-trivial coupling between the translational and rotational motions of the particles. Exploring the behavior of these models at low shear strain rates $\dot{\gamma}$, we map out the phase diagram as a function of packing fraction $\phi$, and a parameter $Q$ that measures the elasticity of collisions (see Figs. 3 and 17). We find for both models that the $(\phi,Q)$ phase diagram includes a sharp first-order phase transition separating a region with Newtonian rheology from a region with Bagnoldian rheology. In both cases, the system transitions to Newtonian behavior as one gets sufficiently close to the jamming transition. The models differ in their behavior at small $Q$, corresponding to strongly inelastic collisions. At small $Q$, model CD is Newtonian at all $\phi$, while CDrot has a first-order Bagnoldian to Newtonian transition as $\phi$ increases.

For CDrot we explore in greater depth the consequences of this rheological transition at small $Q=1$, as one varies the applied shear strain rate $\dot{\gamma}$. We find that in the $(\phi,\dot{\gamma})$ plane, the first-order rheological transition becomes a region of coexisting shear bands of Bagnoldian and Newtonian rheology. The width of this coexistence region appears to vanish as $\dot{\gamma}\to 0$ at a critical $\phi_{c}$ that lies below the jamming $\phi_{J}$. Upon increasing $\dot{\gamma}$ the coexistence region appears to vanish at a critical end-point. Thus we find that model CDrot displays shear banding within a well defined region of the $(\phi,\dot{\gamma})$ phase space, and that this shear banding is a consequence of an underlying Bagnoldian to Newtonian first-order phase transition.

We next consider the behavior of the shear stress $\sigma$ as one crosses through the shear banded coexistence region by varying $\dot{\gamma}$ at fixed $\phi$. We show that if $\dot{\gamma}$ is varied too rapidly for the system to relax to the shear-banded steady state corresponding to the instantaneous value of $\dot{\gamma}$, then there can be large jumps in stress, as well as hysteresis upon cycling $\dot{\gamma}$ up and down through the coexistence region, as the system switches from Bagnoldian to Newtonian rheology or vice versa. Thus the first-order shear banding transition becomes a possible mechanism for DST.

Our object in this work is not to claim that our model is a realistic description of any particular physical system. Nor do we argue that friction is not an important element for the observed DST in real non-Brownian suspensions. Our goal is rather to study a model, built with interactions commonly considered in the literature, yet still sufficiently simple so that we may systematically investigate wide ranges of the rheological phase diagram. We show that this phase diagram depends sensitively on the nature of the dissipative interaction, and that, even for frictionless particles, behavior may be richer than previously believed.

The remainder of our paper is organized as follows. In Sec. II we define our numerical models, the quantities we measure, and our method of simulation. In Sec III we present our numerical results. In Sec. III.1 we give results and map out the $(\phi,Q)$ phase diagram for low strain rates $\dot{\gamma}$ in model CD. We characterize the varying nature of particle collisions as one crosses the Bagnoldian to Newtonian transition line. In Sec. III.2 we give results and map out the $(\phi,Q)$ phase diagram at low $\dot{\gamma}$ for model CDrot. We pay particular attention to how the rotational motion of the particles effects dissipation and leads to a different phase structure at small $Q$, as compared to model CD. In Sec. III.3 we focus on model CDrot in the limit of strongly inelastic collisions at $Q=1$, and show how the Bagnoldian to Newtonian rheological transition becomes a coexistence region of Bagnoldian and Newtonian shear bands. In Sec. III.4 we show how crossing the coexistence region by ramping the strain rate $\dot{\gamma}$ up and down at fixed volume can lead to stress jumps and hysteresis. Finally in Sec. IV we summarize our conclusions and discuss directions for future investigation. Appendixes A and B provide some technical details concerning the location of phase boundaries in model CDrot.

II Model and Simulation Method

We use a well studied model OHern of frictionless, bidisperse, soft-core circular disks in two dimensions, with equal numbers of big and small particles with diameter ratio $d_{b}/d_{s}=1.4$. Particles interact only when they come into contact, in which case they repel with an elastic potential,

[TABLE]

Here $r_{ij}\equiv|\mathbf{r}_{ij}|$, where $\mathbf{r}_{ij}\equiv\mathbf{r}_{i}-\mathbf{r}_{j}$ is the center to center displacement from particle $j$ at position $\mathbf{r}_{j}$ to particle $i$ at $\mathbf{r}_{i}$, and $d_{ij}\equiv(d_{i}+d_{j})/2$ is the average of their diameters. In this work we will use the value $\alpha=2$, corresponding to a harmonic repulsion. The resulting elastic force on particle $i$ from particle $j$ is,

[TABLE]

where $\mathbf{\hat{r}}_{ij}\equiv\mathbf{r}_{ij}/r_{ij}$ is the inward pointing normal direction at the surface of particle $i$.

Because shearing injects energy into the system, we also need to have a mechanism of energy dissipation so as to attain a shear driven steady state. In the present work we will take dissipation to be due solely to particle collisions. We will consider two different models for this dissipation. Both models will include a tangential viscous dissipation, but no tangential elastic friction.

II.1 Model CD

The first, simpler, model is to take the dissipative force proportional to the difference in the center of mass velocities of the two colliding particles,

[TABLE]

where $\mathbf{v}_{i}=d\mathbf{r}_{i}/dt$ is the center of mass velocity of particle $i$. With this dissipation our model becomes a massive version of the Durian “bubble model” for foams Durian , and we have earlier denoted this model as “CD” for “contact dissipation” VOT_rheology . This model has also been used in simulations of convection in granular materials Taguchi ; Luding1 ; Luding2 . Note that the constant $k_{d}$ has different physical units from the constant $k_{e}$.

Given the above forces, particle motion is governed by the deterministic Newton’s equation,

[TABLE]

where $m_{i}$ is the mass of particle $i$ and the sum is over all particles $j$ in contact with particle $i$. In this work we take all particles to have equal mass independent of their size, $m_{i}=m$.

To apply a uniform shear strain rate $\dot{\gamma}$ in the $\mathbf{\hat{x}}$ direction, we use periodic Lees-Edwards boundary conditions LeesEdwards , so that a particle at position $\mathbf{r}=(r_{x},r_{y})$ has images at positions $(r_{x}+mL+n\gamma L,r_{y}+nL)$, with $n$, $m$ integer and $\gamma=\dot{\gamma}t$ the total shear strain at time $t$. We then numerically integrate the equations of motion (4) to determine the particle trajectories, and from that the stress tensor of the system.

If we take the collisional forces to act on a particle at the point of contact, then $\mathbf{f}^{\mathrm{el}}_{ij}$, being normally directed, exerts no torque on the particle. However $\mathbf{f}^{\mathrm{dis}}_{ij}$ can have a non-zero component tangent to the particle’s surface at the point of contact, and so can exert a torque and cause the particles to rotate. However, since such rotational motion does not feed back into the force $\mathbf{f}_{ij}^{\mathrm{dis}}$ of Eq. (3), in model CD the translational motion remains completely decoupled from the rotational motion and so rotational motion may be ignored when computing stresses.

II.2 Model CDrot

This leads us to introduce our next model “CDrot” in which rotational and translational motion couple in a non-trivial fashion. In this model we take the dissipation to be proportional to the velocity difference of colliding particles at the point of contact,

[TABLE]

where $\mathbf{v}_{i}^{C_{j}}$ is the local velocity of particle $i$ at its point of contact with particle $j$. $\mathbf{v}_{i}^{C_{j}}$ thus includes a term arising from the rotation of the particle,

[TABLE]

Here $\theta_{i}$ gives the angular orientation of the disk and $\mathbf{s}_{ji}$ is the moment arm pointing from the center of particle $i$ to the point of contact with particle $j$.

We now need to supplement Eq. (4) for the translational motion with an equation for rotational motion,

[TABLE]

where $I_{i}=\tfrac{1}{2}m(d_{i}/2)^{2}$ is the moment of inertia of particle $i$, and the sum is over all particles $j$ in contact with $i$. We note that such a coupling of rotational to translational motion is necessarily a feature of models that include inter-granular friction, although the tangential force in a frictional model is different from what we use here.

To generalize further, we might wish to split the dissipative force of Eq. (5) into components normal and tangential to the particle’s surface at the point of contact, with separate dissipative constants $k_{d}^{(n)}$ and $k_{d}^{(t)}$ respectively. If we took $k_{d}^{(n)}>0$ but $k_{d}^{(t)}=0$, then again the translational motion decouples from rotational motion, there being no tangential force. We have previously denoted VOT_rheology this model as “CDn” (for normal contact dissipation) and it is well known to have Bagnoldian rheology, with stress $\sim\dot{\gamma}^{2}$, for all values of $\phi$ and $Q$ VOT_CDn ; VOT_CDnQ . In the present work, however, we will keep $k_{d}^{(n)}=k_{d}^{(t)}\equiv k_{d}$ for simplicity.

II.3 Times Scales and Dimensionless Quantities

The above microscopic dynamics possesses two important microscopic time scales VOT_rheology , the elastic and dissipative relaxation times,

[TABLE]

We will measure time in units of $\tau_{e}$. The parameter

[TABLE]

measures the degree of elasticity of the collisions. For the harmonic interaction that we use, if we regarded the elastic potential of Eq. (1) as a spring which did not break when particles lose contact, then $2\pi\tau_{e}$ would give the undamped natural period of oscillation, $2\tau_{d}$ would be the decay time, and $Q$ would be the quality factor. We will denote small $Q$, where the head-on collision of two particles is strongly inelastic, as the “strongly inelastic” region; large $Q$ will then be denoted as the “weakly inelastic” region.

It is also useful to define the time scale VOT_rheology ,

[TABLE]

which has a well defined value in the limit that the particle mass $m\to 0$. This time scale will be useful for describing overdamped systems where inertial effects are unimportant.

Our system consists of a fixed total number particles $N$ in a square box of fixed length $L$. $L$ is chosen to set the particle packing fraction $\phi$,

[TABLE]

The behavior of our system is thus controlled by the three dimensionless parameters, $\phi$, $Q$, and the dimensionless strain rate $\dot{\gamma}\tau_{e}$.

To determine the global rheology of the system we measure the pressure tensor of each configuration, which consists of three pieces LeesEdwards : the elastic part $\mathbf{p}^{\mathrm{el}}$, arising from the repulsive elastic forces of Eq. (2),

[TABLE]

the dissipative part $\mathbf{p}^{\mathrm{dis}}$, arising from the dissipative forces of Eqs. (3) or (5),

[TABLE]

and the kinetic part $\mathbf{p}^{\mathrm{kin}}$ (sometimes called the streaming part),

[TABLE]

where $\delta\mathbf{v}_{i}$ is the fluctuation away from the average velocity profile that characterizes the shear flow. The total pressure tensor is then,

[TABLE]

The average pressure $p$ and shear stress $\sigma$ in the system are then,

[TABLE]

where $\langle\dots\rangle$ represents an ensemble average over configurations in the sheared steady state.

It will be convenient to work in terms of dimensionless quantities. We take the diameter of the small particles $d_{s}$, the mass $m$, and the time $\tau_{e}$ as our units of length, mass, and time respectively. With these choices, stress in two dimensions is measured in units of $m/\tau_{e}^{2}$, and so we can define a dimensionless pressure tensor,

[TABLE]

with similarly,

[TABLE]

as the dimensionless pressure and shear stress.

For a system with Bagnoldian rheology (i.e. $p,\sigma\sim\dot{\gamma}^{2}$, as $\dot{\gamma}\to 0$) we define the dimensionless Bagnold transport coefficients,

[TABLE]

For a system with Newtonian rheology (i.e. $p,\sigma\sim\dot{\gamma}$, as $\dot{\gamma}\to 0$) we define the dimensionless transport coefficients in term of the components of $Q\mathbf{P}/(\dot{\gamma}\tau_{e})=\mathbf{P}/(\dot{\gamma}\tau_{0})$,

[TABLE]

We define $\eta_{p}$ and $\eta_{\sigma}$ in terms of $\mathbf{P}/(\dot{\gamma}\tau_{0})$ and not $\mathbf{P}/(\dot{\gamma}\tau_{e})$ because we have previously shown VOT_rheology that, for systems with Newtonian rheology, the first form becomes independent of $Q$ in the hard-core limit of sufficiently small $\dot{\gamma}\tau_{0}$. $\eta_{\sigma}$ is just the dimensionless shear viscosity of the system, while $\eta_{p}$ is the pressure analog of viscosity.

In much of this work we will focus specifically on the elastic part of the pressure tensor, $\mathbf{p}^{\mathrm{el}}$. In general we find that, at the densities and strain rates we are considering, the elastic part gives the largest contribution to the total $\mathbf{p}$, and it is also the part that is easiest to compute accurately in simulations. Furthermore, any signature of a transition that we find in the elastic part must necessarily have a corresponding signature in the total stress. We will use the superscript “el” when referring to a quantity derived from only the elastic part of the pressure tensor, i.e., $P^{\mathrm{el}}$ for the elastic part of the pressure and $B_{p}^{\mathrm{el}}$ for the elastic part of the Bagnold coefficient.

II.4 Simulation Method

In our numerical simulations, we choose the diameter of the small particles to be $d_{s}=1$, the mass $m=1$, and take the unit of time $\tau_{e}=1$ (which implies the elastic coupling $k_{e}=1$). Since the net force $\mathbf{f}_{i}$ on particle $i$ depends on particle velocities, we use a modified velocity Verlet algorithm: With the position and velocity of particle $i$ at times $t$ and $t+\Delta t$ denoted by $(\mathbf{r}_{i0},\mathbf{v}_{i0})$ and $(\mathbf{r}_{i1},\mathbf{v}_{i1})$, respectively, and the net force by $\mathbf{f}_{i}(\{\mathbf{r}_{j},\mathbf{v}_{j}\})$, our equations are

[TABLE]

The equivalent of Eqs. (21) and (23) are used for the rotational velocity $d\theta_{i}/dt\equiv\dot{\theta}_{i}$. In general we use an integration time step $\Delta t/\tau_{e}=0.05$. The exception is for $Q<0.5$ (large $k_{d}$) where too large a time step may cause particles to unphysically reverse direction rather than just slow down; we then use time steps down to $\Delta t/\tau_{e}=0.01$.

We simulate with two different system sizes. For the smaller size we have $N=1024$ particles, and most of our presented results are for $\dot{\gamma}\tau_{e}=10^{-5}$ and $10^{-6}$. In these cases we simulate to a total shear strain $\gamma=\dot{\gamma}t$ of typically 100. For our larger system of $N=65536$ particles we simulate for a range of strain rates from $\dot{\gamma}\tau_{e}=2\times 10^{-7}$ up to $10^{-3}$. The length of our runs, as measured by the total shear strain $\gamma$, varies substantially for different shear rates. However, even for our lowest shear strain rate the run length is always $\gamma>1$.

III Results

III.1 Model CD in the $\dot{\gamma}\tau_{e}\to 0$ limit

We first discuss our results for model CD. In this section (as well as the next section for model CDrot) we are concerned with inferring the qualitative behavior of our system in the limit of small strain rate $\dot{\gamma}\tau_{e}\to 0$. Because we are simulating over a wide range of parameters $(\phi,Q)$, we use a relatively small system with $N=1024$ particles. Most of our results will be for the strain rate $\dot{\gamma}\tau_{e}=10^{-5}$, although at some $\phi$ we will also consider strain rates of $10^{-6}$ and $10^{-7}$ so as to infer the $\dot{\gamma}\tau_{e}\to 0$ limit. But we caution that the phase diagram we find (see Fig. 3 for CD and Fig. 17 for CDrot) may be only approximate, with possible effects due to the finite values of $N$ and $\dot{\gamma}\tau_{e}$, especially as we get closer to critical points.

III.1.1 Rheological transitions

In Fig. 1 we plot results for the elastic part of the dimensionless pressure $P^{\mathrm{el}}$ vs the collision elasticity parameter $Q$, for several different packing fractions $\phi$. We show results for the strain rate $\dot{\gamma}\tau_{e}=10^{-5}$. We have previously found for this model that the shear driven jamming transition takes place at $\phi_{J}\approx 0.8433$ VOT_RDCD . As $Q$ increases we see that $P^{\mathrm{el}}$ decreases roughly as $1/Q$, but then at a certain $Q_{c}(\phi)$ takes a sharp discontinuous drop. Above $Q_{c}(\phi)$, $P^{\mathrm{el}}$ then increases roughly linearly in $Q$. As $\phi$ increases, the magnitude of this drop decreases, and the width in $Q$ over which the drop occurs broadens somewhat. As will be shown below, this broadening as $\phi$ increases is a consequence of the finite strain rate $\dot{\gamma}\tau_{e}$; the drop will sharpen if one goes to lower values of $\dot{\gamma}\tau_{e}$. Very similar behavior is seen for the shear stress $\sigma^{\mathrm{el}}$, with a discontinuous drop at the same $Q_{c}(\phi)$ as found for pressure.

We now show that for $Q<Q_{c}(\phi)$ the system has Newtonian rheology, while for $Q>Q_{c}(\phi)$ the rheology is Bagnoldian. In Figs. 2(a) and (b) we show the pressure analog of viscosity, $\eta_{p}^{\mathrm{el}}=P^{\mathrm{el}}/(\dot{\gamma}\tau_{0})$, and the Bagnold coefficient for pressure, $B_{p}^{\mathrm{el}}=P^{\mathrm{el}}/(\dot{\gamma}\tau_{e})^{2}$, vs $Q$ for several different strain rates $\dot{\gamma}\tau_{e}$ at the packing fraction $\phi=0.82$. For $Q<Q_{c}$ we see that $\eta_{p}^{\mathrm{el}}$ becomes independent of $\dot{\gamma}\tau_{e}$ for $\dot{\gamma}\tau_{e}$ sufficiently small, thus confirming Newtonian rheology. Moreover, below $Q_{c}$ we find $\eta_{p}^{\mathrm{el}}$ to be independent of $Q$, as we have found previously VOT_rheology . This independence of $\eta_{p}^{\mathrm{el}}$ on $Q$ thus explains the $1/Q$ dependence of $P^{\mathrm{el}}$ seen in Fig. 1 for $Q<Q_{c}$ (since $P^{\mathrm{el}}=\eta_{p}^{\mathrm{el}}\dot{\gamma}\tau_{e}/Q$). For $Q>Q_{c}$ we see that $B_{p}^{\mathrm{el}}$ becomes independent of $\dot{\gamma}\tau_{e}$, thus confirming Bagnoldian rheology. The roughly linear rise in $B_{p}^{\mathrm{el}}$ as $Q$ increases above $Q_{c}$ has been observed by us in other Bagnoldian systems at large $Q$ VOT_CDnQ , and explains the roughly linear increase in $P^{\mathrm{el}}$ seen in Fig. 1 for $Q>Q_{c}$. In both cases we see that the drop at $Q_{c}$ becomes sharper as $\dot{\gamma}\tau_{e}$ decreases. Similar results are shown for the lower packing fraction $\phi=0.76$ in Figs. 2(c) and (d).

III.1.2 Phase diagram

We thus conclude that $Q_{c}(\phi)$ marks a sharp discontinuous transition from a region of Newtonian rheology ($Q<Q_{c}$) to a region of Bagnoldian rheology ($Q>Q_{c}$), in the $\dot{\gamma}\tau_{e}\to 0$ limit. Estimating $Q_{c}$ from the midpoint of the sharp drop in $P^{\mathrm{el}}$ of Fig. 1, in Fig. 3 we plot $Q_{c}$ vs $\phi$, showing results from two different strain rates, $\dot{\gamma}\tau_{e}=10^{-5}$ and $10^{-6}$. We find complete agreement between the two different strain rates, suggesting that our results reasonably give the $\dot{\gamma}\tau_{e}\to 0$ limit. We see that $Q_{c}(\phi)$ appears to diverge as $\phi$ approaches the jamming $\phi_{J}$ from below. Thus in model CD, at any fixed $Q$, one eventually enters the Newtonian region sufficiently close to $\phi_{J}$, and so the jamming criticality is always that of the Newtonian system. Moreover, since we have seen in Figs. 2(a) and (c) that $\eta_{p}^{\mathrm{el}}$ is independent of $Q$ in the Newtonian region, the divergence of $\eta_{p}^{\mathrm{el}}$ as $\phi\to\phi_{J}$ must be the same for all values of $Q$. Thus the value of $\phi_{J}$ and the exponent characterizing the divergence of $\eta_{p}^{\mathrm{el}}$ should be independent of $Q$; the only thing that depends on $Q$ is how close one must get to $\phi_{J}$ in order to enter the $Q$-independent Newtonian region.

Making a power law fit, $Q_{c}(\phi)=Q_{0}+c(\phi_{J}-\phi)^{-x}$, we find good agreement with $x\approx 1.1\pm 0.1$ and $Q_{0}\approx 1.2$. However our data is not sufficiently close to the jamming $\phi_{J}\approx 0.8433$ for us to know whether this exponent represents the true asymptotic behavior as $\phi\to\phi_{J}$, or is only an effective value for the range of $\phi$ we have considered.

We next ask if there is any obvious physical signature that distinguishes the Newtonian phase of our system from the Bagnoldian phase. We find that in the Newtonian phase the particles tend to cluster together, with force chains that percolate throughout the system, while in the Bagnoldian phase particles have few instantaneous contacts VOT_rheology . This is seen by computing the average instantaneous number of contacts per particle $Z$. In Fig. 4 we plot $Z$ vs $Q$ for several different packing fractions $\phi$, for the shear strain rate $\dot{\gamma}\tau_{e}=10^{-5}$. We see that $Z$ takes a sharp drop from values of $O(1)$ to values several orders of magnitude smaller, at the same $Q_{c}$ where $P^{\mathrm{el}}$ has it’s drop.

In Fig. 5 we plot $Z$ vs $Q$ for several different strain rates $\dot{\gamma}\tau_{e}$, at $\phi=0.82$ and 0.76. We see that in the Newtonian phase $Q<Q_{c}$, as $\dot{\gamma}\tau_{e}$ decreases $Z$ becomes independent of $\dot{\gamma}\tau_{e}$ and independent of $Q$. Thus the Newtonian phase is characterized by a finite particle connectivity even as $\dot{\gamma}\tau_{e}\to 0$. In the Bagnoldian phase $Q>Q_{c}$, we see that $Z$ steadily decreases as $\dot{\gamma}\tau_{e}$ decreases. Below we will show that in the Bagnoldian phase the average contact number vanishes linearly, $Z\sim\dot{\gamma}\tau_{e}$, as $\dot{\gamma}\tau_{e}\to 0$.

To understand better the behavior of $Z$, we consider two time scales associated with particle collisions. We define the collision duration time $\tau_{\mathrm{dur}}$ as the average time from the initiation of a particle contact until the breaking of that contact. We define the collision rate $\nu_{\mathrm{coll}}$ as the average number of collisions per unit time divided by the number of particles. As has been argued in Ref. OHL , and as we have shown previously for model CDn VOT_CDnQ , $\tau_{\mathrm{dur}}$ and $\nu_{\mathrm{coll}}$ are related to the average contact number $Z$ by,

[TABLE]

The product $\tau_{\mathrm{dur}}\nu_{\mathrm{coll}}$ measures the fraction of the time between two collisions that a contact on average persists, while the factor 2 is because each contact is shared by two particles. In Fig. 5 we show $2\tau_{\mathrm{dur}}\nu_{\mathrm{coll}}$ vs $Q$ on the same plot as $Z$. Although these two quantities are measured by completely different methods within the simulation, we find perfect agreement, thus providing a consistency check on our calculations.

We now consider $\tau_{\mathrm{dur}}$ and $\nu_{\mathrm{coll}}$ separately. In Figs. 6(a) and (b) we plot $\tau_{\mathrm{dur}}/\tau_{e}$ and $\nu_{\mathrm{coll}}/\dot{\gamma}$ vs $Q$ for several different strain rates $\dot{\gamma}\tau_{e}$ at $\phi=0.82$. We see that for $Q>Q_{c}$ in the Bagnoldian region, $\tau_{\mathrm{dur}}/\tau_{e}$ and $\nu_{\mathrm{coll}}/\dot{\gamma}$ become independent of $\dot{\gamma}\tau_{e}$ at the smallest $\dot{\gamma}\tau_{e}$, thus implying that $Z=2\tau_{\mathrm{dur}}\nu_{\mathrm{coll}}\sim\dot{\gamma}\tau_{e}$ vanishes as $\dot{\gamma}\tau_{e}\to 0$. In the Bagnoldian region, $\tau_{\mathrm{dur}}/\tau_{e}\sim O(1)$, as would be the case for an isolated binary collision. In the Newtonian region $Q<Q_{c}$, $\tau_{\mathrm{dur}}/\tau_{e}$ rapidly increases by several orders of magnitude as $Q$ decreases, while $\nu_{\mathrm{coll}}$ decreases less rapidly than linearly in $\dot{\gamma}$ (so that $\nu_{\mathrm{coll}}/\dot{\gamma}$ increases as $\dot{\gamma}$ decreases), so as to keep $Z=2\tau_{\mathrm{dur}}\nu_{\mathrm{coll}}$ constant. Similar results are shown in Figs. 6(c) and (d) for the lower $\phi=0.76$.

III.1.3 Geometry of collisions

Next, we look more directly at the geometry of individual collisions in the phases above and below $Q_{c}$. To measure this, let us define,

[TABLE]

as the position and velocity of particle $i$ with respect to particle $j$. We then define the angle $\vartheta$ as the angle by which one must rotate $\mathbf{v}_{ij}$ to align it parallel with $\mathbf{r}_{ij}$,

[TABLE]

For two particles just initiating a contact, we must have $\cos\vartheta<0$, so that the particles are driven into each other, as illustrated in Fig. 7(a). In this case we must have $|\vartheta|>90^{\circ}$. For two particles just breaking a contact, we must have $\cos\vartheta>0$, so that the particles are driven away from each other, as illustrated in Fig. 7(b). In this case we must have $|\vartheta|<90^{\circ}$.

Measuring the value of $\vartheta$ each time a contact is initiated and each time a contact is broken, we construct a histogram ${\cal P}(\vartheta)$ which combines both contact initiation and contact breaking events. In Fig. 8 we plot ${\cal P}(\vartheta)$ vs $\vartheta$ for the two packing fractions $\phi=0.82$ and 0.76. We use a bin width of $\Delta\vartheta=3.6^{\circ}$ in constructing this histogram. In each case we show ${\cal P}(\vartheta)$ for three different values of $Q$: one well above $Q_{c}$, one well below $Q_{c}$, and one approximately at $Q_{c}$. We see that for $Q>Q_{c}$, the histogram ${\cal P}(\vartheta)$ has minima at $\vartheta=\pm 90^{\circ}$ and is broadly distributed between; collisions take place with a broad range of impact parameter, $b=|\mathbf{r}_{ij}\times\mathbf{\hat{v}}_{ij}|$, including many normally oriented (i.e. head-on) collisions. For $Q<Q_{c}$, however, ${\cal P}(\vartheta)$ is strongly peaked at $\vartheta=\pm 90^{\circ}$, indicating that collision are primarily glancing, with particles approaching each other close to tangentially. In Fig. 9 we plot the value of the histogram ${\cal P}(\vartheta=90^{\circ})$ vs $Q$, as a measure of the propensity for tangential collisions. We show results at $\phi=0.82$ and $0.76$ at $\dot{\gamma}\tau_{e}=10^{-5}$. We see that ${\cal P}(90^{\circ})$ is essentially zero for $Q>Q_{c}$ in the Bagnoldian region, and then takes a sharp jump upwards at $Q=Q_{c}$ with the transition to the Newtonian region.

Looking at behavior more closely at the level of individual collisions, we find the following. For $Q>Q_{c}$ in the Bagnoldian region, collisions are mostly isolated binary collisions, where a pair of particles collide with a wide range of impact parameter, and then separate before coming into contact with another particle. If we define $v_{ij,N}$ and $v_{ij,T}$ as the root mean square average of the normal and tangential components of $\mathbf{v}_{ij}$ respectively, then for $Q>Q_{c}$ we find $v_{ij,N}\approx v_{ij,T}$, with both decreasing as $Q$ decreases. As $Q$ decreases below $Q_{c}$, however, $v_{ij,N}$ takes a sharp drop to small values, while $v_{ij,T}$ plateaus to a constant independent of $Q$. The drop in $v_{ij,N}$ results from the increasing inelasticity of the collisions as $Q$ gets small, combined with the steric interactions of excluded volume at large particle densities; the decrease in free volume as $\phi$ increases, results in the increase of $Q_{c}(\phi)$. The constant value of $v_{ij,T}$ results from the driving shear strain rate $\dot{\gamma}$, which forces particles to move relative to one another in shear flow, even for the smallest $Q$. Thus below $Q_{c}$ the relative motion between particles is largely tangential. For such tangential motion, the tangential component of the dissipative forces resists relative particle motion, while the applied $\dot{\gamma}$ requires such motion. This results in the build up of extended force chains of elastic force, that drive the particles forward with the required $\dot{\gamma}$. Such force chains result in a sharp increase in the average particle contact number $Z$ and the elastic pressure $P^{\mathrm{el}}$. We leave further investigation of the detailed connection between particle collisions and the transition from Newtonian to Bagnoldian rheology to future work.

III.2 Model CDrot in the $\dot{\gamma}\tau_{e}\to 0$ limit

We now consider behavior in model CDrot, in which the rotational and translational degrees of freedom are coupled. Our results in this section are again for a model with $N=1024$ particles. Our main results are for a strain rate $\dot{\gamma}\tau_{e}=10^{-5}$, although we also use $10^{-6}$ and $10^{-7}$ at several values of $\phi$ so as to better illustrate and locate the rheological phase transitions.

III.2.1 Rheological transitions

In Fig. 10 we plot results for the elastic part of the dimensionless pressure $P^{\mathrm{el}}$ vs the collision elasticity parameter $Q$, for different packing fractions $\phi$ and strain rate $\dot{\gamma}\tau_{e}=10^{-5}$. For the higher packing fractions, $\phi\gtrsim 0.80$, our results for CDrot appear qualitatively the same as found in Fig. 1 for model CD; as $Q$ increases, $P^{\mathrm{el}}$ initially decreases roughly as $\sim 1/Q$, until a $Q_{c}$ is reached at which $P^{\mathrm{el}}$ takes a sharp drop, and then $P^{\mathrm{el}}$ increases roughly linearly as $Q$ increases further. We will see below that this drop at $Q_{c}$ corresponds to a sharp transition from Newtonian to Bagnoldian rheology, as in model CD. For smaller $\phi<0.80$, however, the drop becomes broadened over a wider range of $\phi$, and $P^{\mathrm{el}}$ appears to be approaching a constant as $Q$ further decreases. For the smallest $\phi$, the jump almost disappears, leaving only a kink between the small $Q$ and large $Q$ behaviors. We will see below that this drop/kink corresponds to a non-singular crossover $Q^{*}$ between two different regions of Bagnoldian rheology.

In Figs. 11(a) and (b) we plot $\eta_{p}^{\mathrm{el}}$ and $B_{p}^{\mathrm{el}}$ for strain rates $\dot{\gamma}\tau_{e}=10^{-4}$ to $10^{-7}$, at $\phi=0.82$. We see that the behavior is qualitatively the same as in Figs. 2(a) and (b) for model CD; $\eta_{p}^{\mathrm{el}}$ is independent of $\dot{\gamma}\tau_{e}$ for $Q<Q_{c}$, while $B_{p}^{\mathrm{el}}$ is independent of $\dot{\gamma}\tau_{e}$ for $Q>Q_{c}$. We thus conclude that the rheology is Newtonian for $Q<Q_{c}$, but Bagnoldian for $Q>Q_{c}$. In Figs. 11(c) and (d) we similarly plot $\eta_{p}^{\mathrm{el}}$ and $B_{p}^{\mathrm{el}}$ for the lower packing fraction $\phi=0.76$. Here we see rather different behavior; $\eta_{p}^{\mathrm{el}}$ keeps decreasing as $\dot{\gamma}\tau_{e}$ decreases, while $B_{p}^{\mathrm{el}}$ is essentially independent of $\dot{\gamma}\tau_{e}$ at all $Q$, indicating that the rheology is Bagnoldian both above and below $Q^{*}$.

The same conclusion is obtained by looking at the average instantaneous contact number $Z$, which we plot vs $Q$ in Fig. 12 for several different $\phi$ at $\dot{\gamma}\tau_{e}=10^{-5}$. For $\phi\gtrsim 0.80$ we see $Z\sim O(1)$ at low $Q$, characteristic of the Newtonian phase; as $Q$ increases, $Z$ sharply drops roughly two orders of magnitude, characteristic of the Bagnoldian phase. For $\phi<0.80$, we see $Z<1$ at small $Q$ for all $\phi$, decreasing rapidly as $\phi$ decreases, consistent with a Bagnoldian phase at all $Q$. This becomes clearer in Fig. 13 where we show $Z$ vs $Q$ for strain rates $\dot{\gamma}\tau_{e}=10^{-4}$ to $10^{-7}$ at $\phi=0.82$ and $0.76$. For $\phi=0.82$ in panel (a) we see that, for $Q<Q_{c}$, $Z$ is a constant $O(1)$ independent of $\dot{\gamma}\tau_{e}$, while for $Q>Q_{c}$, $Z$ decreases to zero as $\dot{\gamma}\tau_{e}$ decreases. This is the same behavior seen in Fig. 5(a) for model CD. For $\phi=0.76$ in panel (b), however, we see that $Z$ decreases with decreasing $\dot{\gamma}\tau_{e}$ for all $Q$. We also show in Fig. 13 the quantity $2\tau_{\mathrm{dur}}\nu_{\mathrm{coll}}$ and find everywhere excellent agreement with $Z$, as expected according to Eq. (25).

In Figs. 14(a) and (b) we plot $\tau_{\mathrm{dur}}/\tau_{e}$ and $\nu_{\mathrm{coll}}/\dot{\gamma}$ at $\phi=0.82$ for strain rates $\dot{\gamma}\tau_{e}=10^{-4}$ to $10^{-7}$. We see behavior similar to that of model CD in Figs. 6(a) and (b). For $Q>Q_{c}$ we see that both $\tau_{\mathrm{dur}}/\tau_{e}$ and $\nu_{\mathrm{coll}}/\dot{\gamma}$ are independent of $\dot{\gamma}\tau_{e}$, demonstrating that $Z=2\tau_{\mathrm{dur}}\nu_{\mathrm{coll}}\sim\dot{\gamma}\tau_{e}$ vanishes linearly in $\dot{\gamma}\tau_{e}$ as $\dot{\gamma}\tau_{e}\to 0$. As $Q$ decreases below $Q_{c}$, $\tau_{\mathrm{dur}}$ increases several orders of magnitude while $\nu_{\mathrm{coll}}/\dot{\gamma}$ decreases correspondingly, so that $Z=2\tau_{\mathrm{dur}}\nu_{\mathrm{coll}}$ remains constant, as explicitly seen in Fig. 13(a) for $Q<Q_{c}$.

In Figs. 14(c) and (d) we plot these same quantities, but at the lower packing fraction $\phi=0.76$. Again $\tau_{\mathrm{dur}}/\tau_{e}$ and $\nu_{\mathrm{coll}}/\dot{\gamma}$ are independent of $\dot{\gamma}\tau_{e}$ for $Q>Q^{*}$, demonstrating that in this large $Q$ region $Z\sim\dot{\gamma}\tau_{e}$ vanishes as $\dot{\gamma}\tau_{e}\to 0$. However, for $Q<Q^{*}$, behavior is now different from what is seen in Figs. 14(a) and (b) at the higher $\phi=0.82$. We still find that $\tau_{\mathrm{dur}}/\tau_{e}$ increases as $Q$ decreases below $Q^{*}$, but now that increase is much less dramatic (note the linear vertical scale). We see that $\nu_{\mathrm{coll}}/\dot{\gamma}$ increases with decreasing $\dot{\gamma}\tau_{e}$ for $Q<Q^{*}$, but since from Fig. 13(b) we see that $Z$ decreases with decreasing $\dot{\gamma}\tau_{e}$, and $\tau_{\mathrm{dur}}/\tau_{e}$ appears to be approaching a constant as $\dot{\gamma}\tau_{e}$ decreases, we conclude from Eq. (25) that $\nu_{\mathrm{coll}}$ (and hence $Z$) decreases with decreasing $\dot{\gamma}\tau_{e}$, but more slowly than linearly in $\dot{\gamma}\tau_{e}$. This behavior at $\phi=0.76$ is similar to what we have seen previously in model CDn VOT_CDnQ , where dissipation depends only on the normal component of the velocity difference of the colliding particles.

III.2.2 Phase diagram

From the results in the preceding section we thus conclude that, for $\phi=0.82$, upon increasing $Q$ there is a sharp transition at $Q_{c}$ from a phase where the $\dot{\gamma}\tau_{e}\to 0$ limiting rheology is Newtonian, to a phase where the rheology is Bagnoldian. For $\phi=0.76$, as $Q$ increases above $Q^{*}$, there is a non-singular crossover from a strongly inelastic Bagnoldian rheology to a weakly inelastic Bagnoldian rheology; these two different Bagnoldian regions differ primarily in how $P^{\mathrm{el}}$ varies with $Q$ and how $Z$ varies with $\dot{\gamma}\tau_{e}$. We have previously found similar behavior in model CDn VOT_CDnQ ; Lois . Since we find that the small $Q$ behavior is Newtonian at large $\phi$, and Bagnoldian at small $\phi$, it therefore follows that at small $Q$ there must be a sharp transition from Bagnoldian to Newtonian rheology as $\phi$ increases.

To search for this small $Q$ transition from Bagnoldian to Newtonian rheology, in Fig. 15 we plot the elastic part of the dimensionless pressure $P^{\mathrm{el}}$ vs $\phi$ for several different values of $Q$ at fixed strain rate $\dot{\gamma}\tau_{e}=10^{-5}$. We see a clear step upwards in the vicinity of $\phi\approx 0.79$ – 0.81, depending on the value of $Q$. To confirm that these steps do indeed represent a sharp Bagnoldian-to-Newtonian rheological transition, in Figs. 16(a) and (b) we plot $\eta_{p}^{\mathrm{el}}$ and $B_{p}^{\mathrm{el}}$ vs $\phi$ at several different strain rates $\dot{\gamma}\tau_{e}$, for the value $Q=1$. We see that there is a $\phi_{c}$ such that for $\phi>\phi_{c}$ the curves of $\eta_{p}^{\mathrm{el}}$ become independent of $\dot{\gamma}\tau_{e}$, thus demonstrating Newtonian rheology. For $\phi<\phi_{c}$ the curves of $B_{p}^{\mathrm{el}}$ become independent of $\dot{\gamma}\tau_{e}$, thus demonstrating Bagnoldian rheology. In Figs. 16(c) and (d) we show similar results, with similar conclusions, for $Q=20$. We find that this $\phi_{c}$ tends to lie near the upper shoulder of the steps in Fig. 15. From Fig. 16 we note that there appears to be a small interval in $\phi$ between the point where the system ceases to be Bagnoldian and the point where it starts to be Newtonian. We will address this in greater detail in the next section.

Using the specific methods detailed in Appendix A, we proceed to determine the location of the crossover line $Q^{*}(\phi)$ between the regions of weak inelastic and strong inelastic Bagnoldian rheology at smaller $Q$ (as indicated by Fig. 11(c) and (d)), the transition line $Q_{c}(\phi)$ between the regions of Bagnoldian and Newtonian rheology at larger $Q$ (as indicated by Fig. 11(a) and (b)), and the transition line $\phi_{c}(Q)$ between regions of Bagnoldian and Newtonian rheology at small $Q$ (as indicated by Fig. 16). The resulting phase diagram is shown in Fig. 17. We see that the transition lines $Q_{c}(\phi)$ and $\phi_{c}(Q)$ form one continuous curve that separates the Newtonian phase from the Bagnoldian phase. As in our earlier model CD, our results suggest that the system is always Newtonian at the jamming transition $\phi_{J}$; the independence of $\eta_{p}^{\mathrm{el}}$ on $Q$ in this region, as seen in Fig. 11(a), suggests that the critical parameters of the jamming transition must also be independent of $Q$. Unlike model CD, however, we see that the Newtonian region is confined to relatively large packing fractions, $\phi>\phi_{c}\approx 0.80$. For $\phi<\phi_{c}$, even at small $Q$ where the collisions are strongly inelastic, the rheology remains Bagnoldian; upon increasing $Q$ at fixed $\phi<\phi_{c}$ there is a non-singular crossover $Q^{*}$ between two different regions of Bagnoldian rheology, as we have previously found in model CDn VOT_CDnQ .

III.2.3 Geometry of collisions

We next look at the geometry of collisions, computing the histogram ${\cal P}(\vartheta)$ of the collision angle $\vartheta$ at the time of initiation of a contact and the time of breaking of a contact (see Fig. 7 for the definition of $\vartheta$). In Fig. 18(a) we plot ${\cal P}(\vartheta)$ vs $\vartheta$ at $\phi=0.82>\phi_{c}$, for three values of $Q$: well above, well below, and near $Q_{c}$. In Fig. 18(b) we plot similarly at $\phi=0.76<\phi_{c}$, for $Q$ well above, well below, and near $Q^{*}$. In Fig. 19 we plot ${\cal P}(\vartheta=90^{\circ})$ vs $Q$ for $\phi=0.82$ and $0.76$. Both above and below $\phi_{c}$ we see the same behavior as was found previously in model CD (compare with Figs. 8 and 9). For $Q>Q_{c}$ (at $\phi>\phi_{c}$) or $Q>Q^{*}$ (at $\phi<\phi_{c}$), we see that collisions take place over a broad range of $\vartheta$; tangentially directed collisions at $\vartheta=90^{\circ}$ occur a negligible fraction of the time. For $Q<Q_{c}$ or $Q<Q^{*}$, normally directed collisions are greatly suppressed, and the relative motion of particles is primarily tangential.

III.2.4 Effect of particle rotations

It therefore remains to explain, at low $Q$ where relative particle motion is primarily tangential at any of the $\phi$ values we have considered here, what is the cause of the different rheology for $\phi<\phi_{c}$ (Bagnoldian) compared to $\phi>\phi_{c}$ (Newtonian). Since CDrot couples rotational and translational motion while CD does not, it is natural to look to the rotational motion of the particles for the answer. We therefore consider the rotational velocity of the particles, $\dot{\theta}_{i}\equiv d\theta_{i}/dt$.

In Fig. 20(a) we plot the average $\langle\dot{\theta}_{i}/\dot{\gamma}\rangle$ vs $\phi$ for several of the smaller values of $Q\leq 50$, at $\dot{\gamma}\tau_{e}=10^{-6}$. We see that $\langle\dot{\theta}_{i}/\dot{\gamma}\rangle\approx-1/2$ is roughly independent of $\phi$, as would be expected for circular disks in a uniform shear flow (the minus sign indicates a clockwise rotation). There is no visible feature as $\phi$ passes through $\phi_{c}$. Next we consider the variance of the single particle rotational velocity,

[TABLE]

where $\langle\dots\rangle$ denotes the time average over the steady state, and

[TABLE]

In Fig. 20(b) we plot $\mathrm{var}[\dot{\theta}_{i}/\dot{\gamma}]$ vs $\phi$ for different $Q$, at $\dot{\gamma}\tau_{e}=10^{-6}$. We now see a clear signature of the transition at $\phi_{c}$. As $\phi$ increases within the Bagnoldian phase, $\mathrm{var}[\dot{\theta}_{i}/\dot{\gamma}]$ steadily increases, reaches a peak at $\phi\approx\phi_{c}(Q)$, and then drops abruptly as the Newtonian phase is entered, suggesting that fluctuations of rotation are more constrained in the Newtonian phase. Note, although $\mathrm{var}[\dot{\theta}_{i}/\dot{\gamma}]$ drops in the Newtonian phase, fluctuations remain large compared to the average, $\langle\dot{\theta}_{i}/\dot{\gamma}\rangle\approx-1/2$.

For a clearer understanding of the effects of the rotational motion on the rheology, we now consider behavior at the level of individual collisions. Defining the dimensionless velocity $\mathbf{V}\equiv\mathbf{v}/(d_{s}\dot{\gamma})$, we can write the difference in velocities at the point of contact of particle $i$ colliding with particle $j$ as,

[TABLE]

where $\mathbf{V}_{ij}$ is the difference in center of mass velocities, and from Eq. (6),

[TABLE]

is the contribution from the rotational motion of the particles (we have used $\mathbf{r}_{ij}=-\mathbf{r}_{ji}$). With $\mathbf{\hat{r}}_{ji}$ as the outward unit normal vector at the point of contact on particle $i$, and $\mathbf{\hat{t}}_{ji}\equiv\mathbf{\hat{z}}\times\mathbf{\hat{r}}_{ji}$ as the unit tangent vector at the point of contact, we see that $\mathbf{V}^{R}_{ij}$ is purely tangential. We define,

[TABLE]

as the rotational part of the contact velocity difference,

[TABLE]

as the tangential and normal components of the contact velocity difference, and

[TABLE]

as the tangential and normal components of the center of mass velocity difference. Since $\mathbf{V}^{R}_{ij}$ is purely tangential, $V^{C}_{ij,N}=V_{ij,N}$.

In Fig. 21 we show correlation scatter plots of different components of the contact velocity difference at $\phi=0.76<\phi_{c}$ (left column) and $\phi=0.82>\phi_{c}$ (right column) for $Q=1$ and $\dot{\gamma}\tau_{e}=10^{-5}$. The data points represent $5000$ independent collisions sampled at random times during the duration of each collision. In each panel the white dot in the center of the data locates the average of the distribution, while the black dotted lines give the directions of the principal axes of the covariance matrix.

We consider first behavior at $\phi=0.76<\phi_{c}$. In panel (a) we show the correlation of $V_{ij,T}$ with $V^{R}_{ij}$. We see almost a perfect anticorrellation, with the result from Eq. (30) that the tangential component of the contact velocity difference $V^{C}_{ij,T}\approx 0$. This is confirmed in panel (c), where we show the correlation of $V^{C}_{ij,T}$ with $V^{R}_{ij}$. Despite a wide spread in values for $V^{R}_{ij}$, we find $V^{C}_{ij,T}$ is always close to zero. In panel (e) we show the correlation of the normal component of the contact velocity difference, $V_{ij,N}=V^{C}_{ij,N}$, with the tangential component $V^{C}_{ij,T}$. We see that both components are comparably small. Thus at small $Q$ in the Bagnoldian phase below $\phi_{c}$, the normally directed relative motion between particles is suppressed, relative motion is primarily tangential, but the particles are able to adjust their rotation to approximately cancel out the difference in center of mass motion, so that the tangential component of the contact velocity difference is comparably as small as the normal component. Thus both normal and tangential components of the dissipative force remain small. Behavior here is thus similar to what we have previously found at low $Q$ in model CDn VOT_CDnQ , where the tangential component of the dissipative force is explicitly set to zero.

We next consider the behavior at $\phi=0.82>\phi_{c}$. In panel (b) we show $V_{ij,T}$ vs $V^{R}_{ij}$. Compared to panel (a) we see that there is still correlation, however the scatter of the data is broader and the slope of the major axis of the covariance matrix axis is no longer $-1$, showing that the sum of these two, i.e. $V^{C}_{ij,T}$, remains finite. This is confirmed in panel (d) where we show $V^{C}_{ij,T}$ vs $V^{R}_{ij}$. We see that the spread in both quantities remains comparable, and the nearly vertical and horizontal orientations of the covariance matrix axes show that these quantities are at best very weakly correlated. Thus the tangential component of the contact velocity difference is now largely independent of the rotational motion. This observation is sufficient to explain the correlations seen in panel (b). In panel (f) we show $V_{ij,N}$ vs $V^{C}_{ij,T}$. We see that the spread in the normal component of the contact velocity difference remains small, comparable to that of panel (e), but that the tangential component is comparatively broadly distributed. Thus at small $Q$ in the Newtonian phase above $\phi_{c}$, the normally directed relative motion between particles is suppressed, relative motion is primarily tangential, the particles fluctuate in their rotational motion, but these fluctuations are largely uncorrelated with the tangential component of the contact velocity difference. The normal component of the dissipative force is small, but the tangential component is comparatively large; both are largely unaffected by the rotational motion of the particles. Behavior here is thus similar to what we find in model CD.

We thus see that the transition from Bagnoldian to Newtonian rheology, upon increasing $\phi$ at small $Q$, is due to the changing nature of correlations between rotational and translational motion as the system gets denser. Below $\phi_{c}$ the rotations are almost perfectly anticorrelated with the differences in center of mass motion, so that the tangential part of the contact velocity difference $V^{C}_{ij,T}$ becomes essentially zero, and dissipative forces are small (see Figs 21(a) and (c)). Above $\phi_{c}$, rotational motion becomes less correlated with the translational motion, so that $V^{C}_{ij,T}$ has sizable fluctuations which are uncorrelated with the particle rotations; rotational motion thus has little effect on the dissipative forces, which become sizable (see Figs. 21(b) and (d)). We find that, as $\phi$ increases, the number of particles in mutual contact during any particular collision is increasing. We therefore speculate that this change in the correlation between rotational and translational motion is due to the increasing constraint on rotational motion that arises when the number of contacts per collision increases. We leave further exploration of this effect to future work.

III.3 Model CDrot: Shear banding at finite $\dot{\gamma}\tau_{e}$

In the previous sections we focused on sketching out the behavior in the $(\phi,Q)$ plane as $\dot{\gamma}\tau_{e}\to 0$. In this section we focus more carefully on behavior as a function of $\dot{\gamma}\tau_{e}$. We consider here the specific case of crossing the Bagnoldian to Newtonian transition line at $\phi_{c}$ for small $Q=1$ in model CDrot. We will find that upon increasing $\dot{\gamma}\tau_{e}$, the sharp discontinuous transition becomes a coexistence region of finite width, where a shear band of Bagnoldian rheology coexists with a shear band of Newtonian rheology. Preliminary results suggest that similar behavior exists when crossing the Bagnoldian to Newtonian transition elsewhere at $Q_{c}(\phi)$ in both the models CD and CDrot. To investigate the spatially inhomogeneous configurations associated with this shear banding, it is necessary to consider a much larger system than in the previous sections. Here we use a system with $N=65536$ particles.

III.3.1 Indicators for a coexistence region

In Fig. 22 we plot results for the elastic part of the dimensionless pressure $P^{\mathrm{el}}$ vs packing fraction $\phi$ for different fixed values of the strain rate $\dot{\gamma}\tau_{e}$ at $Q=1$. At the smallest $\dot{\gamma}\tau_{e}$ we see a sharp jump in $P^{\mathrm{el}}$ that signals the transition at $\phi_{c}$ from Bagnoldian to Newtonian rheology. The jump here is indeed sharp, compared to the smoother behavior seen previously in Fig. 15, because the strain rate is much lower and our system is much larger; this $\dot{\gamma}\tau_{e}\to 0$ limiting value $\phi_{c}\approx 0.807$ is also slightly larger than shown in Fig. 17 for the same reason. However as $\dot{\gamma}\tau_{e}$ increases, we see that this jump broadens and the $P^{\mathrm{el}}$ vs $\phi$ curve develops two noticeable kinks rather than a single sharp jump. These kinks are denoted in the figure by the large open squares. At sufficiently large $\dot{\gamma}\tau_{e}$ these kinks disappear, and the $P^{\mathrm{el}}$ vs $\phi$ curve looks smooth. For a given $\dot{\gamma}\tau_{e}$, we denote the kink at the smaller $\phi$ by $\phi_{B}(\dot{\gamma}\tau_{e})$, and the kink at the larger $\phi$ by $\phi_{N}(\dot{\gamma}\tau_{e})$. The dashed lines in Fig. 22 connect these kinks so as to delimit a region in the $(\phi,P^{\mathrm{el}})$ plane. We will soon show that this region denotes a coexistence region where a band of Bagnoldian rheology coexists in mechanical equilibrium with a band of Newtonian rheology.

We next consider the macroscopic friction $\mu\equiv\sigma/p$. Since coexistence between bands of Bagnoldian and Newtonian rheology must result from mechanical equilibrium between states of equal stress, we will consider here the total pressure $p$ and total shear stress $\sigma$, rather than just the elastic parts. We note however, that for the parameters considered here, the dissipative and kinetic parts of the pressure are negligible, as is the kinetic part of the shear stress; the dissipative $\sigma^{\mathrm{dis}}$ contributes non-negligibly, roughly 10% – 20%, to the total $\sigma$. The quantity $\mu$ is of interest to consider since, unlike $p$ or $\sigma$, $\mu$ approaches a finite value as $\dot{\gamma}\tau_{e}\to 0$ in both Bagnoldian and Newtonian systems. Thus for sufficiently small $\dot{\gamma}\tau_{e}$ we might expect that, to a large extent, the curves of $\mu$ vs $\phi$ for different $\dot{\gamma}\tau_{e}$ will overlap.

In Fig. 23 we plot $\mu$ vs $\phi$ for different values of $\dot{\gamma}\tau_{e}$ at $Q=1$. We see that the data for different $\dot{\gamma}\tau_{e}$ do indeed seem to collapse to a common curve, but that this curve has two distinct branches. For a given $\dot{\gamma}\tau_{e}$, the curve switches between the two branches at packing fractions $\phi_{B}(\dot{\gamma}\tau_{e})$ and $\phi_{N}(\dot{\gamma}\tau_{e})$, with $\phi_{B}(\dot{\gamma}\tau_{e})<\phi_{N}(\dot{\gamma}\tau_{e})$. These switching points, as judged by eyeball, are denoted in the figure by the large open circles. As $\dot{\gamma}\tau_{e}$ decreases, $\phi_{B}(\dot{\gamma}\tau_{e})$ and $\phi_{N}(\dot{\gamma}\tau_{e})$ increase towards a common $\phi_{c}$. The two branches we see in $\mu$ suggest a coexistence region, where $\phi_{B}(\dot{\gamma}\tau_{e})$ locates the Bagnoldian side of the coexistence region while $\phi_{N}(\dot{\gamma}\tau_{e})$ locates the Newtonian side.

We now consider a more direct global measure of shear banding, by looking how the velocity profile compares to the linear profile $\langle v_{x}(\mathbf{r})\rangle=\dot{\gamma}\mathbf{r}\cdot\mathbf{\hat{y}}$ , $\langle v_{y}(\mathbf{r})\rangle=0$, that is expected for a uniform homogeneous shear flow. We measure the fluctuation away from this presumed average by,

[TABLE]

For a homogeneously sheared state we expect $\langle\delta v_{x}^{2}\rangle$ to be roughly the same order of magnitude as $\langle v_{y}^{2}\rangle$. However in a shear banded state, where the local strain rate is nonuniform, we would expect $\langle\delta v_{x}^{2}\rangle$, as defined above with respect to the fixed average $\dot{\gamma}$, to become anomalously large; this follows since the velocity profile is no longer linear over the entire width of the system, but only piecewise linear, with different slopes, in each shear band (as we will soon see in Fig. 28).

Defining the ratio of velocity fluctuations defined in Eqs. (35) and (36) as,

[TABLE]

in Fig. 24 we plot $R_{v}$ vs $\phi$ for different fixed values of $\dot{\gamma}\tau_{e}$. We see that at low and high $\phi$, $R_{v}$ is indeed of $O(1)$ as in a homogeneously sheared state. But for intermediate $\phi$ near $\phi_{c}$, we see a large jump in $R_{v}$ as expected for a shear banded state. For a given $\dot{\gamma}\tau_{e}$, we will denote as $\phi_{B}(\dot{\gamma}\tau_{e})$ the value of $\phi$ at which $R_{v}$ takes the sharp jump upwards from $O(1)$, and the higher $\phi_{N}(\dot{\gamma}\tau_{e})$ as the value of $\phi$ at which $R_{v}$ takes the sharp jump back down to $O(1)$.

III.3.2 Phase diagram and shear banding

We now have from $P^{\mathrm{el}}$, $\mu$ and $R_{v}$, three different determinations of the curves $\phi_{B}(\dot{\gamma}\tau_{e})$ and $\phi_{N}(\dot{\gamma}\tau_{e})$, which we interpret as the Bagnoldian and Newtonian boundaries of a coexistence region in the $(\phi,\dot{\gamma}\tau_{e})$ plane. In Fig. 25 we plot these data and see that the results from $P^{\mathrm{el}}$, $\mu$ and $R_{v}$ reasonably agree with one another. Since we regard $R_{v}$ as the most objective of our measurements of these boundaries, we connect the data points from $R_{v}$ in the figure by solid lines, and denote the enclosed shaded region as the coexistence region of a first-order Bagnoldian-to-Newtonian phase transition. As $\dot{\gamma}\tau_{e}$ decreases, we see that $\phi_{B}$ and $\phi_{N}$ appear to come together at a common $\phi_{c}\approx 0.807$. As $\dot{\gamma}\tau_{e}$ increases, the boundary curves $\phi_{B}(\dot{\gamma}\tau_{e})$ and $\phi_{N}(\dot{\gamma}\tau_{e})$ merge together at a critical end-point (denoted schematically in the figure by the “$*$”). While we have not determined the precise location of this critical end-point, we know that it must occur at a $\dot{\gamma}^{*}\tau_{e}<5\times 10^{-4}$, since the curve of $R_{v}$ vs $\phi$ at $\dot{\gamma}\tau_{e}=5\times 10^{-4}$ gives no evidence of shear banding.

It is also interesting to construct the corresponding phase diagram in the $(P,\dot{\gamma}\tau_{e})$ plane, which we show in Fig. 26. It is interesting to note that, as one increases $\dot{\gamma}$ at fixed $\phi$, one crosses from a region of Bagnoldian rheology to a region of Newtonian rheology. In contrast, as one increases $\dot{\gamma}$ at fixed $P$, it is the reverse; i.e., one crosses from a region of Newtonian rheology to a region of Bagnoldian rheology. The vertical dashed line cutting through the coexistence region at constant $P=2\times 10^{-5}$ in the $(P,\dot{\gamma}\tau_{e})$ plane of Fig. 26, maps onto the dashed line in the $(\phi,\dot{\gamma}\tau_{e})$ plane of Fig. 25.

For a direct demonstration that the shaded region in Figs. 25 and 26 is indeed a region of coexisting Bagnoldian and Newtonian shear bands, we can look explicitly at the spatial variation of packing fraction and strain rate within configurations, as we cross this region. If the interpretation of a coexistence region is correct, we expect that the system will phase separate into shear bands of Bagnoldian and Newtonian rheology in mechanical equilibrium with each other, i.e. the two bands will be at equal values of pressure $p$ and shear stress $\sigma$.

We consider crossing the coexistence region along a contour of constant pressure. Such a contour at fixed dimensionless pressure $P=2\times 10^{-5}$ is shown as the dashed line in Figs. 25 and 26. In Fig. 27 we show a gray-scale encoded snapshot of a sample configuration roughly mid way through the coexistence region on this contour, at the global values $(\phi,\dot{\gamma}\tau_{e})=(0.7984,3\times 10^{-5})$. By “global values” we mean the values averaged over the entire system, or equivalently, the values of $\phi$ and $\dot{\gamma}$ that enter as the control parameters of our simulations. We compute the local value of packing fraction $\phi$ and strain rate $\dot{\gamma}$ of each particle by averaging over a circle of radius $4d_{s}$ centered on each particle. Figure 27(a) shows the resulting local $\phi(x,y)$ throughout the configuration, while Fig. 27(b) shows the local $\dot{\gamma}(x,y)\tau_{e}$. Because our system is so large, with $N=65536$, it is difficult to distinguish individual particles in Fig. 27. However the variation in the gray-scale allows one to visualize both the local average $\phi$ and $\dot{\gamma}\tau_{e}$, as well as fluctuations about the local average, as one varies position within the configuration. We see clearly the Newtonian band in the center of the system, with a larger local $\phi$ and smaller local $\dot{\gamma}\tau_{e}$ as compared to the Bagnoldian band at the edges of the system. For $\phi$, the fluctuations about the local average appear to be roughly uniform throughout both bands. For $\dot{\gamma}\tau_{e}$, the fluctuations seem more pronounced in the Bagnoldian band. Presumbaly this is because in the Newtonian band particles stay in contact for long periods of time, thus reducing fluctuations in local velocity and hence the local strain rate.

To quantify the results illustrated in Fig. 27 we average the behavior of the system along the flow direction $\mathbf{\hat{x}}$, and over many different configurations in our shearing ensemble, to plot quantities as a function of the position $y$ transverse to the flow direction. In Fig. 28 we show results at three different points along the contour of constant $P=2\times 10^{-5}$ within the coexistence region, corresponding to the global values of $(\phi,\dot{\gamma}\tau_{e})=(0.792,5\times 10^{-5})$ (left column), $(0.7984,3\times 10^{-5})$ (middle column), and $(0.802,2\times 10^{-5})$ (right column). The top row shows the average local velocity $v_{x}(y)$. We clearly see that the velocity profile is not linear across the entire width of the system, but is rather piecewise linear with regions of two different slopes (recall, the system obeys Lees-Edwards boundary conditions, in the $y$ direction). The middle row shows the average local strain rate $\dot{\gamma}(y)=dv_{x}/dy$. The bottom row shows the average local packing fraction $\phi(y)$. We see that in each case the system has clearly separated into two bands, one with larger local $\dot{\gamma}\tau_{e}$ and smaller local $\phi$, which we identify as the Bagnoldian band, and one with smaller local $\dot{\gamma}\tau_{e}$ and larger local $\phi$, which we identify as the Newtonian band. As the global $\phi$ increases, the Newtonian band increases in width while the Bagnoldian band shrinks. We have explicitly checked that these shear bands are not a transient effect, but rather represent the steady state; computing the profiles as in Fig. 28 gives essentially the same results whether we average over a finite strain interval near the beginning of our run, or over a non-overlapping strain interval near the end of our run.

Our results thus show the behavior commonly expected at a first order phase transition. Within the coexistence region the system separates into two phases in mutual equilibrium, with one phase growing and the other shrinking as one crosses from one side of the coexistence region to the other. However there is one feature of our results that is not quite the same as for a first order transition. If we denote as $(\phi^{*}_{B},\dot{\gamma}^{*}_{B})$ and $(\phi^{*}_{N},\dot{\gamma}^{*}_{N})$ the local values of the packing fraction and strain rate in the center of the Bagnoldian and Newtonian bands (i.e. away from the interface between the bands) then we see that these values are essentially constant everywhere along the contour of constant $P=2\times 10^{-5}$, with $(\phi^{*}_{B},\dot{\gamma}^{*}_{B})=(0.7880,6.2\times 10^{-5})$ and $(\phi^{*}_{N},\dot{\gamma}^{*}_{N})=(0.8090,4.2\times 10^{-6})$. This is expected because the condition of mechanical stability of the coexisting bands requires that $P=P(\phi^{*}_{B},\dot{\gamma}^{*}_{B}\tau_{e})=P(\phi^{*}_{N},\dot{\gamma}^{*}_{N}\tau_{e})$ and $\Sigma=\Sigma(\phi^{*}_{B},\dot{\gamma}^{*}_{B}\tau_{e})=\Sigma(\phi^{*}_{N},\dot{\gamma}^{*}_{N}\tau_{e})$ must be constant along this contour. For a standard first order phase transition, we would expect that these two values $(\phi^{*}_{B},\dot{\gamma}^{*}_{B}\tau_{e})$ and $(\phi^{*}_{N},\dot{\gamma}^{*}_{N}\tau_{e})$ would correspond to the two points $(\phi_{B},\dot{\gamma}_{B}\tau_{e})$ and $(\phi_{N},\dot{\gamma}_{N}\tau_{e})$ where the contour of constant $P=2\times 10^{-5}$ first enters the coexistence region from the Bagnoldian side and where it exits the coexistence region on the Newtonian side. However we find that this is not quite so; we find that $(\phi_{B},\dot{\gamma}_{B}\tau_{e})=(0.790,5.06\times 10^{-5})$ and $(\phi_{N},\dot{\gamma}_{N}\tau_{e})=(0.8062,5.2\times 10^{-6})$, noticeably different from $(\phi^{*}_{B},\dot{\gamma}^{*}_{B}\tau_{e})$ and $(\phi^{*}_{N},\dot{\gamma}^{*}_{N}\tau_{e})$. In Appendix B we consider this difference further and argue that it is a finite-size-effect resulting from the finite length of our system $L$ as compared to the finite width of the interface between the Bagnoldian and Newtonian shear bands seen in Fig. 28.

Note, the fact that the shear bands in Fig. 28 have constant values of $(\phi^{*}_{B},\dot{\gamma}^{*}_{B}\tau_{e})$ and $(\phi^{*}_{N},\dot{\gamma}^{*}_{N}\tau_{e})$, is only because we are crossing the coexistence region on a contour of constant pressure. If we crossed the coexistence region along any other trajectory, for example at constant $\dot{\gamma}$ or constant $\phi$, the values of $\sigma$ and $p$ within the two coexisting shear bands at any global value of $(\phi,\dot{\gamma}\tau_{e})$ would always be equal; but these values would be varying as we moved through the coexistence region along the trajectory, and hence the values of $(\phi^{*}_{B},\dot{\gamma}^{*}_{B}\tau_{e})$ and $(\phi^{*}_{N},\dot{\gamma}^{*}_{N}\tau_{e})$ in the bands would be varying.

Finally, as is the case for an ordinary first order equilibrium phase transition, we expect that the coexistence region of finite width in the $(\phi,\dot{\gamma}\tau_{e})$ phase diagram of Fig. 25 will collapse to a sharp line when plotted in a phase diagram expressed in terms of the appropriate variables that define the stability condition between the coexisting phases. In the present case this is equality of stress. We therefore consider the phase diagram expressed in the $(P,\Sigma)$ plane, where we may regard the dimensionless shear stress $\Sigma$ as the conjugate variable to $\dot{\gamma}\tau_{e}$, and the dimensionless pressure $P$ as the conjugate variable to $\phi$.

In Fig. 29 we plot the values of $(P,\Sigma)$ along the Bagnoldian and Newtonian boundaries of the coexistence region of Fig. 25, specifically the pairs of points $P(\phi_{B}(\dot{\gamma}\tau_{e}),\dot{\gamma}\tau_{e})$ and $\Sigma(\phi_{B}(\dot{\gamma}\tau_{e}),\dot{\gamma}\tau_{e})$, and the pairs of points $P(\phi_{N}(\dot{\gamma}\tau_{e}),\dot{\gamma}\tau_{e})$ and $\Sigma(\phi_{N}(\dot{\gamma}\tau_{e}),\dot{\gamma}\tau_{e})$. We plot such points for the three different determinations of $\phi_{B}(\dot{\gamma}\tau_{e})$ and $\phi_{N}(\dot{\gamma}\tau_{e})$ obtained from $P^{\mathrm{el}}$, $\mu$, and $R_{v}$. We see that all points indeed collapse to a single sharp line, separating a region of Newtonian rheology from a region of Bagnoldian rheology. Upon continuously crossing this line under conditions of controlled $P$ and $\Sigma$, the system will show a discontinuous jump in shear strain rate $\dot{\gamma}\tau_{e}$ and packing fraction $\phi$. This first-order transition line terminates at a critical end-point, indicated schematically by the “$*$” in the figure. For reference, we also indicate in the figure the yield stress line, defined by $\Sigma=\mu_{J}P$, with $\mu_{J}\approx 0.1$ the macroscopic friction at the jamming transition VOT_RDCD . Below the yield stress line lies the region of mechanically stable, static, jammed solid states. As one approaches the yield stress line from above, the Newtonian rheology just below the first-order transition line will merge continuously with the Herschel-Bulkley rheology at finite $P$ just above the yield line.

III.4 Model CDrot: Discontinuous shear thickening

Discontinuous shear thickening (DST) is the observed phenomenon in which there is a discontinuous jump in the shear stress upon continuously increasing the shear strain rate. In this section we show how the first order rheological transition that gives rise to the shear banding also can provide a possible mechanism for DST.

Consider sitting at a fixed packing fraction (i.e. fixed volume) $\phi<\phi_{c}$ and crossing the coexistence region by continuously increasing the strain rate $\dot{\gamma}$. If $\dot{\gamma}$ is varied sufficiently slowly, so that the system has time to reach the true shear banded steady state at each instantaneous value of $\dot{\gamma}$, then the resulting $\Sigma(\dot{\gamma}\tau_{e})$ will be continuous, with kinks where $\dot{\gamma}$ enters and leaves the coexistence region at $\dot{\gamma}_{B}(\phi)$ and $\dot{\gamma}_{N}(\phi)$ (these are the inverse functions of $\phi_{B}(\dot{\gamma}\tau_{e})$ and $\phi_{N}(\dot{\gamma}\tau_{e}))$. We show such a plot of the steady state dimensionless $\Sigma$ vs $\dot{\gamma}\tau_{e}$ at fixed $\phi=0.802$ in Fig. 30.

However if $\dot{\gamma}$ is varied too rapidly for the system to achieve steady state, then we expect to find the phenomenon generally characteristic of first-order transitions. Upon increasing $\dot{\gamma}$, the system will fall out of the true steady state, overshooting the coexistence region phase boundary to persist in the Bagnoldian phase for some range of $\dot{\gamma}>\dot{\gamma}_{B}(\phi)$. Increasing $\dot{\gamma}$ further, the system will at some point become too far out of the proper steady state and will switch to the Newtonian phase with a resulting jump in $\Sigma$. Decreasing $\dot{\gamma}$, one will see a similar phenomenon, but with hysteresis in the location of the stress jump.

To test this expectation we perform simulations at fixed $\phi$ in which the strain rate $\dot{\gamma}$ is ramped up and down at a fixed rate $\kappa$,

[TABLE]

where $(+)$ is for ramping up, and $(-)$ is for ramping down. We implement this in our simulations by changing $\dot{\gamma}$ at time intervals $\Delta t=\tau_{e}$ according to, $\dot{\gamma}(t+\tau_{e})=\dot{\gamma}(t){\rm e}^{\pm\kappa\tau_{e}}$. For convenience we define the dimensionless rate,

[TABLE]

so that $\dot{\gamma}(t+\tau_{e})=\dot{\gamma}(t)10^{\pm\tilde{\kappa}}$. Each time $\dot{\gamma}$ is updated we also rescale particle velocities,

[TABLE]

and angular velocities,

[TABLE]

so that the change in strain rate is imposed throughout the bulk of the system, rather than only at the boundary. We then average $P$, $\Sigma$ and $\dot{\gamma}\tau_{e}$ over a time interval $10^{3}\tau_{e}$, so as to reduce statistical fluctuations, and thus obtain our plotted results.

We start at an initial low strain rate well below $\dot{\gamma}_{B}$, and increase to a largest strain rate well above $\dot{\gamma}_{N}$. We then ramp back down to the initial strain rate. In Fig. 31 we show results for the dimensionless shear stress $\Sigma$ vs $\dot{\gamma}$ as obtained from this process at the fixed $\phi=0.802$. In Fig. 31(a) we show results for the ramping rate $\tilde{\kappa}=10^{-6}$, while in 31(b) we show results for $\tilde{\kappa}=10^{-5}$. In each case we show results for four independent initial configurations.

We see behavior as described above. Upon increasing $\dot{\gamma}\tau_{e}$ the system stays on the Bagnoldian branch, with $\Sigma\sim\dot{\gamma}^{2}\tau_{e}^{2}$, as $\dot{\gamma}\tau_{e}$ increases above $\dot{\gamma}_{B}\tau_{e}$, before transitioning to the Newtonian branch, with $\Sigma\sim\dot{\gamma}\tau_{e}$, inside the coexistence region. Upon decreasing $\dot{\gamma}\tau_{e}$ the system stays on the Newtonian branch throughout the coexistence region, transitioning to the Bagnoldian branch only below $\dot{\gamma}_{B}\tau_{e}$. The jump upwards in $\Sigma$ upon increasing $\dot{\gamma}\tau_{e}$ is smaller and more gradual than the jump downwards upon decreasing $\dot{\gamma}\tau_{e}$. The location of the jumps are found to vary, without any obvious trend, as the initial configuration is varied, while the width of the hysteresis region shrinks as the ramping rate $\tilde{\kappa}$ decreases. We note that we have observed similar stress jumps in model CD when ramping the strain rate near the Bagnoldian-to-Newtonian transition line $Q_{c}(\phi)$. We thus find that the region of coexisting shear bands in the $(\phi,\dot{\gamma}\tau)$ phase diagram provides a mechanism for the system to manifest DST, which is seen to result when the strain rate varies too rapidly for the system to relax to the shear-banded steady state corresponding to the instantaneous value of $\dot{\gamma}$. Our results in Fig. 31 appear qualitatively similar to numerical results obtained previously for particles with intergranular elastic friction OH_friction ; Grob , except in those works the large $\dot{\gamma}$ branch displayed behavior characteristic of the jammed solid while in our case it is a Newtonian liquid.

IV Conclusions

We have demonstrated that in simple two dimensional models of athermal, soft-core, bidisperse, frictionless massive disks, we have a sharp first-order rheological transition from Bagnoldian rheology to Newtonian rheology as one varies the packing fraction $\phi$ and inelasticity of the collisions $Q$. Key to the existence of this transition is the presence of a viscous dissipation acting in the tangential direction at the point of contact of two particles in collision. Within this framework, we find that the system always becomes Newtonian sufficiently close to the jamming transition $\phi_{J}$, however the behavior of the system at small $Q$ (strongly inelastic) depends dramatically on the precise form of the tangential dissipation. When the dissipation couples only to the translational motion (our model CD) then the system at small $Q$ is always Newtonian for any packing fraction $\phi$. But when the rotational motion is coupled to the translational motion (our model CDrot), then rotational motion can compensate for relative translational motion, and act to effectively eliminate the tangential dissipation. In this case we find that our system has Bagnoldian rheology at small $\phi<\phi_{c}$, and Newtonian rheology at large $\phi>\phi_{c}$.

We find that the main distinguishing feature between systems with Newtonian vs Bagnoldian rheology is that in Newtonian rheology the average number of contacts per particle $Z$ stays finite as the strain rate $\dot{\gamma}\to 0$, while in Bagnoldian rheology $Z\to 0$ as $\dot{\gamma}\to 0$. Thus Newtonian rheology is associated with the formation of large connected clusters. The transition from Bagnoldian to Newtonian rheology upon decreasing $Q$ below $Q_{c}(\phi)$ in either model CD or CDrot results from the increasing inelasticity of the collisions, which act to damp out relative motion; below the critical $Q_{c}(\phi)$, the duration of a collision $\tau_{\mathrm{dur}}$ grows dramatically and particles remain in contact. The nature of the transition from Bagnoldian to Newtonian rheology upon increasing $\phi$ at small $Q$ in model CDrot is somewhat more subtle. At low $\phi$, the particles are able to adjust their rotational motion so that the tangential component of the relative contact velocity at a collision $v_{ij,T}^{C}$ is very small or vanishing, greatly reducing or eliminating the tangential dissipation. As a consequence, the tangential relative motion of the particles’ centers of mass $v_{ij,T}$ is not damped, particles may separate after a collision, so no large persistent clusters form and the rheology stays Bagnoldian. As $\phi$ increases however, we find that the average number of particles participating in any given collision increases, and so we speculate that above a critical $\phi_{c}(Q)$, such multi-particle collisions serve to constrain the rotational motion of the particles in contact, tangential dissipation rises, relative motion is damped, and large persistent clusters form.

We have further studied this Bagnoldian-to-Newtonian rheological transition at small $Q$ in model CDrot to explore the behavior of the system as a function of increasing shear strain rate $\dot{\gamma}$. We find that as $\dot{\gamma}$ increases, the discontinuous first-order rheological transition opens up into a region of finite width in the $(\phi,\dot{\gamma})$ plane consisting of coexisting shear bands of Bagnoldian and Newtonian rheology in mechanical equilibrium with each other. Thus the interplay of inertia and tangential dissipation gives rise to a mechanism for shear banding, even in our simple homogeneous system of frictionless, repulsive, spherical particles.

We have also considered the behavior of our model CDrot as one crosses the above coexistence region by ramping the shear strain rate $\dot{\gamma}$ up and down, holding the system at constant volume. When $\dot{\gamma}$ is ramped at a rate too fast for the system to settle into the true steady-state shear-banded configurations, then one can have hysteretic, discontinuous, jumps in stress as the system transitions between the Bagnoldian and Newtonian phases, similar to such phenomena in first-order equilibrium phase transitions. Thus the shear banding coexistence region can provide a possible mechanism for discontinuous shear thickening, even for particles with no microscopic inter-granular friction Hayakawa2 .

It will be interesting to investigate how robust is the shear-banding phase diagram we have found to the inclusion of other physical effects commonly found in experimental systems of interest. In our models the energy dissipation is taken to occur only via the collisions between particles. In non-Brownian suspensions, emulsions, and colloids, however, the particles are embedded within some fluid host matrix, and so one expects there to be a viscous damping, via Stokes drag, between the particles and the host fluid. There should thus be an additional viscous force $-\zeta[\mathbf{v}_{i}-\mathbf{u}(\mathbf{r}_{i})]$, where $\mathbf{v}_{i}$ is the velocity of particle $i$, and $\mathbf{u}(\mathbf{r}_{i})$ is the local fluid velocity at particle $i$; a similar damping of rotational motion should also be included. Such a viscous drag will result in a rheology that is always Newtonian at sufficiently small $\dot{\gamma}$. However, depending on the value of the viscous drag coefficient $\zeta$, upon increasing $\dot{\gamma}$ inertial effects can become important again and there can be a crossover $\dot{\gamma}_{\times}$ where the system crosses smoothly from Newtonian to Bagnoldian-like behavior Bagnold . Indeed, it has been argued that in granular suspensions $\dot{\gamma}_{\times}$ vanishes linearly as the packing fraction increases towards the jamming $\phi_{J}$ Fall . It thus may be that $\dot{\gamma}_{\times}$ is sufficiently small near our Bagnoldian-to-Newtonian transition at $\phi_{c}$, and so the coexistence region may survive for some finite range of $\dot{\gamma}$.

Another interesting modification of our model would be to include random fluctuating forces to model finite temperature effects in colloids. Such thermal forces should also result in a Newtonian rheology at sufficiently small $\dot{\gamma}$. However, if $\tau_{\mathrm{th}}$ is the thermal relaxation time, one might expect such thermal effects to become unimportant when the shear strain rate is sufficiently large, $\dot{\gamma}>\tau_{\mathrm{th}}^{-1}$ OT-thermal ; Ikeda2 ; Kawasaki . Indeed, the shear bands observed in the colloids of Ref. Chikkadi were exactly in this region of $\dot{\gamma}>\tau_{\mathrm{th}}^{-1}$. Thus, if $\tau_{\mathrm{th}}^{-1}$ near $\phi_{c}$ is sufficiently small, again our coexistence region may survive.

Finally it will be interesting to add inter-particle microscopic friction to our model, and study the interplay between the shear banding we observe in this work, and the DST of the frictional models studied in Refs. OH_friction ; Heussinger ; Grob ; Grob2 . Aside from the possible DST in our present model, as discussed in Sec. III.4, it is worth noting that our shear banding transition shares some of the physical phenomena believed to be at work in the DST of the frictional models. In particular, our shear banding transition marks the transition from a Bagnoldian phase, where there are few instantaneous particle contacts and such contacts are short lived, to a Newtonian phase, where there are $O(1)$ contacts per particle and the contacts are persistent. The DST of frictional models is similarly believed to result from a transition between a region where there are few direct particle contacts, to a region where there are many. It thus may be that our shear banding transition will serve as a trigger for the onset of strong frictional contacts that leads to enhanced DST. We must leave these further investigations, outlined above, to future work.

Acknowledgements

This work was supported in part by the European Research Council under the European Union s Seventh Framework Programme (FP7/2007-2013), ERC Grant Agreement No. 306845. Early stages of this work were supported by National Science Foundation Grant No. DMR-1205800 and the Swedish Research Council Grant No. 2010-3725. Simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at PDC and HPC2N.

Appendix A

In this Appendix we describe the specific methods we use to map out the phase diagram of model CDrot, shown in Fig. 17. To determine the transition line $Q_{c}(\phi)$ at large $\phi\gtrsim 0.8$, we use our result, as discussed in connection with Figs. 14(a) and (b), that for the Bagnoldian region $Q>Q_{c}$ we have $Z\sim\dot{\gamma}\tau_{e}$ as $\dot{\gamma}\tau_{e}\to 0$, while for the Newtonian region $Q<Q_{c}$ we have $Z\to$ constant as $\dot{\gamma}\tau_{e}\to 0$. Thus, for sufficiently small $\dot{\gamma}\tau_{e}$ at fixed $Q$ and $\phi$, we expect $Z/\dot{\gamma}\tau_{e}$ will be independent of $\dot{\gamma}\tau_{e}$ for $Q>Q_{c}$, while $Z/\dot{\gamma}\tau_{e}$ will increase with decreasing $\dot{\gamma}\tau_{e}$ for $Q<Q_{c}$. We thus consider two small strain rates $\dot{\gamma}_{1}\tau_{e}=10^{-6}$ and $\dot{\gamma}_{2}\tau_{e}=10^{-5}$, and compute the difference $(Z/\dot{\gamma}_{1}\tau_{e})-(Z/\dot{\gamma}_{2}\tau_{e})$ as a function of $Q$ at fixed $\phi$. By the preceding observations, this quantity should vanish for $Q>Q_{c}$ and be non-zero for $Q<Q_{c}$. In Fig. 32(a) we plot this quantity vs $Q$ for the case with packing fraction $\phi=0.82$. The point where $(Z/\dot{\gamma}_{1}\tau_{e})-(Z/\dot{\gamma}_{2}\tau_{e})$ increases from zero upon decreasing $Q$ thus locates the transition $Q_{c}$.

To locate the crossover line $Q^{*}(\phi)$ at smaller $\phi\lesssim 0.8$, that separates the weakly inelastic Bagnoldian region $Q>Q^{*}$ from the strongly inelastic Bagnoldian region $Q<Q^{*}$, we use the same approach. As discussed in connection with Figs. 14(c) and (d), for $Q>Q^{*}$ we again have $Z\sim\dot{\gamma}\tau_{e}$. For $Q<Q^{*}$ we found that $Z\to 0$ as $\dot{\gamma}\tau_{e}\to 0$, but that this dependence is slower than linear in $\dot{\gamma}\tau_{e}$. Thus $(Z/\dot{\gamma}_{1}\tau_{e})-(Z/\dot{\gamma}_{2}\tau_{e})$ will again vanish for $Q>Q^{*}$ and be non zero for $Q<Q^{*}$. In Fig. 32(b) we plot this quantity vs $Q$ for the case with $\phi=0.76$. The point where $(Z/\dot{\gamma}_{1}\tau_{e})-(Z/\dot{\gamma}_{2}\tau_{e})$ increases from zero upon decreasing $Q$ thus locates the crossover $Q^{*}$. Proceeding in this manner at other values of fixed $\phi$ we thus map out the curves $Q_{c}(\phi)$ and $Q^{*}(\phi)$ in the phase diagram of Fig. 17.

To determine the transition $\phi_{c}(Q)$ at small $Q$, that separates the Bagnoldian region $\phi<\phi_{c}$ from the Newtonian region at $\phi>\phi_{c}$, we consider the Bagnold coefficient $B_{p}^{\mathrm{el}}$. In the Bagnoldian region at sufficiently small $\dot{\gamma}\tau_{e}$ we have $B_{p}^{\mathrm{el}}$ is independent of $\dot{\gamma}\tau_{e}$. In the Newtonian region, where $\eta_{p}^{\mathrm{el}}$ is independent of $\dot{\gamma}\tau_{e}$, we have from Eqs. (19) and (20) that $B_{p}^{\mathrm{el}}\equiv P^{\mathrm{el}}/(\dot{\gamma}\tau_{e})^{2}=\eta_{p}^{\mathrm{el}}/(Q\dot{\gamma}\tau_{e})\sim 1/\dot{\gamma}\tau_{e}$. Thus if we compute the difference $B_{p1}^{\mathrm{el}}-B_{p2}^{\mathrm{el}}$, where the two Bagnold coefficients are computed at $\dot{\gamma}_{1}\tau_{e}=10^{-6}$ and $\dot{\gamma}_{2}\tau_{e}=10^{-5}$ respectively, then this quantity will vanish in the Bagnoldian region $\phi<\phi_{c}$, and become non-zero in the Newtonian region $\phi>\phi_{c}$. In Fig. 33 we show examples of this analysis for the cases $Q=1$ and $Q=20$. The point where $B_{p1}^{\mathrm{el}}-B_{p2}^{\mathrm{el}}$ increases from zero upon increasing $\phi$ thus locates the transition $\phi_{c}$. Proceeding in this manner at other values of fixed $Q$ we thus map out the curve $\phi_{c}(Q)$ to complete the phase diagram of Fig. 17. We caution that the precise location of the transition lines in this phase diagram may be expected to shift slightly if one considers smaller strain rates $\dot{\gamma}\tau_{e}$ and larger system sizes.

Appendix B

In this Appendix we consider the contour of constant pressure $P=2\times 10^{-5}$ that passes through the coexistence region of model CDrot (the dashed line in Fig. 25), and address the issue of why the values of $(\phi^{*}_{B},\dot{\gamma}^{*}_{B}\tau_{e})$ and $(\phi^{*}_{N},\dot{\gamma}^{*}_{N}\tau_{e})$, that describe the local packing fraction and strain rate of the two shear bands in Fig. 28, are found to be slightly different from the values $(\phi_{B},\dot{\gamma}_{B}\tau_{e})$ and $(\phi_{N},\dot{\gamma}_{N}\tau_{e})$ at which the contour enters and exits the coexistence region. We will argue that this difference is due to a finite-size-effect.

To highlight the difference between $(\phi^{*}_{B},\dot{\gamma}^{*}_{B}\tau_{e})$ and $(\phi_{B},\dot{\gamma}_{B}\tau_{e})$, and between $(\phi^{*}_{N},\dot{\gamma}^{*}_{N}\tau_{e})$ and $(\phi_{N},\dot{\gamma}_{N}\tau_{e})$, we consider the behavior of the macroscopic friction $\mu=\sigma/p=\Sigma/P$ as we pass through the coexistence region on the contour of constant pressure. In Fig. 34(a) we plot $\mu$ vs the global $\phi$ along this contour of constant $P=2\times 10^{-5}$. In this figure circles represent points in the homogeneous Bagnoldian phase, triangles are points in the homogeneous Newtonian phase, while squares are points in the shear banded coexistence region. We define $\phi_{B}$ as the highest packing fraction in the homogeneous Bagnoldian phase, and $\phi_{N}$ as the lowest packing fraction in the homogeneous Newtonian phase, as indicated in the figure.

Since we are varying $\phi$ at constant $P$, and since in the coexistence region the shear bands must be in mechanical equilibrium with equal $P$ and $\Sigma$, we thus expect that $\mu=P/\Sigma$ must also remain constant as we pass through the coexistence region at constant $P$. We indeed see in Fig. 34(a) that this is so. However this constant $\mu=0.346$ is not the value of $\mu$ at either $\phi_{B}$ or $\phi_{N}$, but rather we have to enter deeper into the homogeneous phase regions to $\phi^{*}_{B}$ and $\phi^{*}_{N}$, as indicated in the figure, to find homogeneous states with this value of $\mu=0.346$. In Fig. 34(b) we make a similar plot of $\mu$ vs the global $\dot{\gamma}\tau_{e}$, along this same contour of constant $P$. Fig. 34 thus allows us to read off the above defined values of $(\phi^{*}_{B},\dot{\gamma}^{*}_{B}\tau_{e})$ and $(\phi^{*}_{N},\dot{\gamma}^{*}_{N}\tau_{e})$, and we find to these to agree well with the values found in Fig. 28 for the local packing fractions and strain rates within the two coexisting bands. We can also read off from Fig. 34 the values of $(\phi_{B},\dot{\gamma}_{B}\tau_{e})$ and $(\phi_{N},\dot{\gamma}_{N}\tau_{e})$ which give the coexistence region boundaries. We thus see that $\mu(\phi_{B},\dot{\gamma}_{B}\tau_{e})<\mu(\phi^{*}_{B},\dot{\gamma}^{*}_{B}\tau_{e})=\mu(\phi^{*}_{N},\dot{\gamma}^{*}_{N}\tau_{e})<\mu(\phi_{N},\dot{\gamma}_{N}\tau_{e})$, and there appears to be small but sharp jumps in $\mu$ (and hence $\Sigma$) as we cross from the homogeneous regions into the coexistence region.

By sampling other points throughout the shear banded coexistence region, and measuring the packing fraction and strain rate in the two coexisting phases, we can determine the values of $(\phi^{*}_{B},\dot{\gamma}^{*}_{B}\tau_{e})$ and $(\phi^{*}_{N},\dot{\gamma}^{*}_{N}\tau_{e})$ corresponding to other fixed pressures. We find that these values all lie on common curves $\phi^{*}_{B}(\dot{\gamma}\tau_{e})$ and $\phi^{*}_{N}(\dot{\gamma}\tau_{e})$ within the $(\phi,\dot{\gamma}\tau_{e})$ plane, as we show in Fig. 35. We see that the difference between $\phi_{N}(\dot{\gamma}\tau_{e})$ and $\phi_{N}^{*}(\dot{\gamma}\tau_{e})$ is greatest at low $\dot{\gamma}\tau_{e}$, while the difference between $\phi_{B}(\dot{\gamma}\tau_{e})$ and $\phi_{B}^{*}(\dot{\gamma}\tau_{e})$ is greatest near the critical end point.

Returning to Fig. 34, we note that the constant $\mu$ within the shear banding coexistence region is analogous to the Maxwell construction of a mean-field first-order equilibrium phase transition Max . However for a first-order equilibrium transition, we would expect this curve to be continuous and monotonic. Defining $f$ as the fraction of the system in the Bagnoldian phase, we would have expected that $f$ varies smoothly from $f=1$ to $f=0$ as the system crosses from the Bagnoldian side to the Newtonian side. However the difference between $\phi_{N}$ and $\phi_{N}^{*}$, and between $\phi_{B}$ and $\phi_{B}^{*}$ shows that this does not happen. Since we must have for the global average $\phi=f\phi_{B}^{*}+(1-f)\phi_{N}^{*}$, and since $\phi_{B}^{*}<\phi_{B}$, it must be the case that $f<1$ as the system enters the coexistence region at $\phi_{B}$, and similarly since $\phi_{N}^{*}>\phi_{N}$, we must have $f>0$ when the system exists the coexistence region at $\phi_{N}$.

To interpret this unexpected result, we speculate that $\phi^{*}_{B}(\dot{\gamma}\tau_{e})$ and $\phi^{*}_{N}(\dot{\gamma}\tau_{e})$ may be the true boundaries of the coexistence region, and that the stability of the homogeneous phase within the region between $\phi^{*}_{B}(\dot{\gamma}\tau_{e})$ and $\phi_{B}(\dot{\gamma}\tau_{e})$ and between $\phi_{N}(\dot{\gamma}\tau_{e})$ and $\phi_{N}^{*}(\dot{\gamma}\tau_{e})$ is a finite size effect, explained as follows. In the coexistence region, the width of the interface over which the system transitions from the Bagnoldian band to the Newtonian band should have some characteristic size set by a correlation length $\xi$ that remains finite as long as one is not at the critical end point. The widths of the bands themselves are set by the system width $L_{y}$. As $L_{y}\to\infty$, the width of the interface always becomes negligible compared to the widths of the bands. However, for finite $L_{y}$, as one approaches close to the coexistence line from within the shear banded region, the width of one of the bands is decreasing and will ultimately become of order $\xi$. At that point the decrease in free energy (assuming our athermal, non-equilibrium, system has something analogous to a free energy) obtained by phase separating may be smaller than the increase in free energy needed to create the interface between the bands. The homogeneous phase may thus become stabilized. Looking at the profiles of $\dot{\gamma}(y)$ and $\phi(y)$ in Fig. 28, we see that the width of the interface is generally non-negligible compared to the width $L_{y}$ of our system, and so such a finite size effect in our system may be large and noticeable.

We now test this notion by checking whether the boundaries $\phi_{B}(\dot{\gamma}\tau_{e})$ and $\phi_{N}(\dot{\gamma}\tau_{e})$ show any noticeable dependence on the system size. In Fig. 36 we plot the shear banding parameter $R_{v}$ of Eq. (37) vs $\phi$, comparing results for our system with $N=65536$ particles against a system of half the length, with $N=16384$ particles. We show results for two different strain rates, $\dot{\gamma}\tau_{e}=10^{-5}$ and $10^{-4}$. For $\dot{\gamma}\tau_{e}=10^{-5}$, for which Fig. 35 shows essentially no difference between $\phi_{B}$ and $\phi^{*}_{B}$ but a noticeable difference between $\phi_{N}$ and $\phi^{*}_{N}$, we see that there is little finite size effect for $\phi_{B}$, but a strong finite size effect for $\phi_{N}$. For $\dot{\gamma}\tau_{e}=10^{-4}$, for which Fig. 35 shows a noticeable difference between both $\phi_{B}$ and $\phi^{*}_{B}$ and between $\phi_{N}$ and $\phi^{*}_{N}$, we see a noticeable finite size effect in the location of both $\phi_{B}$ and $\phi_{N}$. For both strain rates, we find that the width of the coexistence region shrinks as the system size gets smaller. We thus find, in support of our above speculation, that finite size effects appear in the locations of $\phi_{B}$ and $\phi_{N}$ wherever $\phi_{B}>\phi^{*}_{B}$ and $\phi_{N}<\phi^{*}_{N}$. We thus suppose that we would find $\phi_{N}(\dot{\gamma}\tau_{e})\to\phi^{*}_{N}(\dot{\gamma}\tau_{e})$, and $\phi_{B}(\dot{\gamma}\tau_{e})\to\phi^{*}_{B}(\dot{\gamma}\tau_{e})$, if we could look at systems with much larger $L_{y}$.

We note that another manifestation of such a finite size effect may be directly seen in Fig. 28. We can determine the width $w_{B}$ of the Bagnoldian band by the intersection of the packing fraction profile $\phi(y)$ with the global average $\phi$, so that the fraction of the system in the Bagnoldian band is $f=w_{B}/L_{y}$. Such a definition gives the correct relation $\phi=f\phi^{*}_{B}+(1-f)\phi^{*}_{N}$. We can then proceed similarly using the strain rate profile $\dot{\gamma}(y)$, to define an $f^{\prime}$ so that $\dot{\gamma}=f^{\prime}\dot{\gamma}^{*}_{B}+(1-f^{\prime})\dot{\gamma}^{*}_{N}$. We then find that, although $f$ and $f^{\prime}$ are close, they are not equal. This can be verified by comparing the location of the vertical dashed lines in Figs. 28 (d) – (f) with those in (g) – (i). The difference between $f$ and $f^{\prime}$ is a consequence of the fact that the width of the interface is a non-negligible fraction of the system width $L_{y}$.