2.1 Microscopic basic equations

We consider an assembly of $N$ frictionless monodisperse spherical particles of diameter $d$ contained in a box of volume $V$ and immersed in a liquid of viscosity $\unicode[STIX]{x1D702}_{0}$ . A simple steady shear with shear rate $\dot{\unicode[STIX]{x1D6FE}}$ is applied to the system. The coordinate is chosen such that the flow is in the $x$ -direction and the velocity gradient is in the $y$ -direction. We consider the overdamped equation of motion

(2.1) $$\begin{eqnarray}\displaystyle \mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D701}_{ij}^{(N)}(\dot{\boldsymbol{r}}_{j}-\dot{\unicode[STIX]{x1D6FE}}y_{j}\boldsymbol{e}_{x})=\boldsymbol{F}_{i}^{(p)}\quad (i=1,\ldots ,N), & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{r}_{i}$ and $\dot{\boldsymbol{r}}_{i}$ are the position and velocity of particle $i$ , respectively, $\boldsymbol{e}_{x}$ is the unit vector in the $x$ -direction, $\boldsymbol{F}_{i}^{(p)}$ is the interparticle force exerted on particle $i$ from other particles and $\{\unicode[STIX]{x1D701}_{ij}^{(N)}\}\text{}_{i,j=1}^{N}$ is the resistance matrix of the suspension, which depends on the configuration of the particles, $\{\boldsymbol{r}_{i}\}\text{}_{i=1}^{N}$ . Note that $\{\unicode[STIX]{x1D701}_{ij}^{(N)}\}\text{}_{i,j=1}^{N}$ is a $3N\times 3N$ matrix, where each component $\unicode[STIX]{x1D701}_{ij}^{(N)}$ is a $3\times 3$ matrix. In particle suspensions, the inertia of the particles is absorbed by the background fluid and hence insignificant. Thus we neglect it in (2.1). We also neglect the rotation of the particles and the thermal fluctuating force exerted on the particles from the solvent in (2.1).

The time evolution of an arbitrary observable $A(\unicode[STIX]{x1D71E})$ is determined by the Liouville equation

(2.2) $$\begin{eqnarray}\displaystyle {\dot{A}}(\unicode[STIX]{x1D71E}(t))=\dot{\unicode[STIX]{x1D71E}}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D71E}}A(\unicode[STIX]{x1D71E}(t)):=\text{i}{\mathcal{L}}A(\unicode[STIX]{x1D71E}(t)), & & \displaystyle\end{eqnarray}$$

where $\text{i}{\mathcal{L}}$ is the Liouvillian. In simple shear flows, $\unicode[STIX]{x1D71E}$ is given by $\unicode[STIX]{x1D71E}=\{\boldsymbol{r}_{i},\boldsymbol{v}_{i}\}\text{}_{i=1}^{N}$ , where

(2.3) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{i}:=\dot{\boldsymbol{r}}_{i}-\dot{\unicode[STIX]{x1D6FE}}y_{i}\boldsymbol{e}_{x}=\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D701}_{ij}^{(N)-1}\boldsymbol{F}_{j}^{(p)} & & \displaystyle\end{eqnarray}$$

is the peculiar velocity, which is the velocity in the sheared frame. For (2.1), $\text{i}{\mathcal{L}}$ is given by

(2.4) $$\begin{eqnarray}\displaystyle \text{i}{\mathcal{L}}:=\mathop{\sum }_{i=1}^{N}\left(\boldsymbol{v}_{i}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{r}_{i}}+\dot{\boldsymbol{v}}_{i}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{v}_{i}}\right)=\mathop{\sum }_{i=1}^{N}\left(\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D701}_{ij}^{(N)-1}\boldsymbol{F}_{j}^{(p)}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{r}_{i}}+\dot{\boldsymbol{v}}_{i}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{v}_{i}}\right), & & \displaystyle\end{eqnarray}$$

where $\dot{\boldsymbol{v}}_{i}$ is evaluated as

(2.5) $$\begin{eqnarray}\displaystyle \dot{\boldsymbol{v}}_{i}=\ddot{\boldsymbol{r}}_{i}-\dot{\unicode[STIX]{x1D6FE}}{\dot{y}}_{i}\boldsymbol{e}_{x}=\ddot{\boldsymbol{r}}_{i}-\dot{\unicode[STIX]{x1D6FE}}\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D701}_{ij}^{(N)-1}F_{j,y}^{(p)}\boldsymbol{e}_{x}. & & \displaystyle\end{eqnarray}$$

Note that (2.3) is utilized in the second equality of (2.5). In the formulation of overdamped Liouville equation, we neglect $\ddot{\boldsymbol{r}}_{i}$ in (2.5), because $\boldsymbol{v}_{i}$ rather than $\ddot{\boldsymbol{r}}_{i}$ resides in the equation of motion, equation (2.3). This leads us to the expression

(2.6) $$\begin{eqnarray}\displaystyle \text{i}{\mathcal{L}}=\mathop{\sum }_{i=1}^{N}\left(\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D701}_{ij}^{(N)-1}\boldsymbol{F}_{j}^{(p)}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{r}_{i}}-\dot{\unicode[STIX]{x1D6FE}}F_{i,y}^{(p)}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}F_{i,x}^{(p)}}\right). & & \displaystyle\end{eqnarray}$$

Note, however, that neglecting $\ddot{\boldsymbol{r}}_{i}$ in (2.5) does not mean the absence of $\dot{\boldsymbol{v}}_{i}$ as can be seen in (2.5). In fact, because of the existence of $\boldsymbol{v}_{i}$ , trajectories of the particles can be curved.

The Liouville equation of the microscopic stress tensor $\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\unicode[STIX]{x1D71E})$ ( $\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}=x,y,z$ ) reads

(2.7) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\unicode[STIX]{x1D71E})=\mathop{\sum }_{i=1}^{N}\left(\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D701}_{ij}^{(N)-1}\boldsymbol{F}_{j}^{(p)}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\unicode[STIX]{x1D71E})}{\unicode[STIX]{x2202}\boldsymbol{r}_{i}}-\dot{\unicode[STIX]{x1D6FE}}F_{i,y}^{(p)}\frac{\unicode[STIX]{x2202}\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\unicode[STIX]{x1D71E})}{\unicode[STIX]{x2202}F_{i,x}^{(p)}}\right), & & \displaystyle\end{eqnarray}$$

where $\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\unicode[STIX]{x1D71E})$ is given by

(2.8) $$\begin{eqnarray}\displaystyle \tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\unicode[STIX]{x1D71E}):=-\frac{1}{2V}\mathop{\sum }_{i=1}^{N}(r_{i,\unicode[STIX]{x1D6FD}}F_{i,\unicode[STIX]{x1D6FC}}^{(p)}+r_{i,\unicode[STIX]{x1D6FC}}F_{i,\unicode[STIX]{x1D6FD}}^{(p)}). & & \displaystyle\end{eqnarray}$$

In simple shear flows, the only non-zero components of $\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\unicode[STIX]{x1D71E})$ are $\tilde{\unicode[STIX]{x1D70E}}_{xy}(\unicode[STIX]{x1D71E})$ and the diagonal ones, from which the microscopic pressure is given by $\tilde{P}(\unicode[STIX]{x1D71E})=-(\tilde{\unicode[STIX]{x1D70E}}_{xx}(\unicode[STIX]{x1D71E})+\tilde{\unicode[STIX]{x1D70E}}_{yy}(\unicode[STIX]{x1D71E})+\tilde{\unicode[STIX]{x1D70E}}_{zz}(\unicode[STIX]{x1D71E}))/3$ . Note that we mainly consider the Cauchy stress, which contributes to the divergence at $\unicode[STIX]{x1D711}\approx \unicode[STIX]{x1D711}_{J}$ , but it is possible to define the kinetic stress by the peculiar velocity.

2.2 Exact equations for the stress

Macroscopic equation of continuity of the stress tensor is obtained by multiplying (2.7) by the non-equilibrium distribution function $f(\unicode[STIX]{x1D71E},t)$ and integrating over $\unicode[STIX]{x1D71E}$ ,

(2.9) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\text{d}}{\text{d}t}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}+\frac{1}{2}\dot{\unicode[STIX]{x1D6FE}}(\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FC}x}\unicode[STIX]{x1D70E}_{y\unicode[STIX]{x1D6FD}}+\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FD}x}\unicode[STIX]{x1D70E}_{y\unicode[STIX]{x1D6FC}})=-\frac{1}{2V}\mathop{\sum }_{i}\left\langle \mathop{\sum }_{j}\unicode[STIX]{x1D701}_{ij}^{(N)-1}F_{j,\unicode[STIX]{x1D6FD}}^{(p)}F_{i,\unicode[STIX]{x1D6FC}}^{(p)}\right\rangle +(\unicode[STIX]{x1D6FC}\leftrightarrow \unicode[STIX]{x1D6FD})\nonumber\\ \displaystyle & & \displaystyle \quad -\,\frac{1}{2V}\sum _{i,j}\left\langle \mathop{\sum }_{k}\unicode[STIX]{x1D701}_{ik}^{(N)-1}F_{k,\unicode[STIX]{x1D706}}^{(p)}r_{j,\unicode[STIX]{x1D6FD}}\frac{\unicode[STIX]{x2202}F_{j,\unicode[STIX]{x1D6FC}}^{(p)}}{\unicode[STIX]{x2202}r_{i,\unicode[STIX]{x1D706}}}\right\rangle +(\unicode[STIX]{x1D6FC}\leftrightarrow \unicode[STIX]{x1D6FD}),\end{eqnarray}$$

where

(2.10) $$\begin{eqnarray}\displaystyle \langle \cdots \rangle :=\int \text{d}\unicode[STIX]{x1D71E}f(\unicode[STIX]{x1D71E},t)\cdots & & \displaystyle\end{eqnarray}$$

is the macroscopic average, $\sum _{i,j}^{\prime }$ denotes the summation over $i$ and $j$ with $i\neq j$ and the macroscopic stress tensor denotes

(2.11) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}:=\int \text{d}\unicode[STIX]{x1D71E}f(\unicode[STIX]{x1D71E},t)\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\unicode[STIX]{x1D71E}). & & \displaystyle\end{eqnarray}$$

It might be noteworthy that the equation of continuity of the stress tensor, equation (2.9), is consistent with that for the Enskog theory of moderately dense inertial suspensions (Hayakawa, Takada & Garzó Reference Hayakawa, Takada and Garzó2017). In fact, the two terms on the right-hand side of (2.9), which are proportional to $\unicode[STIX]{x1D701}_{ij}^{(N)-1}$ and originate from particle contacts, correspond to the collision integral terms in the Enskog theory.

4.1 Grad’s 13-moment-like expansion

To obtain the stress tensor from (4.5), it is still necessary to have the distribution function $f(\unicode[STIX]{x1D71E},t)$ at hand to evaluate the statistical averages on the right-hand side. However, the exact expression of $f(\unicode[STIX]{x1D71E},t)$ for many-body problems is unknown and thus we should resort to approximations. Here we adopt Grad’s 13-moment-like expansion for $f(\unicode[STIX]{x1D71E},t)$ . This method is well established to approximate the distribution function in the kinetic theory of dilute or moderately dense gases (Grad Reference Grad1949; Herdegen & Hess Reference Herdegen and Hess1982; Jenkins & Richman Reference Jenkins and Richman1985b ,Reference Jenkins and Richman a ; Tsao & Koch Reference Tsao and Koch1995; Sangani et al. Reference Sangani, Mo, Tsao and Koch1996; Garzó Reference Garzó2002; Santos, Garzó & Dufty Reference Santos, Garzó and Dufty2004; Kremer Reference Kremer2010; Garzó Reference Garzó2013; Chamorro, Reyes & Garzó Reference Chamorro, Reyes and Garzó2015; Hayakawa & Takada Reference Hayakawa and Takada2016, Reference Hayakawa and Takada2017; Hayakawa et al. Reference Hayakawa, Takada and Garzó2017). It is an expansion in terms of the heat and stress currents in addition to the five conserved currents for the collisional invariants. For simple shear without spatial inhomogeneity, the current for the heat and the conserved quantities are negligible, and hence the velocity distribution function is dominated by the stress current,

(4.6) $$\begin{eqnarray}\displaystyle f(\boldsymbol{v})\approx f_{eq}(\boldsymbol{v})\left[1+\frac{V}{2T_{K}}\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(K)}\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(K)}(\boldsymbol{v})\right], & & \displaystyle\end{eqnarray}$$

where summation over repeated indices $\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}$ is taken, e.g.  $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}}:=\text{tr}(\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}})=\unicode[STIX]{x1D70E}_{xx}+\unicode[STIX]{x1D70E}_{yy}+\unicode[STIX]{x1D70E}_{zz}$ . Here, $\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(K)}(\boldsymbol{v}):=-mv_{\unicode[STIX]{x1D6FC}}v_{\unicode[STIX]{x1D6FD}}/V$ is the microscopic kinetic stress, where $m$ and $\boldsymbol{v}$ are the mass and velocity of the particle, $\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(K)}:=\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(K)}/P^{(K)}+\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}$ is the normalized deviatoric stress, where $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(K)}:=\int \text{d}\boldsymbol{v}f(\boldsymbol{v})\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(K)}(\boldsymbol{v})$ and $P^{(K)}:=-\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}}^{(K)}/3$ are the macroscopic kinetic stress and the pressure, $T_{K}:=-\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}}^{(K)}/(3n)$ is the kinetic temperature, where $n$ is the average number density, and $f_{eq}(\boldsymbol{v})$ is the equilibrium distribution function. Note that the kinetic pressure satisfies the relation $P^{(K)}=nT_{K}$ . This distribution function gives reasonably precise description of non-equilibrium gases, e.g. continuous as well as discontinuous shear thickening (Chamorro et al. Reference Chamorro, Reyes and Garzó2015; Hayakawa & Takada Reference Hayakawa and Takada2016, Reference Hayakawa and Takada2017; Hayakawa et al. Reference Hayakawa, Takada and Garzó2017). It is also notable that (4.6) satisfies the Green–Kubo formula within the Bhatnagar–Gross–Krook (BGK) approximation (Hayakawa & Takada Reference Hayakawa and Takada2016), while it further incorporates the normal stress differences, which is not the case for the Green–Kubo formula.

The distribution function (4.6) cannot be directly applied to dense non-Brownian suspensions, where the Cauchy stress dominates the kinetic stress. A possible extension of this expansion for non-Brownian suspensions would be

(4.7) $$\begin{eqnarray}\displaystyle f(\unicode[STIX]{x1D71E},t)\approx f_{eq}(\unicode[STIX]{x1D71E}_{v})f_{eq}(\unicode[STIX]{x1D71E}_{r})\left[1+\frac{V}{2T}\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(t)\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\unicode[STIX]{x1D71E}_{r})\right], & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D71E}:=\{\unicode[STIX]{x1D71E}_{r},\unicode[STIX]{x1D71E}_{v}\}$ with $\unicode[STIX]{x1D71E}_{r}:=\{\boldsymbol{r}_{i}\}\text{}_{i=1}^{N}$ and $\unicode[STIX]{x1D71E}_{v}:=\{\boldsymbol{v}_{i}\}\text{}_{i=1}^{N}$ , and the kinetic stress is replaced by the Cauchy stress. Here,

(4.8) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}:=\frac{\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}}{P}+\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}} & & \displaystyle\end{eqnarray}$$

is the normalized deviatoric stress and the appropriate definition of the temperature $T$ for non-Brownian suspensions will be discussed in § 4.4. Here it is postulated that the distribution function is factorized into peculiar velocity-dependent and position-dependent parts, where the peculiar velocity-dependent part can be approximated by Gaussian, $f_{eq}(\unicode[STIX]{x1D71E}_{v})$ , and the position-dependent part is approximated by an expansion around equilibrium, $f_{eq}(\unicode[STIX]{x1D71E}_{r})$ , with the stress current. This expansion around equilibrium is non-trivial for non-Brownian suspensions, where the equilibrium state is absent. Nonetheless, we will show in § 6 that the velocity distribution is nearly Gaussian and thus the factorization of (4.7) seems to be valid, and the expansion of the position distribution with the stress current is applicable for the evaluation of the stress.

From (4.7), the macroscopic average of an arbitrary observable $A(\unicode[STIX]{x1D71E}_{r})$ which depends on $\unicode[STIX]{x1D71E}_{r}$ is given by

(4.9) $$\begin{eqnarray}\langle A(\unicode[STIX]{x1D71E}_{r}(t))\rangle \approx \langle A(\unicode[STIX]{x1D71E}_{r})\rangle _{eq}+\frac{V}{2T}\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(t)\langle A(\unicode[STIX]{x1D71E}_{r})\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\unicode[STIX]{x1D71E}_{r})\rangle _{eq},\end{eqnarray}$$

where the first term on the right-hand side is the canonical term with

(4.10) $$\begin{eqnarray}\displaystyle \langle \cdots \rangle _{eq}:=\int \text{d}\unicode[STIX]{x1D71E}_{r}f_{eq}(\unicode[STIX]{x1D71E}_{r})\cdots & & \displaystyle\end{eqnarray}$$

and the second term is its non-canonical correction. Note that (4.9) is formally equivalent to the multiple relaxation-time approximation of the Green–Kubo formula (Suzuki & Hayakawa Reference Suzuki and Hayakawa2015), where the dimensionless tensor $\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}$ plays the role of the multiple relaxation times.

4.2 Approximate expressions

Let us evaluate the two averages on the right-hand side of (4.5) via the approximate formula (4.9). Detailed evaluation of the averages is shown in § B.1. These two averages are written as sums of terms of the form

(4.11) $$\begin{eqnarray}\displaystyle \mathop{\sum }_{i=1}^{N}\mathop{\sum }_{j\neq i}\mathop{\sum }_{k\neq i}\langle X(\boldsymbol{r}_{ij})Y(\boldsymbol{r}_{ik})\rangle , & & \displaystyle\end{eqnarray}$$

where $X(\boldsymbol{r}_{ij})$ or $Y(\boldsymbol{r}_{ik})$ abbreviates term which depends on $\boldsymbol{r}_{ij}$ or $\boldsymbol{r}_{ik}$ , respectively. For instance, for the first term on the right-hand side (4.5), $X(\boldsymbol{r}_{ij})$ and $Y(\boldsymbol{r}_{ik})$ are given by $X(\boldsymbol{r}_{ij})=\hat{r}_{ij,\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D6E5}_{ij}^{xy},Y(\boldsymbol{r}_{ik})=\hat{r}_{ik,\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D6E5}_{ik}^{xy}$ , or $X(\boldsymbol{r}_{ij})=\hat{r}_{ij,\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D6E5}_{ij}^{xy},Y(\boldsymbol{r}_{ik})=\hat{r}_{ik,\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D6E5}_{ik}^{xy}$ . Terms of the form (4.11) are evaluated by (4.9) as

(4.12) $$\begin{eqnarray}\displaystyle & & \displaystyle \mathop{\sum }_{i=1}^{N}\mathop{\sum }_{j\neq i}\mathop{\sum }_{k\neq i}\langle X(\boldsymbol{r}_{ij})Y(\boldsymbol{r}_{ik})\rangle \nonumber\\ \displaystyle & & \displaystyle \quad =\mathop{\sum }_{i=1}^{N}\mathop{\sum }_{j\neq i}\mathop{\sum }_{k\neq i}\left\{\langle X(\boldsymbol{r}_{ij})Y(\boldsymbol{r}_{ik})\rangle _{eq}+\frac{V}{2T}\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70E}}\langle X(\boldsymbol{r}_{ij})Y(\boldsymbol{r}_{ik})\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70E}}\rangle _{eq}\right\}.\end{eqnarray}$$

Note that $\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70E}}$ depends on the relative coordinate, e.g.  $\boldsymbol{r}_{lm}$ , but either $l$ or $m$ must be identical to e.g.  $i$ ; otherwise the correlation decouples and vanishes. Hence, the non-equilibrium term is given by

(4.13) $$\begin{eqnarray}\displaystyle \mathop{\sum }_{i=1}^{N}\mathop{\sum }_{j\neq i}\mathop{\sum }_{k\neq i}\mathop{\sum }_{l\neq i}\langle X(\boldsymbol{r}_{ij})Y(\boldsymbol{r}_{ik})\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70E}}(\boldsymbol{r}_{il})\rangle _{eq}, & & \displaystyle\end{eqnarray}$$

which can be decomposed into four-, three- and two-body correlations as

(4.14) $$\begin{eqnarray}\displaystyle & & \displaystyle \mathop{\sum }_{i=1}^{N}\mathop{\sum }_{j\neq i}\mathop{\sum }_{k\neq i}\mathop{\sum }_{l\neq i}\langle X(\boldsymbol{r}_{ij})Y(\boldsymbol{r}_{ik})\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70E}}(\boldsymbol{r}_{il})\rangle _{eq}=\sum _{i,j,k,l}\langle X(\boldsymbol{r}_{ij})Y(\boldsymbol{r}_{ik})\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70E}}(\boldsymbol{r}_{il})\rangle _{eq}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\sum _{i,j,k}\langle X(\boldsymbol{r}_{ij})Y(\boldsymbol{r}_{ik})\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70E}}(\boldsymbol{r}_{ik})\rangle _{eq}+\sum _{i,j}\langle X(\boldsymbol{r}_{ij})Y(\boldsymbol{r}_{ij})\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70E}}(\boldsymbol{r}_{ij})\rangle _{eq}.\end{eqnarray}$$

Here, the four-, three- and two-body terms are given by

(4.15) $$\begin{eqnarray}\displaystyle & \displaystyle \sum _{i,j,k,l}\langle X(\boldsymbol{r}_{ij})Y(\boldsymbol{r}_{ik})\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70E}}(\boldsymbol{r}_{il})\rangle _{eq}=Nn^{3}\int \text{d}^{3}\boldsymbol{r}\int \text{d}^{3}\boldsymbol{r}^{\prime }\int \text{d}^{3}\boldsymbol{r}^{\prime \prime }g^{(4)}(\boldsymbol{r},\boldsymbol{r}^{\prime },\boldsymbol{r}^{\prime \prime })X(\boldsymbol{r})Y(\boldsymbol{r}^{\prime })\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70E}}(\boldsymbol{r}^{\prime \prime }), & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(4.16) $$\begin{eqnarray}\displaystyle & \displaystyle \sum _{i,j,k}\langle X(\boldsymbol{r}_{ij})Y(\boldsymbol{r}_{ik})\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70E}}(\boldsymbol{r}_{ik})\rangle _{eq}=Nn^{3}\int \text{d}^{3}\boldsymbol{r}\int \text{d}^{3}\boldsymbol{r}^{\prime }g^{(3)}(\boldsymbol{r},\boldsymbol{r}^{\prime })X(\boldsymbol{r})Y(\boldsymbol{r}^{\prime })\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70E}}(\boldsymbol{r}^{\prime }), & \displaystyle\end{eqnarray}$$

(4.17) $$\begin{eqnarray}\displaystyle & \displaystyle \sum _{i,j(i\neq j)}\langle X(\boldsymbol{r}_{ij})Y(\boldsymbol{r}_{ij})\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70E}}(\boldsymbol{r}_{ij})\rangle _{eq}=Nn^{3}\int \text{d}^{3}\boldsymbol{r}g(r)X(\boldsymbol{r})Y(\boldsymbol{r})\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70E}}(\boldsymbol{r}), & \displaystyle\end{eqnarray}$$

respectively, where $g^{(4)}(\boldsymbol{r},\boldsymbol{r}^{\prime },\boldsymbol{r}^{\prime \prime })$ and $g^{(3)}(\boldsymbol{r},\boldsymbol{r}^{\prime })$ are the quadruplet-correlation function (Hansen & McDonald Reference Hansen and McDonald2006)

(4.18) $$\begin{eqnarray}g^{(4)}(\boldsymbol{r},\boldsymbol{r}^{\prime },\boldsymbol{r}^{\prime \prime }):=\frac{1}{Nn^{3}}\sum _{i,j,k,l}\langle \unicode[STIX]{x1D6FF}(\boldsymbol{r}-\boldsymbol{r}_{ij})\unicode[STIX]{x1D6FF}(\boldsymbol{r}^{\prime }-\boldsymbol{r}_{ik})\unicode[STIX]{x1D6FF}(\boldsymbol{r}^{\prime \prime }-\boldsymbol{r}_{il})\rangle _{eq}\end{eqnarray}$$

and the triplet-correlation function (Hansen & McDonald Reference Hansen and McDonald2006)

(4.19) $$\begin{eqnarray}g^{(3)}(\boldsymbol{r},\boldsymbol{r}^{\prime }):=\frac{1}{Nn^{2}}\sum _{i,j,k}\langle \unicode[STIX]{x1D6FF}(\boldsymbol{r}-\boldsymbol{r}_{ij})\unicode[STIX]{x1D6FF}(\boldsymbol{r}^{\prime }-\boldsymbol{r}_{ik})\rangle _{eq},\end{eqnarray}$$

and the summations $\sum _{i,j,k}^{\prime \prime }$ , and $\sum _{i,j,k,l}^{\prime \prime \prime }$ are performed over different particles. For instance, $\sum _{i,j,k}^{\prime \prime }$ is performed for $i,j,k$ with $i\neq j$ , $j\neq k$ and $k\neq i$ . In the spatial integrations over $g(r)$ , $g^{(3)}(\boldsymbol{r},\boldsymbol{r}^{\prime })$ , or $g^{(4)}(\boldsymbol{r},\boldsymbol{r}^{\prime },\boldsymbol{r}^{\prime \prime })$ , it is crucial that these correlation functions are accompanied by the delta functions $\unicode[STIX]{x1D6FF}(r-d)$ , $\unicode[STIX]{x1D6FF}(r-d)\unicode[STIX]{x1D6FF}(r^{\prime }-d)$ , or $\unicode[STIX]{x1D6FF}(r-d)\unicode[STIX]{x1D6FF}(r^{\prime }-d)\unicode[STIX]{x1D6FF}(r^{\prime \prime }-d)$ , respectively. This feature can be explicitly traced back in (4.5), and is a consequence of the hard-core collision of the particles, equation (3.15). This implies that only the contact values of the correlation functions contribute. By virtue of this feature, we can conveniently approximate $g^{(3)}(\boldsymbol{r},\boldsymbol{r}^{\prime })$ and $g^{(4)}(\boldsymbol{r},\boldsymbol{r}^{\prime },\boldsymbol{r}^{\prime \prime })$ .

Let us consider $g^{(3)}(\boldsymbol{r},\boldsymbol{r}^{\prime })$ for illustration. First of all, we adopt the factorization approximation (Kirkwood Reference Kirkwood1935), $g^{(3)}(\boldsymbol{r},\boldsymbol{r}^{\prime })\approx g(r)g(r^{\prime })g(|\boldsymbol{r}-\boldsymbol{r}^{\prime }|)$ . Although this approximation is not accurate in general, it has been argued that it is valid at contacts, where $r,r^{\prime },|\boldsymbol{r}-\boldsymbol{r}^{\prime }|\approx d$ (Alder Reference Alder1964; Grouba, Zorin & Sevastianov Reference Grouba, Zorin and Sevastianov2004). Furthermore, we have shown in Suzuki & Hayakawa (Reference Suzuki and Hayakawa2015) that only the radial contacts contribute to the divergence in the vicinity of the jamming point, i.e.

(4.20) $$\begin{eqnarray}\displaystyle g^{(3)}(\boldsymbol{r},\boldsymbol{r}^{\prime })\approx g(r)g(r^{\prime }) & & \displaystyle\end{eqnarray}$$

for $r,r^{\prime }\approx d$ . Similarly, $g^{(4)}(\boldsymbol{r},\boldsymbol{r}^{\prime },\boldsymbol{r}^{\prime \prime })$ can be approximated as

(4.21) $$\begin{eqnarray}\displaystyle g^{(4)}(\boldsymbol{r},\boldsymbol{r}^{\prime },\boldsymbol{r}^{\prime \prime })\approx g(r)g(r^{\prime })g(r^{\prime \prime })g(|\boldsymbol{r}-\boldsymbol{r}^{\prime }|)g(|\boldsymbol{r}^{\prime }-\boldsymbol{r}^{\prime \prime }|)g(|\boldsymbol{r}^{\prime \prime }-\boldsymbol{r}|)\approx g(r)g(r^{\prime })g(r^{\prime \prime })\qquad & & \displaystyle\end{eqnarray}$$

for $r,r^{\prime },r^{\prime \prime }\approx d$ , as far as divergence is concerned. We will examine the validity of the factorization approximation for the evaluation of the stress in § 6.

These approximations, together with $g(r)\unicode[STIX]{x1D6FF}(r-d)=g_{0}(\unicode[STIX]{x1D711})\unicode[STIX]{x1D6FF}(r-d)$ , enable us to express the two averages on the right-hand side of (4.5) in terms of polynomials of the radial distribution function at contact, $g_{0}(\unicode[STIX]{x1D711})$ . From tedious but straightforward calculation as shown in §§ B.1 and B.2, we reach the approximate equation of the stress

(4.22) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\text{d}}{\text{d}t}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}+\frac{1}{2}\dot{\unicode[STIX]{x1D6FE}}(\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FC}x}\unicode[STIX]{x1D70E}_{y\unicode[STIX]{x1D6FD}}+\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FD}x}\unicode[STIX]{x1D70E}_{y\unicode[STIX]{x1D6FC}})\approx \frac{\unicode[STIX]{x1D701}\dot{\unicode[STIX]{x1D6FE}}^{2}}{4d}\mathop{\sum }_{\ell =1}^{2}\{\!-\unicode[STIX]{x1D711}^{\ast 3}g_{0}(\unicode[STIX]{x1D711})^{2}{\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(\ell :c2)}-\unicode[STIX]{x1D711}^{\ast 2}g_{0}(\unicode[STIX]{x1D711}){\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(\ell :c1)}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\unicode[STIX]{x1D6EC}\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}}[\unicode[STIX]{x1D711}^{\ast 4}g_{0}(\unicode[STIX]{x1D711})^{3}{\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}}^{(\ell :nc3)}+\unicode[STIX]{x1D711}^{\ast 3}g_{0}(\unicode[STIX]{x1D711})^{2}{\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}}^{(\ell :nc2)}+\unicode[STIX]{x1D711}^{\ast 2}g_{0}(\unicode[STIX]{x1D711}){\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}}^{(\ell :nc1)}]\!\},\end{eqnarray}$$

where the two terms with coefficients ${\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(\ell :c2)}$ and ${\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(\ell :c1)}$ are the canonical contributions, and the three terms with coefficients ${\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}}^{(\ell :nc3)}$ , ${\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}}^{(\ell :nc2)}$ and ${\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}}^{(\ell :nc1)}$ are the non-equilibrium corrections. The coefficients ${\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(\ell :c1)}$ , ${\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(\ell :c2)}$ and ${\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}}^{(\ell :nc3)}$ , ${\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}}^{(\ell :nc2)}$ , ${\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}}^{(\ell :nc1)}$ are numbers which arise from angular integrals. Here, we have introduced a dimensionless scalar

(4.23) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6EC}:=\frac{\unicode[STIX]{x1D701}\dot{\unicode[STIX]{x1D6FE}}d^{2}}{4T}, & & \displaystyle\end{eqnarray}$$

and $\unicode[STIX]{x1D711}^{\ast }=6\unicode[STIX]{x1D711}/\unicode[STIX]{x03C0}=nd^{3}$ denotes the dimensionless number density.

4.3 Implications of the symmetry

A specific feature of non-Brownian suspensions under simple shear is the symmetry under the parity ‘ $\hat{x}\rightarrow -\hat{x}$ and ${\hat{y}}\rightarrow -{\hat{y}}$ ’, which follows from ${\mathcal{P}}(\hat{x},{\hat{y}})$ introduced in (3.14). This implies that only the parity-even terms in (4.22) survive. Furthermore, even though parity even, terms odd with respect to $\hat{z}$ vanish. These features can be summarized as follows,

(4.24) $$\begin{eqnarray}\displaystyle \int \text{d}{\mathcal{S}}\,\unicode[STIX]{x1D6E9}(-\hat{x}{\hat{y}})\hat{x}^{i}{\hat{y}}^{j}\hat{z}^{k}\neq 0\quad \text{if and only if }i+j=\text{even and }k=\text{even}, & & \displaystyle\end{eqnarray}$$

where $\int \text{d}{\mathcal{S}}\cdots \,$ expresses an angular integral with respect to $\hat{\boldsymbol{r}}$ . As a consequence of this property we have

(4.25) $$\begin{eqnarray}\displaystyle \int \text{d}{\mathcal{S}}\unicode[STIX]{x1D6E9}(-\hat{x}{\hat{y}})\hat{x}^{i}{\hat{y}}^{j}\hat{z}^{k}=0\quad \text{if }i+j+k=\text{odd}. & & \displaystyle\end{eqnarray}$$

The terms proportional to $g_{0}(\unicode[STIX]{x1D711})^{3}$ in (4.22), which are cubic with respect to $\hat{\boldsymbol{r}}$ , vanish because of (4.25),

(4.26) $$\begin{eqnarray}\displaystyle {\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}}^{(1:nc3)}={\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}}^{(2:nc3)}=0. & & \displaystyle\end{eqnarray}$$

The same observation holds for the canonical terms proportional to $g_{0}(\unicode[STIX]{x1D711})^{2}$ , which are also cubic in $\hat{\boldsymbol{r}}$ ,

(4.27) $$\begin{eqnarray}\displaystyle {\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(1:c2)}={\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(2:c2)}=0. & & \displaystyle\end{eqnarray}$$

From (4.26) and (4.27), (4.22) reduces to

(4.28) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\text{d}}{\text{d}t}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}+\frac{1}{2}\dot{\unicode[STIX]{x1D6FE}}(\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FC}x}\unicode[STIX]{x1D70E}_{y\unicode[STIX]{x1D6FD}}+\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FD}x}\unicode[STIX]{x1D70E}_{y\unicode[STIX]{x1D6FC}})\approx \frac{\unicode[STIX]{x1D701}\dot{\unicode[STIX]{x1D6FE}}^{2}}{4d}\mathop{\sum }_{\ell =1}^{2}\{\!-\unicode[STIX]{x1D711}^{\ast 2}g_{0}(\unicode[STIX]{x1D711}){\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(\ell :c1)}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\unicode[STIX]{x1D6EC}\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}}[\unicode[STIX]{x1D711}^{\ast 3}g_{0}(\unicode[STIX]{x1D711})^{2}{\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}}^{(\ell :nc2)}+\unicode[STIX]{x1D711}^{\ast 2}g_{0}(\unicode[STIX]{x1D711}){\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}}^{(\ell :nc1)}]\!\},\end{eqnarray}$$

or, explicitly in components,

(4.29) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D70E}_{xy}\\ -3P\\ \unicode[STIX]{x1D70E}_{xx}\\ \unicode[STIX]{x1D70E}_{yy}\end{array}\right)+\left(\begin{array}{@{}c@{}}\frac{1}{2}\unicode[STIX]{x1D70E}_{yy}\\ \unicode[STIX]{x1D70E}_{xy}\\ \unicode[STIX]{x1D70E}_{xy}\\ 0\end{array}\right)\approx \frac{\unicode[STIX]{x1D701}\dot{\unicode[STIX]{x1D6FE}}}{4d}\unicode[STIX]{x1D711}^{\ast 2}g_{0}(\unicode[STIX]{x1D711})\left\{-\boldsymbol{{\mathcal{B}}}+\unicode[STIX]{x1D6EC}(\unicode[STIX]{x1D711}^{\ast }g_{0}(\unicode[STIX]{x1D711})\boldsymbol{{\mathcal{A}}}^{(2)}+\boldsymbol{{\mathcal{A}}}^{(1)})\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D6F1}_{xy}\\ \unicode[STIX]{x1D6F1}_{xx}\\ \unicode[STIX]{x1D6F1}_{yy}\end{array}\right)\right\}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where the matrices $\boldsymbol{{\mathcal{A}}}^{(m)}$ ( $m=1,2$ ) and the vector $\boldsymbol{{\mathcal{B}}}$ are given by (C 64)–(C 66). Note that $\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}$ is given by (4.8) and hence (4.29) is a closed set of equations for $P$ , $\unicode[STIX]{x1D70E}_{xy}$ , $\unicode[STIX]{x1D70E}_{xx}$ and $\unicode[STIX]{x1D70E}_{yy}$ .

Note that (4.28) includes terms proportional to $g_{0}(\unicode[STIX]{x1D711})^{2}$ or $g_{0}(\unicode[STIX]{x1D711})$ . Hence, we decompose the stress tensor as $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}=\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(2)}+\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(1)}$ , where $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(m)}$ is $O(g_{0}(\unicode[STIX]{x1D711})^{m})$ ( $m=1,2$ ), because (3.13) suggests that it is separable near the jamming point. Then, equation (4.28) is cast into two equations, each for $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(2)}$ and $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(1)}$ , respectively,

(4.30) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}t}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(2)}+\frac{1}{2}\dot{\unicode[STIX]{x1D6FE}}(\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FC}x}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FD}y}^{(2)}+\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FD}x}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}y}^{(2)})=\frac{\unicode[STIX]{x1D701}\dot{\unicode[STIX]{x1D6FE}}^{2}}{4d}\unicode[STIX]{x1D6EC}\unicode[STIX]{x1D711}^{\ast 3}g_{0}(\unicode[STIX]{x1D711})^{2}\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}}^{(2)}\mathop{\sum }_{\ell =1}^{2}{\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}}^{(\ell :nc2)}. & \displaystyle\end{eqnarray}$$

(4.31) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}t}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(1)}+\frac{1}{2}\dot{\unicode[STIX]{x1D6FE}}(\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FC}x}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FD}y}^{(1)}+\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FD}x}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}y}^{(1)})=\frac{\unicode[STIX]{x1D701}\dot{\unicode[STIX]{x1D6FE}}^{2}}{4d}\unicode[STIX]{x1D711}^{\ast 2}g_{0}(\unicode[STIX]{x1D711})\mathop{\sum }_{\ell =1}^{2}\{-{\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(\ell :c1)}+\unicode[STIX]{x1D6EC}\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}}^{(1)}{\mathcal{S}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D706}}^{(\ell :nc1)}\}.\qquad & \displaystyle\end{eqnarray}$$

Here, we have introduced

(4.32) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(m)}:=\frac{\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(m)}}{P^{(m)}}+\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\quad (m=1,2), & & \displaystyle\end{eqnarray}$$

where $P^{(m)}$ is the component of the pressure of $O(g_{0}(\unicode[STIX]{x1D711})^{m})$ ( $m=1,2$ ). The valid range of (4.30) and (4.31) from (4.28) is discussed in § D.2.

Another implication of the symmetry is that the uniform profile of the velocity distribution is maintained. This issue is discussed in appendix G.

4.4 Interpretation of the temperature

In contrast to the case of the kinetic theory of dilute and moderately dense gases, or even dense inertial suspensions near the jamming point, the determination method of the temperature $T$ is not clear for non-Brownian suspensions. As discussed in § 4.1, in the kinetic theory, the temperature is determined by the equation of state $P^{(K)}=nT_{K}$ , where $P^{(K)}$ and $T_{K}$ are the kinetic pressure and the kinetic temperature, respectively. We attempt to introduce the temperature by the Cauchy stress, which dominates the kinetic stress in dense non-Brownian suspensions. Because $T$ appears in the position-dependent part in (4.7), $T$ should be defined by the position-dependent Cauchy stress. Note that $T$ determines the magnitude of the non-equilibrium correction of the distribution function, equation (4.7). Accordingly, $\unicode[STIX]{x1D6EC}$ introduced in (4.23) determines the non-equilibrium correction of the stress, equations (4.28) or (4.29). In this paper, let us introduce $T$ by the equation of state

(4.33) $$\begin{eqnarray}\displaystyle P^{(eq)}=nT[1+2\unicode[STIX]{x1D711}g_{0}(\unicode[STIX]{x1D711})], & & \displaystyle\end{eqnarray}$$

where $P^{(eq)}$ is the equilibrium part of the pressure determined from (4.28), (4.29) or (4.31) with $\unicode[STIX]{x1D6EC}=0$ . In other words, $P^{(eq)}$ can be estimated only by $f_{eq}(\unicode[STIX]{x1D71E}_{r})$ in (4.7). The solution of (4.31) is given in (D 2), (D 3) in § D.1, from which we obtain $P^{(eq)}$ as ( $P^{(1)}$ with $\unicode[STIX]{x1D6EC}=0$ )

(4.34) $$\begin{eqnarray}\displaystyle P^{(eq)}=0.0060\times \frac{\unicode[STIX]{x1D701}\dot{\unicode[STIX]{x1D6FE}}}{4d}\unicode[STIX]{x1D711}^{\ast 2}g_{0}(\unicode[STIX]{x1D711})=0.022\times \frac{\unicode[STIX]{x1D701}\dot{\unicode[STIX]{x1D6FE}}}{4d}\unicode[STIX]{x1D711}^{2}g_{0}(\unicode[STIX]{x1D711}), & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D711}^{\ast }=6\unicode[STIX]{x1D711}/\unicode[STIX]{x03C0}$ . Thus, from (4.33) and (4.34), we obtain

(4.35) $$\begin{eqnarray}\displaystyle T=0.022\times \frac{\unicode[STIX]{x1D701}\dot{\unicode[STIX]{x1D6FE}}}{4nd}\frac{\unicode[STIX]{x1D711}^{2}g_{0}(\unicode[STIX]{x1D711})}{1+2\unicode[STIX]{x1D711}g_{0}(\unicode[STIX]{x1D711})}\approx 0.022\times \frac{\unicode[STIX]{x1D701}\dot{\unicode[STIX]{x1D6FE}}d^{2}}{8nd^{3}}\unicode[STIX]{x1D711}=0.0014\unicode[STIX]{x1D701}\dot{\unicode[STIX]{x1D6FE}}d^{2}, & & \displaystyle\end{eqnarray}$$

where we have approximated $1+2\unicode[STIX]{x1D711}g_{0}(\unicode[STIX]{x1D711})\approx 2\unicode[STIX]{x1D711}g_{0}(\unicode[STIX]{x1D711})$ in the second equality. This determines $\unicode[STIX]{x1D6EC}$ as

(4.36) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6EC}=\frac{\unicode[STIX]{x1D701}\dot{\unicode[STIX]{x1D6FE}}d^{2}}{4T}\approx 174, & & \displaystyle\end{eqnarray}$$

which is independent of dimensional physical variables such as $\dot{\unicode[STIX]{x1D6FE}}$ , $d$ or $\unicode[STIX]{x1D702}_{0}$ , and is merely a number. This value is larger than the value $\unicode[STIX]{x1D6EC}=0.04$ determined by fitting the absolute values of the shear and pressure viscosities to the result of molecular dynamics (MD) simulation (cf. §§ 4.5 and 5). This might be due to the fact that $\unicode[STIX]{x1D701}$ in (4.36) is not equivalent to $\unicode[STIX]{x1D701}_{0}$ which sets the magnitude of the viscosities in MD simulation. In fact, $\unicode[STIX]{x1D701}$ is given by (4.3) as $\unicode[STIX]{x1D701}=3\unicode[STIX]{x1D701}_{0}/(128\unicode[STIX]{x1D716}^{2})$ , where $\unicode[STIX]{x1D716}\ll 1$ is the magnitude of the separation of contacting particles. The ratio $\unicode[STIX]{x1D701}/\unicode[STIX]{x1D701}_{0}=174/0.04$ corresponds to $\unicode[STIX]{x1D716}\approx 0.002$ , but since we cannot determine the magnitude of $\unicode[STIX]{x1D716}$ in our framework, we will leave $\unicode[STIX]{x1D6EC}$ as a fitting parameter.

The important implication of (4.35) is

(4.37) $$\begin{eqnarray}\displaystyle T\propto \unicode[STIX]{x1D701}\dot{\unicode[STIX]{x1D6FE}}d^{2}\propto 3\unicode[STIX]{x03C0}d^{3}\unicode[STIX]{x1D702}_{0}\dot{\unicode[STIX]{x1D6FE}}, & & \displaystyle\end{eqnarray}$$

which is consistent with the ‘effective temperature’ of suspensions (Ono et al. Reference Ono, O’Hern, Durian, Langer, Liu and Nagel2002; Eisenmann et al. Reference Eisenmann, Kim, Mattsson and Weitz2010) defined via the Stokes–Einstein relation, $D=T_{eff}/(3\unicode[STIX]{x03C0}\,\text{d}\unicode[STIX]{x1D702}_{0})$ , where $D$ is the diffusion coefficient. It has been shown by experiment (Leighton & Acrivos Reference Leighton and Acrivos1987a ,Reference Leighton and Acrivos b ; Breedveld et al. Reference Breedveld, van den Ende, Tripathi and Acrivos1998, Reference Breedveld, van den Ende, Bosscher, Jongschaap and Mellema2002; Eisenmann et al. Reference Eisenmann, Kim, Mattsson and Weitz2010) and simulation (Foss & Brady Reference Foss and Brady1999; Heussinger et al. Reference Heussinger, Berthier and Barrat2010; Olsson Reference Olsson2010) that $D\sim \text{d}^{2}\dot{\unicode[STIX]{x1D6FE}}$ holds below the jamming point and $D$ exhibits only a weak dependence on the density, which suggests $T_{eff}=3\unicode[STIX]{x03C0}d\unicode[STIX]{x1D702}_{0}D\sim 3\unicode[STIX]{x03C0}d^{3}\unicode[STIX]{x1D702}_{0}\dot{\unicode[STIX]{x1D6FE}}$ .

4.5 Viscosities and $\unicode[STIX]{x1D707}$ – $J$ rheology in the steady state

Let us analyse the rheology in the steady state, i.e. consider (4.28) with $\text{d}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}/\text{d}t=0$ . Then, equations (4.30) and (4.31) can be solved analytically. The analytic solutions for $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(2)}$ and $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(1)}$ are given by (D 1)–(D 4) in § D.1, together with (3.13):

(4.38) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(2)}\sim \frac{\unicode[STIX]{x1D701}\dot{\unicode[STIX]{x1D6FE}}}{d}g_{0}(\unicode[STIX]{x1D711})^{2}\sim \frac{\unicode[STIX]{x1D701}\dot{\unicode[STIX]{x1D6FE}}}{d}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}^{-2}, & \displaystyle\end{eqnarray}$$

(4.39) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}^{(1)}\sim \frac{\unicode[STIX]{x1D701}\dot{\unicode[STIX]{x1D6FE}}}{d}g_{0}(\unicode[STIX]{x1D711})\sim \frac{\unicode[STIX]{x1D701}\dot{\unicode[STIX]{x1D6FE}}}{d}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}^{-1}. & \displaystyle\end{eqnarray}$$

In particular, the normalized shear viscosity $\unicode[STIX]{x1D702}_{s}^{\ast }:=\unicode[STIX]{x1D702}_{s}/\unicode[STIX]{x1D702}_{0}$ and pressure viscosity $\unicode[STIX]{x1D702}_{n}^{\ast }:=\unicode[STIX]{x1D702}_{n}/\unicode[STIX]{x1D702}_{0}$ are given by

(4.40) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}_{s}^{\ast }\approx 0.0073\unicode[STIX]{x1D711}^{\ast 3}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}^{-2}+0.080\unicode[STIX]{x1D711}^{\ast 2}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}^{-1}, & \displaystyle\end{eqnarray}$$

(4.41) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}_{n}^{\ast }\approx 0.045\unicode[STIX]{x1D711}^{\ast 3}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}^{-2}+0.0005\unicode[STIX]{x1D711}^{\ast 2}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}^{-1}\approx 0.045\unicode[STIX]{x1D711}^{\ast 3}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}^{-2}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6EC}=0.04$ is adopted, which is determined by fitting the absolute values of $\unicode[STIX]{x1D702}_{s}^{\ast }$ and $\unicode[STIX]{x1D702}_{n}^{\ast }$ to the results of the MD simulation shown in the next section. The empirical formula $g_{0}(\unicode[STIX]{x1D711})=g_{CS}(\unicode[STIX]{x1D711}_{f})(\unicode[STIX]{x1D711}_{J}-\unicode[STIX]{x1D711}_{f})/(\unicode[STIX]{x1D711}_{J}-\unicode[STIX]{x1D711})$ with $\unicode[STIX]{x1D711}_{f}=0.49$ and $g_{CS}(\unicode[STIX]{x1D711}):=(1-\unicode[STIX]{x1D711}/2)/(1-\unicode[STIX]{x1D711})^{3}$ valid for $\unicode[STIX]{x1D711}_{f}<\unicode[STIX]{x1D711}<\unicode[STIX]{x1D711}_{J}$ (Torquato Reference Torquato1995) is adopted, from which we obtain $g_{0}(\unicode[STIX]{x1D711})\approx 0.848\,\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}^{-1}$ for $\unicode[STIX]{x1D711}_{J}=0.639$ . Although it is widely recognized that $\unicode[STIX]{x1D711}_{J}\approx 0.64$ for monodisperse frictionless hard spheres without any solvent, it has been reported that the value of $\unicode[STIX]{x1D711}_{J}$ is not uniquely identified and depends on the protocols used to generate the jammed configurations (Ciamarra, Coniglio & de Candia Reference Ciamarra, Coniglio and de Candia2010). In this work, we adopt $\unicode[STIX]{x1D711}_{J}=0.639$ , which is obtained by the conjugate-gradient protocol (O’Hern et al. Reference O’Hern, Silbert, Liu and Nagel2003). The numerical coefficients in (4.40) and (4.41) are determined analytically in rational forms, but we only show their approximate values in decimals for brevity.

It is evident from (4.40) and (4.41) that $\unicode[STIX]{x1D702}_{s}^{\ast }$ and $\unicode[STIX]{x1D702}_{n}^{\ast }$ exhibit

(4.42) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}_{s}^{\ast }\sim \unicode[STIX]{x1D702}_{n}^{\ast }\sim \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}^{-2} & & \displaystyle\end{eqnarray}$$

near the jamming point. On the other hand, $O(\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}^{-1})$ term in $\unicode[STIX]{x1D702}_{n}^{\ast }$ is small for finite $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}$ , while the corresponding term in $\unicode[STIX]{x1D702}_{s}^{\ast }$ is significant. Noting $J=1/\unicode[STIX]{x1D702}_{n}^{\ast }$ , $J$ and $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}$ are uniquely invertible via (4.41) as

(4.43) $$\begin{eqnarray}\displaystyle J^{1/2}\approx 4.74\unicode[STIX]{x1D711}^{\ast -3/2}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}, & & \displaystyle\end{eqnarray}$$

or, equivalently,

(4.44) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D711}=\unicode[STIX]{x1D711}_{J}-0.211\unicode[STIX]{x1D711}^{\ast 3/2}J^{1/2}. & & \displaystyle\end{eqnarray}$$

From (4.40) and (4.41), we obtain the stress ratio $\unicode[STIX]{x1D707}$ as

(4.45) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}\approx 0.163+1.78\unicode[STIX]{x1D711}^{\ast -1}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D711}=0.163+0.377\unicode[STIX]{x1D711}^{\ast 1/2}J^{1/2}. & & \displaystyle\end{eqnarray}$$

Note that $\unicode[STIX]{x1D707}$ in the jamming limit, $\unicode[STIX]{x1D707}(J\rightarrow 0)=0.163$ , is independent of $\unicode[STIX]{x1D6EC}$ .

Figure 3. Comparison of theory with MD simulation. (a) The normalized shear viscosity $\unicode[STIX]{x1D702}_{s}^{\ast }=\unicode[STIX]{x1D702}_{s}/\unicode[STIX]{x1D702}_{0}$ and the pressure viscosity $\unicode[STIX]{x1D702}_{n}^{\ast }=\unicode[STIX]{x1D702}_{n}/\unicode[STIX]{x1D702}_{0}$ , and (b) the stress ratio $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D70E}_{xy}/P$ . The results of the theory and MD for the shear (pressure) viscosity are shown in dashed (solid) lines and open circles (squares) in (a), and those for the stress ratio are shown in solid line and open diamonds in (b), respectively. The theoretical result adopts the fitting parameter $\unicode[STIX]{x1D6EC}=0.04$ . The result of MD is for $\unicode[STIX]{x1D711}=$ 0.632, 0.631, 0.63, 0.626, 0.62, 0.61, 0.6, 0.58, 0.54 and 0.5.

6.1 Peculiar velocity distribution function

The distribution of the peculiar velocity $\boldsymbol{v}_{i}=\dot{\boldsymbol{r}}_{i}-\dot{\unicode[STIX]{x1D6FE}}y_{i}\boldsymbol{e}_{x}$ is measured by MD simulation. As shown in figure 6, the peculiar velocity distribution function is nearly Gaussian in spite of the absence of any inertial effect in the dynamics. This result implies that the factorization of the distribution function assumed in (4.7) is valid. The result of our simulation supports another theoretical assumption that the base state is nearly equilibrium.

Figure 6. Distribution of the peculiar velocity measured by MD simulation (circles). The conditions are $\unicode[STIX]{x1D711}=0.63$ , $N=1000$ , $\dot{\unicode[STIX]{x1D6FE}}^{\ast }=1\times 10^{-5}$ . The number of data averaged for each point is $10^{5}$ . The figure shows the $x$ - and $y$ -components. Also shown is the Gaussian distribution fitted to the measured data (solid line).

6.2 Grad’s 13-moment-like expansion and factorization of multi-body correlations

We validate the expansion of the distribution of the position with the stress current, equation (4.7). For this purpose, we consider the radial distribution function at the steady state, $g_{\dot{\unicode[STIX]{x1D6FE}}}(\boldsymbol{r})$ , which is anisotropic under shear. It is given by

(6.1) $$\begin{eqnarray}\displaystyle g_{\dot{\unicode[STIX]{x1D6FE}}}(\boldsymbol{r}):=\frac{1}{Nn}\left\langle \mathop{\sum }_{i}\mathop{\sum }_{j\neq i}\unicode[STIX]{x1D6FF}(\boldsymbol{r}-\boldsymbol{r}_{ij})\right\rangle =\frac{1}{Nn}\int \text{d}\unicode[STIX]{x1D71E}f(\unicode[STIX]{x1D71E})\mathop{\sum }_{i}\mathop{\sum }_{j\neq i}\unicode[STIX]{x1D6FF}(\boldsymbol{r}-\boldsymbol{r}_{ij}). & & \displaystyle\end{eqnarray}$$

From (4.7), after trivially integrating out $f_{eq}(\unicode[STIX]{x1D71E}_{v})$ , we obtain

(6.2) $$\begin{eqnarray}\displaystyle g_{\dot{\unicode[STIX]{x1D6FE}}}(\boldsymbol{r})=g(r)+\frac{1}{Nn}\frac{V}{2T}\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\left\langle \mathop{\sum }_{i}\mathop{\sum }_{j\neq i}\unicode[STIX]{x1D6FF}(\boldsymbol{r}-\boldsymbol{r}_{ij})\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\right\rangle _{\!\!eq}, & & \displaystyle\end{eqnarray}$$

where $g(r)$ is the radial distribution function at equilibrium. By substituting the expression for the stress $\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}$ , equation (2.8) and that for the interparticle force $\boldsymbol{F}_{ij}^{(p)}$ therein, equation (3.8), we obtain

(6.3) $$\begin{eqnarray}\displaystyle & & \displaystyle g_{\dot{\unicode[STIX]{x1D6FE}}}(\boldsymbol{r})=g(r)-\frac{1}{Nn}\frac{\unicode[STIX]{x1D701}_{e}\dot{\unicode[STIX]{x1D6FE}}d^{2}}{8T}\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\left\langle \mathop{\sum }_{i}\mathop{\sum }_{j\neq i}\mathop{\sum }_{k\neq i}\unicode[STIX]{x1D6FF}(\boldsymbol{r}-\boldsymbol{r}_{ij})r_{ik,\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D6FF}(r_{ik}-d){\mathcal{P}}(\hat{x}_{ik},{\hat{y}}_{ik})\hat{r}_{ik,\unicode[STIX]{x1D6FC}}\right\rangle _{\!\!eq}\nonumber\\ \displaystyle & & \displaystyle \quad =g(r)-\frac{1}{Nn}\frac{\unicode[STIX]{x1D701}_{e}\dot{\unicode[STIX]{x1D6FE}}d^{3}}{8T}\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\left\langle \mathop{\sum }_{i}\mathop{\sum }_{j\neq i}\mathop{\sum }_{k\neq i}\unicode[STIX]{x1D6FF}(\boldsymbol{r}-\boldsymbol{r}_{ij})\unicode[STIX]{x1D6FF}(r_{ik}-d){\mathcal{P}}(\hat{x}_{ik},{\hat{y}}_{ik})\hat{r}_{ik,\unicode[STIX]{x1D6FC}}\hat{r}_{ik,\unicode[STIX]{x1D6FD}}\right\rangle _{\!\!eq}.\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

In the second equality, we have utilized $r_{ik,\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D6FF}(r_{ik}-d)=d\hat{r}_{ik,\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D6FF}(r_{ik}-d)$ . Equation (6.3) is expressed in terms of the triplet-correlation function $g^{(3)}(\boldsymbol{r},\boldsymbol{r}^{\prime })$ , equation (4.19), as

(6.4) $$\begin{eqnarray}\displaystyle g_{\dot{\unicode[STIX]{x1D6FE}}}(\boldsymbol{r})=g(r)-n\frac{\unicode[STIX]{x1D701}_{e}\dot{\unicode[STIX]{x1D6FE}}d^{3}}{8T}\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\int \text{d}^{3}\boldsymbol{r}^{\prime }g^{(3)}(\boldsymbol{r},\boldsymbol{r}^{\prime })\unicode[STIX]{x1D6FF}(r^{\prime }-d){\mathcal{P}}(\hat{x}^{\prime },{\hat{y}}^{\prime })\hat{r}_{\unicode[STIX]{x1D6FC}}^{\prime }\hat{r}_{\unicode[STIX]{x1D6FD}}^{\prime }. & & \displaystyle\end{eqnarray}$$

We apply the factorization approximation to $g^{(3)}(\boldsymbol{r},\boldsymbol{r}^{\prime })$ ,

(6.5) $$\begin{eqnarray}\displaystyle g^{(3)}(\boldsymbol{r},\boldsymbol{r}^{\prime })\approx g(r)g(r^{\prime })[1+h(|\boldsymbol{r}-\boldsymbol{r}^{\prime }|)], & & \displaystyle\end{eqnarray}$$

which leads to

(6.6) $$\begin{eqnarray}\displaystyle g_{\dot{\unicode[STIX]{x1D6FE}}}(\boldsymbol{r}) & {\approx} & \displaystyle g(r)\left\{1-n\frac{\unicode[STIX]{x1D701}_{e}\dot{\unicode[STIX]{x1D6FE}}d^{3}}{8T}\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\int \text{d}^{3}\boldsymbol{r}^{\prime }g(r^{\prime })[1+h(|\boldsymbol{r}-\boldsymbol{r}^{\prime }|)]\unicode[STIX]{x1D6FF}(r^{\prime }-d){\mathcal{P}}(\hat{x}^{\prime },{\hat{y}}^{\prime })\hat{r}_{\unicode[STIX]{x1D6FC}}^{\prime }\hat{r}_{\unicode[STIX]{x1D6FD}}^{\prime }\right\}\nonumber\\ \displaystyle & = & \displaystyle g(r)\left\{1-n\frac{\unicode[STIX]{x1D701}_{e}\dot{\unicode[STIX]{x1D6FE}}d^{5}}{8T}g(d)\,\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\int \text{d}{\mathcal{S}}^{\prime }[1+h(|\boldsymbol{r}-\boldsymbol{r}^{\prime }|)]{\mathcal{P}}(\hat{x}^{\prime },{\hat{y}}^{\prime })\hat{r}_{\unicode[STIX]{x1D6FC}}^{\prime }\hat{r}_{\unicode[STIX]{x1D6FD}}^{\prime }\right\}.\end{eqnarray}$$

Here, $h(r)=g(r)-1$ is the pair-correlation function, and $\int \text{d}{\mathcal{S}}^{\prime }\cdots \,$ denotes angular integral with respect to $\hat{\boldsymbol{r}}^{\prime }$ .

Figure 7. Definition of the angle of the separation vector. $\unicode[STIX]{x1D703}=135^{\circ }$ and $45^{\circ }$ are the compressional and extensional directions of the shear, respectively: $\unicode[STIX]{x1D703}=0^{\circ }$ is the ‘flow direction’, which is the direction of the flow of the solvent.

Let us consider the angular dependence of the contact value of the radial distribution function in the $(x,y)$ -plane. From (6.6), this is given by

(6.7) $$\begin{eqnarray}\displaystyle g_{\dot{\unicode[STIX]{x1D6FE}}}(d,\unicode[STIX]{x1D703})\approx g(d)\left\{1-n\frac{\unicode[STIX]{x1D701}_{e}\dot{\unicode[STIX]{x1D6FE}}d^{5}}{8T}g(d)\,\unicode[STIX]{x1D6F1}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\int \text{d}{\mathcal{S}}^{\prime }[1+h(|\boldsymbol{r}-\boldsymbol{r}^{\prime }|)]{\mathcal{P}}(\hat{x}^{\prime },{\hat{y}}^{\prime })\hat{r}_{\unicode[STIX]{x1D6FC}}^{\prime }\hat{r}_{\unicode[STIX]{x1D6FD}}^{\prime }\right\},\qquad & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D703}$ is the angle of $\hat{\boldsymbol{r}}$ with respect to the $x$ -axis, and $\unicode[STIX]{x1D703}=135^{\circ }$ and $45^{\circ }$ are the directions of compression and extension, respectively (cf. figure 7). We adopt the approximate formula for the delta-function contribution of $h(r)$ (Donev et al. Reference Donev, Torquato and Stillinger2005),

(6.8) $$\begin{eqnarray}\displaystyle h(r)\approx \frac{1}{4\unicode[STIX]{x1D711}\unicode[STIX]{x1D6FF}}\left[\frac{6A}{\displaystyle \left(\frac{r/d-1}{\unicode[STIX]{x1D6FF}}+C\right)^{4}}+\frac{B}{\displaystyle \left(\frac{r/d-1}{\unicode[STIX]{x1D6FF}}+C\right)^{2}}\right], & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}\approx (\unicode[STIX]{x1D711}_{J}-\unicode[STIX]{x1D711})/(3\unicode[STIX]{x1D711}_{J})$ and $A=3.43,B=1.45,C=2.25$ . For the evaluation of the right-hand side of (6.7), we perform numerical integration.

The comparison of the both sides of (6.7) is shown in figure 8. The left-hand side is the radial distribution function at the steady state and measured by MD simulation (circles). The conditions for the simulation are described in the caption. The angular integration on the right-hand side has been evaluated as a double integral with respect to the two angles $(\unicode[STIX]{x1D719}_{1},\unicode[STIX]{x1D719}_{2})$ , where $\cos \unicode[STIX]{x1D719}_{1}=\hat{\boldsymbol{r}}\boldsymbol{\cdot }\hat{\boldsymbol{r}}^{\prime }$ and $\unicode[STIX]{x1D719}_{2}$ is the azimuthal angle of $\hat{\boldsymbol{r}}^{\prime }$ around $\hat{\boldsymbol{r}}$ . These angles are bounded in the range $\unicode[STIX]{x03C0}/3<\unicode[STIX]{x1D719}_{1}<\unicode[STIX]{x03C0}$ and $0<\unicode[STIX]{x1D719}_{2}<2\unicode[STIX]{x03C0}$ , respectively. The region $0<\unicode[STIX]{x1D719}_{1}<\unicode[STIX]{x03C0}/3$ is excluded to avoid overlap. The contact value of the equilibrium radial distribution function, $g(d)$ , is given by $g_{0}(\unicode[STIX]{x1D711})\approx 0.848/(\unicode[STIX]{x1D711}_{J}-\unicode[STIX]{x1D711})$ for $\unicode[STIX]{x1D711}_{J}=0.639$ , which is approximately 100 for $\unicode[STIX]{x1D711}=0.63$ . The result for the right-hand side shown in solid line in figure 8 is reduced to $1/3$ to fit the amplitude. We see that the peak in the compression direction ( $\unicode[STIX]{x1D703}=135^{\circ }$ ) is captured by the theory, although the peak in the direction $\unicode[STIX]{x1D703}=0^{\circ }$ is not. The reason why the theory still reasonably predicts the stress is that the particles in contact with $\unicode[STIX]{x1D703}=0^{\circ }$ does not contribute to the stress, because they are driven with the same velocity by the uniform shear.

Figure 8. Angular distribution of the contact value of the radial distribution function at the steady state measured by MD simulation (circle). The conditions are $\unicode[STIX]{x1D711}=0.63$ , $N=1000$ , $\dot{\unicode[STIX]{x1D6FE}}^{\ast }=1\times 10^{-5}$ . The number of data averaged for each point is $10^{4}$ . Also shown are the theoretical result of Grad’s expansion (solid line) and the direct measurement of Grad’s expansion (squares).

For further validation, we have measured the right-hand side of (4.7) directly by MD simulation. We measure the distribution of the microscopic stress $\tilde{\unicode[STIX]{x1D70E}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}$ with respect to the angle of the separation vector $\boldsymbol{r}_{ij}$ , which resides in the $(x,y)$ -plane. To extract the contribution from the contacting pairs, we have assigned 1 to $f_{eq}(\boldsymbol{r}_{ij})$ when the two particles are in contact, and 0 otherwise. The factor of $(\unicode[STIX]{x03C0}/180)^{2}/Nn\times \unicode[STIX]{x0394}r\unicode[STIX]{x0394}\unicode[STIX]{x1D703}\unicode[STIX]{x0394}\unicode[STIX]{x1D719}$ is multiplied to obtain values which correspond to the radial distribution function. Here $\unicode[STIX]{x1D719}$ is the angle of $\boldsymbol{r}_{ij}$ around the $x$ -axis and $\unicode[STIX]{x0394}r$ , $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}$ , $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}$ are the size of the bins in each direction. The result is shown in squares in figure 8. We see that the directly measured distribution has the same tendency with the theoretically evaluated distribution, i.e. both exhibit a peak in the compression direction, $\unicode[STIX]{x1D703}=135^{\circ }$ , but the peak in $\unicode[STIX]{x1D703}=0^{\circ }$ is not visible.

It should be noticed that the directly measured distribution does not reflect the factorization approximation, equation (6.5). Hence, this result implies the validity of (6.5). As a check for the four-body correlation, we have measured the four-point susceptibility, $\unicode[STIX]{x1D712}_{4}$ , by MD simulation (cf. appendix I). The result is shown in figure 9. We see that no divergence is found in $\unicode[STIX]{x1D712}_{4}$ near the jamming point. This implies that the divergence of the stress is determined by the divergence of the radial distribution function at contact, $g_{0}$ . Thus, our theoretical treatment based on the factorization approximation is sufficient to discuss the singularities of the stress at the jamming point.

Figure 9. Four-point susceptibility measured by MD simulation. The conditions are $\unicode[STIX]{x1D711}=0.60-0.635$ , $N=1000$ , $\dot{\unicode[STIX]{x1D6FE}}^{\ast }=1\times 10^{-5}$ . The number of data averaged for each point is 1000. The results are for $\unicode[STIX]{x1D711}=0.60$ (thin solid line), 0.62 (dashed line), 0.63 (thick solid line) and 0.635 (dot-dashed line).