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

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.

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 [1,2]. In the vicinity of the jamming transition point, various physical quantities, namely, the stress, the interparticle contact number, and the viscosity, behave critically [3][4][5][6][7][8]. 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 [9,10].
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 [11], and the validity of the scaling has been widely discussed to date [5,[12][13][14]. 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 [12,[15][16][17][18]. 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 [2]. The force between the particles, i and j, in contact is modeled by a "linear spring-dashpot" [19], i.e., f ij = (kξ ij − ηξ ij )n ij , with the stiffness k and viscosity coefficient η. The force is parallel to the normal unit vector n ij = r ij /|r ij |, where r ij ≡ r i − r j , with the particle positions, r i and r j denoting the relative positions. In addition, ξ ij = R i + R j − |r ij | > 0 is the overlap between the particles, andξ 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 ≡ η/k and are adjusted such that the normal restitution coefficient of the particles is exactly zero, i.e., e = exp(−π/ 2mk/η 2 − 1) = 0 [19].
We randomly distribute the N particles in an L×L×L cubic periodic box and relax the system to a mechanically stable state [20]. Then, we apply simple shear deformations to the system under the Lees-Edwards boundary conditions [21]. In each time step, we apply affine deformation to the system by replacing every particle position (x i , y i , z i ) with r i = (x i +∆γy i , y i , z i ) (i = 1, . . . , N ) and then numerically integrate the equations of motion, mr i = j f ij , with a small time increment ∆t [22,23]. Here, ∆γ is the strain increment; hence, the shear rate is defined asγ ≡ ∆γ/∆t.
In our MD simulations, we control the volume fraction of the particles ϕ and the shear rateγ. To control the shear rate, we change both ∆γ and ∆t within the constraints ∆γ ≤ 10 −6 and ∆t ≤ 0.1t 0 . In addition, we measure the mechanical responses of the system to simple shear deformations by the shear stress Here, f el ijx = kξ ij n ijx is the x-component of the elastic force, and r ijy is the y-component of the relative position r ij between the particles i and j, which are in contact. For each ϕ andγ, we compute the mean value σ and fluctuations of the shear stress in a steady state, where the applied strain is in the range 1 < γ < 5. We also take ensemble averages of σ and χ σ (the definitions of which are given in Sec. III B) over at least 20 different initial configurations.

A. Average shear stress
We first present the dependence of the average shear stress σ on the volume fraction ϕ and the shear ratė γ in Fig. 1a. Specifically, the values of σ under different combinations of the parameters as functions of ϕ are shown. In the low ϕ regime, σ plateaus for allγ. We can also tell that σ low scales linearly withγ. This Newtonianlike shear rate dependence is considered the consequence of the effective overdamped dynamics due to the zero restitution coefficient.
Between these two qualitatively different volume fraction regimes, we observe a steep growth in σ . As intuitively expected, this sharp increase in σ is observed in the vicinity of the jamming point (ϕ J ≈ 0.6461; see Appendix B for the determination of ϕ J under shear [5,29]). However, the stress growth is most prominent at a volume fraction that is clearly smaller than ϕ J at finiteγ. Furthermore, asγ increases, the growth becomes less steep, and the onset volume fraction of the stress growth shifts towards the low ϕ side.

B. 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 σ by the susceptibility χ σ , defined as: where σ local is the time-and particle-averaged value of the particle-based local stress and j∈contact is the sum over the neighbors ( σ 2 local is the corresponding second-order moment) 1 . With this definition, the average of σ i over the particles is identical to the macroscopic value σ, σ = 1 N N i σ i . This susceptibility χ σ quantifies the degree of collectivity in the stress fluctuations: χ σ 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 χ σ as a function of the volume fraction ϕ.
In the low ϕ regime, χ σ increases with increasingγ. However, interestingly, for low rates (γ ≤ 10 −5 ), χ σ hardly depends onγ. This behavior is in contrast to that of σ , which depends linearly onγ for all shear ratesγ in the low ϕ regime. Regarding the volume fraction dependence in this regime, χ σ grows weakly with increasing volume fraction.
In the high ϕ regime, the opposite trend is observed: χ σ becomes smaller when eitherγ or ϕ increases. Still, theγ dependence disappears for low values ofγ (in this case,γ ≤ 10 −6 ), in accordance with the behavior in the low ϕ regime.
At an intermediate value of ϕ between these two regimes, χ σ exhibits a clear peak. As the shear rateγ increases, the height of the peak decreases, and the position shifts towards the low ϕ direction. Note that if we further increase the shear rate toγ = 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 ≤ ϕ ≤ 0.65. In accordance with the convergence of χ σ both in the high and low ϕ regimes, the height and position of the peak become almost constant forγ ≤ 10 −6 . This total convergence of the susceptibility χ σ in the low rate regime over all values of ϕ 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 B, although we leave the precise identification for a future study.
Hereafter, we call the height and position of this peak χ max σ (γ) and ϕ χ max σ (γ), respectively (we omit the explicit notation for theγ dependence below). 1 The only difference between the definitions of σ 2 and σ 2 local is the order in which the averages are taken over particles and time.

C. Stress-strain curves
To further obtain an intuitive understanding of the parameter dependence of the susceptibility χ σ , we plot typical stress-strain curves for the systems under various combinations of the volume fraction ϕ and the shear ratė γ (γ = 10 −6 , 10 −5 and 10 −4 ) in Fig. 2. For the whole parameter space investigated here, the average stress σ becomes larger with both increasing ϕ and increasingγ, as presented in Fig. 1a. However, the dependence oṅ γ changes significantly depending on ϕ: while the order of σ remains the same regardless of the value oḟ γ at a high volume fraction (ϕ = 0.65 > ϕ χ max σ , Fig. 2 light gray curves) 2 , it scales linearly withγ at a low volume fraction (ϕ = 0.62 < ϕ χ 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 ϕ χ max σ under a slow shear rate (γ = 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 ϕ (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 ϕ = 0.65. Importantly, the probability distribution of σ, P (σ), exhibits a power-law-like shape for ϕ = ϕ χ max σ andγ = 10 −6 , 10 −5 , indicating that this susceptibility peak reflects the criticality expected forγ → 0 (see Appendix C). As the shear rate increases, the spikes become less sharp and less frequent (γ = 10 −5 ; Fig. 7b dark gray curve), and finally, the whole stress-strain curve becomes almost detached from that for a low ϕ atγ = 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 (γ = 0) [15]. Notably, under a high shear rate (γ = 10 −4 ), P (σ) exhibits a clear unimodal shape without power-law tails at either end (see Appendix C). This observation suggests that the increase in σ becomes more similar to a cross-over rather than a phase transition because of the effect of the strong external field (see Sec. IV A 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 σ can be viewed as a natural "conjugate" variable to the strength of the external field (namely, the shear rateγ). However, σ changes its value by orders of magnitude depending onγ 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.

A. Correspondence to conventional criticality in equilibrium systems
To further explore the parameter dependence of the shear stress σ 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 β, which is the control parameter of the criticality in the Ising model, corresponds to the volume fraction ϕ in our system. Similarly, the external magnetic field h and the magnetization m correspond to the shear rateγ and the mean stress σ , respectively. Moreover, in both systems, as the external field (h orγ) increases, the change in the order parameter (m or σ ) becomes less steep, and the whole plot shifts towards the less-ordered side. Regarding the susceptibility (χ or χ σ ), we observe peaks at a value of the control parameter (β or ϕ) 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 γ but the shear rateγ, 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 [30,31]. In this sense, the shear stress σ can be viewed as a natural conjugate variable to the external field and thus as an order parameter. However, since σ 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 oṅ γ even in the "disordered", dilute state. In the next sec-tion, we instead measure the average and susceptibility of the interparticle contact number, the values of which are expected to exhibit lessγ dependence.

B. Contact number
The interparticle contact number z characterizes the jamming transition most directly in terms of the microscopic structures [4]. For the jamming transition in quiescent systems without external fields, z changes discontinuously from zero to approximately z C at the critical point ϕ 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 [2]. According to Maxwell's condition, z 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 ϕ in Fig. 3. For the definition of the susceptibility χ z , we employ a definition similar to Eq. 2.
The dependence of the average contact number z on the volume fraction ϕ is qualitatively very similar to that of the average stress σ : it is almost constant in the low ϕ regime and then shows sudden growth around ϕ J , after which the growth rate decreases in the high ϕ regime. However, the dependence onγ is significantly different from that of σ : in the low ϕ regime, the plateau disappears for highγ, and the shear rate dependence is not linear. Furthermore, the values of z at the highest ϕ hardly depend onγ.
The susceptibility of the contact number χ z behaves qualitatively very similarly to that of χ σ : it exhibits a clear peak near the jamming point ϕ J , and the peak height and position change in the same way as χ σ wheṅ γ increases. One major difference from χ σ is that the peak position and height of χ z obviously change even in the low rate limitγ ≤ 10 −6 , where χ σ becomes constant.
This qualitative difference intriguingly suggests that the characteristic lengths that govern σ and z (ξ σ and ξ 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 [32][33][34][35]. This length scale, often referred to as l c , is a candidate for ξ 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 (γ = 0) [36,37]. This knowledge implies that l c and ξ z can be different in nature, since l c is measured in the absence of an external field (γ = 0), while ξ z should be measured in the steady state γ > O(1). As an example of a correlation length measured in a dynamic situation, Refs. [14,22,38] reported that the correlation length of the nonaffine velocities of particles diverges in the limit of ϕ → ϕ J andγ → 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 [39]. Instead, in ref. [39], 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γ → 0 [26,28]. However, this length scale can be well defined only in the high ϕ regime, where the Herschel-Bulkley law is valid and cannot describe the total convergence of χ σ over the whole ϕ 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 oḟ γ in Fig. 4. In this plot, we compare the results for χ σ and χ z . These ridges can be regarded as the dissipativesystem counterpart of the Widom lines by definition. Both Widom lines seem to converge to ϕ J in the limit ofγ → 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.

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 χ z are qualitatively consistent with those for χ σ , 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γ → 0.
As a future direction, an investigation of whether modification of the physical dimension [38][39][40], the damping  Fig. 1(b). The dotted lines depict the location of the jamming point ϕJ.
In Fig. 5, we plot the average shear stress σ in the system as a function of the shear rateγ with ϕ = 0.65.
To roughly estimate the yield stress σ 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 [15,25,27]. We employ the FIRE algorithm [20] and terminate the iteration when the maximum magnitude of the force exerted on one particle meets f max < 10 −9 . The strain increment is ∆γ = 5 × 10 −5 [15]. Fig. 5 shows that the average stress converges to the AQS value at very slow shear rates (γ ≤ 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, σ ∼ σ Y +γ n . This simple estimation provides n ∼ 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 [26,28]. tion. 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 (γ ≥ 0.5), we locate the jamming point as ϕ J ≈ 0.6461 (Fig. 6). This value of ϕ J is consistent with the one estimated by directly fitting the diverging trend of the viscosity [5]. We stress that the data in Fig. 6 is the average over 60 samples.

Appendix C: Stress-strain curves in a normal plot
In this section, we present a normal plot of the stressstrain 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 ϕ = ϕ χ 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 σ for the same combinations of parameters. At a low volume fraction ϕ = 0.62, the PDF is almost Dirac's delta function for all shear rates (the width is very narrow). At a high volume fraction ϕ = 0.65, the PDF is unimodal, with a large width for all shear rates. At ϕ = ϕ χ max σ , however, we observe shear rate dependence. Although the PDF exhibits a powerlaw-like shape for slow shear rates (γ ≤ 10 −5 ), it becomes rather regular unimodal shape for a high shear rate (γ = 10 −4 ).

Appendix D: Ising model
In this section, we recapitulate the famous selfconsistent 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 ddimensional Ising-type spin system on a regular lattice whose Hamiltonian H is written as: where S i ∈ {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 1 2 magnetization m as: where z c = 2d is the spin coordination number. If we employ d = 3 and J = 1, the critical inverse temperature β c is obtained as β c = 1/Jz ≈ 0.167. We plot the values of the magnetization m and the susceptibility χ = dm/dh as functions of β for various values of h in Fig. 9.