Critical Scaling of Diffusion Coefficients and Size of Rigid Clusters of Soft Athermal Particles Under Shear

We numerically investigate the self-diffusion coefficient and correlation length (i.e., the typical size of the collective motions) in sheared soft athermal particles. Here we find that the rheological flow curves on the self-diffusion coefficient are collapsed by the proximity to the jamming transition density. This feature is in common with the well-established critical scaling of flow curves on shear stress or viscosity. We furthermore reveal that the divergence of the correlation length governs the critical behavior of the diffusion coefficient, where the diffusion coefficient is proportional to the correlation length and the strain rate for a wide range of the strain rate and packing fraction across the jamming transition density.


INTRODUCTION
Transport properties of soft athermal particles, e.g., emulsions, foams, colloidal suspensions, and granular materials, are important in science and engineering technology [1]. In many manufacturing processes, these particles are forced to flow through pipes, containers, etc. Therefore, the transportation of "flowing particles" is of central importance for industrial applications [2] and thus there is a need to understand how the transport properties are affected by the rheology of the particles.
On the other hand, the mass transport or self-diffusion of soft athermal particles seems to be controversial. As is the rheological flow behavior on shear stress or viscosity, the diffusivity of the particles under shear is also dependent on bothγ and φ. Its dependence onγ is weakened with the increase ofγ , i.e., the diffusivity D exhibits a crossover from a linear scaling D ∼γ to the sub-linear scaling D ∼γ q at a characteristic shear rateγ c , where the exponent is smaller than unity, q < 1 [20][21][22][23][24][25]. For example, in molecular dynamics (MD) simulations of Durian's bubble model in two dimensions [20,21] and frictionless granular particles in three dimensions [22], the diffusivity varies from D ∼γ (γ <γ c ) to D ∼γ 0.8 (γ >γ c ). These results also imply that the diffusivity does not depend on spatial dimensions. However, another crossover from D ∼γ to D ∼γ 1/2 was suggested by the studies of amorphous solids (though the scaling D ∼γ 1/2 is the asymptotic behavior in rapid flowsγ ≫γ c ) [23][24][25]. In addition, it was found in MD simulations of soft athermal disks that, in a sufficiently small flow rate range, the diffusivity changes from D ∼γ (φ < φ J ) toγ 0.78 (φ ≃ φ J ) [26], implying that the crossover shear rateγ c vanishes as the system approaches jamming from below φ → φ J .
Note that the self-diffusion of soft athermal particles shows a clear difference from the diffusion in glass; no plateau is observed in (transverse) mean square displacements (MSDs) [20,[22][23][24][25][26]. The absence of sub-diffusion can also be seen in quasistatic simulations (γ → 0) of soft athermal disks [27] and MD simulations of granular materials sheared under constant pressure [28].
Because the self-diffusion can be associated with collective motions of soft athermal particles, physicists have analyzed spatial correlations of velocity fluctuations [3] or non-affine displacements [29] of the particles. Characteristic sizes of collectively moving regions, i.e., rigid clusters, are then extracted as a function ofγ and φ. However, there is a lack of consensus on the scaling of the sizes. For example, the size of rigid clusters ξ diverges as the shear rate goes to zeroγ → 0 so that the power-law scaling ξ ∼γ −s was suggested, where the exponent varies from s = 0.23 to 0.5 depending on numerical models and flow conditions [22,28]. The dependence of ξ on φ is also controversial. If the system is below jamming, the critical scaling is given by ξ ∼ | φ| −w , where different exponents in the range between 0.5 ≤ w ≤ 1.0 have been reported by various simulations [3,[29][30][31]. In contrast, if the system is above jamming, the size becomes insensitive to φ (or exceeds the system size L) as only L is the relevant length scale, i.e., ξ ∼ L, in a quasi-static regime [4,[23][24][25].
From a scaling argument [27,28], a relation between the diffusivity and size of rigid clusters was proposed as where d 0 is the particle diameter. It seems that previous results above jamming, i.e., asγ increases, D/γ changes from constant tȯ γ −1/2 and corresponding ξ undergoes from L toγ −1/2 , support this argument [23][24][25]. However, the link between the diffusivity and rigid clusters below jamming is still not clear.
In this paper, we study the self-diffusion of soft athermal particles and the size of rigid clusters. The particles are driven by simple shear flows and their fluctuating motions around a mean velocity field are numerically calculated. From numerical results, we extract the diffusivity of the particles and explain its dependence on the control parameters, i.e.,γ and φ. We investigate wide ranges of the control parameters in order to unify our understanding of the diffusivity in both slow and fast flows (γ <γ c andγ >γ c ) and both below and above jamming (φ < φ J and φ > φ J ). Our main result is the critical scaling of the diffusivity D, which parallels the critical scaling of the size of rigid clusters ξ . We find that the linear relation between D/γ and ξ (Equation 1) holds over the whole ranges ofγ and φ if finite-size effects are not important. In the following, we show our numerical method in section 2 and numerical results in section 3. In section 4, we discuss and conclude our results and future works.

METHODS
We perform MD simulations of two-dimensional disks. In order to avoid crystallization of the system, we randomly distribute an equal number of small and large disks (with diameters d S and d L = 1.4d S ) in a L×L square periodic box [32]. The total number of disks is N = 8192 and the packing fraction of the disks φ is controlled around the jamming transition density φ J ≃ 0.8433 [3]. We introduce an elastic force between the disks, i and j, in contact as f e ij = kδ ij n ij , where k is the stiffness and n ij ≡ r ij /|r ij | with the relative position r ij ≡ r i − r j is the normal unit vector. The elastic force is linear in the overlap δ ij ≡ R i + R j − |r ij | > 0, where R i (R j ) is the radius of the disk i (j). We also add a damping force to every disk as where η, v i , and u(r) are the damping coefficient, velocity of the disk i, and external flow field, respectively. The stiffness and damping coefficient determine a time scale as t 0 ≡ η/k.
To simulate simple shear flows of the system, we impose the external flow field u(r) = (γ y, 0) under the Lees-Edwards boundary condition [33], whereγ is the shear rate. Motions of the disks are described by overdamped dynamics [3,4,11], i.e., j =i f e ij + f d i = 0, so that we numerically integrate the disk velocity v i = u(r i ) + η −1 j =i f el ij with a time increment t = 0.1t 0 . Figure 1 displays snapshots of our MD simulations, where the system is sheared along the horizontal arrows in ( Figure 1A). Here, the shear rate is fixed toγ = 10 −7 t −1 0 and the packing fraction increases from φ = 0.8 ( Figure 1A) to 0.9 ( Figure 1D). We observe that, as φ increases, the magnitude of non-affine velocities, |v i − u(r i )|, is pronounced and mobile particles (satisfying |v i − u(r i )| > 10d 0γ with the mean disk diameter d 0 ≡ (d S + d L )/2) form clusters (filled circles). The size of these clusters, or rigid clusters, increases with φ ( Figures 1B,C) and reaches the box size L if the system is far above jamming ( Figure 1D).
In the following, we analyze the data in a steady state, where shear strain applied to the system is larger than unity, and scale every time and length by t 0 and d 0 , respectively.

RESULTS
In this section, we show our numerical results of the selfdiffusion of soft athermal particles (section 3.1). We also extract rigid clusters (Figure 1) from numerical data in order to relate their sizes to the diffusivity (section 3.2). We explain additional data of the rheology and non-affine displacements in Supplemental Materials.

Diffusion
We analyze the self-diffusion of soft athermal particles by the transverse component of mean squared displacement (MSD) [24][25][26]28], Here, y i (τ ) is the y-component of particle displacement and the ensemble average . . . is taken over different choices of the initial time (see Supplemental Materials for the detail) 1 .  (Figure 2B)γ . The horizontal axes are the time interval scaled by the shear rate, γ ≡γ τ , i.e., the shear strain applied to the system for the duration τ . As can be seen, every MSD exhibits a crossover to the normal diffusive behavior, (τ ) 2 ∼γ τ (dashed lines), around a crossover strain γ = γ c ≃ 1 regardless of φ andγ . The MSDs below jamming (φ < φ J ) monotonously increase with the increase of packing fraction, while they (almost) stop increasing once the packing fraction exceeds the jamming point (φ > φ J ) (Figure 2A). The dependence of MSDs on the shear rate is monotonous; their heights decrease with the increase ofγ ( Figure 2B). These trends well-correspond with the fact that the non-affine displacements are amplified in slow flows of dense systems, i.e.,γ t 0 ≪ 1 and φ > φ J [34]. In addition, different from thermal systems under shear [14,35,36], any plateaus are not observed in the MSDs. Therefore, neither "caging" nor "sub-diffusion" of the particles exists in our sheared athermal systems [22,27,37].
To quantify the normal diffusion of the disks, we introduce the diffusivity (or diffusion coefficient) as 2 .
(3) Figure 3A shows double logarithmic plots of the diffusivity (Equation 3) over the shear rate D/γ , where symbols represent the packing fraction φ (as listed in the legend). The diffusivity over the shear rate increases with φ. If the system is above jamming φ > φ J , it is a monotonously decreasing function oḟ γ . On the other hand, if the system is below jamming φ < φ J , it exhibits a crossover from plateau to a monotonous decrease around a characteristic shear rate, e.g.,γ 0 t 0 ≃ 10 −3 for φ = 0.80 [24,25].
In summary, the diffusivity of the disks scales as in the quasi-static regime (γ <γ c ) and D ∼γ 1−ν/λ (5) in the plastic flow regime (γ >γ c ), where the critical exponents are estimated as ν = 1.0 and λ = 4.0. From Equations (4) and (5), we find that the diffusivity below jamming (φ < φ J ) is linear in the shear rate D ∼γ in slow flows, whereas its dependence on the shear rate is algebraic D ∼γ 3/4 in fast flows. A similar trend has been found in molecular dynamics studies of simple shear flows below jamming [22,26,37]. In addition, the proportionality for the diffusivity below jamming diverges at the transition as | φ| −1 (Equation 4), which we will relate to a length scale diverging as the system approaches jamming from below (section 3.2). The diffusivity above jamming (φ > φ J ) implies the crossover from D ∼ | φ| 0.2γ 0.7 tȯ γ 3/4 =γ 0.75 , which reasonably agrees with the prior work on soft athermal disks under shear [26]. Interestingly, the crossover shear rate vanishes at the transition asγ c ∼ | φ| 4.0 , which is reminiscent of the fact that the crossover from the Newtonian or yield stress to the plastic flow vanishes at the onset of jamming (see Supplemental Materials).

Rigid Clusters
We now relate the diffusivity D to rigid clusters of soft athermal particles under shear (Figure 1). The rigid clusters represent collective motions of the particles which tend to move in the same direction [34]. According to the literature of jamming [3,4,29], we quantify the collective motions by a spatial correlation function where v y (x, y) is the transverse velocity field and the ensemble average . . . is taken over disk positions and time (in a steady state). Figure 4 shows the normalized correlation function C(x)/C(0), where the horizontal axis (x-axis) is scaled by the mean disk diameter d 0 .
As can be seen, the correlation function exhibits a well-defined minimum at a characteristic length scale x = ξ (as indicated by the vertical arrow for the case of φ = 0.84 in Figure 4A). Because the minimum is negative C(ξ ) < 0, the transverse velocities are most "anti-correlated" at x = ξ [3]. Therefore, if we assume that the rigid clusters are circular 4 , their mean diameter is comparable in size with ξ . The length scale ξ increases with the increase of φ ( Figure 4A) but decreases with the increase ofγ ( Figure 4B). These results are consistent with the fact that the collective behavior is most enhanced in slow flows of dense systems [34]. As reported in Olsson and Teitel [3], we examine critical scaling of the length scale. Figure 5A displays scaling collapses of the data of ξ , where the critical exponents, ν = 1.0 and 4 In a short range of x, non-affine velocities tend to be aligned in the same direction, i.e., v y (x i , y i )v y (x i + x, y i ) > 0, such that the correlation function C(x) is a positive decreasing function of x. If the non-affine velocities start to align in the opposite direction, i.e., v y (x i , y i )v y (x i + x, y i ) < 0, the correlation function becomes negative, C(x) < 0. Because the correlation function is minimum when the non-affine velocities are located on either side of a vortex-like structure, a typical size of rigid clusters can be defined as the distance at which C(x) becomes minimum. In a long distance of x, due to the randomness of nonaffine velocities, the correlation function eventually decays to zero with or without oscillations [34].  Figure 3B. If the shear rate is smaller than the characteristic value, i.e.,γ <γ c ≃ 10 4 | φ| λ , the data below jamming (φ < φ J ) exhibit plateau, whereas those above jamming (φ > φ J ) diverge with the decrease of shear rate. Therefore, if we assume that the data above jamming follow the power-law (solid line) with the slope −0.40 ± 0.02 (approximately −0.4), the length scale in the quasi-static regime (γ <γ c ) can be described as | φ| ν ξ ∼ J ± (γ /| φ| λ ) with the scaling functions, J − (x) ∼ const. for φ < φ J and J + (x) ∼ x −0.4 otherwise. Note that, however, the length scale is limited to the system size L (shaded region in Figure 5A) and should be scaled as ξ ∼ L above jamming in the quasi-static limitγ → 0 [4,24,29]. This means that the system size is the only relevant length scale [38] and thus we conclude ξ ∼ L in slow flows of jammed systems. On the other hand, ifγ >γ c , all the data are collapsed onto a single power law (dotted line in Figure 5A). Therefore, the scaling functions are given by J ± (x) ∼ x −z such that the length scale scales as ξ ∼γ −z | φ| λz−ν . Because this relation is independent of | φ|, the exponent should be z = ν/λ as confirmed in Figure 5A.
In summary, the length scale, or the mean size of rigid clusters, scales as in the quasi-static regime (γ <γ c ) and in the plastic flow regime (γ >γ c ), where the critical exponents, ν and λ, are the same with those for the diffusivity (Equations 4 and 5). The critical divergence below jamming in the quasi-static regime, i.e., ξ ∼ | φ| −1 (Equation 6), is consistent with the results of quasi-static simulations (γ → 0) of sheared athermal disks [27,29]. In addition, the scaling ξ ∼γ −1/4 in the plastic flow regime (Equation 7) is very close to the prior work on athermal particles under shear [22]. From the results of the diffusivity and length scale, we discuss how the rigid clusters contribute to the diffusion of the particles. From Equations (4) to (7), the linear relation D ∼ d 0 ξγ (Equation 1) holds below jamming (regardless oḟ γ ) and in the plastic flow regime (regardless of φ). We stress that the divergence of the diffusivity over the shear rate in the quasi-static regime, i.e., D/γ ∼ | φ| −1 (Equation 4), is caused by the diverging length scale below jamming, i.e., ξ ∼ | φ| −1 (Equation 6). As shown in Figure 5B, the linear relation (dotted line) well explains our results below jamming, where the length scale ξ is smaller than 10d 0 . The linear relation (dotted line) also well captures our results above jamming 5 until the length scale eventually reaches the system size L/2 ≃ 44d 0 (shaded region).

DISCUSSIONS
In this study, we have numerically investigated rheological and transport properties of soft athermal particles under shear. Employing MD simulations of two-dimensional disks, we have clarified how the rheology, self-diffusion, and size of rigid clusters vary with the control parameters, i.e., the externally imposed shear rateγ and packing fraction of the disks φ.
Our main result is the critical scaling of the diffusivity (section 3.1), where their dependence on bothγ and φ is reported (Equations 4-7). The diffusivity has been calculated on both sides of jamming (by a single numerical protocol) to unify the understanding of self-diffusion of soft particulate systems. We found that (i) the diffusivity below jamming exhibits a crossover from the linear scaling D ∼γ to the power-law D ∼γ 3/4 . Such a crossover can also be seen in previous simulations [20][21][22]. In addition, (ii) the diffusivity below jamming diverges as D ∼ | φ| −1 if the system is in the quasi-static regime (γ <γ c ), whereas (iii) the diffusivity (both below and above jamming) is insensitive to φ if the system is in the plastic flow regime (γ >γ c ). Note that (iv) the crossover shear rate vanishes at the onset of jamming asγ c ∼ | φ| 4.0 . These results (ii)-(iv) are the new findings of this study. On the other hand, we found that (v) the diffusivity above jamming is weakly dependent on φ (as D ∼ | φ| 0.2 ) in the quasistatic regime and (vi) shows a slight change from D ∼γ 0.7 tȯ γ 3/4 . Though the result (v) is the new finding, the result (vi) contrasts with the prior studies, where the diffusivity exhibits a crossover from D ∼γ toγ 1/2 [24,25,28]. Because our scaling D ∼γ 0.7 in the quasi-static regime is consistent with Olsson [26], where the same overdamped dynamics are used, we suppose that the discrepancy is caused by numerical models and flow conditions, where the Lennard-Jones potential is used for the interaction between the particles in Lemaître and Caroli [24] and Chattoraj et al. [25] and frictionless/frictional disks are sheared under constant pressure in Kharel and Rognon [28]. 5 If the system is above jamming (φ > φ J ) and in the quasi-static regime (γ <γ c ), the diffusivity over the shear rate scales as D/γ ∼ | φ| 0.2γ −0.3 (see Equation 4). However, if φ > φ J andγ <γ c , the length scale is ξ ∼ | φ| 0.6γ −0.4 as long as ξ does not exceed the system size L. Though the linear relation D/γ ∼ ξ is not exact in this regime, their exponents are close to each other such that the dotted line in Figure 5B fairly well explains the data. To figure out the relation between D/γ and ξ in this regime more carefully, further studies of different system sizes may be necessary because there is a possibility of the coexistence of different length scales, e.g., not only the rigid clusters but also string-like structures with different fractal dimensions as can be seen in Figure 1.
We have also examined the relation between the diffusivity and typical size of rigid clusters ξ (section 3.2). Below jamming, we found the critical divergence ξ ∼ | φ| −1 in the quasistatic regime as previously observed in quasi-static simulations (γ → 0) of sheared athermal disks [29]. Note that our exponent is larger than that obtained from previous studies of soft athermal particles [3,30], where the critical divergence ξ ∼ | φ| −ν is characterized by the exponent in the range between 0.6 < ν < 0.7. Though the meaning of ξ is similar, in the sense that it quantifies the size of collective motions of the particles under shear, we found that the quality of data collapse of ξ is even worse if we use the exponents ν = 0.6 and 0.7. So far, we cannot figure out the cause of this discrepancy and postpone this problem as a future work. In the plastic flow regime, the size becomes independent of φ and scales as ξ ∼γ −1/4 . This is consistent with the previous result of sheared athermal particles [22] (and is also close to the result of thermal glasses under shear [36]). Above jamming, however, the size exhibits a crossover from ξ ∼ L toγ −1/4 which contrasts with the crossover from ξ ∼ L toγ −1/2 previously reported in simulations of amorphous solids [24,25,28]. From our scaling analyses, we found that the linear relation D ∼ d 0 ξγ (Equation 1) holds below jamming (for ∀γ ) and in the plastic flow regime (for ∀φ), indicating that the self-diffusion is enhanced by the rotation of rigid clusters [3,28].
In our MD simulations, we fixed the system size to L ≃ 88d 0 . However, systematic studies of different system sizes are needed to clarify the relation between D and ξ ∼ L above jamming, especially in the quasi-static limitγ → 0 [24,25]. In addition, our analyses are limited to two dimensions. Though previous studies suggest that the diffusivity is independent of the dimensionality [20][21][22], a recent study of soft athermal particles reported that the critical scaling of shear viscosity depends on dimensions [39]. Therefore, it is important to check whether the critical scaling (Equations 4 and 5) is different in three-dimensional systems. In addition, O'Hern et al. found that the critical exponents near jamming are sensitive to interaction potentials [32]. Therefore, it is important to examine whether our exponents, ν and λ, are sensitive to the details of numerical models, e.g., contact forces, particle inertia, microscopic friction, particle size distributions, etc. with the aim of testing universality class of the diffusivity and size of rigid clusters. Moreover, the analysis of local structures of instantaneous configurations of the disks is important to investigate the rigidity of the clusters [40] and the relation between the diffusivity and shear viscosity may be interesting because it gives a Stokes-Einstein like relation for the nonequilibrium systems studied here.
In our recent study of non-Brownian suspensions [15], we found that the viscosity exhibits a crossover when the system approaches the jamming transition. To describe such a crossover, a correction term should be added to the conventional scaling function as F (x) + G(x), where two different exponents are introduced as F (x) ∼ x −z and G(x) ∼ x −w . This means that the scaling collapse of the flow curves [3], where the scaling function is given by a single power law, can be violated if the system is very close to jamming. Similarly, our scaling collapses of D ( Figure 3B) and ξ ( Figure 5A) could break down by the same reason. Therefore, it is necessary to examine whether the crossovers of D and ξ exist in the case that the system is in the vicinity of the jamming transition.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
KS and TK designed the research and wrote the paper. KS performed the research.

FUNDING
This work was supported by KAKENHI Grant Nos. 16H04025 and 18K13464 from JSPS.