#  Dynamic susceptibilities in dense soft athermal spheres under a finite-rate
shear

Norihiro Oyama<EMAIL_ADDRESS>Graduate School of Arts and
Sciences, The University of Tokyo, Tokyo 153-8902, Japan  Mathematics for
Advanced Materials-OIL, AIST, Sendai 980-8577, Japan  Takeshi Kawasaki
<EMAIL_ADDRESS>Department of Physics, Nagoya University,
Nagoya 464-8602, Japan  Kuniyasu Saitoh  Department of Physics, Faculty of
Science, Kyoto Sangyo University, Motoyama, Kamigamo, Kita-ku, Kyoto 603-8555,
Japan

###### Abstract

The mechanical responses of dense packings of soft athermal spheres under a
finite-rate shear are studied by means of molecular dynamics simulations. We
investigate the volume fraction and shear rate dependence of the fluctuations
in the shear stress and the interparticle contact number. In particular, we
quantify them by defining the susceptibility as the ratio of the global to
local fluctuations. The obtained susceptibilities form ridges on the volume
fraction-shear rate plane, which are reminiscent of the Widom lines around the
critical point in an equilibrium phase transition.

##  I Introduction

Soft condensed matters comprising bubbles, emulsions, or powder particles are
generally referred to as “soft athermal particle systems”. Soft athermal
particles are characterized by their (quasi-) elastic interactions, and
thermal motion is negligible since they are large in size. When their density
increases quasistatically, a transition from the liquid state, where the
stress is zero, to the amorphous solid state, where the stress is finite,
occurs. This transition is called the jamming transition Liu and Nagel (1998);
O’Hern  _et al._ (2003). In the vicinity of the jamming transition point,
various physical quantities, namely, the stress, the interparticle contact
number, and the viscosity, behave critically Durian (1995); van Hecke (2010);
Kawasaki  _et al._ (2015); Olsson (2019); Ikeda  _et al._ (2020); Saitoh  _et
al._ (2020). The jamming transition is similar to the glass transition
observed in thermal particle systems such as atomic, molecular, and colloidal
systems; recently, however, they have been revealed to be distinct Ikeda  _et
al._ (2012, 2013).

The rheology of athermal particles with shear flow also exhibits critical
behaviors caused by the jamming transition. In particular, a scaling function
for the flow curve regarding the volume fraction and shear rate has been
proposed Olsson and Teitel (2007), and the validity of the scaling has been
widely discussed to date Kawasaki  _et al._ (2015); Vågberg  _et al._ (2016);
Bonn  _et al._ (2017); Saitoh and Kawasaki (2020). First, the jamming
transition can be strictly defined in the athermal quasistatic limit; thus,
under a finite-rate shear, the existence of a jamming transition is not
obvious. Most conventional jamming transition studies are concerned with the
criticality of macroscopic mean quantities, whereas with the finite-rate
shear, physical quantities such as the shear stress continuously increase with
increasing volume fraction, and no remarkable singularity is observed
Heussinger and Barrat (2009); Vågberg  _et al._ (2014, 2016); Vescovi and
Luding (2016); Nagasawa  _et al._ (2019). In the statistical mechanics of
thermal equilibrium systems, a naive phase transition picture is often
captured by the fluctuation of physical quantities. In previous studies on the
jamming transition, little discussion on the fluctuation has been made,
although it is potentially significant. Accordingly, this work focuses on the
fluctuation of the physical quantities and clarifies the jamming transition
behavior under a finite shear rate.

In this work, we investigate the stress response of soft athermal particles
using molecular dynamics simulations with a finite-rate shear flow. We measure
the volume fraction dependence of the shear stress under a constant shear
rate, and then, near the jamming transition point, which is characterized by
the athermal quasistatic (AQS) limit, we find that the fluctuation of the
stress exhibits a peak. We also find that the peak height diverges and the
peak position converges to the jamming transition point when we decrease the
shear rate towards the AQS limit, which is reminiscent of the Widom line near
the critical point in an equilibrium phase transition. Despite this
similarity, the mechanism of these fluctuations in dense athermal particles is
still not apparent due to their strong nonequilibriumness. Hence, to clarify
the mechanism, we investigate the time evolution of the stress when the stress
fluctuation is enhanced, and we reveal that under a wide range of finite
rates, the system transiently acquires rigidity intermittently. We furthermore
obtain the Widom line from the contact number fluctuations, which converge to
the jamming transition point in the AQS limit, yet its trace is not identical
to that of the stress fluctuation. These findings deepen our understanding of
the jamming transition under a finite-rate shear and provide us with
extensible knowledge for various phase transition phenomena under an external
field.

This paper is constructed as follows. First, we introduce the numerical
simulation method. Next, we discuss the average shear stress and its
fluctuation. Then, we examine the stress-strain curve and contact number
fluctuations. Afterward, we draw the Widom lines obtained from the stress and
contact number fluctuations. Finally, we summarize the results and give our
perspectives.

##  II Numerical methods

We employ molecular dynamics (MD) simulations of soft athermal particles in
three dimensions. To avoid crystallization of the system, we prepare a 50:50
binary mixture of $N$ particles, where different kinds of particles have the
same mass $m$ and different diameters, $d$ and $1.4d$ O’Hern  _et al._ (2003).
The force between the particles, $i$ and $j$, in contact is modeled by a
“linear spring-dashpot” Luding (2005), i.e.,
$\bm{f}_{ij}=(k\xi_{ij}-\eta\dot{\xi}_{ij})\bm{n}_{ij}$, with the stiffness
$k$ and viscosity coefficient $\eta$. The force is parallel to the normal unit
vector $\bm{n}_{ij}=\bm{r}_{ij}/|\bm{r}_{ij}|$, where
$\bm{r}_{ij}\equiv\bm{r}_{i}-\bm{r}_{j}$, with the particle positions,
$\bm{r}_{i}$ and $\bm{r}_{j}$ denoting the relative positions. In addition,
$\xi_{ij}=R_{i}+R_{j}-|\bm{r}_{ij}|>0$ is the overlap between the particles,
and $\dot{\xi}_{ij}$ is its time derivative, where $R_{i}$ ($R_{j}$) is the
radius of particle $i$ ($j$). The stiffness and viscosity coefficient
determine the time scale as $t_{0}\equiv\eta/k$ and are adjusted such that the
_normal restitution coefficient_ of the particles is exactly zero, i.e.,
$e=\text{exp}(-\pi/\sqrt{2mk/\eta^{2}-1})=0$ Luding (2005).

We randomly distribute the $N$ particles in an $L\times L\times L$ cubic
periodic box and relax the system to a mechanically stable state Bitzek  _et
al._ (2006). Then, we apply simple shear deformations to the system under the
Lees-Edwards boundary conditions Lees and Edwards (1972). In each time step,
we apply affine deformation to the system by replacing every particle position
$(x_{i},y_{i},z_{i})$ with $\bm{r}_{i}=(x_{i}+\Delta\gamma y_{i},y_{i},z_{i})$
($i=1,\dots,N$) and then numerically integrate the equations of motion,
$m\ddot{\bm{r}}_{i}=\sum_{j}\bm{f}_{ij}$, with a small time increment $\Delta
t$ Saitoh and Mizuno (2016a, 2017). Here, $\Delta\gamma$ is the strain
increment; hence, the shear rate is defined as
$\dot{\gamma}\equiv\Delta\gamma/\Delta t$.

In our MD simulations, we control the volume fraction of the particles
$\varphi$ and the shear rate $\dot{\gamma}$. To control the shear rate, we
change both $\Delta\gamma$ and $\Delta t$ within the constraints
$\Delta\gamma\leq 10^{-6}$ and $\Delta t\leq 0.1t_{0}$. In addition, we
measure the mechanical responses of the system to simple shear deformations by
the shear stress

$\sigma=-\frac{1}{L^{3}}\sum_{i,j}f_{ijx}^{\mathrm{el}}r_{ijy}~{}.$ (1)

Here, $f_{ijx}^{\mathrm{el}}=k\xi_{ij}n_{ijx}$ is the $x$-component of the
elastic force, and $r_{ijy}$ is the $y$-component of the relative position
$\bm{r}_{ij}$ between the particles $i$ and $j$, which are in contact. For
each $\varphi$ and $\dot{\gamma}$, we compute the mean value
$\langle\sigma\rangle$ and fluctuations of the shear stress in a steady state,
where the applied strain is in the range $1<\gamma<5$. We also take ensemble
averages of $\langle\sigma\rangle$ and $\chi_{\sigma}$ (the definitions of
which are given in Sec. III.2) over at least $20$ different initial
configurations.

##  III Results

Figure 1:  Simulation results of (a) the average stress $\langle\sigma\rangle$
and (b) the susceptibility $\chi_{\sigma}$ as functions of the volume fraction
$\varphi$. Different markers are used to distinguish the different shear rate
$\dot{\gamma}$ values, as shown in the legend in (b). The dotted lines depict
the location of the jamming point $\varphi_{\rm J}$.

###  III.1 Average shear stress

We first present the dependence of the average shear stress
$\langle\sigma\rangle$ on the volume fraction $\varphi$ and the shear rate
$\dot{\gamma}$ in Fig. 1a. Specifically, the values of $\langle\sigma\rangle$
under different combinations of the parameters as functions of $\varphi$ are
shown.

In the low $\varphi$ regime, $\langle\sigma\rangle$ plateaus for all
$\dot{\gamma}$. We can also tell that $\sigma_{\rm low}$ scales linearly with
$\dot{\gamma}$. This Newtonian-like shear rate dependence is considered the
consequence of the effective overdamped dynamics due to the zero restitution
coefficient.

In the high $\varphi$ regime, $\langle\sigma\rangle$ increases with increasing
$\varphi$. In particular, when $\varphi$ is high enough and the system
exhibits a clear yielding behavior, the $\dot{\gamma}$ dependence of
$\langle\sigma\rangle$ follows the famous Herschel-Bulkley law Herschel and
Bulkley (1926): $\langle\sigma\rangle\sim\sigma_{\rm Y}+\dot{\gamma}^{n}$ (see
Appendix A for the flow curve when $\varphi=0.65$ Maeda and Takeuchi (1978);
Heussinger and Barrat (2009); Lin  _et al._ (2014); Saitoh  _et al._ (2019);
Oyama  _et al._ (2020a)).

Between these two qualitatively different volume fraction regimes, we observe
a steep growth in $\langle\sigma\rangle$. As intuitively expected, this sharp
increase in $\langle\sigma\rangle$ is observed in the vicinity of the jamming
point ($\varphi_{\rm J}\approx 0.6461$; see Appendix B for the determination
of $\varphi_{\rm J}$ under shear Kawasaki  _et al._ (2015); Kawasaki and
Miyazaki (2020)). However, the stress growth is most prominent at a volume
fraction that is clearly smaller than $\varphi_{\rm J}$ at finite
$\dot{\gamma}$. Furthermore, as $\dot{\gamma}$ increases, the growth becomes
less steep, and the onset volume fraction of the stress growth shifts towards
the low $\varphi$ side.

###  III.2 Susceptibility of the shear stress

We next focus on the fluctuation of the shear stress. In particular, we
quantify the enhancement of the collectivity in the fluctuations that
accompanies the rapid increase in $\langle\sigma\rangle$ by the susceptibility
$\chi_{\sigma}$, defined as:

$\displaystyle\chi_{\sigma}\equiv
N(\langle\sigma^{2}\rangle-\langle\sigma\rangle^{2})/(\langle\sigma_{\rm
local}^{2}\rangle-\langle\sigma_{\rm local}\rangle^{2}),$ (2)

where $\langle\sigma_{\rm local}\rangle$ is the time- and particle-averaged
value of the particle-based local stress
$\sigma_{i}\equiv-\frac{N}{2V}\sum_{j\in{\rm contact}}f_{ijx}^{\rm
el}(t)r_{ijy}(t)$ and $\sum_{j\in{\rm contact}}$ is the sum over the neighbors
($\langle\sigma_{\rm local}^{2}\rangle$ is the corresponding second-order
moment)111 The only difference between the definitions of
$\langle\sigma^{2}\rangle$ and $\langle\sigma_{\rm local}^{2}\rangle$ is the
order in which the averages are taken over particles and time.. With this
definition, the average of $\sigma_{i}$ over the particles is identical to the
macroscopic value $\sigma$, $\sigma=\frac{1}{N}\sum_{i}^{N}\sigma_{i}$. This
susceptibility $\chi_{\sigma}$ quantifies the degree of collectivity in the
stress fluctuations: $\chi_{\sigma}$ is expected to diverge with increasing
system size $N$ when the whole system behaves collectively, as in a system
located near a critical point. In Fig. 1b, we plot the measurement results of
$\chi_{\sigma}$ as a function of the volume fraction $\varphi$.

In the low $\varphi$ regime, $\chi_{\sigma}$ increases with increasing
$\dot{\gamma}$. However, interestingly, for low rates ($\dot{\gamma}\leq
10^{-5}$), $\chi_{\sigma}$ hardly depends on $\dot{\gamma}$. This behavior is
in contrast to that of $\langle\sigma\rangle$, which depends linearly on
$\dot{\gamma}$ for all shear rates $\dot{\gamma}$ in the low $\varphi$ regime.
Regarding the volume fraction dependence in this regime, $\chi_{\sigma}$ grows
weakly with increasing volume fraction.

In the high $\varphi$ regime, the opposite trend is observed: $\chi_{\sigma}$
becomes smaller when either $\dot{\gamma}$ or $\varphi$ increases. Still, the
$\dot{\gamma}$ dependence disappears for low values of $\dot{\gamma}$ (in this
case, $\dot{\gamma}\leq 10^{-6}$), in accordance with the behavior in the low
$\varphi$ regime.

At an intermediate value of $\varphi$ between these two regimes,
$\chi_{\sigma}$ exhibits a clear peak. As the shear rate $\dot{\gamma}$
increases, the height of the peak decreases, and the position shifts towards
the low $\varphi$ direction. Note that if we further increase the shear rate
to $\dot{\gamma}=10^{-2}$, we no longer observe a peak, at least in the range
of the volume fraction that we have investigated, i.e., $0.62\leq\varphi\leq
0.65$. In accordance with the convergence of $\chi_{\sigma}$ both in the high
and low $\varphi$ regimes, the height and position of the peak become almost
constant for $\dot{\gamma}\leq 10^{-6}$. This total convergence of the
susceptibility $\chi_{\sigma}$ in the low rate regime over all values of
$\varphi$ suggests that the length scale that governs the stress fluctuation
spans the whole system in this regime. We discuss the possible candidates for
this length scale in Sec. IV.2, although we leave the precise identification
for a future study.

Hereafter, we call the height and position of this peak $\chi_{\sigma}^{\rm
max}(\dot{\gamma})$ and $\varphi_{\chi_{\sigma}^{\rm max}}(\dot{\gamma})$,
respectively (we omit the explicit notation for the $\dot{\gamma}$ dependence
below).

Figure 2:  Stress-strain curves for various combinations of the shear rate
$\dot{\gamma}$ and the volume fraction $\varphi$. The vertical axis follows a
log-scale. Results for (a) $\dot{\gamma}=10^{-6}$, (b) $\dot{\gamma}=10^{-5}$,
and (c) $\dot{\gamma}=10^{-4}$ are shown. In all panels, the results for three
volume fractions are compared: $\varphi=0.62$ (black), $\varphi=0.65$ (light
gray), and $\varphi=\varphi_{\chi_{\sigma}^{\rm max}}$ (dark gray).

###  III.3 Stress-strain curves

To further obtain an intuitive understanding of the parameter dependence of
the susceptibility $\chi_{\sigma}$, we plot typical stress-strain curves for
the systems under various combinations of the volume fraction $\varphi$ and
the shear rate $\dot{\gamma}$ ($\dot{\gamma}=10^{-6},10^{-5}$ and $10^{-4}$)
in Fig. 2. For the whole parameter space investigated here, the average stress
$\langle\sigma\rangle$ becomes larger with both increasing $\varphi$ and
increasing $\dot{\gamma}$, as presented in Fig. 1a. However, the dependence on
$\dot{\gamma}$ changes significantly depending on $\varphi$: while the order
of $\langle\sigma\rangle$ remains the same regardless of the value of
$\dot{\gamma}$ at a high volume fraction
($\varphi=0.65>\varphi_{\chi_{\sigma}^{\rm max}}$, Fig. 2 light gray
curves)222Note that since $\langle\sigma\rangle$ obeys the Herschel-Bulkley
law at this high volume fraction ($\varphi=0.65$), the order of the stress
becomes larger if we apply a much faster shear., it scales linearly with
$\dot{\gamma}$ at a low volume fraction
($\varphi=0.62<\varphi_{\chi_{\sigma}^{\rm max}}$, Fig. 2 black curves).
However, the shapes of the stress-strain curves in these different regimes are
similar in that the fluctuations are suppressed.

By contrast, the shape of the stress-strain curves dramatically changes in the
vicinity of $\varphi_{\chi_{\sigma}^{\rm max}}$ under a slow shear rate
($\dot{\gamma}=10^{-6}$; Fig. 7a dark gray curve): we observe spiky peaks,
with the height of the baseline being on the order of the stress at low
$\varphi$ (see Appendix C for normal plots of the stress-strain curves where
the spiky shapes are more appreciable). The heights of the spikes are larger
than the baseline by at most two orders of magnitude and barely reach the
curve for $\varphi=0.65$. Importantly, the probability distribution of
$\sigma$, $P(\sigma)$, exhibits a power-law-like shape for
$\varphi=\varphi_{\chi_{\sigma}^{\rm max}}$ and
$\dot{\gamma}=10^{-6},10^{-5}$, indicating that this susceptibility peak
reflects the criticality expected for $\dot{\gamma}\to 0$ (see Appendix C). As
the shear rate increases, the spikes become less sharp and less frequent
($\dot{\gamma}=10^{-5}$; Fig. 7b dark gray curve), and finally, the whole
stress-strain curve becomes almost detached from that for a low $\varphi$ at
$\dot{\gamma}=10^{-4}$ (Fig. 7c). Since the magnitudes of the stress at the
baseline and the peak top are comparable to those for low and high volume
fractions respectively, we consider that these spikes are formed because the
system goes back and forth between fluid-like low-stress states and solid-like
large-stress states. That is, the whole system collectively changes its
“state” during the time evolution, as indicated by the susceptibility peak. We
mention that similar repetitive transitions between fluid-like and solid-like
states have also been observed under the AQS shear ($\dot{\gamma}=0$)
Heussinger and Barrat (2009). Notably, under a high shear rate
($\dot{\gamma}=10^{-4}$), $P(\sigma)$ exhibits a clear unimodal shape without
power-law tails at either end (see Appendix C). This observation suggests that
the increase in $\sigma$ becomes more similar to a cross-over rather than a
phase transition because of the effect of the strong external field (see Sec.
IV.1 for the qualitative similarity between our system and the conventional
critical phenomena).

##  IV Discussion

In this section, we discuss the similarity between our system and the
conventional critical phenomenon: the ferromagnetic transition in the Ising
model under an external field. Based on this analogy, we can tell that the
shear stress $\sigma$ can be viewed as a natural “conjugate” variable to the
strength of the external field (namely, the shear rate $\dot{\gamma}$).
However, $\sigma$ changes its value by orders of magnitude depending on
$\dot{\gamma}$ even in the “disordered”, low-stress phase. In this sense, it
is qualitatively different from conventional standard order parameters that
are normalized to be between zero and one in most cases. Therefore, we further
conduct the same analysis for an alternative candidate for an order parameter,
i.e., the interparticle contact number $z$.

###  IV.1 Correspondence to conventional criticality in equilibrium systems

To further explore the parameter dependence of the shear stress $\sigma$ and
its fluctuations, we rely on an analogy with a well-understood phase
transition. Here, in particular, we discuss an analogy with one of the most
famous examples: the Ising model under a magnetic field (see Appendix D for a
brief recapitulation of the mean-field solution). As shown in Figs. 1(a,b),
the average and the susceptibility of the stress exhibit qualitative
similarities with the magnetization and the susceptibility in the Ising model
(Appendix D): the inverse temperature $\beta$, which is the control parameter
of the criticality in the Ising model, corresponds to the volume fraction
$\varphi$ in our system. Similarly, the external magnetic field $h$ and the
magnetization $m$ correspond to the shear rate $\dot{\gamma}$ and the mean
stress $\langle\sigma\rangle$, respectively. Moreover, in both systems, as the
external field ($h$ or $\dot{\gamma}$) increases, the change in the order
parameter ($m$ or $\langle\sigma\rangle$) becomes less steep, and the whole
plot shifts towards the less-ordered side. Regarding the susceptibility
($\chi$ or $\chi_{\sigma}$), we observe peaks at a value of the control
parameter ($\beta$ or $\varphi$) that is shifted from the critical point when
an external field is present. The height of these peaks decreases with
increasing external field, and the position shifts towards the small-order
side. We emphasize that the counterpart of the magnetic field in our system is
not the strain $\gamma$ but the shear rate $\dot{\gamma}$, which is the
conjugate of the stress in effective energy dissipation. Hence, the free
energy of the Ising model corresponds to the dissipation function in our
system and is consistent with the empirical knowledge that the dissipation
system takes precedence over the dynamics of the minimum energy dissipation
Unger  _et al._ (2004); Török  _et al._ (2007). In this sense, the shear
stress $\sigma$ can be viewed as a natural conjugate variable to the external
field and thus as an order parameter. However, since $\sigma$ is dependent not
only on the existence of contacts but also on the degree of overlapping of
each contact, it changes its value by orders of magnitude depending on
$\dot{\gamma}$ even in the “disordered”, dilute state. In the next section, we
instead measure the average and susceptibility of the interparticle contact
number, the values of which are expected to exhibit less $\dot{\gamma}$
dependence.

###  IV.2 Contact number

The interparticle contact number $z$ characterizes the jamming transition most
directly in terms of the microscopic structures van Hecke (2010). For the
jamming transition in quiescent systems without external fields, $z$ changes
discontinuously from zero to approximately $z_{\rm C}$ at the critical point
$\varphi_{\rm J}$, above which physical quantities such as the pressure or the
shear modulus change in a power-law manner, as in the case of the conventional
second-order phase transitions O’Hern  _et al._ (2003). According to Maxwell’s
condition, $z_{\rm C}=2d$ holds for frictionless soft athermal spheres, where
$d$ is the spatial dimension of the system. Here, we plot the average and the
susceptibility of the interparticle contact number $z$ (we do not exclude
rattlers to compute $z$) under a finite-rate shear as functions of $\varphi$
in Fig. 3. For the definition of the susceptibility $\chi_{z}$, we employ a
definition similar to Eq. 2.

The dependence of the average contact number $\langle z\rangle$ on the volume
fraction $\varphi$ is qualitatively very similar to that of the average stress
$\langle\sigma\rangle$: it is almost constant in the low $\varphi$ regime and
then shows sudden growth around $\varphi_{\rm J}$, after which the growth rate
decreases in the high $\varphi$ regime. However, the dependence on
$\dot{\gamma}$ is significantly different from that of $\langle\sigma\rangle$:
in the low $\varphi$ regime, the plateau disappears for high $\dot{\gamma}$,
and the shear rate dependence is not linear. Furthermore, the values of
$\langle z\rangle$ at the highest $\varphi$ hardly depend on $\dot{\gamma}$.

The susceptibility of the contact number $\chi_{z}$ behaves qualitatively very
similarly to that of $\chi_{\sigma}$: it exhibits a clear peak near the
jamming point $\varphi_{\rm J}$, and the peak height and position change in
the same way as $\chi_{\sigma}$ when $\dot{\gamma}$ increases. One major
difference from $\chi_{\sigma}$ is that the peak position and height of
$\chi_{z}$ obviously change even in the low rate limit $\dot{\gamma}\leq
10^{-6}$, where $\chi_{\sigma}$ becomes constant.

This qualitative difference intriguingly suggests that the characteristic
lengths that govern $\sigma$ and $z$ ($\xi_{\sigma}$ and $\xi_{z}$,
respectively) are different. Let us enumerate several candidates from previous
studies. For example, it is known that the correlation length of the deviation
from the continuum description diverges at the jamming point Ellenbroek  _et
al._ (2006, 2009); Lerner  _et al._ (2014); Mizuno  _et al._ (2016). This
length scale, often referred to as $l_{c}$, is a candidate for $\xi_{z}$. On
the other hand, the isotropic as-quenched state has recently been shown to be
qualitatively different from the sheared nonequilibrium steady state in terms
of the stability against perturbation, even in the AQS limit
($\dot{\gamma}=0$) Karmakar  _et al._ (2010); Oyama  _et al._ (2020b). This
knowledge implies that $l_{c}$ and $\xi_{z}$ can be different in nature, since
$l_{c}$ is measured in the absence of an external field ($\gamma=0$), while
$\xi_{z}$ should be measured in the steady state $\gamma>{\cal O}(1)$. As an
example of a correlation length measured in a dynamic situation, Refs. Saitoh
and Mizuno (2016a, b); Saitoh and Kawasaki (2020) reported that the
correlation length of the nonaffine velocities of particles diverges in the
limit of $\varphi\to\varphi_{\rm J}$ and $\dot{\gamma}\to 0$ in two-
dimensional packings of soft frictionless disks. However, this correlation
length has been shown to remain finite even in the same limit in three
dimensions Oyama  _et al._ (2019). Instead, in ref. Oyama  _et al._ (2019),
the authors introduced the correlation length of the vortex clusters, which
diverges in that limit. As another example of a dynamical correlation length,
the one associated to the yielding criticality is also known to diverge in the
limit of $\dot{\gamma}\to 0$ Lin  _et al._ (2014); Oyama  _et al._ (2020a).
However, this length scale can be well defined only in the high $\varphi$
regime, where the Herschel-Bulkley law is valid and cannot describe the total
convergence of $\chi_{\sigma}$ over the whole $\varphi$ regime. As discussed
here, multiple candidates exist, with the possibility that none of them is the
desired one. Although identifying the governing length scale by comparing all
these candidates is an important issue, we leave it as a future problem.

Finally, we present the ridges obtained by connecting the peaks of the
susceptibilities under different values of $\dot{\gamma}$ in Fig. 4. In this
plot, we compare the results for $\chi_{\sigma}$ and $\chi_{z}$. These ridges
can be regarded as the dissipative-system counterpart of the Widom lines by
definition. Both Widom lines seem to converge to $\varphi_{\rm J}$ in the
limit of $\dot{\gamma}\to 0$, as expected. Moreover, these two lines follow
different paths, as is the case for the conventional equilibrium systems,
e.g., the Widom lines around the liquid-gas critical point.

Figure 3:  (a) Average $\langle z\rangle$ and (b) susceptibility $\chi_{z}$ of
the interparticle contact number as functions of the volume fraction
$\varphi$. Different markers represent different shear rates, as shown in the
legend in Fig. 1(b). The dotted lines depict the location of the jamming point
$\varphi_{\rm J}$.  Figure 4:  Widom lines, or a plot of the locations of the
susceptibility peaks $\varphi_{\chi_{\alpha}^{\rm max}}$ on a
$\varphi-\dot{\gamma}$ plane, where $\alpha\in\\{\sigma,z\\}$. Different
symbols represent different definitions of the susceptibility, as shown in the
legend. The dotted line shows the location of the jamming point $\varphi_{\rm
J}\approx 0.6461$ estimated at $\dot{\gamma}=0$ (see Appendix B). Error bars
indicate the range of $\varphi$ for which the values of $\chi_{\alpha}$ are
greater than 90% of $\chi_{\alpha}^{\rm max}$.

##  V Summary and overview

In this work, we conducted MD simulations for dense packings of soft athermal
spheres under a finite-rate shear and investigated the dependence of the
statistics of the shear stress on the shear rate and the volume fraction. The
average stress changes largely in the vicinity of the jamming point; moreover,
the onset volume fraction for the stress growth becomes smaller when the shear
rate increases. Interestingly, this sudden stress growth is accompanied by the
formation of a peak of the susceptibility. To further understand this
susceptibility peak, we investigated the time evolution of the stress. We
found that the stress-strain curve exhibits spiky peaks at the volume fraction
where the susceptibility peak is observed. These peaks are formed since the
system can temporally gain solidity with the aid of the external shear, while
it is fluidic otherwise. We furthermore measured the average and
susceptibility of the interparticle contact number as an example of a
normalized order parameter in our system. The results for $\chi_{z}$ are
qualitatively consistent with those for $\chi_{\sigma}$, although the length
scales that govern these two fluctuations seem different. We furthermore
visualized the Widom lines in our system, or the ridges of the susceptibility
peaks for both the stress and contact number. As the equilibrium phase diagram
shows, two Widom lines follow different paths, although both seem converge to
a critical point in the limit $\dot{\gamma}\to 0$.

As a future direction, an investigation of whether modification of the
physical dimension Radjai and Roux (2002); Saitoh and Mizuno (2016b); Oyama
_et al._ (2019), the damping coefficient Andreotti  _et al._ (2012); Kawasaki
_et al._ (2014); Vågberg  _et al._ (2017), or the local dissipation mechanisms
(e.g., introduction of the tangential friction Otsuki and Hayakawa (2009,
2011)) leads to any qualitative changes should be carried out.

###### Acknowledgements.

We thank Atsushi Ikeda, Kota Mitsumoto, and Yusuke Hara for the fruitful
discussions. This work was financially supported by JSPS KAKENHI Grant Numbers
18H01188, 18K13464, 19K03767, 20H05157, 20H00128, 20H01868, 20J00802, and
20K14436.

##  Appendix A Flow curve for a dense system

In Fig. 5, we plot the average shear stress $\langle\sigma\rangle$ in the
system as a function of the shear rate $\dot{\gamma}$ with $\varphi=0.65$.

To roughly estimate the yield stress $\sigma_{\rm Y}$, we also conduct an AQS
simulation. In the AQS simulation, instead of integrating the equation of
motion, we minimize the potential energy of the system Maeda and Takeuchi
(1978); Heussinger and Barrat (2009); Saitoh  _et al._ (2019). We employ the
FIRE algorithm Bitzek  _et al._ (2006) and terminate the iteration when the
maximum magnitude of the force exerted on one particle meets $f_{\rm
max}<10^{-9}$. The strain increment is $\Delta\gamma=5\times 10^{-5}$
Heussinger and Barrat (2009).

Fig. 5 shows that the average stress converges to the AQS value at very slow
shear rates ($\dot{\gamma}\leq 10^{-6}$). In other words, these shear rates
can be considered as in the so-called _quasistatic regime_. We regard the
average stress under the AQS shear as the yield stress and further fit the
numerical results to the Herschel-Bulkley law,
$\langle\sigma\rangle\sim\sigma_{\rm Y}+\dot{\gamma}^{n}$. This simple
estimation provides $n\sim 0.62$, and the obtained curve captures the
numerical data very well. Note, however, that we must take into account the
finite size effects to accurately evaluate the Herschel-Bulkley parameters,
namely, the yield stress and the critical exponent Lin  _et al._ (2014); Oyama
_et al._ (2020a).

Figure 5:  The dependence of the average stress $\langle\sigma\rangle$ on the
shear rate $\dot{\gamma}$ at $\varphi=0.65$. The markers represent the results
of the MD simulation, and the dashed line is the fitting to the Herschel-
Bulkley law. The Herschel-Bulkley exponent here is $n\sim 0.62$. The dotted
line depicts the value of $\langle\sigma\rangle$ measured under the AQS shear.

##  Appendix B Determination of the jamming point

We determine the precise location of the jamming point $\varphi_{\rm J}$
following ref. Kawasaki and Miyazaki (2020), the protocol of which we briefly
explain here. We first prepare a random particle configuration in a fixed
volume system with a linear dimension $L$. Then, we minimize the total
potential energy of the system to obtain a mechanically equilibrated
configuration, with the pressure $p$ being controlled such that $p\approx 0$.
We furthermore apply shear in an AQS manner until the system reaches a steady
state ($\gamma=1$) with a strain increment $\Delta\gamma=10^{-3}$. The volume
fraction at the steady state ($\gamma>0.5$) can be used as a well-defined
jamming point Kawasaki and Miyazaki (2020). With this protocol, the value of
the jamming point does not depend on the initial configuration. Because
achieving numerically the exact mechanically equilibrated configuration with
zero pressure is almost impossible, we set the target pressure $P=10^{-5}$.
Following this protocol and averaging over the values in the steady state
($\gamma\geq 0.5$), we locate the jamming point as $\varphi_{\rm J}\approx
0.6461$ (Fig. 6). This value of $\varphi_{\rm J}$ is consistent with the one
estimated by directly fitting the diverging trend of the viscosity Kawasaki
_et al._ (2015). We stress that the data in Fig. 6 is the average over 60
samples.

Figure 6:  Average over 60 samples. Dashed line represents the estimated value
of $\varphi_{\rm J}=0.6461$ (see the text regarding how to locate it).

##  Appendix C Stress-strain curves in a normal plot

In this section, we present a normal plot of the stress-strain curves for
various combinations of parameters in Fig. 7 (the ones used for Fig. 2 in the
main text are employed). In the plots for $\varphi=\varphi_{\chi_{\sigma}^{\rm
max}}$, we observe sharp spikes, especially for a slow shear rate.

In Fig. 8, we plot the probability distribution function (PDF) of the shear
stress $\sigma$ for the same combinations of parameters. At a low volume
fraction $\varphi=0.62$, the PDF is almost Dirac’s delta function for all
shear rates (the width is very narrow). At a high volume fraction
$\varphi=0.65$, the PDF is unimodal, with a large width for all shear rates.
At $\varphi=\varphi_{\chi_{\sigma}^{\rm max}}$, however, we observe shear rate
dependence. Although the PDF exhibits a power-law-like shape for slow shear
rates ($\dot{\gamma}\leq 10^{-5}$), it becomes rather regular unimodal shape
for a high shear rate ($\dot{\gamma}=10^{-4}$).

Figure 7:  Stress-strain curves for systems with different combinations of the
shear rate $\dot{\gamma}$ and the volume fraction $\varphi$. Results for (top
row) $\varphi=0.65>\varphi_{\chi_{\sigma}^{\rm max}}$, (middle row)
$\varphi=\varphi_{\chi_{\sigma}^{\rm max}}$, and (bottom row)
$\varphi=0.62<\varphi_{\chi_{\sigma}^{\rm max}}$ are shown. From left to
right, the shear rate increases as follows: $\dot{\gamma}=10^{-6},10^{-5}$ and
$10^{-4}$.  Figure 8:  PDFs of the shear stress $\sigma$ for different shear
rates $\dot{\gamma}=10^{-6},10^{-5}$ and $10^{-4}$. Different symbols indicate
different volume fractions $\varphi$, as shown in the legend.

##  Appendix D Ising model

In this section, we recapitulate the famous self-consistent equation for the
magnetization of the Ising model under an external magnetic field, which is
derived with a mean-field approximation. Assume we have a $d$-dimensional
Ising-type spin system on a regular lattice whose Hamiltonian ${\cal H}$ is
written as:

$\displaystyle{\cal H}=-J\sum_{\langle
ij\rangle}S_{i}S_{j}-h\sum_{i=1}^{N}S_{i},$ (3)

where $S_{i}\in\\{1,-1\\}$ is the spin variable at the site $i$, $h$ stands
for the strength of the external field, and $J$ represents the coupling
constant. Then, with a mean-field approximation, we can derive a self-
consistent equation for the spin $\frac{1}{2}$ magnetization $m$ as:

$\displaystyle m=\text{tanh}(\beta Jz_{c}m+\beta h),$ (4)

where $z_{c}=2d$ is the spin coordination number. If we employ $d=3$ and
$J=1$, the critical inverse temperature $\beta_{c}$ is obtained as
$\beta_{c}=1/Jz\approx 0.167$. We plot the values of the magnetization $m$ and
the susceptibility $\chi=dm/dh$ as functions of $\beta$ for various values of
$h$ in Fig. 9.

Figure 9:  (a) Magnetization $h$ and (b) susceptibility $\chi$ obtained from
the self-consistent Eq. 4 as functions of the inverse temperature $\beta$.
Different colors represent different values of $h$, as shown in the legend.

## References

  * Liu and Nagel (1998) A. J. Liu and S. R. Nagel, Nature 396, 21 (1998).
  * O’Hern  _et al._ (2003) C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. E 68, 011306 (2003), arXiv:0304421 [cond-mat] .
  * Durian (1995) D. J. Durian, Phys. Rev. Lett. 75, 4780 (1995).
  * van Hecke (2010) M. van Hecke, J. Phys. Condens. Matter 22, 033101 (2010), arXiv:0911.1384 .
  * Kawasaki  _et al._ (2015) T. Kawasaki, D. Coslovich, A. Ikeda, and L. Berthier, Phys. Rev. E 91, 012203 (2015).
  * Olsson (2019) P. Olsson, Phys. Rev. Lett. 122, 108003 (2019).
  * Ikeda  _et al._ (2020) A. Ikeda, T. Kawasaki, L. Berthier, K. Saitoh, and T. Hatano, Phys. Rev. Lett. 124, 058001 (2020).
  * Saitoh  _et al._ (2020) K. Saitoh, T. Hatano, A. Ikeda, and B. P. Tighe, Phys. Rev. Lett. 124, 118001 (2020), arXiv:1909.10142 .
  * Ikeda  _et al._ (2012) A. Ikeda, L. Berthier, and P. Sollich, Phys. Rev. Lett. 109, 018301 (2012), arXiv:1203.0825 .
  * Ikeda  _et al._ (2013) A. Ikeda, L. Berthier, and P. Sollich, Soft Matter 9, 7669 (2013), arXiv:1302.4271 .
  * Olsson and Teitel (2007) P. Olsson and S. Teitel, Phys. Rev. Lett. 99, 178001 (2007).
  * Vågberg  _et al._ (2016) D. Vågberg, P. Olsson, and S. Teitel, Phys. Rev. E 93, 052902 (2016), arXiv:1510.03312 .
  * Bonn  _et al._ (2017) D. Bonn, M. M. Denn, L. Berthier, T. Divoux, and S. Manneville, Rev. Mod. Phys. 89, 035005 (2017).
  * Saitoh and Kawasaki (2020) K. Saitoh and T. Kawasaki, Front. Phys. 8, 1 (2020).
  * Heussinger and Barrat (2009) C. Heussinger and J.-L. Barrat, Phys. Rev. Lett. 102, 218303 (2009).
  * Vågberg  _et al._ (2014) D. Vågberg, P. Olsson, and S. Teitel, Phys. Rev. Lett. 112, 208303 (2014).
  * Vescovi and Luding (2016) D. Vescovi and S. Luding, Soft Matter 12, 8616 (2016), arXiv:1609.07414 .
  * Nagasawa  _et al._ (2019) K. Nagasawa, K. Miyazaki, and T. Kawasaki, Soft Matter 15, 7557 (2019), arXiv:1907.08383 .
  * Luding (2005) S. Luding, J. Phys. Condens. Matter 17, S2623 (2005).
  * Bitzek  _et al._ (2006) E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, and P. Gumbsch, Phys. Rev. Lett. 97, 170201 (2006).
  * Lees and Edwards (1972) A. W. Lees and S. F. Edwards, J. Phys. C Solid State Phys. 5, 1921 (1972).
  * Saitoh and Mizuno (2016a) K. Saitoh and H. Mizuno, Phys. Rev. E 94, 022908 (2016a).
  * Saitoh and Mizuno (2017) K. Saitoh and H. Mizuno, Phys. Rev. E 96, 012903 (2017).
  * Herschel and Bulkley (1926) W. H. Herschel and R. Bulkley, Kolloid-Zeitschrift 39, 291 (1926).
  * Maeda and Takeuchi (1978) K. Maeda and S. Takeuchi, Phys. Status Solidi 49, 685 (1978).
  * Lin  _et al._ (2014) J. Lin, E. Lerner, A. Rosso, and M. Wyart, Proc. Natl. Acad. Sci. 111, 14382 (2014), arXiv:1403.6735 .
  * Saitoh  _et al._ (2019) K. Saitoh, N. Oyama, F. Ogushi, and S. Luding, Soft Matter 15, 3487 (2019).
  * Oyama  _et al._ (2020a) N. Oyama, H. Mizuno, and A. Ikeda, arXiv 1, 1 (2020a), arXiv:2011.12568 .
  * Kawasaki and Miyazaki (2020) T. Kawasaki and K. Miyazaki, arXiv , 1 (2020), arXiv:2003.10716 .
  * Unger  _et al._ (2004) T. Unger, J. Török, J. Kertész, and D. E. Wolf, Phys. Rev. Lett. 92, 214301 (2004).
  * Török  _et al._ (2007) J. Török, T. Unger, J. Kertész, and D. E. Wolf, Phys. Rev. E 75, 011305 (2007), arXiv:0607162 [cond-mat] .
  * Ellenbroek  _et al._ (2006) W. G. Ellenbroek, E. Somfai, M. van Hecke, and W. van Saarloos, Phys. Rev. Lett. 97, 258001 (2006), arXiv:0604157 [cond-mat] .
  * Ellenbroek  _et al._ (2009) W. G. Ellenbroek, M. van Hecke, and W. van Saarloos, Phys. Rev. E 80, 061307 (2009), arXiv:0911.0944 .
  * Lerner  _et al._ (2014) E. Lerner, E. DeGiuli, G. Düring, and M. Wyart, Soft Matter 10, 5085 (2014), arXiv:1312.2146 .
  * Mizuno  _et al._ (2016) H. Mizuno, L. E. Silbert, and M. Sperl, Phys. Rev. Lett. 116, 068302 (2016).
  * Karmakar  _et al._ (2010) S. Karmakar, E. Lerner, and I. Procaccia, Phys. Rev. E 82, 055103 (2010), arXiv:1008.3967 .
  * Oyama  _et al._ (2020b) N. Oyama, H. Mizuno, and A. Ikeda,  (2020b), arXiv:2009.02635 .
  * Saitoh and Mizuno (2016b) K. Saitoh and H. Mizuno, Soft Matter 12, 1360 (2016b).
  * Oyama  _et al._ (2019) N. Oyama, H. Mizuno, and K. Saitoh, Phys. Rev. Lett. 122, 188004 (2019).
  * Radjai and Roux (2002) F. Radjai and S. Roux, Phys. Rev. Lett. 89, 064302 (2002).
  * Andreotti  _et al._ (2012) B. Andreotti, J.-L. Barrat, and C. Heussinger, Phys. Rev. Lett. 109, 105901 (2012).
  * Kawasaki  _et al._ (2014) T. Kawasaki, A. Ikeda, and L. Berthier, EPL (Europhysics Lett. 107, 28009 (2014), arXiv:1404.4778 .
  * Vågberg  _et al._ (2017) D. Vågberg, P. Olsson, and S. Teitel, Phys. Rev. E 95, 052903 (2017).
  * Otsuki and Hayakawa (2009) M. Otsuki and H. Hayakawa, Phys. Rev. E 80, 011308 (2009), arXiv:0904.3267 .
  * Otsuki and Hayakawa (2011) M. Otsuki and H. Hayakawa, Phys. Rev. E 83, 051301 (2011), arXiv:1006.3597 .