A microscopic theory for discontinuous shear thickening of frictional granular materials

We extend a recent theory for the rheology of frictionless granular materials [K. Suzuki and H. Hayakawa, Phys. Rev. Lett. 2015, 115, 098001] to the case of frictional disks in two dimensions. Employing a frictional contact model for molecular dynamics simulations, we derive difference equations of the shear stress, the granular temperature, and the spin temperature from the generalized Green-Kubo formula, where all the terms are given by microscopic expressions. The numerical solutions of the difference equations not only describe the flow curve, but also reproduce the hysteresis of shear stress, which can be the signature of discontinuous shear thickening of frictional disks.


Introduction
Understanding of the rheology of granular materials is important for industrial applications so that its microscopic descriptions have been desired in science and technology [1]. It is now known that the rheology drastically changes if the fraction exceeds the jamming point, i.e. the shear stress changes from the liquid-to solid-branches of flow curves. Besides the strong dependence on the fraction, the friction between constituent particles also triggers the anomalous flow behavior, i.e. the discontinuous shear thickening (DST), where the shear stress discontinuously jumps from liquid-to solid-branches, and vice versa [2][3][4]. Recently, a microscopic theory for the rheology of frictionless granular materials have been developed by Suzuki and one of authors [5], where the critical scaling of shear viscosity and that of granular temperature below jamming were well predicted without any fitting parameters.
In this paper, we apply their theory to the case of frictional disks in two-dimension such that the DST observed in molecular dynamics (MD) simulations [2,4] can be explained.

Theory
First, we describe translational motions of frictional disks by the SLLOD equation [6]. We assume that the system is homogeneously sheared along the x-axis under the Lees-Edwards boundary condition. Then, time derivatives of the position, r i = (x i , y i ), and peculiar momentum, p i = (p ix , p iy ), of i-th particle (i = 1, . . . , N) are given bẏ respectively, where m i and F i = (F ix , F iy ) are the mass of i-th particle and the force acts on the i-th particle, respectively. On the right-hand-sides of Eq. (1), each element of the two-rank tensor,γ, is defined asγ αβ =γδ xα δ yβ (α, β = x, y) with the shear rate,γ, such thatγ : r i = (γy i , 0). Next, rotational motions of the disks are described by the Euler equation of motions, where time derivatives of the angular position, ϕ i , and angular velocity, ω i , of i-th particle are given bẏ respectively. In Eq. (2), I i and M i represent the moment of inertia and the torque acts on the i-th particle, respectively. From Eqs. (1) and (2), the phase variables are defined as

Frictional contact model
To study the DST of frictional disks, we employ a frictional contact model [7], where the total force acts on the i-th particle is divided into the normal and tangential directions as F i = F in + F it (the subscripts, n and t, are, respectively, used for the normal and the tangential components). Here, both the normal and tangential forces are modeled by linear spring-dashpot and are further decomposed into (reversible) elastic and (irreversible) viscous parts as it , respectively. For later use, we also introduce the elastic and viscous parts of the total force as F i = F (el) In the MD simulations of frictional disks, the DST is enhanced for strong friction (the friction coefficient, µ = 2, is used for Coulomb's friction [2,4]). Therefore, in this study, we neglect the switch from static friction, F it , to dynamical one, µ|F in |, where all the disks in contact do not slip with each other.
In the normal force, F in , the elastic and viscous parts are given by F (el) in = k n j i Θ(ξ i j )ξ i j n i j and F (vis) in = −η n j i Θ(ξ i j )(ṙ i j · n i j )n i j , respectively, where k n and η n are introduced as a spring constant and viscosity coefficient, respectively. Here, Θ(ξ i j ) is the Heaviside step function (which is 1 for ξ i j > 0 and 0 otherwise), where ξ i j ≡ R i +R j −r i j , R i (or R j ), r i j ≡ |r i −r j |, n i j ≡ (r i −r j )/r i j , andṙ i j ≡ṙ i −ṙ j represent an overlap between the disks (i and j), the radius of i-th (or j-th) disk, the distance between the disks, the unit vector parallel to the relative position, and the relative velocity, respectively.
In the tangential force, F it , the elastic and viscous parts are written as respectively, where k t and η t are the spring constant and the viscosity coefficient, respectively, and the disks, i and j, get into contact at timẽ t i j . Here, v i jt = n i j × (R i ω i + R j ω j ) +ṙ i j − (ṙ i j · n i j )n i j represents the tangential relative velocity at the contact point [7], which is perpendicular to the normal unit vector between the disks, i.e. v i jt · n i j = 0. Then, the elastic part, F (el) it , is proportional to the relative displacement, The total torque acts on the i-th particle also consists of elastic and viscous parts as

Hamiltonian and Liouville equation
In our system, the Hamiltonian is introduced as H = where U n and U t are elastic potentials stored in the normal and tangential directions, respectively. Because the Hamiltonian depends on time through the phase variable, Γ, its time derivative is given byḢ (1) and (2), the time derivative is rewritten asḢ = −γS σ xy − 2R + J. Here, are introduced as the shear stress 1 , the dissipation function, and the rate of energy production by the couple stress [8], respectively, where S represents the system area.
If we introduce a dynamical variable, A(t), as a function of phase variable, Γ, its time derivative is given by 1 It is also written as σ xy = i m −1 where iL ≡Γ · (∂/∂Γ) is defined as the Liouville operator [6]. From Eqs. (1) and (2), the Liouville operator is decomposed as iL = iL el + iLγ + iL vis 2 . On the other hand, the time development of N-body distribution function, f (t), is described by the Liouville equation as where Λ ≡ (∂/∂Γ) ·Γ is defined as the compression factor [6]. The Liouville equation is also written as where iL † ≡ iL + Λ is the adjoint of Liouville operator so that the Liouville operator is non-Hermitian, i.e. iL † iL [6]. From Eqs. (1) and (2), we find that the compression factor is given by which only consists of irreversible viscous parts, F (vis) i and M (vis) i . Then, it is readily found that the first and second terms on the right-hand-side of Eq. (8) are reduced to is introduced as the coordination number of the i-th particle. If we use the moment of inertia for a two-dimensional disk, I i = m i R 2 i /2, and assume that every mass is identical, m i = m, we find that the compression factor is proportional to the total number of contacts, i.e. Λ = −(η n + 3η t )N c /m with N c ≡ i z i .

Generalized Green-Kubo formula
To describe the shear stress of frictional disks, let us remind the measurement of flow curves [2]: After the system is equilibrated under shear, the shear rate is increased fromγ toγ + ∆γ (at time t = 0) such that the shear stress relaxes from σ xy γ to σ xy γ+∆γ , where . . . γ represents the average in a non-equilibrium steady state under shear withγ. The generalized Green-Kubo formula around the steady state characterized byγ is expressed as where Ω ≡ βḢ − Λ = −γβS σ xy − 2βR + βJ − Λ with the inverse of granular temperature, β ≡ 1/ T γ , is the work function [6,9] and σ xy γ 0. Assuming an exponential decay of the time correlation function, σ xy (s)Ω γ ≃ σ xy Ω γ e −s/τ rel , with a relaxation time, τ rel , and taking a long time limit (t → ∞), we reduce Eq. (9) to σ xy γ+∆γ ≃ σ xy γ + τ rel σ xy Ω γ , where we have rewritten the steady state value as σ xy γ+∆γ ≡ σ xy (∞) γ+∆γ (see Ref. [5] for more rigorous treatment). If we assume that the non-equilibrium steady state is not far from the equilibrium state, we can rewrite the correlation function as 2 Each part is explicitly written as iL el = i {m −1 a 0 = 2τ rel k 2 n i, j,k,l (i j,k l) Θ(ξ i j )Θ(ξ kl )ξ i j ξ kl r i j r kl n i jx n klx n i jy n kly γ a 1 = (mS ) −1 τ rel i, j,k (i j,i k) 3η 2 n Θ(ξ i j ) 2 r 2 i j n 2 i jx n 2 i jy + 3η 2 t Θ(ξ i j ) 2 r 2 i j n 4 i jy + η 2 n Θ(ξ i j )Θ(ξ ik ) n i jx n ikx + n i jy n iky r i j r ik n i jx n ikx n i jy n iky +η n η t Θ(ξ i j )Θ(ξ ik )r i j r ik n i jx n i jy n 2 iky n i jx n iky − n i jy n ikx + η 2 t Θ(ξ i j )Θ(ξ ik )r i j r ik n 2 i jy n 2 iky n i jy n iky + n i jx n ikx γ a 2 = (mS ) −1 τ rel i, j,k (i j,i k) 6η 2 t Θ(ξ i j ) 2 r 2 i j n 4 i jy + 2η 2 t Θ(ξ i j )Θ(ξ ik )r i j r ik n 2 i jy n 2 iky γ b 0 = 2τ rel k 2 n i, j,k,l (i j,k l) Θ(ξ i j )Θ(ξ kl )ξ i j ξ kl r i j r kl n i jx n klx n i jy n kly γ n Θ(ξ i j ) 2 r 2 i j n 2 i jx n 2 i jy + 3η 2 t Θ(ξ i j ) 2 r 2 i j n 4 i jy + η 2 n Θ(ξ i j )Θ(ξ ik ) n i jx n ikx + n i jy n iky r i j r ik n i jx n ikx n i jy n iky +η n η t Θ(ξ i j )Θ(ξ ik )r i j r ik n i jx n i jy n 2 iky n i jx n iky − n i jy n ikx + η 2 t Θ(ξ i j )Θ(ξ ik )r i j r ik n 2 i jy n 2 iky n i jy n iky + n i jx n ikx Θ(ξ i j )Θ(ξ kl )ξ i j ξ kl r i j r kl n i jx n klx n i jy n kly γ n Θ(ξ i j ) 2 r 2 i j n 2 i jx n 2 i jy + 3η 2 t Θ(ξ i j ) 2 r 2 i j n 4 i jy + η 2 n Θ(ξ i j )Θ(ξ ik ) n i jx n ikx + n i jy n iky r i j r ik n i jx n ikx n i jy n iky +η n η t Θ(ξ i j )Θ(ξ ik )r i j r ik n i jx n i jy n 2 iky n i jx n iky − n i jy n ikx + η 2 t Θ(ξ i j )Θ(ξ ik )r i j r ik n 2 i jy n 2 iky n i jy n iky + n i jx n ikx γ iky γ σ xy Ω γ ≃ σ xy Ω eq + τ rel σ xy Ω 2 eq , where . . . eq represents the ensemble average in the equilibrium state. Then, we find that the shear stress is expanded into the series of relaxation time as σ xy γ+∆γ ≃ σ xy γ + τ rel σ xy Ω eq + τ 2 rel σ xy Ω 2 eq . The parity of phase variables requires σ xy Ω 2 eq = 0 [5] and finally, the shear stress is given by σ xy γ+∆γ ≃ σ xy γ + τ rel σ xy Ω eq , where the equilibrium averages, σ xy Ω eq = −γβS σ 2 xy eq − 2β σ xy R eq + β σ xy J eq − σ xy Λ eq , can be explicitly calculated from Eqs. (3)-(5) and (8). Applying the parallel procedure to the granular temperature, T γ , and rotational temperature, T t γ 3 , we find difference equations as where the coefficients (a 0 , a 1 , a 2 , b 0 , b 1 , b 2 , c 0 , c 1 , and c 2 ) are listed in Table 1 4 . To explicitly calculate the coefficients in Table 1, we need to calculate the three-and four-body correlations [5] and have to introduce the radial distribution function for frictional disks. In addition, it is necessary to determine the two time scales, i.e. the relaxation time, τ rel , and duration of contact, τ i j , which is far beyond the scope of this work. Therefore, in the following analysis, we simply regard the coefficients as fitting parameters. 3 The granular temperature and rotational temperature are introduced as T ≡ N −1 i p 2 i /2m and T t ≡ N −1 i I i ω 2 i , respectively. 4 To calculate the correlation functions, the elastic part of tangential

Summary
In this study, we have extended the previous microscopic theory for frictionless granular particles [5] to the case of two-dimensional frictional disks, where all the flow behavior, i.e. Bagnold's scaling, σ * xy γ ∼γ * 2 , the plastic flow, σ * xy γ ∼γ * , the rate-independent yielding stress, σ * xy γ ∼ σ Y , and the hysteresis of shear stress as the signature of DST, is well described by the difference equations (10)-(12). All the terms in Eqs. (10)-(12) including the coefficients have been derived from the generalized Green-Kubo formula, Eq. (9), where the canonical ensembles are assumed for the system under shear and the time correlation function between the shear stress and work function, i.e. σ xy (s)Ω γ , is approximated by the exponential function with the relaxation time, τ rel . Although it is needed to introduce the radial distribution function for frictional disks and determine the two time scales, i.e. the relaxation time and duration of contact, the coefficients in the difference equations have been provided by the microscopic expressions as listed in Table 1. However, we have used the coefficients as numerical parameters to draw the flow curves of frictional disks (Fig. 1). From our results, we can expect that the coupling with rotational temperature, which is absent in the case of frictionless particles [5], triggers both the transition from the liquid-like branch to solid-like one and the hysteresis behavior of shear stress.
Nevertheless, the transition from the liquid-like behavior (unjammed state) to the solid-like branch (jammed state) during the increase of shear rate (the open circles in Fig. 1) is much smoother than that observed in MD simulations [2,4]. In addition, we have not observed a sudden drop from the solid-like branch to the liquid-like one (the open squares in Fig. 1). Therefore, we will need further investigations of our microscopic approach to the DST, where more rigorous treatments of the coefficients, time scales, and the radial distribution function might be crucial. In addition, it is important to implement the recent theory of non-Brownian suspensions, where the competing divergences of shear viscosity with the smooth shift of jamming point explain the discontinuous shear thickening [11].
In conclusion, we have developed a microscopic theory for the DST of two-dimensional frictional disks, where the difference equations of the shear stress, the granular temperature, and the rotational temperature derived from the generalized Green-Kubo formula well explain all the flow behavior including the hysteresis of the shear stress.