Fluctuations and like-torque clusters at the onset of the discontinuous shear thickening transition in granular materials

The main mechanism driving rheological transitions is usually mechanical perturbation by shear unjamming mechanism. Investigating discontinuous shear thickening is challenging because the shear counterintuitively acts as a jamming mechanism. Moreover, at the brink of this transition, a thickening material exhibits fluctuations that extend both spatially and temporally. Despite recent extensive research, the origins of such spatiotemporal fluctuations remain unidentified. Here, we numerically investigate the fluctuations in injected power in discontinuous shear thickening in granular materials. We show that a simple fluctuation relation governs the statistics of power fluctuations. Furthermore, we reveal the formation of like-torque clusters near thickening and identify an unexpected relation between the spatiotemporal fluctuations and the collective behavior due to the formation of like-torque clusters. We expect that our general approach should pave the way to unmasking the origin of spatiotemporal fluctuations in discontinuous shear thickening. https://doi.org/10.1038/s42005-021-00574-8 OPEN

A s proposed by Liu and Nagel 1 , the dynamics of frictionless particles with elasticity, such as granular media, gels, and emulsions, seizes at high-volume fractions above a jamming density ϕ > ϕ J . The dynamics in a jammed state can be rejuvenated by shearing, and rheology describes the response of matter to shearing. At low flow rates, _ γ ! 0, and below jamming, particulate matter exhibits two different rheologies: (1) inertia-dominated flows, such as granular media [2][3][4][5] , show quadratic-Bagnold rheology 6 , σ xy $ _ γ 2 , and (2) in overdamped systems such as emulsions, in which inertia plays no role, the Newtonian regime is recovered 7 . In the Newtonian regime, shear stress is linearly related to the flow rate, σ xy $ _ γ. Above jamming, the dynamics is arrested, and a threshold stress, σ Y , must be overcome to invoke the flow. In this phase, the relationship between σ xy and _ γ is σ xy ¼ σ Y þ k_ γ n , where k and n are constants. As a result, the viscosity, η ¼ σ xy =_ γ, decreases as a function of rate. This behavior, known as shear thinning, is crucial in technological applications, because shear thinning facilitates the pumping of gels and emulsions.
Counterintuitively, shear may act as a jamming mechanism for granular media. Jamming by shear, or shear jamming (SJ), was originally observed in experiments on granular materials under quasi-static shear, in which jammed states are induced when the shear strain exceeds a critical value 8,9 . This situation changes when the system is subjected to a steady shear: for sufficiently dense granular materials, there exists a critical rate at which the viscosity drastically increases 10 . This abrupt increase in viscosity at a critical rate is generally called discontinuous shear thickening (DST). Whether SJ and DST are related remains unclear 11 .
Likewise, DST occurs in a wide range of complex fluids, such as Brownian and non-Brownian suspensions 12 . In suspensions, the transition results from an interplay between a stabilizing mechanism and frictional contacts due to the roughness of particles. A stabilizing mechanism provides a stress scale for the onset of the transition: when the shear force exceeds this threshold stress, the system undergoes DST. Moreover, at low flow speeds, the stabilizing mechanism keeps particles apart. At high flow speeds, when the threshold stress is overcome by the shear forces, a proliferation of frictional contacts results in an abrupt increase in viscosity by more than an order of magnitude-the hallmark of DST 13 . This behavior is well captured by the Wyart-Cates model, which describes DST as a transition from frictionless to frictional rheologies, owing to the proliferation of frictional contacts 14 . The stabilizing mechanism is system dependent; however, the frictional contacts are essential for DST 15 .
The stabilizing mechanism is well understood in suspensions: (1) in non-Brownian suspensions, this mechanism is usually a repulsive interaction, such as electrostatic repulsion of double layers [15][16][17] , and (2) in Brownian suspensions, this role is played by thermal forces 18 . However, the origin of the onset stress of DST in granular materials is elusive, and consequently DST is not well understood in granular media. Our results shed light on this issue.
Despite the different mechanisms underlying DST in all the aforementioned soft material systems, DST has one consistent aspect: near the thickening transition, a system becomes unstable with spatiotemporal fluctuations. In suspensions, temporal fluctuations appear as oscillations and chaotic time series akin to turbulence 19 , an effect dubbed rheochaos 20 . Moreover, spatial fluctuations result in intermittent stress heterogeneities. These stress anomalies propagate along the vorticity direction and are referred to as vorticity bands 21,22 . The strength of these fluctuations is enhanced with increasing stress, as confirmed by recent experiments using advanced techniques for the measurement of local rheology 19,23,24 . When confined to a micro-channel, dense suspensions of hard colloids exhibit oscillations in the density and velocity of particles 25 . Remarkably, these oscillations can be reproduced in a numerical model when frictional contacts are incorporated 26 . Similar spatiotemporal fluctuations have also been reported in extensive detail for frictional granular materials undergoing DST, in a series of publications by Grob et al. 27,28 and Saw et al. 29 .
These spatiotemporal fluctuations are known to increase dramatically as the onset of DST is approached, yet the origin of these fluctuations remains a mystery. One reason for this lack of understanding is that, similarly to the glass transition 30 , usual measures, such as pair correlation functions, show no dramatic change signatures in the micro-physics across DST 31 . Here, we shed light on the spatiotemporal fluctuations by performing a fluctuation analysis of the injected power. We show that the statistics of power fluctuations obeys a simple fluctuation relation. We specifically focus on rare power fluctuations caused by frictional forces, and we systematically discriminate different contributions resulting in negative power fluctuations. This discriminative approach results in finding a mirror-image decomposition graph. One of the co-authors of this work has previously reported a decomposition graph for frictionless particles undergoing jamming transition 32 . It was shown that a decomposition graph is a blueprint of collective behavior giving rise to rare power fluctuations on both sides of the transition point. Similarly, here we find the formation of like-torque clusters that underly fluctuations in total torque of the system.

Results
Rheology across the thickening transition. We perform two dimensional molecular dynamics simulations of bidisperse frictional disks in a simple shear flow in the underdamped regime (inertial). We neglect thermal and hydrodynamic forces to focus on the central role of frictional forces in DST. However, a generalization of our approach to include hydrodynamic and thermal effects in the overdamped regime would be straightforward. Details on the simulations are given in the "Methods." A typical flow curve of the system is given in Supplementary Figs. 2 and 3 and their corresponding description in Supplementary Note 2, as well as the coordination number in Supplementary Fig. 4, where the onset of thickening is given by a critical shear rate, _ γ c , below which the system is fluid and above which the stress, as well as the viscosity, shows a sudden increase by more than an order of magnitude. For _ γ < _ γ c , the rheology is Bagnold, and for _ γ > _ γ c , the shear stress has a weak dependence on the shear rate. Which mechanism determines the onset stress in DST in granular media is unclear. We adopt a similar approach to that proposed by Kawasaki and Berthier 18 to describe the onset stress in DST in Brownian suspensions. In Fig. 1, we depict the shear stress, σ xy , rescaled by thermal stress given by , where d is the average diameter, T g is the granular temperature (defined in Methods 4), and k n is the elastic constant. DST can be seen to occur when the shear stress in the system reaches the stress scale given by the thermal stress, σ T . We demonstrate that most of the physical measures, such as pressure ( Supplementary  Fig. 2) and coordination number ( Supplementary Fig. 3), represent the inherent state of the system in a manner similar to the rheology. Next, we investigate negative fluctuations in power.
Injected power. We define the power injected by the local shear as (proof in Supplementary Information): where _ γ l is the local shear rate, σ xy is the local shear stress, and δσ xy = σ xy,t − σ yx,t is the local stress difference, which is equal to difference in the off-diagonal components of the tangential part of the shear stress. The derivation of Eq. (1) is given in Supplementary Note 1, where we use the elasticity theory developed specially for micro-polar fluids 33,34 . The stress is computed locally in the rectangle of area A = A bin equal to the length of the system along the shearing direction and the width of w = 4R, with R = 0.7 the radius of larger particles. Therefore, the total power injected into the bin is equal to p bin = A bin p. We cross-checked our results for a wider bin of width w = 4R and obtained similar results. A typical probability distribution function (PDF) of power is shown in the inset in Fig. 2a. The distribution is bell shaped, with exponential tails. Notably, such a PDF has been reported for various far-from equilibrium systems, for example, work on frictional particles 35 , turbulence 36 , frictionless disks 32 , and colloidal suspensions 37 . Thus, this PDF has some common features across various nonequilibrium systems that have been overlooked to date. A power equal to p bin given by Eq. (1) is dissipated in a subsystem, in a manner akin to entropy production in a thermodynamic system. However, owing to a negative fluctuation, the subsystem can give up the power, in a manner akin to entropy consumption 38 . The PDF of power p bin is given by Pðp bin Þ. We conjecture an instantaneous detailed fluctuation relation comparing the ratio of the PDF of the entropy production and consumption rates via: in which 1/β = T e /τ e is the ratio of the effective temperature to the elastic timescale. Similarly to the other fluctuation relations, Eq.
(2) is based on the notion that the probability for entropy production must be exponentially larger than that for the entropy consumption. In the main panel of Fig. 2a, we plot this ratio for various shear rates. A linear dependence is clearly recovered, and thus the fluctuation relation is verified by our data. To clearly show how the effective temperature T e changes across the transition, in Fig. 2b we display the effective temperature T e versus the shear rate. The effective temperature is computed via a direct linear regression by using Eq. (2). Whereas a Bagnold dependence is obtained in the fluid branch, T e resembles the behavior of the shear stress in the thickened branch. Therefore, our proposed fluctuation relation gives rise to an effective temperature that behaves similarly to the rheology. An important aspect affecting fluctuation relations is the time-reversal symmetry, which, in our case, is broken because of the dissipation. Therefore, the existence of a simple fluctuation relation for DST is significant and should motivate further application of thermodynamic principles to DST. Equation (2) gives rise to a trivial result hpi > 0, where h:::i refers to a time averaging operation. This result can be easily verified by numerical simulations and can be interpreted as an equivalent to the second law of thermodynamics for this system. Yet, the PDF of power, as depicted in the inset of Fig. 2a, has a substantial negative part. These negative events, p < 0, correspond to rare fluctuations whose origin is unknown to date. Below, we investigate these rare events to gain insights into fluctuations in DST.
Rare fluctuations due to frictional forces. We identify a rare fluctuation when the second term of the right-hand-side (rhs) of Eq. (1) becomes negative, i.e., p t ¼ δσ xy _ γ l < 0. We note that the stress difference, δσ xy , is related to the torque caused by frictional forces. The first term of the rhs of Eq. (1) can also cause a negative fluctuation. However, in this work, we focus on only the second term to investigate the role played by torque due to friction. In the inset of Fig. 3, the probability for such a negative fluctuation Pðδσ xy _ γ l < 0Þ is plotted versus the shear rate. Pðδσ xy _ γ l < 0Þ is almost constant (within error bars) in the fluid branch, and it decreases monotonically in the thickened branch. Because p t is a product of two terms, the local stress difference δσ xy and the local shear rate _ γ l , a negative p t can be due to either (1) δσ xy < 0 with _ γ l > 0, or (2) _ γ l < 0 with δσ xy > 0. This can be mathematically expressed as a decomposition relation: in which P(,) is a joint probability. In the main panel of Fig. 3, we display Pðδσ À xy ; _ γ þ l Þ (left Y-axis) and Pð_ γ À l ; δσ þ xy Þ (right Y-axis) with blue circles and red squares, respectively. A mirror-image pattern emerges, indicating a decomposition of Pðδσ xy _ γ l < 0Þ into two distinct branches for Pðδσ À xy ; _ γ þ l Þ and Pð_ γ À l ; δσ þ xy Þ. Moreover, this decomposition has two unique features. First, Pðδσ À xy ; _ γ þ l Þ and Pð_ γ À l ; δσ þ xy Þ are approximately mirror images of each other, and second, in both thickened and fluid states, their dependence on _ γ l is opposite, meaning that when Pðδσ À xy ; _ γ þ l Þ increases by _ γ l , Pð_ γ À l ; δσ þ xy Þ inversely decreases, and vice versa. This finding provides evidence that rare fluctuations in DST are governed by simple decomposition relations. We now focus on this relationship between different contributions giving rise to negative fluctuations to determine what lessons might be learned about the nature of DST.
We emphasize again that the local shear rate, _ γ l , and local stress difference, δσ xy , when averaged over long time periods, are definitely positive. However, a rare fluctuation results in an instantaneous negative _ γ l or δσ xy . We start with the interpretation of a rare fluctuation due to _ γ l < 0. A negative local shear rate occurs when a given layer of particles is slower than the one below it. Therefore, a negative fluctuation of _ γ l < 0 corresponds to a non-monotonic change in local drift velocity due to a retarded layer. In Fig. 3 for _ γ < _ γ c , Pð_ γ À l ; δσ þ xy Þ increases as _ γ c is approached. Thus, the flow becomes non-monotonic as the transition point is approached from below. This increasing nonmonotonicity might be related to the well-known instability near _ γ c 27-29 . In the thickened branch for _ γ > _ γ c , the instability is washed out, and Pð_ γ À l ; δσ þ xy Þ decreases by increasing the shear rate. This finding is consistent with conventional wisdom, Fig. 1 Onset stress of discontinuous shear thickening (DST). Shear stress rescaled by a stress scale given by the granular temperature, versus the shear rate _ γ, provides the onset of DST. The system undergoes discontinuous shear thickening when the rescaled stress exceeds unity, and the red dashed line is a guide for eye to mark the threshold. This threshold behavior is akin to that of a Brownian suspension near DST 18 . Each point is an average over a total strain of 20, with the number of configurations equal to the total strain. The packing fraction is ϕ = 0.81, the number of particles is N = 16,384, and the system size is L ¼ ffiffiffi ffi N p . Error bars are smaller than symbols.
because the shear is a bias, and in the absence of instability, it removes all the retarded layers at large _ γ, thus explaining the decreasing trend in Pð_ γ À l ; δσ þ xy Þ in the thickened branch. Pðδσ À xy ; _ γ þ l Þ has opposite behavior with respect to Pð _ γ À l ; δσ þ xy Þ across the transition region. To interpret the behavior of Pðδσ À xy ; _ γ þ l Þ, we first describe the relation of δσ xy to the microphysics. The stress difference is related to the micro-physics via the total torque τ acting on particles according to σ xy,t − σ yx,t ∝ ∑ i τ i , in which i runs through all particles in the (sub-)system 34,39 (a formal definition of torque is given in Supplementary Information). Thus, p t ∝ τ. As a result, Pðδσ À xy ; _ γ þ l Þ corresponds to a rare fluctuation due to the negative total torque of a subsystem. Fig. 3 reveals that Pðδσ À xy ; _ γ þ l Þ decreases as the transition point in the fluid branch is approached. This decrease in the probability of negative torque fluctuations may be interpreted as a uniformity of torque across subsystems upon approaching the thickening point. We investigate this possibility in the next section.

Reinspection of
Formation of torque clusters. To gain a better understanding of the enhancement of the uniformity of torque, we display subsequent snapshots below thickening as the system is sheared from left to right in Fig. 4a-d. In snapshot a, the system is homogeneous except for anomalies appearing as tiny clusters of large negative (blue) and positive (red) like-torque particles. These The inset shows the probability distribution function (PDF) of injected power into a bin versus the power divided by the mean power for various global shear rates. The system is divided into narrow rectangular bins along the shearing direction. It displays a bell-type distribution with exponential tails. In the main panel, we examine the fluctuation relation proposed in Eq. (2), which is successfully verified by our data. In this equation, PDF of power, P(p bin ), is divided by that of negative power, P(−p bin ). Blue and red symbols correspond to regimes below and above thickening, respectively. The dashed line is a guide for eye to show the expected linear dependence. b The effective temperature T e , derived from the fluctuation relation, as a function of the shear rate _ γ is plotted. Bagnold behavior is found in the fluid branch, which crosses over to a rate-independent behavior in the thickened branch, in accordance with the rheology (Supplementary Fig. 1a). The effective temperature increases almost two orders of magnitude after thickening. The dashed line is a guide for eye to show the Bagnold scaling, $ _ γ 2 , in the phase below the thickening transition. Error bars are smaller than symbols. The packing fraction is ϕ = 0.81, and the number of particles is N = 16,384. We create an ensemble by storing configuration of the system at every strain difference of one particle diameter. In total, we average over 2560 configurations.
domains of large torque particles show that large torque is spatially localized. These clusters grow spatially in snapshot b, with more red (positive) clusters. In contrast, in snapshot c, the blue (negative) clusters appear, and finally in snapshot d, the red clusters nearly percolate in the system. To show how the total torque of the system changes in these snapshots, we display the mean torque in the system as a function of strain in panel e. In this figure, the corresponding mean torque of the snapshots is marked by red letters/arrows. The mean torque exhibits an oscillatory behavior whose amplitude is first enhanced, then decays, and finally fades to 0 at the end of the interval. Moreover, the torque subsequently undergoes another oscillatory behavior for γ > 1470. This is a typical pattern that repeats itself throughout the simulations for _ γ < _ γ c . Here, we reach an important conclusion: the oscillatory behavior of the torque (panel e) originates from competition between clusters of large negative and positive like-torque particles. Interestingly, similar spontaneous oscillations have been reported in other shear thickening systems 20,40,41 . A simple mechanism describes the formation and growth of the like-torque clusters. In the fluid branch, particles do not participate in enduring contacts, because the coordination number z < z J = 3, and most of the collisions are binary akin to "spinning tops" on a plane. As dictated by the equations of motion, a nonzero frictional force results in parallel torques, and successive binary collisions between neighboring particles underlie the growth of like-torque clusters.
A typical configuration of the system is displayed after thickening in the thickened branch in Fig. 4f, where the total absolute value of torque in the system, τ abs ¼ ∑ N i¼1 jτ i j, decreases more than two orders of magnitude with respect to that in the fluid branch. For consistency, we use a similar range for color coding. In snapshot f, there are a few anomalous domains of opposite-torque particles. If we decrease the range of torque for color coding, these clusters will span the entire system. Interestingly, in these clusters, particles with anomalously large positive torque co-exist with those of anomalously large negative torque. A close-up of such a cluster is given in panel g. A question arises as to whether the growth of like-torque clusters might correlate with the rheology. To answer this question, in Fig. 5a, d, the color code of each particle corresponds to total torque per particle, and in panels b and e that corresponds to total shear stress per particle. In panel a, the instability has just begun, and there are small anomalies of like-torque clusters. Panel b depicts the same snapshot, in which the color code corresponds to total shear stress per particle. It can be seen that the stress in the system is homogeneous (a close-up is given in panel c). In panel d, the instability has set in and one can see spatially larger anomalous like-torque clusters, and at the same time, in panel e the stress field is quite heterogeneous. Moreover, stress-bearing chains of positive magnitude exist along the compression direction, and similar chains of negative magnitude stress exist along the dilation direction (a close-up is shown in panel f). We emphasize that by correlation between stress and torque, we do not refer to the spatial correlation between liketorque clusters and stress-bearing chains. Instead, the stress in the system becomes heterogeneous (in a patterned manner) when like-torque clusters grow. In a quantitative demonstration, in panel i, the second moments of shear stress per particle (blue circles) and torque per particle (red squares) are shown as a function of the strain. One can see that <σ 2 xy > changes proportionately with <τ 2 >. However, this proportionality is broken in panel j, which is above thickening. A corresponding snapshot can be seen in panels g and h. In Supplementary Fig. 5, we demonstrate this phenomenon with more snapshots. Therefore, our results imply a previously unnoticed link between the collective behavior of torque with the stress heterogeneities, the latter have been a subject of recent debate 19,23,[27][28][29] .

Discussion
We investigated power fluctuations in a model system undergoing DST. We showed that different contributions giving rise to rare fluctuations caused by frictional forces are governed by a simple mirror-image decomposition relation that underlies the collective behaviors across the thickening transition. The joint probability of rare fluctuations due to frictional forces Pðδσ À xy ; _ γ þ l Þ decreases as the thickening transition is approached, thus indicating a collective behavior in the torque. Consequently, we discovered clusters of like-torque particles. We showed that (1) the growth of these clusters directly correlates with the rheology, (2) the formation of the clusters results in spatially heterogeneous structures of stress-bearing chains below the thickening, and (3) a competition between negative and positive like-torque clusters underlies the temporal fluctuations in total torque in the system. After thickening, we observe clusters of opposite-torque particles. Moreover, in this regime, particles grind against each other, as they would be obliged to do by frictional forces that glue particles to one another. This opposite motion is similar to that of "gears" in a mechanical watch. The correlation between the growth of like-torque clusters and the rheology may motivate theoretical work to explore possible connections. Investigating whether this correlation results in a causality of stress heterogeneity as a result of like-torque clusters, or vice versa, should prove interesting.
In most models for suspensions, dissipation occurs through the fluid's viscosity, η, of the host, and as a result, the inertia of particles Fig. 3 Decomposition of probabilities of negative power. Inset: probability of negative power, Pðδσ xy _ γ l < 0Þ where P is the net probability, δσ xy is the stress difference, and _ γ l is the local shear rate defined by Eq. (12), as a function of the global shear rate, _ γ. A power equal to δσ xy _ γ l > 0 must be dissipated; however, in a negative event this power is given up by the subsystem, δσ xy _ γ l < 0. In the thickened phase, Pðδσ xy _ γ l < 0Þ decreases monotonically with increasing the shear rate. No discontinuity is observed at the critical shear rate. Main panel: because (3), we calculate Pðδσ À xy ; _ γ þ l Þ, the probability for negative power due to negative stress difference, and Pð _ γ À l ; δσ þ xy Þ, the probability for negative power due to negative local strain rate, directly from our raw data to see if they can be decomposed into two separate curves. Joint probabilities Pðδσ À xy ; _ γ þ l Þ and Pð _ γ À l ; δσ þ xy Þ are displayed by the blue circles and red squares, respectively. Notably, a decomposition of the joint probabilities into a mirror-image curve has been found. Whereas Pðδσ À xy ; _ γ þ l Þ is reduced almost twice at _ γ c ¼ 10 À5 , and then increases as a function of the shear rate, Pð_ γ À l ; δσ þ xy Þ is enhanced by the same rate, and it decreases with the shear rate in the fluid phase. is neglected. The rheology of this regime is characterized by constitutive relations depending on the viscous number, J. For suspensions with large inertia and for granular materials, the dissipation mainly occurs in the collisions between particles, and the rheology is characterized by an inertia number, I. For a suspension undergoing DST at an intermediate regime, the rheology is characterized by a combination of I 2 and J, as shown recently by Dong and Trulsson 42 . Using this combination, the authors have obtained unified constitutive laws of rheology. Determining how fluctuations arise in suspensions with inertial particles, and whether torque clusters exist and persist in this regime, should prove interesting.
Although the well-known Wyart-Cates model successfully describes the transition from a frictionless to a frictional rheology in DST of overdamped suspensions, growing experimental 24,43 and theoretical 44 evidence supports the existence of correlations and collective behavior across DST. A recent study has suggested that DST might be a critical phenomenon 45 . Most of these studies stem from two notable reports by Lootens et al. 46,47 , in which the authors show giant fluctuation of stress, power-law distribution, and periodicity near the onset of transition of the rheology in a dense suspension. This finding urges revisiting DST beyond mean-field theory, as in a recent such study by Thomas et al. 44 . Our work provides a general framework for investigation of the collective behavior in DST for suspensions. Application of our method to overdamped systems would require care, because the overdamped limit is described by a balance of forces and torques on each particle individually, and as a result, the total torque on all particles is set to 0. Thus, we do not expect formation of like-torque clusters to occur in the overdamped regime. Yet, fluctuations in stress and rate can arise in an overdamped system as a result of flow fluctuations, and our method can be easily adapted for the investigation of negative fluctuations. We believe that the investigation of collective motion/ Fig. 4 Torque across discontinuous shear thickening transition. Panels (a-d) show subsequent snapshots of the system as it undergoes instability at global shear rate _ γ ¼ 4:467 10 À6 . The color-coding corresponds to the total torque of each particle, in which the blue particles have τ ≤ −5 × 10 −5 , the red particles have τ ≥ 5 × 10 −5 , and the green particles have nearly 0 torque. As the system is sheared, domains of like-torque particles nucleate and grow, thus resulting in an enhancement in the rotational degrees of freedom, which underlies the well-known instability near the thickening transition. The mean torque as a function of strain is shown in panel (e), where the position of each snapshot as a function of strain is marked. Error bars are smaller than symbols. In panel (f), a snapshot of the system is shown above thickening for _ γ ¼ 1:122 10 À5 , where one can see localized clusters of opposite-torque particles. The binning of the color bar is similar to panels (a-d), and corresponds to total torque in the range τ ≤ 5 × 10 −5 and τ ≥ 5 × − 10 −5 . A close-up of such a cluster is given in panel (g).  5 Like-torque clusters and rheology across the discontinuous shear thickening transition. A snapshot of the system is depicted by two different measures: panel (a) total torque per particle and panel (b) total shear stress per particle. The global shear rate is equal to _ γ ¼ 4:467 10 À6 . Small clusters of very large positive (red) and negative (blue) torque particles can be seen; at the same time, the stress in the system is homogeneous. This phenomenon can be better seen in a close-up in panel (c). In panel (d), a configuration of the system with a strain difference of three particle diameter with respect to panel (a) is given; the global shear rate is not changed. In this panel, the instability has set in, and the like-torque clusters are grown. Panel (e) shows that the stress becomes highly heterogeneous. In panel (f), a close-up is given. Stress-bearing chains of large negative stress (blue) can be seen along the dilation direction, while large positive chains (red) elongate along the compression direction. The population of like-torque clusters and the stress-bearing chains are qualitatively correlated. In panels (g, h), we show that such a correlation does not exist above thickening for _ γ ¼ 1:122 10 À5 . Panel (i) The second moments of stress per particle (blue circles) and torque per particle (red squares) are depicted versus the strain. This panel shows that fluctuations in the stress and torque are correlated below thickening, at _ γ ¼ 4:467 10 À6 , supporting the hypothesis of a qualitative correlation between the population of like-torque clusters and the stress-bearing chains. By correlation, we mean that the stress and torque change proportionately as a function of strain. Such a correlation cannot be seen above thickening (panel (j)). Panels (i, f) show instantaneous averaging over particles and error bars are smaller than symbol size. rotation in an overdamped system will be essential for the understanding of intermittent fluctuations in these systems. In an elegant study, Vagberg et al. 48 have shown that Newtonian rheology arises as a result of clustering. Therefore, in a system with Newtonian rheology, clustering is expected to be much more pronounced than that in a system with Bagnold rheology.
We note that Chattoraj et al. 49,50 have recently demonstrated that a deformed very dense frictional system exhibits oscillatory instability, as a result of a pair of complex eigenvalues of the Hessian. Chattoraj et al. have also discussed a possible relation between the oscillatory amplifications in a frictional system with a long-standing problem in earthquake physics, remote triggering 51 . These results again emphasize the roles of frictional forces, and suggest that a more thorough investigation is warranted.
A fundamental question about particulate matter is whether a temperature-like quantity exists that describes fluctuations in these systems. A primary development was made by the so-called Edwards' entropy 52 , which gives rise to angoricity-a granular temperature. Several methods have been developed to test Edwards' ideas. Probably the most renowned method is the overlapping histogram method by Dean and Lefevre 53 , which has been used in many recent studies [54][55][56] . One large obstacle to applying Edwards' theory is that the density of states Ω(ϕ) and the corresponding partition function may be unknown for a given granular ensemble 57 . More modern treatments of this problem have been described by Makse and Kurchan 58 and Zheng et al. 59 . Our fluctuation analysis approach provides an alternative method to compute an effective stress-temperature that is similar to the rheology. We expect that this method will pave the way to better understanding of fluctuations in many different rheological phase transitions.

Methods
We use a linear dashpot spring to model both normal and tangential forces 60 . Two particles at positions r i and r j , with radii a i and a j , respectively, interact when they overlap, δ = |r i − r j | − (a i + a j ) < 0. A spring whose force is proportional to the overlap δ acts as a repulsive mechanism between two colliding particles. The interaction force along the normal direction is then given by: where k n and η n are the elastic and damping constants, respectively, and n ij is a unit vector along the line connecting the centers of two particles n ij = (r i − r j )/|r i − r j |.
With ω i and ω j , the angular velocities, the total tangential velocity at the contact point can be written as: Integrating the tangential velocity v ij,t from the initiation of contact to the current time gives the tangential overlap as ξ ¼ R t coll 0 jv ij;t jdt 0 . A spring proportional to ξ acts in the tangential direction along the contact plane to model the static friction where k t and η t are the spring and damping coefficients, and t ij is a unit vector along the contact plane, t ij × n ij = 0. A torque proportional to the tangential force acts on each particle. Accordingly, the total force is equal to: from which translational and rotational degrees of freedom are coupled in this model. We use a 50:50 bidisperse mixture of particles whose ratio of radii is 1.4. The diameter of small particles is chosen to be the unit of length. The mass is equal to the area of each particle. The spring constants are k n = 1 and k t = 0.5k n , and the damping coefficients are η n = η t = 1. Time is measured in units of the elastic timescale τ e ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi md=k e p , where d = 1 is the diameter of small particles, and m = πd 2 /4. The magnitude of the tangential force is bound by Coulomb's frictional law |f ij,t | ≤ μ|f ij,n |, where μ = 1 is the coefficient of friction.
We use Large-scale Atomic/Molecular Massively Parallel Simulator to integrate equations of motion of particles. This process is performed by using pair style "gran/hooke/history" to model the interactions plus Lees-Edwards boundary conditions by using "deform." Granular temperature T g is defined as the kinetic temperature of particles as: where i runs over particles, v i,x is the horizontal component of the velocity of particle i, and v(y) is the drift velocity at altitude y.
Shear stress is defined as: σ xy ¼ σ xy;n þ σ xy;t ð9Þ where the first and second terms on the rhs refer to the shear stress caused by normal forces f n and tangential forces f t between particles, respectively. Moreover, the normal and tangential parts are equal to: where A = wL is the area of the bin, subscript i runs over all particles inside a bin, and j runs over the neighbors of particle i. A visual depiction is given in Supplementary Fig. 1. In addition, δr ij is the distance between particles i and j with two components δx ij and δy ij along the x-and y-directions, respectively. Similarly, f t ij;x and f t ij;y are the components of f t ij . If both particles i and j are inside a given bin, then their contribution to the local stress is given by the above equation. If i and j belong to different bins, a factor of 1/2 of the above terms will contribute to the shear stress of their host bins. We also compute the stress per particle. In that case, the subscript j in the above equations runs over the neighbors of particle i. We emphasize that there are three different types of stress in this study: (1) local shear stress, (2) global shear stress, which is the sum over all local stresses, and (3) the stress per particle. For simplicity, we do not use subscripts to distinguish these three different cases, instead, we specifically mention them wherever applicable.
To compute the instantaneous local shear rate _ γ i at a given time t, we compute the average drift velocity v i,x inside a given bin i according to: where j runs over all particles inside bin i, v j,x is the horizontal component of the velocity of particle j along the streaming direction, and n i is the number of particles inside bin i. The local shear rate at bin i is then given by a mid-point derivative of the velocity profile: where w is the bin width. We refer to the local shear rate as _ γ l ; however, we do not use any subscript for the global shear rate _ γ.

Data availability
All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.