# Discontinuous Shear Thickening of a Moderately Dense Inertial Suspension of
Hydrodynamically Interacting Frictionless Soft Particles: Some New Findings

Satoshi Takada<EMAIL_ADDRESS>Department of Mechanical Systems
Engineering and Institute of Engineering, Tokyo University of Agriculture and
Technology, 2–24–16 Naka-cho, Koganei, Tokyo 184–8588, Japan  Kazuhiro Hara
<EMAIL_ADDRESS>Department of Industrial Technology and Innovation,
Tokyo University of Agriculture and Technology, 2–24–16 Naka-cho, Koganei,
Tokyo 184–8588, Japan  Hisao Hayakawa<EMAIL_ADDRESS>Yukawa
Institute for Theoretical Physics, Kyoto University, Kitashirakawa Oiwakecho,
Sakyo-ku, Kyoto 606–8502, Japan

###### Abstract

We demonstrate that discontinuous shear thickening (DST) can take place even
in a moderately dense inertial suspension of hydrodynamically interacting
frictionless soft particles. The results which demonstrate this fact are
obtained using the Lubrication-Friction Discrete Element Method. Our
simulation indicates that DST can be observed for lower densities if the
inertia of suspended particles and their softness are both of a marked nature.
We also confirm that the DST behavior is effectively approximated by the
kinetic theory under these conditions without consideration of hydrodynamic
interactions.

Introduction.— Discontinuous shear thickening (DST) is a process, whereby the
viscosity of a dense suspension changes abruptly at a certain critical shear
rate [1, 2, 3, 4]. The normal stress difference also changes abruptly at the
same time [5, 6]. DST is observed in suspensions of solid particles in liquid
media, and also in frictional dry granular materials [7], where it is closely
related to the phenomenon of shear jamming [8, 9, 10, 11, 12]. DST is thus of
industrial importance, for the design and performance of such apparatus such
as sporting equipment [13], traction controls [14], and protective vests [15]
as well as materials handling in general. DST also attracts much interest
among physicists as an instructive example of a nonequilibrium phase
transition. Therefore, the correct physical understanding of DST is both of
theoretical interest and a key requirement of the improved design of
industrial applications based on densely packed particles.

The physical origin of DST is still the subject of debate. At present, it
seems as though the most likely origin of DST lies in inter-particle
frictional forces [7, 4, 16, 17]. Other potential causal mechanisms for DST
such as order-disorder transition [18, 19, 20], and hydrodynamic clusters [22,
23, 24, 1, 25, 26, 21] have also been proposed.

It is known that a DST-like phenomenon caused by the ignited-quenched
transition exists in inertial suspensions; a model of aerosols has been
developed in which collisions between particles play an important role [27,
28, 29]. There are also several theoretical studies of inertial suspensions
consisting of hard-core frictionless particles based on the kinetic theory,
that is, without hydrodynamic interactions between particles [30, 31, 32, 33,
34, 35, 36, 37, 38, 39, 40]. These theoretical studies indicate that DST or
the ignited-quenched transition becomes the continuous shear thickening if the
volume fraction $\varphi$ is larger than a few percent [32, 36, 39, 40]. This
behavior is completely different from the DST commonly observed in colloidal
suspensions, where DST can be observed only in dense suspensions.

In this Letter, we demonstrate the occurrence of DST and the ignited-quenched
transition when hydrodynamically interacting particles are soft, frictionless,
and in moderately dense concentrations. This finding may be regarded as a new
mechanism for DST in moderately dense suspensions. In earlier work, we have
reported a complementary analysis based on the kinetic theory of moderately
dense inertial suspensions for frictionless soft particles without
hydrodynamic interaction [41]; this work is based on an extension of the
kinetic theory of dilute inertial suspensions [42]. Thus, clarifying the role
of hydrodynamic interactions among particles is also a motivation of our
present study.

Langevin model.— We consider $N$ monodisperse frictionless soft particles
(each particle of mass $m$ and diameter $d$), which are suspended in a fluid
and are confined in a three-dimensional cubic box with the linear size $L$. We
assume that the contact force between particles is described by the harmonic
potential

$U(r)=\frac{\varepsilon}{2}\left(1-\frac{r}{d}\right)^{2}\Theta\left(1-\frac{r}{d}\right),$
(1)

where $r$ is the inter-particle distance and $\varepsilon$ is the energy scale
to characterize the repulsive interaction, and $\Theta(x)$, defined as
$\Theta(x)=1$ ($x\geq 0$) and $0$ ($x<0$) is the step function. Although
clustering effects caused by attractive interactions between particles cannot
be ignored in realistic situations, such effects are suppressed if particles
are charged [43, 44, 45]. Moreover, particles are also prevented from
clustering if the temperature is high enough [46, 30, 47, 48].

The equation of motion of the suspended particle $i$ (its position
$\bm{r}_{i}$) under simple shear with the shear rate $\dot{\gamma}$ is given
by

$\frac{d\bm{p}_{i}}{dt}=\sum_{j\neq
i}\bm{F}_{ij}-\sum_{j}\overleftrightarrow{\zeta_{ij}}\bm{p}_{j}+\bm{\zeta}^{\mathrm{Sh}}_{i}+\bm{\xi}_{i},$
(2)

where we have introduced the peculiar momentum $\bm{p}_{i}\equiv
m(\bm{v}_{i}-\dot{\gamma}y_{i}\hat{\bm{e}}_{x})$ with the unit vector
$\hat{\bm{e}}_{x}$ parallel to the $x$ direction and the velocity $\bm{v}_{i}$
of the $i-$th particle, $\bm{F}_{ij}\equiv-\partial
U(r_{ij})/\partial\bm{r}_{ij}$ is the inter-particle force between the $i$–th
and $j$–th particles with $\bm{r}_{ij}\equiv\bm{r}_{i}-\bm{r}_{j}$ and
$r_{ij}\equiv|\bm{r}_{ij}|$.

The terms containing $\overleftrightarrow{\zeta_{ij}}$ and $\bm{\zeta}^{\rm
Sh}$ in Eq. (2) originate from the hydrodynamic force acting on the $i$–th
particle. If smooth particles are embedded in a Stokes fluid, the particles
cannot collide with each other because the lubrication force between
contacting particles pushes them apart, hydrodynamically, before they make
contact [49, 50, 51, 52, 53]. We note that previous studies of colloidal
suspensions have ignored the effects of inertia in the equation of motion.
Therefore, it is important to clarify the role of both hydrodynamic
interactions and inertia effects on the rheology of the system we are modeling
[54]. For this purpose, we focus on situations in which the fluid motion can
be described by the Stokes equation, and in which, therefore, the resistance
matrix $\overleftrightarrow{\zeta_{ij}}$ depends only on the positions of the
$i-$th and $j-$th particles. Then, we adopt the Lubrication-Friction Discrete
Element Method (LF-DEM) which considers only the short-range lubrication force
between particles [4, 48]. Such a simplification, whereby we ignore long-range
hydrodynamic interactions, may be justified in relatively dense particle
concentrations.

Within the framework of the LF-DEM the resistance matrix
$\overleftrightarrow{\mathsf{\zeta}_{ij}}$ in the Stokes flow is expressed as
[52, 48, 55]

$\zeta_{ij,\alpha\beta}=\begin{cases}\displaystyle
3\pi\frac{\eta_{0}d_{\mathrm{H}}}{m}\delta_{\alpha\beta}+\displaystyle\sum_{k\neq
i}\frac{1}{m}A_{ik,\alpha\beta}^{(11)}\Theta(r_{\mathrm{c}}-r_{ik})&(i=j)\\\
\displaystyle-\frac{1}{m}A_{ij,\alpha\beta}^{(11)}\Theta(r_{\mathrm{c}}-r_{ij})&(i\neq
j)\end{cases},$ (3)

where $\eta_{0}$ is the viscosity of the solvent, $r_{\mathrm{c}}\equiv
d_{\mathrm{H}}(1+\delta)$ is the cutoff length of the lubrication force with
the hydrodynamic diameter $d_{\mathrm{H}}$ (which is $d_{\rm H}\leq d$), and
$A^{(11)}_{ij,\alpha\beta}$ is given by Refs. [55, 56]. Although smooth hard-
core particles are not allowed to contact each other [49, 50, 51, 52, 53], the
LF-DEM does allow contact between rough particles, and so we introduce the
roughness parameter $\delta\equiv(d-d_{\mathrm{H}})/d_{\mathrm{H}}$ [4, 48,
12]. This parameter equates to a simplified description of dimples located on
the surface of particles. Previous papers have adopted $\delta=0.05$ [4, 48,
12], but we are free to choose larger or smaller values of $\delta$.
Similarly, the vector $\bm{\zeta}^{\mathrm{Sh}}_{i}$ is the contribution from
the shear flow due to the lubrication force:

$\zeta^{\mathrm{Sh}}_{i,\alpha}=2\eta_{0}\dot{\gamma}\sum_{j}\widetilde{G}^{(11)}_{ij,xy\alpha},$
(4)

where $\widetilde{G}^{(11)}_{ij,xy\alpha}$ is presented in Refs. [55, 57].

The noise parameter $\bm{\xi}_{i}(t)=\xi_{i,\alpha}(t)\hat{\bm{e}}_{\alpha}$
satisfies the fluctuation-dissipation relation:

$\langle\bm{\xi}_{i}(t)\rangle=\bm{0},\quad\langle\xi_{i,\alpha}(t)\xi_{j,\beta}(t^{\prime})\rangle=2mT_{\mathrm{env}}\zeta_{ij,\alpha\beta}\delta(t-t^{\prime}),$
(5)

where $\langle\cdot\rangle$ expresses the average over the noise and
$T_{\mathrm{env}}$ is the environmental (solvent) temperature. Note that the
resistance matrix $\overleftrightarrow{\zeta_{ij}}$ is not diagonal if
hydrodynamic interaction exists between the particles. The implementation of
the noise parameter $\bm{\xi}_{i}(t)$ satisfying Eq. (5) is explained in Ref.
[55]. We also note that $\bm{\zeta}^{\mathrm{Sh}}_{i}$ in Eq. (2) expresses
the hydrodynamic force in the presence of the shear flow, which is
proportional to the shear rate $\dot{\gamma}$ [49, 50, 51, 52, 53]. Now, let
us introduce two dimensionless parameters for later convenience:

$\varepsilon^{*}\equiv\frac{\varepsilon}{md^{2}\zeta^{2}},\quad\xi_{\mathrm{env}}\equiv\sqrt{\frac{T_{\mathrm{env}}}{m}}\frac{1}{d\zeta},$
(6)

where $\zeta\equiv 3\pi\eta_{0}d/m$ is the Stokesian drag acting on a sphere.

The assumptions behind Eqs. (2) and (5) are summarized as follows. (i)
Particles are suspended in a fluid in which the motion of the particles is
agitated by white Gaussian noise as in Eq. (5), (ii) the environmental
temperature $T_{\mathrm{env}}$ is independent of the motion of the suspended
particles, (iii) the hydrodynamic interactions among the particles are only
described by the lubrication force, and (iv) the effects of gravity can be
ignored. Although gravitational sedimentation plays some role in aerosols,
such an effect is negligible within the observation time for small, suspended
particles [28, 29]. In addition, inertial suspension can be regarded as a
model of a colloidal suspension, in which the sedimentation effect is
negligible and the hydrodynamic interactions among particles play important
roles.

We adopt SLLOD dynamics [59, 58] to simulate shear flow under the Lees-Edwards
boundary condition [60]. As far as we have checked, the uniform flow is stable
once the system reaches a steady state.

Rheology.— In this system, the stress tensor is given by

$\displaystyle\sigma_{\alpha\beta}$
$\displaystyle=\frac{1}{NL^{3}}\sum_{i}\Bigg{\langle}\left(-mv_{i,\alpha}v_{i,\beta}-\frac{1}{2}\sum_{j\neq
i}r_{ij,\alpha}F^{(\rm el)}_{ij,\beta}\right.$ $\displaystyle\hskip
50.00008pt\left.+\overline{\sigma}^{\mathrm{St}}_{i,\alpha\beta}+\frac{1}{2}\sum_{j\neq
i}\overline{\sigma}_{ij,\alpha\beta}^{\mathrm{H}}\right)\Bigg{\rangle},$ (7)

$\overline{\sigma}^{\mathrm{St}}_{i,\alpha\beta}$ is given by [52]

$\overline{\sigma}^{\mathrm{St}}_{i,\alpha\beta}=\frac{5}{12}\pi
d_{\mathrm{H}}^{3}\eta_{0}\dot{\gamma}(\delta_{\alpha x}\delta_{\beta
y}+\delta_{\alpha y}\delta_{\beta x}),$ (8)

and the hydrodynamic stress $\overline{\sigma}^{\mathrm{H}}_{ij,\alpha\beta}$
between the $i$–th and the $j$–th particles is given by [52]

$\overline{\sigma}^{\mathrm{H}}_{ij,\alpha\beta}=-2\eta_{0}G^{(11)}_{ij,\alpha\beta\gamma}V_{ij,\gamma}+2\eta_{0}\dot{\gamma}M^{(1)}_{ij,\alpha\beta},$
(9)

with $\bm{V}_{i}\equiv\bm{v}_{i}-\dot{\gamma}y_{i}\hat{\bm{e}}_{x}$ and
$V_{ij,\gamma}$ is $\gamma$ component of $\bm{V}_{i}-\bm{V}_{j}$. The explicit
expressions of $G^{(11)}_{ij,\alpha\beta\gamma}$ and
$M^{(1)}_{ij,\alpha\beta}$ are presented in Refs. [55, 61].

Let us introduce the Peclet number $\mathrm{Pe}$ defined as [62]

$\mathrm{Pe}\equiv\frac{3\pi\eta_{0}d^{3}}{4T_{\mathrm{env}}}\dot{\gamma}.$
(10)

Since $\xi_{\rm env}$ introduced in Eq. (6) is proportional to
$\sqrt{mT_{\mathrm{env}}/\eta_{0}^{2}d^{4}}$, as shown in Eq. (6), $\xi_{\rm
env}$ increases as the inertia becomes important. Using $\dot{\gamma}^{*}$, we
introduce the dimensionless viscosity $\eta^{*}$ as

$\eta^{*}\equiv\frac{\sigma_{xy}}{nT_{\mathrm{env}}\mathrm{Pe}}.$ (11)

Figure 1: Plots of (a) $\tilde{\eta}$ and (b) $\theta$ against $\mathrm{Pe}$
for $\delta=0.02$ (open circles), $0.05$ (open squares), $0.10$ (open
triangles), and $0.25$ (crosses) with fixed values of $\varphi=0.30$,
$\varepsilon^{*}=10^{4}$, and $\xi_{\mathrm{env}}=1.0$. The dashed line
represents the prediction of the kinetic theory in Ref. [41].

We examine $N=1000$ and $10$ ensemble averages in the simulations. We control
$\mathrm{Pe}$, $\xi_{\mathrm{env}}$, $\varepsilon^{*}$, and $\delta$ in the
range $0.25\leq\mathrm{Pe}\leq 25$, $1\leq\xi_{\mathrm{env}}\leq 10^{2}$,
$10^{2}\leq\varepsilon^{*}\leq 10^{8}$, and $0.02\leq\delta\leq 0.25$,
respectively.

Figure 1 shows how the scaled viscosity
$\tilde{\eta}\equiv\eta^{*}/\eta_{\mathrm{a}}$ and the dimensionless
temperature $\theta\equiv T/T_{\mathrm{env}}$ depend on $\mathrm{Pe}$ with
fixed values of $\varphi(=\pi Nd^{3}/(6L^{3}))=0.30$,
$\varepsilon^{*}=10^{4}$, and $\xi_{\mathrm{env}}=1.0$ for various $\delta$,
where we have introduced the kinetic temperature $T\equiv\sum_{i=1}^{N}\langle
m_{i}\bm{V}_{i}^{2}\rangle/(3N)$ and the empirical expression of the apparent
viscosity $\eta_{\mathrm{a}}$ in the low shear limit
$\eta_{\mathrm{a}}=1+(5/2)\varphi+4\varphi^{2}+42\varphi^{3}$ [63]. Note that
we do not have any theoretical result for $\eta_{\rm a}$ until
$O(\varphi^{3})$ as far as we have been able to ascertain, though the
expression until $O(\varphi^{2})$ is well known for hard-core suspensions
[64].

Figure 2: Plots of (a) $\tilde{\eta}$ and (b) $\theta$ against $\mathrm{Pe}$
for $\delta=0.05$ (open squares) and $\delta=0.02$ (open circles) with fixed
values of $\varphi=0.40$, $\varepsilon^{*}=10^{4}$, and
$\xi_{\mathrm{env}}=1.0$. The dashed line represents the prediction of the
kinetic theory in Ref. [41].

Figure 2 is the corresponding result for $\varphi=0.40$ with $\delta=0.05$ and
$\delta=0.02$. Figures 1(a) and 2(a) exhibit remarkable results, in which DST
can be observed around $\mathrm{Pe}=10$. It is noteworthy that DST is observed
even for $\delta=0.02$ (Figs. 1 (a) and 2(a)), though the discontinuous jump
of the viscosity is enhanced for larger values of $\delta$. Although the
quantitative agreement between the kinetic theory without consideration of
hydrodynamic interactions [41] and hydrodynamic simulation is poor for
$\varphi=0.40$ (see Fig. 2(a)), the kinetic theory captures the qualitative
behavior of the DST. Moreover, the agreement between the theory and simulation
is reasonable for $\varphi=0.30$ (Fig. 1 (a)). We note that the results of the
simulation based on the LF-DEM approach the theoretical prediction reported in
Ref. [41] as $\delta$ increases, though we have used the empirical expression
for $\eta_{\rm a}$. Figures 1(b) and 2(b) indicate that these DSTs are caused
by the ignited-quenched transition.

Thus, we have confirmed that DST and the ignited-quenched transition of
inertial suspensions for soft and frictionless particles can survive even if
hydrodynamic interactions between particles exist. We have also confirmed the
relevance of the kinetic theory developed in Ref. [41] to describe the
suspensions with hydrodynamic interactions.

Figure 3(a) is the plot of $\tilde{\eta}$ against ${\rm Pe}$ for various
$\xi_{\rm env}$ with fixed values of $\delta=0.05$ and $\varphi=0.40$. We
verify that DST still persists for $\xi_{\rm env}=10$ but disappears around
$10\lesssim\xi_{\rm env}\lesssim 20$. Figure 3(b) is the plot of
$\tilde{\eta}/\varepsilon^{*4/3}$ against ${\rm Pe}$ for various
$\varepsilon^{*}$ with fixed values of $\delta=0.05$ and $\varphi=0.40$.
Except for the case of $\varepsilon^{*}=10^{2}$, $\tilde{\eta}$ exhibits DST.
As the softness parameter $\varepsilon^{*}$ increases, the upper branch also
increases (see Fig. 3(b)). In accordance with Ref. [41] and as shown in Figs.
3(b) and (c), the viscosity in the upper branch is scaled as
$\tilde{\eta}/\varepsilon^{*4/3}$ except for $\varepsilon^{*}=10^{2}$. Thus,
we conclude that DST with frictionless soft particles can be observed for
systems with $\xi_{\mathrm{env}}\leq 20$ and $\varepsilon^{*}\geq 10^{4}$.

Figure 3: (a) Plot of $\tilde{\eta}$ against $\mathrm{Pe}$ for various
$\xi_{\mathrm{env}}$. (b) Plot of $\tilde{\eta}/\varepsilon^{*4/3}$ for
various $\varepsilon^{*}$. Here, we fix (a) $\varepsilon^{*}=10^{4}$ and (b)
$\xi_{\mathrm{env}}=1$, respectively. (c) A detailed view of the high
$\mathrm{Pe}$ regime of (b). Parameters are fixed as $\varphi=0.40$ and
$\delta=0.05$.

Concluding remarks.— In this Letter, we have successfully demonstrated the
existence of DST-like changes in viscosity and kinetic temperature for
moderately dense inertial suspensions composed of frictionless soft particles
using the LF-DEM. The kinetic theory without consideration of hydrodynamic
interaction effectively describes the discontinuous behaviors of both
viscosity and kinetic temperature. Our results also reveal a new mechanism for
DST caused by the ignited-quenched transition for frictionless soft particles.

It is difficult to observe the simultaneously large and discontinuous changes
of $\eta^{*}$ and $\theta$ obtained in our model if the solvent is a liquid.
Indeed, if the kinetic temperature in the ignited phase becomes $10^{6}$ times
larger than that in the quenched phase, the liquid in the ignited phase might
be evaporated because of the strong stirring effect of suspended particles,
though our model does not include such effects. Even if we consider particles
suspended in a gas, it might be difficult to avoid the melting of solid
particles. Nevertheless, an indication of the existence of DST-like changes of
$\eta^{*}$ and $\theta$ for frictionless soft particles, even among relatively
dense suspensions, is important.

Let us discuss whether this behavior is observable in experiments of colloidal
suspensions. Here, let us assume that the diameter, the mass density, and
Young’s modulus of colloidal particles are given by $d=2\times 10^{-6}\
\mathrm{m}$, $\rho=10^{3}\ \mathrm{kg/m^{3}}$, and $Y=10\ \mathrm{GPa}$, and
thus, the mass $m$ becomes $m\sim 4\times 10^{-15}\ \mathrm{kg}$. When we
consider a case when the solvent is water at room temperature, of which the
viscosity is $\eta_{0}\sim 10^{-3}\ \mathrm{Pa\cdot s}$, the corresponding
drag coefficient becomes $\zeta=3\pi\eta_{0}d/m\sim 5\times 10^{6}\
\mathrm{s}^{-1}$. We can also evaluate $\varepsilon^{*}\sim 6\times
10^{5}\delta\sim 6\times 10^{3}$ with the assumption of $\delta\sim 10^{-2}$
[65]. The value of $\varepsilon^{*}$ is close to that used in Fig. 3. Note
that the behavior of the viscosity can be scaled by ${\varepsilon^{*}}^{4/3}$.
This means that the DST discussed in this paper can be observed in the region
of large Pe. A serious problem exists in the value of $\xi_{\rm env}$. If we
assume that $T_{\rm env}=4\times 10^{-21}$ J corresponding to the room
temperature, $\xi_{\rm env}$ becomes about $2\times 10^{-4}$. This value is
much smaller than the critical value
($10<\xi_{\mathrm{env}}^{\mathrm{c}}<10^{2}$) of the absence of the DST, which
means that we can expect the existence of the DST. We have verified that the
DST still can be observed for realistic values of parameters in the prediction
of the kinetic theory [41]. Thus, the observation of such a DST for
frictionless small grains which are supposed to be Brownian suspensions is
challenging for experimentalists.

Needless to say, it is important to analyze suspensions of frictional grains
for larger $\varepsilon^{*}$ corresponding to the typical experimental setup
for colloidal suspensions. This will be our forthcoming task.

Acknowledgements.— The authors thank Takeshi Kawasaki, Takashi Uneyama, Michio
Otsuki, and Pradipto for their helpful comments. This research was partially
supported by Grants-in-Aid from MEXT for Scientific Research (Grant Nos.
JP20K14428 and JP21H01006).

## References

  * [1] N. J. Wagner and J. F. Brady, Shear thickening in colloidal dispersions, Phys. Today 62, 27 (2009).
  * [2] E. Brown and H. M. Jaeger, Shear thickening in concentrated suspensions: phenomenology, mechanisms and relations to jamming, Rep. Prog. Phys. 77, 046602 (2014).
  * [3] C. Ness, R. Seto, and R. Mari, The Physics of Dense Suspensions, Ann. Rev. Cond. Mat. Phys. 13, 97 (2022).
  * [4] R. Seto, R. Mari, J. F. Morris, and M. M. Denn, Discontinuous Shear Thickening of Frictional Hard-Sphere Suspensions, Phys. Rev. Lett. 111, 218301 (2013).
  * [5] H. M. Laun, Normal stresses in extremely shear thickening polymer dispersions, J. Non-Newton. Fluid Mech. 54, 87 (1994).
  * [6] C. D. Cwalina and N. J. Wagner, Material properties of the shear-thickened state in concentrated near hard-sphere colloidal dispersions, J. Rheol. 58, 949 (2014).
  * [7] M. Otsuki and H. Hayakawa, Critical scaling near jamming transition for frictional granular particles, Phys. Rev. E 83, 051301 (2011).
  * [8] D. Bi, J. Zhang, B. Chakraborty, and R. P. Behringer, Jamming by shear, Nature 480, 355 (2011).
  * [9] A. Fall, F. Bertrand, D. Hautemayou, C. Mezière, P. Moucheront, A. Lemaítre, and G. Ovarlez, Macroscopic Discontinuous Shear Thickening versus Local Shear Jamming in Cornstarch, Phys. Rev. Lett. 114, 098301 (2015).
  * [10] I. R. Peters, S. Majumdar, and H. M. Jaeger, Direct observation of dynamic shear jamming in dense suspensions, Nature 532, 214 (2016).
  * [11] M. Otsuki and H. Hayakawa, Shear jamming, discontinuous shear thickening, and fragile states in dry granular materials under oscillatory shear, Phys. Rev. E 101, 032905 (2020).
  * [12] Pradipto and H. Hayakawa, Simulation of dense non-Brownian suspensions with the lattice Boltzmann method: shear jammed and fragile states, Soft Matter 16, 945 (2020), [Correction] 18, 685 (2022).
  * [13] C. Fischer, S. A. Braun, P.-E. Bourban, V. Michaud, C. J. G. Plummer, and J.-A. E. Månson, Smart Materials and Structures Dynamic properties of sandwich structures with integrated shear-thickening fluids, Smart Mater. Struct. 15, 1467 (2006).
  * [14] H. M. Laun, R. Bung, and F. Schmidt, Rheology of extremely shear thickening polymer dispersions (passively viscosity switching fluids), J. Rheol. 35, 999 (1991).
  * [15] Y. S. Lee, E. D. Wetzel, and N. J. Wagner, The ballistic impact characteristics of Kevlar® woven fabrics impregnated with a colloidal shear thickening fluid, J. Mater. Sci. 38, 2825 (2003).
  * [16] T. Kawasaki and L. Berthier, Discontinuous shear thickening in Brownian suspensions, Phys. Rev. E 98, 012609 (2018).
  * [17] M. Wyart and M. E. Cates, Discontinuous Shear Thickening without Inertia in Dense Non-Brownian Suspensions, Phys. Rev. Lett. 112, 098302 (2014).
  * [18] R. L. Hoffman, Discontinuous and Dilatant Viscosity Behavior in Concentrated Suspensions. I. Observation of a Flow Instability, Trans. Soc. Rheol. 16, 155 (1972).
  * [19] R. L. Hoffman, Discontinuous and dilatant viscosity behavior in concentrated suspensions. II. Theory and experimental tests, J. Colloid Interface Sci. 46, 491 (1974).
  * [20] R. L. Hoffman, Explanations for the cause of shear thickening in concentrated colloidal suspensions, J. Rheol. 42, 111 (1998).
  * [21] M. Wang, S. Jamali, and J. F. Brady, A hydrodynamic model for discontinuous shear-thickening in dense suspensions, J. Rheol. 64, 379 (2020).
  * [22] J. F. Brady and G. Bossis, The rheology of concentrated suspensions of spheres in simple shear flow by numerical simulation, J. Fluid Mech. 155, 105 (1985).
  * [23] J. F. Brady and G. Bossis, Stokesian Dynamics, Ann. Rev. Fluid Mech. 20, 111 (1988).
  * [24] J. Bender and N. J. Wagner, Reversible shear thickening in monodisperse and bidisperse colloidal dispersions, J. Rheol. 40, 899 (1996).
  * [25] J. R. Melrose and R. C. Ball, Continuous shear thickening transitions in model concentrated colloids—The role of interparticle forces, J. Rheol. 48, 937 (2004).
  * [26] R. S. Farr, J. R. Melrose, and R. C. Ball, Kinetic theory of jamming in hard-sphere startup flows, Phys. Rev. E 55, 7203 (1997).
  * [27] S. K. Friedlander, Smoke, Dust and Haze: Fundamentals of Aerosol Dynamics, 2nd ed. (Oxford Univ. Press, Oxford, 2000).
  * [28] D. L. Koch and R. J. Hill, Inertial Effects in Suspension and Porous-Media Flows, Ann. Rev. Fluid Mech. 33, 619 (2001).
  * [29] M. Hu, J. Peng, K. Sun, D. Yue, S. Guo, A. Wiedensohler, and Z. Wu, Estimation of Size-Resolved Ambient Particle Density Based on the Measurement of Aerosol Number, Mass, and Chemical Size Distributions in the Winter in Beijing, Environ. Sci. Technol. 46, 9941 (2012).
  * [30] H. Hayakawa and S. Takada, Kinetic theory of discontinuous shear thickening for a dilute gas-solid suspension, Prog. Theor. Exp. Phys. 2019, 083J01 (2019).
  * [31] H.-W. Tsao and D. L. Koch, Simple shear flows of dilute gas–solid suspensions, J. Fluid Mech. 296, 211 (1995).
  * [32] A. S. Sangani, G. Mo. H.-W. Tsao, and D. L. Koch, Simple shear flows of dense gas-solid suspensions at finite Stokes numbers, J. Fluid Mech. 313, 309 (1996).
  * [33] D. L. Koch and A. S. Sangani, Particle pressure and marginal stability limits for a homogeneous monodisperse gas-fluidized bed: kinetic theory and numerical simulations, J. Fluid Mech. 400, 229 (1999).
  * [34] D. L. Koch and R. J. Hill, Inertial Effects in Suspension and Porous-Media Flows, Ann. Rev. Fluid Mech. 33, 619 (2001).
  * [35] V. Garzó, S. Tennetti, S. Subramaniam, and C. M. Hrenya, Enskog kinetic theory for monodisperse gas–solid flows, J. Fluid Mech. 712, 129 (2012).
  * [36] S. Saha and M. Alam, Revisiting ignited–quenched transition and the non-Newtonian rheology of a sheared dilute gas–solid suspension, J. Fluid. Mech. 833, 206 (2017).
  * [37] M. Alam, S. Saha, and R. Gupta, Unified theory for a sheared gas–solid suspension: from rapid granular suspension to its small-Stokes-number limit, J. Fluid Mech. 870, 1175 (2019).
  * [38] S. Saha and M. Alam, Burnett-order constitutive relations, second moment anisotropy and co-existing states in sheared dense gas–solid suspensions, J. Fluid Mech. 887, 25 (2020).
  * [39] H. Hayakawa, S. Takada, and V. Garzó, Kinetic theory of shear thickening for a moderately dense gas-solid suspension: From discontinuous thickening to continuous thickening, Phys. Rev. E 96, 042903 (2017), [Erratum] 101, 069904(E) (2020).
  * [40] S. Takada, H. Hayakawa, A. Santos, and V. Garzó, Enskog kinetic theory of rheology for a moderately dense inertial suspension, Phys. Rev. E 102, 022907 (2020).
  * [41] S. Takada, K. Hara, and H. Hayakawa, Kinetic theory of discontinuous shear thickening of a moderately dense inertial suspension of frictionless soft particles, arXiv:2207.05752.
  * [42] S. Sugimoto and S. Takada, Two-Step Discontinuous Shear Thickening of Dilute Inertial Suspensions Having Soft-Core Potential, J. Phys. Soc. Jpn. 89, 084803 (2020), [Addendum] 89, 127001 (2020).
  * [43] B. Derjaguin and L. D. Landau, Theory of the stability of strongly charged lyophobic sols and of the adhesion of strongly charged particles in solutions of electrolytes, Acta Physicochim. URSS 14, 633 (1941).
  * [44] E. J. W. Verwey and J. T. G. Overbeek, Theory of the Stability of Lyophobic Colloids: The Interaction of Sol Particles Having an Electric Double Layer (Elsevier, Amsterdam, 1948).
  * [45] J. Israelachvili, Intermolecular and Surface Forces, 3rd ed. (Academic Press, New York, 2011).
  * [46] T. Kawasaki, A. Ikeda, and L. Berthier, Thinning or thickening? Multiple rheological regimes in dense suspensions of soft particles, EPL 107, 28009 (2014).
  * [47] G. Bossis and J. F. Brady, The rheology of Brownian suspensions, J. Chem. Phys. 91, 1866 (1989).
  * [48] R. Mari, R. Seto, J. F. Morris, and M. M. Denn, Discontinuous shear thickening in Brownian suspensions by dynamic simulation, Proc. Natl. Acad. Sci. USA 112, 15326 (2015).
  * [49] J. Happel and H. Brenner, Low Reynolds number hydrodynamics —with special applications to particulate media— (Noordhoff International Publishing, Leiden, 1973).
  * [50] D. J. Jeffrey and Y. Onishi, Calculation of the resistance and mobility functions for two unequal rigid spheres in low-Reynolds-number flow, J. Fluid Mech. 139, 261 (1984).
  * [51] D. J. Jeffrey, The calculation of the low Reynolds number resistance functions for two unequal spheres, Phys. Fluids A 4, 16 (1992).
  * [52] S. Kim and S. J. Karrila, Microhydrodynamics: Principles and Selected Applications (Dover, New York, 2005).
  * [53] K. Ichiki, A. E. Kobryn, and A. Kovalenko, Resistance functions for two unequal spheres in linear flow at low Reynolds number with the Navier slip boundary condition arXiv:1302.0461.
  * [54] Our description of the motion of the particles can have an inertial term inserted via Newton’s equations. Inclusion of the inertial terms has a practical advantage, because the simulation of under-damped particles is much easier than that of over-damped particles.
  * [55] See the Supplemental Material for this paper.
  * [56] The expression of $A^{(11)}_{ij,\alpha\beta}$ is given by $A^{(11)}_{ij,\alpha\beta}=X^{A(11)}_{ij}\hat{k}_{ij,\alpha}\hat{k}_{ij,\beta}+Y^{A(11)}_{ij}\left(\delta_{\alpha\beta}-\hat{k}_{ij,\alpha}\hat{k}_{ij,\beta}\right)$, with $\hat{\bm{k}}_{12}\equiv(\bm{r}_{2}-\bm{r}_{1})/|\bm{r}_{2}-\bm{r}_{1}|$, $X^{A(11)}_{ij}=(3/4)\pi d_{\mathrm{H}}/h_{ij}$, and $Y^{A(k\ell)}_{ij}=(1/2)\pi d_{\mathrm{H}}\ln(1/h_{ij})$, where we have introduced $h_{ij}\equiv 2(r_{ij}/d_{\mathrm{H}}-1)$ (see also the Supplemental Material [55] for the details).
  * [57] The expression of $\tilde{G}^{(11)}_{ij,\alpha\beta\gamma}=G^{(11)}_{\beta\gamma\alpha}$ is given by $G^{(11)}_{ij,\alpha\beta\gamma}=X^{G(11)}_{ij}\left(\hat{k}_{ij,\gamma}\hat{k}_{ij,\alpha}-\frac{1}{3}\delta_{\gamma\alpha}\right)\hat{k}_{ij,\beta}+Y^{G(11)}_{ij}\left(\hat{k}_{ij,\gamma}\delta_{\alpha\beta}+\hat{k}_{ij,\alpha}\delta_{\gamma\beta}-2\hat{k}_{ij,\gamma}\hat{k}_{ij,\alpha}\hat{k}_{ij,\beta}\right)$, with $X^{G(k\ell)}_{ij}=(3/8)\pi d_{\mathrm{H}}^{2}/h_{ij}$ and $Y^{G(k\ell)}_{ij}=(1/8)\pi d_{\mathrm{H}}^{2}\ln(1/h_{ij})$ (see also the Supplemental Material [55] for the details).
  * [58] D. J. Evans and G. Morriss, Statistical Mechanics of Nonequilibrium Liquids 2nd ed. (Cambridge Univ. Press, Cambridge, 2008).
  * [59] D. J. Evans and G. P. Morriss, Nonlinear-response theory for steady planar Couette flow, Phys. Rev. A 30, 1528 (1984).
  * [60] A. W. Lees and S. F. Edwards, The computer study of transport processes under extreme conditions, J. Phys. C: Solid State Phys. 5, 1921 (1972).
  * [61] The expressions of $G^{(11)}_{ij,\alpha\beta\gamma}$ and $M^{(1)}_{ij,\alpha\beta}$ are, respectively, given by $G^{(11)}_{ij,\alpha\beta\gamma}=\tilde{G}^{(11)}_{ij,\gamma\alpha\beta}$ and $M^{(11)}_{ij,\alpha\beta\gamma\delta}=X^{M(11)}_{ij}\hat{k}^{(0)}_{ij,\alpha\beta\gamma\delta}+Y^{M(11)}_{ij}\hat{k}^{(1)}_{ij,\alpha\beta\gamma\delta}+Z^{M(11)}_{ij}\hat{k}^{(2)}_{ij,\alpha\beta\gamma\delta}$, with $X^{M(11)}_{ij}=(1/8)\pi d_{\mathrm{H}}^{3}/h_{ij}$, $Y^{M(11)}_{ij}=(1/10)\pi d_{\mathrm{H}}^{3}\ln(1/h_{ij})$, and $Z^{M(11)}_{ij}\simeq(5/6)\pi d_{\mathrm{H}}^{3}$, where we also have introduced $\hat{k}^{(0)}_{ij,\alpha\beta\gamma\delta}$, $\hat{k}^{(1)}_{ij,\alpha\beta\gamma\delta}$, and $\hat{k}^{(2)}_{ij,\alpha\beta\gamma\delta}$ presented in Ref. [55].
  * [62] M. L. Hunt, R. Zenit, C. S. Campbell, and C. E. Brennen, Revisiting the 1954 suspension experiments of R. A. Bagnold, J. Fluid Mech. 452, 1 (2002).
  * [63] C. G. de Kruif, E. M. F. van Iersel, A. Vrij, and W. B. Russel, Hard sphere colloidal dispersions: Viscosity as a function of shear rate and volume fraction, J. Chem. Phys. 83, 4717 (1985).
  * [64] G. K. Batchelor and J. T. Green, The determination of the bulk stress in a suspension of spherical particles to order $c^{2}$, J. Fluid Mech. 56, 401 (1972).
  * [65] We can estimate the dimensionless softness for this system. Introducing the spring constant $k$ in Eq. (1) as $k\equiv\varepsilon/d^{2}$, we get $\varepsilon^{*}=k/(m\zeta^{2})$. The spring constant is related to Young’s modulus $Y$ as $k\sim(\mathcal{A}/d)Y$, where $\mathcal{A}=\pi r_{A}^{2}$ is the area of the contact surface between two particles with the radius $r_{A}$. Thus, we rewrite $\varepsilon^{*}$ as $\varepsilon^{*}=\mathcal{A}Y/(md\zeta^{2})$. Let us evaluate $\mathcal{A}$. We may assume that the contact distance is evaluated as $d-d_{\rm H}$, which can be approximated as $1-\cos\theta_{0}\approx\theta_{0}^{2}/2$ where $\theta_{0}$ is the angle of the edge of the contact area. Thus, $\theta_{0}$ can be evaluated as $\theta_{0}\approx\sqrt{2\delta}$, where $\delta$ is the roughness parameter. Similarly, the radius of the contact area $r_{A}$ is evaluated as $(d/2)\sin\theta_{0}\approx(d/2)\theta_{0}=\sqrt{2\delta}d/2$. Thus, we obtain $\mathcal{A}\approx\pi d^{2}\delta/2$. Substituting these values into $\varepsilon^{*}\sim\pi dY\delta/(m\zeta^{2})$ we obtain $\varepsilon^{*}\sim 6\times 10^{5}\delta\sim 6\times 10^{3}$. This value is close to that which we used in this paper.

Supplemental Materials for “Discontinuous Shear Thickening of a Moderately
Dense Inertial Suspension of Hydrodynamically Interacting Frictionless Soft
Particles: Some New Findings”

This Supplemental Material explains the details of calculations that are not
included in the main text. In Sect. I, we briefly explain how we generate
random numbers when the resistance matrix in Eq. (5) has non-diagonal
elements. In Sect. II, we summarize the hydrodynamic interaction when we
consider the lubrication force between particles. In Sect. III, we check the
convergence of the results depending on the strain.

##  I Generation of random numbers obeying nonuniform distribution

In this section, let us explain how we generate random numbers obeying
nonuniform distribution based on Ref. [1]. We consider a situation where the
random number $X_{k}$ ($k=1,2,\cdots,3N$) obeys

$\langle X_{k}(t)X_{\ell}(t^{\prime})\rangle=2D_{k\ell}\delta(t-t^{\prime}).$
(S1)

Here, $X_{k}$ and $D_{k\ell}$, respectively, relates to the noise
$\bm{\xi}_{i}$ and the resistance matrix $\overleftrightarrow{\zeta}$ as

$\displaystyle\begin{bmatrix}X_{3i-2}\\\ X_{3i-1}\\\
X_{3i}\end{bmatrix}=\begin{bmatrix}\xi_{i,x}\\\ \xi_{i,y}\\\
\xi_{i,z}\end{bmatrix},$ (S2a)
$\displaystyle\begin{bmatrix}D_{3i-2,3j-2}&D_{3i-2,3j-1}&D_{3i-2,3j}\\\
D_{3i-1,3j-2}&D_{3i-1,3j-1}&D_{3i-1,3j}\\\
D_{3i,3j-2}&D_{3i,3j-1}&D_{3i,3j}\end{bmatrix}=\begin{bmatrix}\zeta_{ij,xx}&\zeta_{ij,xy}&\zeta_{ij,xz}\\\
\zeta_{ij,yx}&\zeta_{ij,yy}&\zeta_{ij,yz}\\\
\zeta_{ij,zx}&\zeta_{ij,zy}&\zeta_{ij,zz}\end{bmatrix},$ (S2b)

for $i,j=1,2,\cdots,N$. It is useful to introduce $Y_{i}$ as

$X_{k}=\sum_{\ell=1}^{k}\alpha_{k\ell}Y_{\ell},$ (S3)

where $Y_{k}$ satisfies

$\langle Y_{k}(t)\rangle=0,\quad\langle
Y_{k}(t)Y_{\ell}(t^{\prime})\rangle=2\delta_{k\ell}\delta(t-t^{\prime}).$ (S4)

Now, the coefficients $\alpha_{k\ell}$ should be given by (see Ref. [1])

$\displaystyle\alpha_{11}$ $\displaystyle=D_{11}^{1/2},$ (S5a)
$\displaystyle\alpha_{k1}$ $\displaystyle=D_{k1}/\alpha_{11},$ (S5b)
$\displaystyle\alpha_{kk}$
$\displaystyle=\left(D_{kk}-\sum_{\ell=1}^{k-1}\alpha_{k\ell}^{2}\right)^{1/2}\quad(k>1),$
(S5c) $\displaystyle\alpha_{k\ell}$
$\displaystyle=\left(D_{k\ell}-\sum_{m=1}^{\ell-1}\alpha_{km}\alpha_{\ell
m}\right)/\alpha_{\ell\ell}\quad(k>\ell>1).$ (S5d)

By solving Eqs. (S5) with Eq. (S3), we can evaluate the random force
$\bm{\xi}_{i}$.

##  II Hydrodynamic interaction

In this section, we summarize the hydrodynamic interactions between suspended
particles when the fluid motion is expressed as the Stokes equation [2, 3, 4,
5, 6]. If two particles ($i$–th and $j$–th particles, $i\neq j$) interact with
each other, the following lubrication forces and stresses act between them:

$\begin{bmatrix}\bm{F}^{\mathrm{H}}_{ij}\\\ \bm{F}^{\mathrm{H}}_{ji}\\\
\overleftrightarrow{\overline{\sigma}^{\mathrm{H}}_{ij}}\\\
\overleftrightarrow{\overline{\sigma}^{\mathrm{H}}_{ji}}\end{bmatrix}=-\eta_{0}\begin{bmatrix}\overleftrightarrow{A_{ij}^{(11)}}&\overleftrightarrow{A_{ij}^{(12)}}&\overleftrightarrow{\widetilde{G}_{ij}^{(11)}}&\overleftrightarrow{\widetilde{G}_{ij}^{(12)}}\\\
\overleftrightarrow{A_{ij}^{(21)}}&\overleftrightarrow{A_{ij}^{(22)}}&\overleftrightarrow{\widetilde{G}_{ij}^{(21)}}&\overleftrightarrow{\widetilde{G}_{ij}^{(22)}}\\\
\overleftrightarrow{G_{ij}^{(11)}}&\overleftrightarrow{G_{ij}^{(12)}}&\overleftrightarrow{M_{ij}^{(11)}}&\overleftrightarrow{M_{ij}^{(12)}}\\\
\overleftrightarrow{G_{ij}^{(21)}}&\overleftrightarrow{G_{ij}^{(22)}}&\overleftrightarrow{M_{ij}^{(21)}}&\overleftrightarrow{M_{ij}^{(22)}}\end{bmatrix}\begin{bmatrix}\bm{v}_{i}-\overleftrightarrow{E}\bm{x}_{i}\\\
\bm{v}_{j}-\overleftrightarrow{E}\bm{x}_{j}\\\ -\overleftrightarrow{E}\\\
-\overleftrightarrow{E}\end{bmatrix},$ (S6)

where $\bm{F}^{\mathrm{H}}_{ij}$ is the force acting on $i$–th particle due to
the lubrication effect between $i$–th and $j$–th particles. Because we are
only interested in the simple shear condition, the matrix
$\overleftrightarrow{E}$ is reduced to

$E_{\alpha\beta}=\dot{\gamma}\delta_{\alpha x}\delta_{\beta y},$ (S7)

where $\dot{\gamma}$ is the shear rate.

It is known that the coefficients $A^{(k\ell)}_{ij,\alpha\beta}$,
$G^{(k\ell)}_{ij,\alpha\beta\gamma}$, and
$M^{(k\ell)}_{ij,\alpha\beta\gamma\delta}$ in Eq. (S6) satisfy the following
symmetries [5]:

$\displaystyle A^{(k\ell)}_{ij,\alpha\beta}$ $\displaystyle=A^{(\ell
k)}_{ij,\beta\alpha},$ (S8a) $\displaystyle
G^{(k\ell)}_{ij,\alpha\beta\gamma}$ $\displaystyle=\widetilde{G}^{(\ell
k)}_{ij,\gamma\alpha\beta},$ (S8b) $\displaystyle
M^{(k\ell)}_{ij,\alpha\beta\gamma\delta}$ $\displaystyle=M^{(\ell
k)}_{ij,\gamma\delta\alpha\beta}.$ (S8c)

The explicit forms of the coefficients are written as

$\displaystyle A^{(k\ell)}_{ij,\alpha\beta}$
$\displaystyle=X^{A(k\ell)}_{ij}\hat{k}_{ij,\alpha}\hat{k}_{ij,\beta}+Y^{A(k\ell)}_{ij}\left(\delta_{\alpha\beta}-\hat{k}_{ij,\alpha}\hat{k}_{ij,\beta}\right),$
(S9a) $\displaystyle G^{(k\ell)}_{ij,\alpha\beta\gamma}$
$\displaystyle=X^{G(k\ell)}_{ij}\left(\hat{k}_{ij,\alpha}\hat{k}_{ij,\beta}-\frac{1}{3}\delta_{\alpha\beta}\right)\hat{k}_{ij,\gamma}+Y^{G(k\ell)}_{ij}\left(\hat{k}_{ij,\alpha}\delta_{\beta\gamma}+\hat{k}_{ij,\beta}\delta_{\alpha\gamma}-2\hat{k}_{ij,\alpha}\hat{k}_{ij,\beta}\hat{k}_{ij,\gamma}\right),$
(S9b) $\displaystyle M^{(k\ell)}_{ij,\alpha\beta\gamma\delta}$
$\displaystyle=X^{M(k\ell)}_{ij}\hat{k}^{(0)}_{ij,\alpha\beta\gamma\delta}+Y^{M(k\ell)}_{ij}\hat{k}^{(1)}_{ij,\alpha\beta\gamma\delta}+Z^{M(k\ell)}_{ij}\hat{k}^{(2)}_{ij,\alpha\beta\gamma\delta},$
(S9c)

respectively, where we have introduced $X^{A(k\ell)}_{ij}$,
$Y^{A(k\ell)}_{ij}$, $X^{G(k\ell)}_{ij}$, $Y^{G(k\ell)}_{ij}$,
$X^{M(k\ell)}_{ij}$, $Y^{M(k\ell)}_{ij}$, and $Z^{M(k\ell)}_{ij}$ which are
functions of the separation length $h_{ij}\equiv
2(r_{ij}-d_{\mathrm{H}})/d_{\mathrm{H}}=2(r_{ij}/d_{\mathrm{H}}-1)$, as well
as $\hat{{k}}^{(0)}_{ij,\alpha\beta\gamma\delta}$,
$\hat{{k}}^{(1)}_{ij,\alpha\beta\gamma\delta}$, and
$\hat{k}^{(2)}_{ij,\alpha\beta\gamma\delta}$ as [5]

$\displaystyle\hat{k}^{(0)}_{ij,\alpha\beta\gamma\delta}$
$\displaystyle\equiv\frac{3}{2}\left(\hat{k}_{ij,\alpha}\hat{k}_{ij,\beta}-\frac{1}{3}\delta_{\alpha\beta}\right)\left(\hat{k}_{ij,\gamma}\hat{k}_{ij,\delta}-\frac{1}{3}\delta_{\gamma\delta}\right),$
(S10a) $\displaystyle\hat{k}^{(1)}_{ij,\alpha\beta\gamma\delta}$
$\displaystyle\equiv\frac{1}{2}\left(\hat{k}_{ij,\alpha}\hat{k}_{ij,\gamma}\delta_{\beta\delta}+\hat{k}_{ij,\beta}\hat{k}_{ij,\gamma}\delta_{\alpha\delta}+\hat{k}_{ij,\alpha}\hat{k}_{ij,\delta}\delta_{\beta\gamma}+\hat{k}_{ij,\beta}\hat{k}_{ij,\delta}\delta_{\alpha\gamma}-4\hat{k}_{ij,\alpha}\hat{k}_{ij,\beta}\hat{k}_{ij,\delta}\hat{k}_{ij,\gamma}\right),$
(S10b) $\displaystyle\hat{k}^{(2)}_{ij,\alpha\beta\gamma\delta}$
$\displaystyle\equiv\frac{1}{2}\left(\delta_{\alpha\gamma}\delta_{\beta\delta}+\delta_{\alpha\delta}\delta_{\beta\gamma}-\delta_{\alpha\beta}\delta_{\gamma\delta}+\hat{k}_{ij,\alpha}\hat{k}_{ij,\beta}\delta_{\gamma\delta}+\hat{k}_{ij,\gamma}\hat{k}_{ij,\delta}\delta_{\alpha\beta}-\hat{k}_{ij,\alpha}\hat{k}_{ij,\gamma}\delta_{\beta\delta}\right.$
$\displaystyle\hskip
30.00005pt\left.-\hat{k}_{ij,\beta}\hat{k}_{ij,\gamma}\delta_{\alpha\delta}-\hat{k}_{ij,\alpha}\hat{k}_{ij,\delta}\delta_{\beta\gamma}-\hat{k}_{ij,\beta}\hat{k}_{ij,\delta}\delta_{\alpha\gamma}+\hat{k}_{ij,\alpha}\hat{k}_{ij,\beta}\hat{k}_{ij,\delta}\hat{k}_{ij,\gamma}\right),$
(S10c)

When we consider their leading terms for $h_{ij}$, these coefficients for the
monodisperse particles are reduced to [3, 4]. Note that these expressions are
dimensional. We should be careful when we compare these results with those
given in Refs. [3, 4, 5].

$\displaystyle X^{A(11)}_{ij}$
$\displaystyle=-X^{A(12)}_{ij}=-X^{A(21)}_{ij}=X^{A(22)}_{ij}=3\pi
d_{\mathrm{H}}\frac{1}{4}\frac{1}{h_{ij}},$ (S11a) $\displaystyle
Y^{A(11)}_{ij}$
$\displaystyle=-Y^{A(12)}_{ij}=-Y^{A(21)}_{ij}=Y^{A(22)}_{ij}=3\pi
d_{\mathrm{H}}\frac{1}{6}\ln\frac{1}{h_{ij}},$ (S11b) $\displaystyle
X^{G(11)}_{ij}$
$\displaystyle=-X^{G(12)}_{ij}=X^{G(21)}_{ij}=-X^{G(22)}_{ij}=\frac{3}{8}\pi
d_{\mathrm{H}}^{2}\frac{1}{h_{ij}},$ (S11c) $\displaystyle Y^{G(11)}_{ij}$
$\displaystyle=-Y^{G(12)}_{ij}=Y^{G(21)}_{ij}=-Y^{G(22)}_{ij}=\frac{1}{8}\pi
d_{\mathrm{H}}^{2}\ln\frac{1}{h_{ij}},$ (S11d) $\displaystyle X^{M(11)}_{ij}$
$\displaystyle=X^{M(12)}_{ij}=X^{M(21)}_{ij}=X^{M(22)}_{ij}=\frac{1}{8}\pi
d_{\mathrm{H}}^{3}\frac{1}{h_{ij}},$ (S11e) $\displaystyle Y^{M(11)}_{ij}$
$\displaystyle=Y^{M(22)}_{ij}=\frac{1}{10}\pi
d_{\mathrm{H}}^{3}\ln\frac{1}{h_{ij}},\quad
Y^{M(12)}_{ij}=Y^{M(21)}_{ij}=\frac{1}{40}\pi
d_{\mathrm{H}}^{3}\ln\frac{1}{h_{ij}},$ (S11f) $\displaystyle Z^{M(11)}_{ij}$
$\displaystyle=Z^{M(22)}_{ij}\simeq\frac{5}{6}\pi d_{\mathrm{H}}^{3},\quad
Z^{M(12)}_{ij}=Z^{M(21)}_{ij}\simeq-\frac{3}{16}\pi d_{\mathrm{H}}^{3},$
(S11g)

respectively.

Let us introduce the peculiar velocity $\bm{V}_{i}$ as

$\bm{V}_{i}\equiv\bm{v}_{i}-\overleftrightarrow{E}\bm{x}_{i}.$ (S12)

Then, the force acting on $i$–th particle from $j$–th particle is given by

$\bm{F}^{\mathrm{H}}_{ij}=-\eta_{0}\left(\overleftrightarrow{A^{(11)}_{ij}}\bm{V}_{i}+\overleftrightarrow{A^{(12)}_{ij}}\bm{V}_{j}\right)+\eta_{0}\left(\overleftrightarrow{\widetilde{G}^{(11)}_{ij}}+\overleftrightarrow{\widetilde{G}^{(12)}_{ij}}\right)\overleftrightarrow{E},$
(S13)

or equivalently,

$\displaystyle F^{\mathrm{H}}_{ij,\alpha}$
$\displaystyle=-\eta_{0}\left(A^{(11)}_{ij,\alpha\beta}V_{i,\beta}+A^{(12)}_{ij,\alpha\beta}V_{j,\beta}\right)+\eta_{0}\left(\widetilde{G}^{(11)}_{ij,\alpha\beta\gamma}+\widetilde{G}^{(12)}_{ij,\alpha\beta\gamma}\right)E_{\beta\gamma}$
$\displaystyle=-\eta_{0}\left(A^{(11)}_{ij,\alpha\beta}V_{i,\beta}+A^{(12)}_{ij,\alpha\beta}V_{j,\beta}\right)+\eta_{0}\left(G^{(11)}_{ij,\beta\gamma\alpha}+G^{(21)}_{ij,\beta\gamma\alpha}\right)E_{\beta\gamma}$
$\displaystyle=-\eta_{0}\left(A^{(11)}_{ij,\alpha\beta}V_{i,\beta}+A^{(12)}_{ij,\alpha\beta}V_{j,\beta}\right)+2\eta_{0}\dot{\gamma}G^{(11)}_{ij,xy\alpha},$
(S14)

where we have used the symmetric property of the tensor $G^{(k\ell)}_{ij}$
(see Eq. (S8b)). We can rewrite $F^{\mathrm{H}}_{ij,\alpha}$ furthermore. With
the aid of the identity

$\displaystyle
A^{(11)}_{ij,\alpha\beta}V_{i,\beta}+A^{(12)}_{ij,\alpha\beta}V_{j,\beta}$
$\displaystyle=A^{(11)}_{ij,\alpha\beta}V_{i,\beta}-A^{(11)}_{ij,\alpha\beta}V_{j,\beta}$
$\displaystyle=A^{(11)}_{ij,\alpha\beta}V_{ij,\beta},$ (S15)

with $\bm{V}_{ij}\equiv\bm{V}_{i}-\bm{V}_{j}$, we can rewrite Eq. (S14) as

$F^{\mathrm{H}}_{ij,\alpha}=-\eta_{0}A^{(11)}_{ij,\alpha\beta}V_{ij,\beta}+2\eta_{0}\dot{\gamma}G^{(11)}_{ij,xy\alpha}.$
(S16)

Similarly, $F_{ji,\alpha}$ is also written as

$\displaystyle F^{\mathrm{H}}_{ji,\alpha}$
$\displaystyle=\eta_{0}A^{(22)}_{ij,\alpha\beta}V_{ij,\beta}-2\eta_{0}\dot{\gamma}G^{(11)}_{ij,xy\alpha}$
$\displaystyle=\eta_{0}A^{(11)}_{ij,\alpha\beta}V_{ij,\beta}-2\eta_{0}\dot{\gamma}G^{(11)}_{ij,xy\alpha}$
$\displaystyle=-F^{\mathrm{H}}_{ij,\alpha}.$ (S17)

Next, let us calculate the hydrodynamic stress tensor
$\overleftrightarrow{\overline{\sigma}_{ij}^{\mathrm{H}}}$:

$\displaystyle\overline{\sigma}^{\mathrm{H}}_{ij,\alpha\beta}$
$\displaystyle=-\eta_{0}\left(G^{(11)}_{ij,\alpha\beta\gamma}V_{i,\gamma}+G^{(12)}_{ij,\alpha\beta\gamma}V_{j,\gamma}\right)+\eta_{0}\left(M^{(11)}_{ij,\alpha\beta\gamma\delta}+M^{(12)}_{ij,\alpha\beta\gamma\delta}\right)E_{\gamma\delta}$
$\displaystyle=-\eta_{0}\left(G^{(11)}_{ij,\alpha\beta\gamma}V_{i,\gamma}+G^{(12)}_{ij,\alpha\beta\gamma}V_{j,\gamma}\right)+\eta_{0}\dot{\gamma}\left(M^{(11)}_{ij,\alpha\beta
xy}+M^{(12)}_{ij,\alpha\beta xy}\right).$ (S18)

Using the similar procedure as that in $F_{i,\alpha}$, we obtain the
followings:

$\displaystyle
G^{(11)}_{ij,\alpha\beta\gamma}V_{i,\gamma}+G^{(12)}_{ij,\alpha\beta\gamma}V_{j,\gamma}$
$\displaystyle=G^{(11)}_{ij,\alpha\beta\gamma}V_{i,\gamma}-G^{(11)}_{ij,\alpha\beta\gamma}V_{j,\gamma}$
$\displaystyle=G^{(11)}_{ij,\alpha\beta\gamma}V_{ij,\gamma}.$ (S19)

and

$\displaystyle M^{(11)}_{ij,\alpha\beta xy}+M^{(12)}_{ij,\alpha\beta xy}$
$\displaystyle=\left(X^{M}_{11}+X^{M}_{12}\right)\hat{k}^{(0)}_{\alpha\beta
xy}+\left(Y^{M}_{11}+Y^{M}_{12}\right)\hat{k}^{(1)}_{\alpha\beta
xy}+\left(Z^{M}_{11}+Z^{M}_{12}\right)\hat{k}^{(2)}_{\alpha\beta xy}$
$\displaystyle=M^{(1)}_{ij,\alpha\beta},$ (S20)

where we have introduced

$M^{(1)}_{ij,\alpha\beta}\equiv X^{M(1)}_{ij}\hat{k}^{(0)}_{ij,\alpha\beta
xy}+Y^{M(1)}_{ij}\hat{k}^{(1)}_{ij,\alpha\beta
xy}+Z^{M(1)}_{ij}\hat{k}^{(2)}_{ij,\alpha\beta xy},$ (S21)

with

$\displaystyle X^{M(1)}_{ij}$ $\displaystyle\equiv
X^{M(11)}_{ij}+X^{M(12)}_{ij}=\frac{1}{4}\pi
d_{\mathrm{H}}^{3}\frac{1}{h_{ij}},$ (S22a) $\displaystyle Y^{M(1)}_{ij}$
$\displaystyle\equiv Y^{M(11)}_{ij}+Y^{M(12)}_{ij}=\frac{1}{8}\pi
d_{\mathrm{H}}^{3}\ln\frac{1}{h_{ij}},$ (S22b) $\displaystyle Z^{M(1)}_{ij}$
$\displaystyle\equiv Z^{M(11)}_{ij}+Z^{M(12)}_{ij}=\frac{31}{48}\pi
d_{\mathrm{H}}^{3}.$ (S22c)

Then, we can evaluate
$\overleftrightarrow{\overline{\sigma}_{ij}^{\mathrm{H}}}$ as

$\overline{\sigma}_{ij,\alpha\beta}^{\mathrm{H}}=-\eta_{0}G^{(11)}_{ij,\alpha\beta\gamma}V_{ij,\gamma}+\eta_{0}\dot{\gamma}M^{(1)}_{ij,\alpha\beta}.$
(S23)

We can also calculate $\overleftrightarrow{\sigma_{ji}^{\mathrm{H}}}$, whose
$\alpha\beta$ component is given by

$\overline{\sigma}_{ji,\alpha\beta}^{\mathrm{H}}=-\eta_{0}\left(G^{(21)}_{ij,\alpha\beta\gamma}V_{i,\gamma}+G^{(22)}_{ij,\alpha\beta\gamma}V_{j,\gamma}\right)+\eta_{0}\dot{\gamma}\left(M^{(21)}_{ij,\alpha\beta
xy}+M^{(22)}_{ij,\alpha\beta xy}\right),$ (S24)

where

$G^{(21)}_{ij,\alpha\beta\gamma}V_{i,\gamma}+G^{(22)}_{ij,\alpha\beta\gamma}V_{j,\gamma}=-G^{(22)}_{ij,\alpha\beta\gamma}V_{ij,\gamma},$
(S25)

and

$M^{(2)}_{ij,\alpha\beta}\equiv X^{M(2)}_{ij}\widehat{k}^{(0)}_{ij,\alpha\beta
xy}+Y^{M(2)}_{ij}\widehat{k}^{(1)}_{ij,\alpha\beta
xy}+Z^{M(2)}_{ij}\widehat{k}^{(2)}_{ij,\alpha\beta xy},$ (S26)

with

$\displaystyle X^{M(2)}_{ij}$ $\displaystyle\equiv
X^{M(21)}_{ij}+X^{M(22)}_{ij}=X^{M(1)}_{ij},$ (S27a) $\displaystyle
Y^{M(2)}_{ij}$ $\displaystyle\equiv
Y^{M(21)}_{ij}+Y^{M(22)}_{ij}=Y^{M(1)}_{ij},$ (S27b) $\displaystyle
Z^{M(2)}_{ij}$ $\displaystyle\equiv
Z^{M(21)}_{ij}+Z^{M(22)}_{ij}=Z^{M(1)}_{ij}.$ (S27c)

There, the hydrodynamic stress due to the hydrodynamic interaction between
$i$–th and $j$–th particles is given by

$\displaystyle\overline{\sigma}^{\mathrm{H}}_{ij,\alpha\beta}$
$\displaystyle=\overline{\sigma}_{ij,\alpha\beta}+\overline{\sigma}_{ji,\alpha\beta}$
$\displaystyle=-\eta_{0}\left(G^{(11)}_{ij,\alpha\beta\gamma}-G^{(22)}_{ij,\alpha\beta\gamma}\right)V_{ij,\gamma}+\eta_{0}\dot{\gamma}\left(M^{(1)}_{ij,\alpha\beta}+M^{(2)}_{ij,\alpha\beta}\right)$
$\displaystyle=-2\eta_{0}G^{(11)}_{ij,\alpha\beta\gamma}V_{ij,\gamma}+2\eta_{0}\dot{\gamma}M^{(1)}_{ij,\alpha\beta}.$
(S28)

In the LF-DEM simulation, we introduce the outer-cutoff of the interaction
range $r_{\mathrm{max}}$ which is given by
$r_{\mathrm{max}}/d_{\mathrm{H}}=1.5$.

##  III Convergence of the results with respect to the strain increment of the
noise

In this section, we verify the convergence of the results regardless of the
update interval of the noise term. Although the noise term depends on the
configuration of particles when the hydrodynamic interaction exists, the
generation of the noise is time-consuming as compared to the case when the
resistance matrix is diagonal. Thus, we need to look for the optimal update
interval of the noise to implement the simulation code.

Let us check how the choice of the update step affects the results. Since we
are interested in the behavior under a simple shear with $\dot{\gamma}$, the
dimensionless time scale in the simulation is the strain
$\gamma\equiv\dot{\gamma}t$. Therefore, the time increment in the simulation
is expressed as the step strain $\Delta\gamma$ in the simulation.

Figure S1 shows that the evolution of the temperature is insensitive to the
choice of the interval, at least, for $\Delta\gamma\lesssim 10^{-1}$.
Therefore, in this paper, we adopt the updated interval of the noise term
$\Delta\gamma=10^{-1}$.

Figure S1: Plot of $\theta$ against the strain $\gamma$ for the step strain
$\Delta\gamma=10^{-1}$ (solid line) and $10^{-2}$ (dashed line) under the
conditions of $\varphi=0.30$, $\varepsilon^{*}=10^{4}$, and
$\xi_{\mathrm{env}}=1.0$.

## References

  * [1] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Oxford Univ. Press, Oxford, 1987).
  * [2] J. Happel and H. Brenner, Low Reynolds number hydrodynamics —with special applications to particulate media— (Noordhoff International Publishing, Leiden, 1973).
  * [3] D. J. Jeffrey and Y. Onishi, Calculation of the resistance and mobility functions for two unequal rigid spheres in low-Reynolds-number flow, J. Fluid Mech. 139, 261 (1984).
  * [4] D. J. Jeffrey, The calculation of the low Reynolds number resistance functions for two unequal spheres, Phys. Fluids A 4, 16 (1992).
  * [5] S. Kim and S. J. Karrila, Microhydrodynamics: Principles and Selected Applications (Dover, New York, 2005).
  * [6] K. Ichiki, A. E. Kobryn, and A. Kovalenko, Resistance functions for two unequal spheres in linear flow at low Reynolds number with the Navier slip boundary condition, arXiv:1302.0461.