Characterizing rare fluctuations in soft particulate flows

Soft particulate media include a wide range of systems involving athermal dissipative particles both in non-living and biological materials. Characterization of flows of particulate media is of great practical and theoretical importance. A fascinating feature of these systems is the existence of a critical rigidity transition in the dense regime dominated by highly intermittent fluctuations that severely affects the flow properties. Here, we unveil the underlying mechanisms of rare fluctuations in soft particulate flows. We find that rare fluctuations have different origins above and below the critical jamming density and become suppressed near the jamming transition. We then conjecture a time-independent local fluctuation relation, which we verify numerically, and that gives rise to an effective temperature. We discuss similarities and differences between our proposed effective temperature with the conventional kinetic temperature in the system by means of a universal scaling collapse.

One can see that the instantaneous velocity profile is non-monotonic and deviates from the averaged linear profile at each point. Furthermore, at four points, marked by the maroon arrows, the velocity gradient becomes negative, i.e., δv i < 0. We refer to this as negative shear velocity or negative velocity gradient. In this simulation the system size is L = 30, the packing fraction is φ = 0.9, and the shear rate isγ = 0.02.   (φ-φ J )/γ q/y P(φ)/γ q   Scaling collapse of T e and T g for highly dissipative system. T g /|δφ| y and T e /|δφ| y vsγ/|φ − φ c | y/q are plotted. Filled symbols refer to T g and hollow symbols correspond to T e . The color code shows different shear rates. Here, we adopt exactly the same exponents and critical densities as we used for the quasi-elastic systems. This demonstrates the robustness of universality against changes in the dissipation in the system. The damping coefficient is γ = 8.95 and the spring constant is Y = 100.

A. Numerical Constants
Numerical constants are as follows: i) the time step is dt = 10 −4 , ii) the elastic and dissipative constants for the linear Dashpot model are Y = 100 and γ = 0.315, respectively.
The latter constants, give rise to a time scale t = γ/Y = 3.15 × 10 −3 which can be used to non-dimensionalize other quantities.

B. Artificial attractive force and coefficient of restitution
As a consequence of the dissipation in the linear Dashpot model, absolute value of the post-collisional normal relative velocity, g , will be smaller than the pre-collisional value, g, i.e., g < g. Coefficient of the restitution is defined as = g /g which is then in the range 0 < < 1 where 0 and 1 correspond to fully dissipative and elastic collisions, respectively.
By integrating Newton's equation of motion, the coefficient of the restitution, , can be easily calculated where force. To prevent this, we enforce zero force between colliding particles when the net force becomes negative. This modification increases the restitution coefficient to 0.96 [1].

C. Measurements
Simulations are done in a large range of parameters. The packing fraction is considered in the range 0.7 < φ < 0.9 where the jamming is φ J = 0.84 [2]. The shear rate is in the range 0.02 <γ < 0.2. Dimensionless shear rate,γ can be defined asγ = t ×γ =γ × γ/Y . Therefore, the dimensionless shear rate will be in the range 10 −5 γ 10 −4 . For the finite size scaling analysis, we adopt different system sizes; L = 10, 20, 30, and 50. As a result, number of particles varies between 100 to 3000.
Rather than fixing the total time for simulations with different shear rates, we fix the total strain to 35 × L where L is the box size. As a consequence, the total simulation time increases for smaller values of shear rate. We start to measure physical quantities when 5 × L strain is reached. This limit ensures balance between the injection of energy and dissipation due to the Dashpot model.

Supplementary Note 2. RELEVANT PHYSICAL QUANTITIES
Let us first define the shear stress σ xy which is often measured and used in our calculations.
At the spatial position x, the local shear stress σ xy (x) is calculated according to where A = 2 · L is the area of of the rectangular bin at position x, v i,x and v i,y are the x and y components of the velocity of the ith particle, respectively, v y (x) is the local drift velocity at x, r ij,x is the x component of the vector which connects the center of particle i to that of particle j, and F ij,y is the y component of the vector force exerted by particle i to j. We only consider the spring (potential) force in the calculation of the stress tensor.
The first sum in the right hand side of Supplementary Eq. (2) runs through all particles at spatial position x, and the second sum runs over all neighboring particles that interact with particle i. If the neighboring particle belongs to a neighboring bin, we only consider half of the second term for each bin. With this local shear stress, it is straightforward to calculate the global shear stress Σ xy in the system: Granular temperature T g is defined as the random kinetic energy of particles. Here, the randomness means that the drift velocity is subtracted from velocity of each particle, so that the temperature can be defined as: where the sum over i runs over N number of particles.

A. Global injected power
In Lees-Edwards Boundary Conditions (LEBC) the energy injected into the system is defined as where Σ xy is the time averaged global shear stress in the system,γ is the shear rate, and L

B. Local injected power
Using Supplementary Eq.5, it is straightforward to calculate the local injected power which is the main physical quantity investigated in our study. We consider the space between centers of two neighboring bins as a thermodynamic subsystem. Each bin is identified by an integer number i, which runs from 0 (the most left bin), to L/2 − 1 (the most right bin).
Then the local power flux in 2D can be defined as with where σ where , denotes the time averaging. Since we have already shown that velocity gradient and shear stress are mutually exclusive, then to a good extent we can write σ (i) xy × δv i = σ (i) xy δv i . Using the fact that the time averaged mean shear stress is independent of the spatial position, we get: Because of the linear velocity profile, the time averaged mean velocity gradient is also the same for each point Since the distance between two neighboring bins is equal to 2, therefore δv /2 is equals to the global shear rate,γ, as a results 2) Reverse flow p i < 0.

Supplementary Note 4. ORIGIN OF THE NEGATIVE POWER
According to Supplementary Eq.6, a reverse flow −p can be either caused by the reversal of the shear force or the reversal of the shear velocity. In order to clarify the origin of negative power in the preceding sections, we illustrate cases where either of these quantities become negative.

A. Negative shear velocities
In LEBCs, the velocity profile v y (x) attains a trivial linear function whereγ is the shear rate defined as the boundary velocity v B divided by the system size L and x is the spatial coordinate. Supplementary Eq. 13 is only recovered when instantaneous velocity profiles over different times are averaged. However, instantaneous profiles possess fluctuations. For the time averaged profile, the difference of the drift velocity, δv i , between two subsequent bins, i, and, i + 1, is simply where d bin = 2 is the distance between centers of two successive bins. and may even become negative, i.e., δv i < 0. We refer to the latter as negative shear velocity or negative velocity gradient. Whereas after averaging over different snapshots we can recover a linear profile given by Supplementary Eq. 13, we observe huge deviations from the linear profile for any instantaneous snapshot. Supplementary Fig. 1 shows the local drift velocity v y (x) as function of the spatial position x. The solid blue line shows the averaged linear velocity profile. However, one can see that the instantaneous velocity profile (the golden bars) deviates dramatically from the linear dependence. Moreover, we find that difference of the drift velocity between neighboring bins can also be negative, i.e., δv i < 0.
Examples of the negative shear velocities are marked by the maroon arrows. The packing fraction for this snapshot is φ = 0.9. For φ > φ J , flow profiles are non-monotonic which is the main reason for negative velocity gradients. However, for φ < φ J , although instantaneous velocity profiles have some fluctuations, the flow profiles respect monotonicity.
To investigate the origin of the negative shear velocities let us consider the following: If distribution functions of the drift velocities (v y ) inside neighboring bins overlap, then the reversal of the shear velocity will be a trivial consequence of the overlapping PDFs. To investigate this, we introduce a characteristic quantity as α = (v y − v y (x)) 2 1/2 /d binγ which is ratio of the width of the distribution function of the drift velocities and difference of the velocity drift between two neighboring bins. If α 1 then the negative velocity gradients are trivial. However for α << 1 probability of overtaking from PDF of the neighboring bin will be negligible and these events will be rather due to interesting physics. To get insights about the nature of the negative shear velocities, within a good approximation one can suppose that distribution of velocities inside the bins are Gaussian; P (v y ) ∝ e −(vy−vy(x)) 2 /2σ 2 where σ 2 = (v y − v y (x)) 2 = 2T g /m, and m = π × 0.5 2 . Therefore, α ≈ T g /γ. According to our computations, we find that for φ >> φ J , T g = 0.22γ. Therefore, for that region α ≈ 0.3.
Furthermore, for φ ≈ φ J , T g = 0.32γ 1.5 which accordingly α ≈ 0.32 √γ which for the probed range of the shear rates we obtain 10 −3 < α < 10 −2 . Also for φ << φ J , T g ∝γ 2 , and thus 10 −4 < α < 10 −3 . Therefore, one can see that for all packing fractions we find α < 1, and thus the negative shear velocity events are not trivially due to the overlap of the PDFs of the drift velocities of the particles.

B. Negative shear stress
The shear stress consists of two independent terms i) kinetic and ii) static (the first and second terms on the rhs of Supplementary Eq. 2). In the range of the shear rates in our study, the kinetic term is an order of magnitude smaller than the static term. Therefore, we only focus on the behavior of the static term. A particle can contribute in the static stress tensor only when it interacts with at least one of its neighbors. A typical example is given in Supplementary Fig. 2 where two colliding particles at distance δr are depicted.
The force between two particles, F , is given by the linear Dashpot model with components F x = F cos θ and F y = F sin θ where cos θ = δx/δr and sin θ = δy/δr. The contributing term to the shear stress is given by In the above equation, strength of the force, F , and distance between particles, dr, are both strictly positive. Therefore, the only term which can become negative is δx · δy. We show in Supplementary Fig. 2 that it is indeed the local orientation of force chains which can result to fluctuations in the sign of the stress tensor. Specifically, when the neighboring particle is either at first or third quadrants (cases a and b in Supplementary Fig. 2), Supplementary Eq.15 gives rise to a positive contribution. However, for the second or fourth quadrants (cases c and d in Supplementary Fig. 2) a negative contribution to the shear stress is obtained. We shall call the former and the latter as the positive and negative orientations, respectively. As a result, a negative stress is related to the local orientations of force chains in the system. Now the question is that why physics of the system depends on the local orientation of the force chains. One should note that the symmetry of the system is broken by the applied external shearing. A positive orientation corresponds to a head-to-head collision which is naturally caused by the shearing. However, the negative orientation can only occur in the case of a backup collision. A backup collision is provided only when two particles move opposite to the direction dictated by the shearing flow. This is the reason why negative stress is a rare event and it is suppressed by the shearing field.
Inspection of Fig. 2 in the main part of the manuscript reveals that for a given packing fraction, φ, probability of the negative shear stress decreases upon increasing the shear rate. This can be justified by the fact that shearing favors the head-to-head collisions which result to positive shear stress. On the other hand, for a given shear rate,γ, probability of the negative shear stress also decreases by increasing the packing fraction, φ. This can be justified by the fact that by increasing the packing fraction, one systematically increases coordination number of particles, z. At jamming density, φ J , the mean coordination number is z c = 4. For the fluid case z < z c , and for the jammed phase z ≥ z c . Therefore, one can see that by increasing the packing fraction, number of contacts per particle will increase. Since the shearing suppresses backup collisions, increasing the coordination number of particles will result to more head-to-head collisions. This results to a decrease of probability of negative events upon increasing the packing fraction.

BALANCE OF PROBABILITIES
Negative power p < 0 can occur for i) negative shear stress σ − xy and positive velocity gradient δv + or ii) negative velocity gradient δv − and σ + xy . Therefore, probability for negative power, P (p < 0), should be equal to the sum of the two joint probabilities of the aforementioned cases through In Supplementary Fig. 3 we compile all our data to show the validity of this equation.
Near jamming, where the negative events are suppressed, we obtain a large scatter due to large error bars.

Supplementary Note 6. TIME SCALE FOR THE FLUCTUATION RELATION
In the main part of the manuscript we conjecture a fluctuation relation which requires a time scale τ . In order to determine this time scale, we need to look again at the PDFs of the injected power ( Fig. 1(a)-Inset in the main manuscript). As it was mentioned, for very small shear rates, the PDF becomes symmetric around p = 0. As a result, Pγ →0 (p) Pγ →0 (−p) which gives zero slope for ln P (p)/P (−p) vs p. This is the lower limit of our fluctuation relation. On the other hand, for very large shear rates, the negative part of the PDFs diminish which results to approximately zero probability of negative power Pγ →∞ (−p) 0.
This determines the upper limit of the shear rate above which the proposed fluctuation relation can not be applied. By a more careful look at our data, we decided the following range, 0.005 <γ < 0.5. In order to determine the relevant time scale in this range of shear rates, we look at the ratio of the elastic and the dissipative parts of the inter-particle forces.
We observe that the former is larger than the latter by an order of magnitude. Therefore, in this range, the only relevant time scale comes from the elastic interactions, viz. τ = τ e .
The elastic time scale represents the typical time scale of a collision which is given by the following relation where m and M are masses of the small and large particles, respectively, and Y is the pre-factor of repulsive force in the linear Dashpot model. Each term in Supplementary Eq. 17 refers to half of the period of oscillation caused by the linear repulsive force.

Supplementary Note 7. CONSISTENCY OF CRITICAL EXPONENTS
A recent estimation of the critical exponents of jamming has been done by Vagberg et al. [3] in which the authors establish a scaling theory for the "average quantities" such as pressure p and shear stress σ in the system where the averages are taken over an ensemble of configurations in the steady state. Since p and σ are both components of a unified tensor, their scaling exponents y are assumed (and numerically shown) to be the same. The authors show that such average quantities follow a universal scaling form which, in the leading scaling form, is given by σ =γ y/zν f (δφ/γ 1/zν ) with f (x) being a universal scaling function and three critical exponents y, z and ν. Such scaling ansatz suggests that when the data for σ/γ q with q = y/zν is plotted versus δφ/γ 1/zν = δφ/γ q/y (or equivalently σ/δφ y versusγ/δφ y/q ), all data at different values of φ andγ will collapse onto a single master curve. As stressed in the above paper, the exponent y is specific to the observable being measured, and indeed is so for q. Therefore, our exponents for T e and T g obtained in data collapse in Fig. 4 of our manuscript, need not necessarily be the same as those reported in the paper [3] for p and σ. On the other hand, the critical exponents should characterize the true algebraic divergence of the observables asymptotically close to the athermal jamming critical point, i.e., φ → φ J andγ → 0 which basically needs simulations of very large system sizes. Since our fluctuation observation is specific to small system sizes we can not naturally approach to the true asymptotic region. In this sense, we can rather compare our findings with other studies in which the chosen parameters are close enough to ours.
Hatano [4][5][6] reported critical exponents of the granular temperature and shear stress in the range of shear rates similar to ours. According to these papers, our critical exponent q = 1.5 for T g is in perfect agreement with that reported by Hatano for the kinetic temperature.
However, in order to establish a connection between the scaling behavior of our data and those observed by Hatano (as an independent test for our data), we have measured "average power flux" P (φ) as a function of the packing fraction φ for various shear ratesγ and system size L = 30-see Main Supplementary Fig. 4. The vertical dashed line refers to the critical point φ J = 0.84335(10) for the jamming transition reported in [3]. As shown in the Inset of Supplementary Fig. 4, the mean power flux at the jamming transition exhibits a scaling behavior withγ, i.e., P (φ = φ J ) ∼γ q with q = 1.72(10). Our exponent q agrees with the critical exponent y γ = 0.63(2) reported in [6] for average shear stress (note that P (φ = φ J ) ∝ σ.γ ∼γ yγ .γ ∼γ 1+yγ ). Since we adapt a more accurate critical point than that used by Hatano, our q exponent is slightly larger than the value reported by Hatano, but still in the same range within the numerical accuracy. Rescaled data for power flux witḣ Although the exponents q and y are specific to the corresponding observable which is measured, their ratio q/y = 1/zν should be independent of the measured observable. In this sense our obtained exponent q/y = 0.75(15) is larger than that obtained by Hatano as q/y = 0.42 (4). The reason is that here we demand a data collapse for a rather wide range of packing fractions far from the critical jamming point. In order to show that our exponents are in perfect agreement with those reported by Hatano, we focused on the critical region and generated more data around the critical point- Supplementary Fig. 7(a). We then estimated the exponent q/y around φ J and find that q/y = 0.49(10) in perfect agreement with that reported by Hatano withing the error bars- Supplementary Fig. 7(b).

Supplementary Note 8. VERIFICATION OF FR IN HIGHLY DISSIPATIVE REGIME
So far we have considered only quasi-elastic particles which correspond to granular materials. The value t = γ/Y = 3.15 × 10 −3 corresponds to the restitution coefficient ≈ 0.98