Uncovering instabilities in the spatiotemporal dynamics of a shear-thickening cornstarch suspension

Recent theories predict that discontinuous shear-thickening (DST) involves an instability, the nature of which remains elusive. Here, we explore unsteady dynamics in a dense cornstarch suspension by coupling long rheological measurements under constant shear stresses to ultrasound imaging. We demonstrate that unsteadiness in DST results from localized bands that travel along the vorticity direction with a specific signature on the global shear rate response. These propagating events coexist with quiescent phases for stresses slightly above DST onset, resulting in intermittent, turbulent-like dynamics. Deeper into DST, events proliferate, leading to simpler, Gaussian dynamics. We interpret our results in terms of unstable vorticity bands as inferred from recent model and numerical simulations.


I. INTRODUCTION
Instabilities are commonly observed in simple, Newtonian fluids when forced to flow under increasingly large Reynolds numbers. Such hydrodynamic instabilities ultimately lead to fully developed turbulence [1,2], yet following multiple pathways in which vortices [3,4] or turbulent puffs and spots [5,6] may arise and mediate unsteady or chaotic large-scale flow dynamics. While inertia is at the heart of flow instabilities in simple fluids, non-Newtonian fluids, such as polymer or self-assembled surfactant solutions, may display instabilities at vanishingly small Reynolds numbers due to elasticity [7,8] or due to a strong coupling between the flow and the fluid microstructure [9]. A typical example is provided by wormlike micellar solutions where shear-induced alignment associated with high viscoelasticity leads to shearbanded flows that transition to elastic turbulence at high Weissenberg numbers [10].
Shear-thickening, the process by which the viscosity of concentrated particulate dispersions dramatically increases above some critical load, is another widespread phenomenon that can be interpreted in terms of an underlying instability. Yet it still remains largely debated and far from fully understood. While the shear-induced growth of particle clusters, referred to as "hydroclusters", has long been invoked to explain shear-thickening, especially in Brownian suspensions of small colloidal particles [11,12], it was recently recognized that shear-thickening involves solid friction activated through shear-induced compressive stresses, at least for non-Brownian particles [13][14][15][16][17].
In particular, the Wyart and Cates model [15] provides a minimal framework where shear-thickening is described as a transition from a low-viscosity, lubricated assembly of nonfrictional particles to a high-viscosity, and possibly fully jammed, frictional contact network. Although supported both by simulations [18] and experiments [19,20], this model only provides a zero-dimensional picture that ignores the spatial and dynamical aspects of the transition. In particular, the prediction of hysteresis associated with S-shaped flow curves hints at an instability that should involve the separation of the system into shear bands oriented along the vorticity direction of the flow and referred to as vorticity bands in the literature [21]. In the specific case of hard non-Brownian particles, stress balance across the interface between bands prevents the existence of steady vorticity bands [20,22]. This could explain the large temporal fluctuations that are ubiquitously observed in shear-thickening systems [19,20,[23][24][25].
The exact nature of the instability and its origin remain however elusive. For instance, it is unclear whether particle migration, free surface instability, vorticity shear-banding and/or gradient shear-banding are at play [26][27][28] and whether the instability involves truly chaotic dynamics [20,29] as contained in early phenomenological models [30,31]. Therefore a full spatiotemporal picture of such unstable flows is required from experiments that rely on sufficiently long and statistically stationary measurements.
Here we focus on the dynamics of a non-Brownian dispersion of cornstarch particles that displays discontinuous shear-thickening (DST) in the hysteretic region of the phase diagram [15]. Dense starch suspensions are a popular system to investigate shear-thickening [32][33][34][35][36][37][38], even though they raise technical challenges linked to particle polydispersity and porosity, to sedimentation and potential migration under shear, or to possible adhesion between particles [26,39,40]. Using a rheometer in a concentric-cylinder geometry, we solicit the suspension at a constant shear stress and simultaneously image the local flow behavior at a "mesoscopic" scale of a few particle sizes with ultrasonic echography. We first provide a detailed statistical analysis of the shear rate fluctuations as a function of the distance to the DST on-)) ) )) ) )) ) )) ) )) ) set. These global dynamics show intermittent, turbulentlike statistics close to DST onset that progressively give way to Gaussian statistics far above DST. We then turn to the spatially and temporally resolved features of the flow to get deeper physical insights in such dynamical regimes. Our major findings are (i) that the bulk suspension remains homogeneously sheared whatever the applied stress in the DST regime and (ii) that the dynamics of the shear rate at short time scales display a diffusive behavior scattered with ballistic events which correspond to unstable macroscopic bands that propagate along the vorticity direction and proliferate as the stress is increased.

II. STATISTICAL ANALYSIS OF THE SHEAR RATE DYNAMICS
The mechanical response of a cornstarch dispersion at 41% wt. in a density-matched mixture of water and cesium chloride (Appendix A 1) is monitored under a cylindrical Taylor-Couette flow generated by a stresscontrolled rheometer (TA Instruments ARG2). A smallgap concentric-cylinder geometry, sketched in Fig. 1, was carefully designed to minimize solvent evaporation, iner-tia, particle sedimentation and migration, and to avoid instability of the free surface and bubbles trapped in the dispersion, allowing for measurements over at least 20 hours (Appendix A 2). Upon an increasing ramp of imposed shear stress, the flow curve, shear stress σ vs shear rateγ, of our cornstarch dispersion displays the hallmark of DST. As seen in Fig. 2a, the shear rate increases smoothly with σ up to a critical stress σ c 12 Pa at whichγ shows a sudden discontinuity. Above σ c ,γ remains on average constant and equal toγ c 2 s −1 i.e. the viscosity η = σ/γ c increases linearly with stress, which is a signature of DST [24]. As already reported in previous studies [19,20], the shear rate exhibits complex fluctuations over the whole vertical part of the flow curve for σ > σ c .
To gain insight into these dynamics, we impose various constant shear stresses ranging from 5 to 200 Pa during long steps of at least 2,500 s each. The power spectral density (PSD) of the shear rate is displayed in Fig. 2c. We may distinguish three different regimes depending on the imposed stress σ. Below σ c , fluctuations ofγ are negligible and the PSD simply reflects the sum of mechanical and electrical noises: it roughly decreases as 1/f at low frequencies before reaching noise level (black curve in Fig. 2c). The large peak at f 0.04 Hz corresponds to the rotation period of the inner cylinder, which always shows in the PSDs below σ c due to slight imperfections of the apparatus. For σ c < σ 55 Pa, the PSD shows a complex series of peaks at low frequencies followed by a strong decrease at higher frequencies (purple and blue curves in Fig. 2c). For even higher stresses σ 55 Pa, the peaks give way to a plateau with a cutoff frequency of about 0.5 Hz above which the spectrum decreases as f −4 . Thus, based on the shape of the PSDs, we hereafter identify a stable regime for σ < σ c , a first unstable and statistically complex DST regime for σ c < σ 55 Pa (regime I), and a second unstable yet statistically simpler DST regime for 55 σ 200 Pa (regime II). Note that we do not claim that regime II represents the ultimate stage of the DST transition as instability of the free surface above 200 Pa prevented higher stress values to be explored while ensuring homogeneous shearing conditions.
To further analyze the fine properties of the shear rate fluctuations in both regimes I and II, we now borrow statistical tools from turbulence physics [41,42] and focus our attention on the statistical properties of the increments of the shear rate δγ(t, ∆t) =γ(t + ∆t) −γ(t). For a given value of the time lag ∆t, we first introduce the probability density function (PDF) P(δγ, ∆t) of δγ(t, ∆t) taken as a random variable in time, see Fig. 3a. We then investigate the second, third and fourth moments of δγ, respectively through the variance δγ 2 , the skewness S and the kurtosis K displayed in Fig. 3b-d as a function of ∆t (see Appendix A 3 for definitions).
In regime II, P(δγ, ∆t) remains Gaussian for all time lags ∆t whereas in regime I, it shows strong asymmetry with exponential tails at short time scales and only assumes a Gaussian shape at long time scales (∆t 10 s in the specific case of σ = 16 Pa shown in Fig. 3). Such a striking difference between regimes I and II shows even more clearly on the moments of the increments. For both regimes I and II, δγ 2 saturates to twice the variance of the shear rate, which indicates statistical independence ofγ(t) andγ(t + ∆t) at large ∆t. Meanwhile S and K both converge to zero, confirming loss of correlation and Gaussian statistics for large time lags. At short time lags, however, the dynamics appear much less trivial and differ significantly. In regime II, the skewness and the kurtosis are featureless and always close to zero, consistently with Gaussian distributions. Still, the variance displays a ballistic growth, i.e. δγ 2 ∝ ∆t 2 , before plateauing to its asymptotic value for ∆t 1 s. In regime I, a similar initial ballistic-like behavior is followed by diffusive dynamics, i.e. δγ 2 ∝ ∆t, for 0.5 ∆t 5 s, coupled to a strong negative peak in the skewness and large decreasing values of the kurtosis. Together with exponential tails in the PDFs, this suggests the presence of intermittent events with large negative amplitude at short time scales in regime I.
To summarize, above σ c , the shear rate first displays non-Gaussian, intermittent fluctuations at short time scales, characterized by a cross-over between ballistic and diffusive dynamics as the time lag increases (regime I).
As the shear stress is increased deeper into DST, the fluctuations become fully Gaussian with only a ballistic behavior at short time scales (regime II). In order to unveil the origin of such complex dynamics, we turn to ultrasonic imaging which allows us to map the local tangential velocity v(r, z, t) of the cornstarch suspension as a function of the radial direction r across the gap, the vorticity direction z and time t [43] (Appendix A 4).

III. LOCAL INSIGHTS FROM ULTRASONIC IMAGING
A direct comparison between the shear rate response and local velocity data recorded just above DST onset (σ = 12 Pa) is shown over long periods in Figs. 13 and 14 in Appendix A 5 and over the two recurring spatiotemporal patterns observed throughout regime I in Fig. 4, see also Supplemental Video 1. A first important observation is that in spite of strong slippage at both walls, the dispersion always remains homogeneously sheared in the bulk i.e. we do not observe any solidlike, shear-jammed state. A more thorough analysis of flow profiles, detailed in Appendix A 6 and displayed in Fig. 14, shows that the effective shear rate in the bulk suspensionγ B (t) correlates very well withγ(t). Therefore the global shear rate most likely reflects the local dynamics at least close to DST onset.
A second, even more striking piece of information that can be extracted from velocity maps in regime I is the presence of propagating events along the vorticity direction that alternate with quiescent phases. Indeed, from the one-hour experiment at σ = 12 Pa, we could extract ∼20 sequences of duration ∼10 s where the fluctuations ofγ(t) and v(r, z, t) remain essentially flat (Fig. 4a) and a similar number of sequences where v(r, z, t) provides clear evidence for the propagation of a band along z that is delineated in time by peaks in both the global and effective shear rates (Fig. 4b). As seen in Fig. 13, such global peaks may be present in bothγ(t) and v(r, z, t) with no sign of propagating events in between. This is of course because ultrasound imaging only captures a vertical slice of the sample and can miss traveling bands that do not span the whole azimuthal direction. Therefore, we shall now restrict the statistical analysis proposed above in Section II to sequences where these two different patterns can be unambiguously identified.
Interestingly, during quiescent phases, the variance δγ 2 q of the shear rate increments exhibits a purely diffusive scaling with time lag ∆t (blue symbols in Fig. 4c), which implies that the suspension applies random kicks on the surface of the rotating cylinder. On the other hand, the variance δγ 2 p corresponding to sequences that include propagating events displays a short-time ballistic scaling up to ∆t 0.5 s (red symbols in Fig. 4c). This characteristic timescale corresponds to the duration of the acceleration phases during global peaks inγ that precede and follow propagating events. Ultrasound images further show that propagating events travel along the vorticity direction with a speed c 10 mm s −1 . They extend over roughly 10 mm in the z-direction while they span the whole gap in the r-direction, see also Supplemental Video 1. Such traveling bands are associated with a peculiar temporal signature in the local velocity as shown in Fig. 4d. Within the band, the suspension flows at a velocity that is about 20% smaller than average. The band is itself preceded by a front where the suspension moves about 10% faster than average. Remarkably, the velocity profiles always show homogeneous shear across the gap, see Fig. 4e. The amount of wall slip is observed to transiently increase (decrease resp.) at the rotating (fixed resp.) cylinder while the effective shear rate in the bulk remains roughly constant.
Finally, coming back to the global shear rate dynamics, we use the two patterns associated with quiescent phases and propagating events as a projection basis to reconstruct the statistics of the wholeγ(t) time series. As shown in Fig. 4c, a linear combination of δγ 2 q and δγ 2 p provides an excellent fit of the global variance δγ 2 . Thus, singling out those two independent patterns based on local velocity maps allows us to recover the full dynamics from only about 10% of the time series. The same analysis was repeated for stresses spanning both regimes I and II. Figure 5a shows that the prop- agation speed of traveling bands stays roughly constant throughout regime I but shows large variability around increasing mean values in regime II. Meanwhile, the occurrence frequency of propagating events weakly increases in regime I before jumping to large values that are comparable to the cutoff frequencies observed in the PSDs in regime II, see Fig. 5b. Strikingly, linear combinations of the variances δγ 2 q and δγ 2 p extracted at σ = 12 Pa remain very efficient to reconstruct δγ 2 for all imposed stresses, see Fig. 12 in Appendix A 3. This suggests that propagating bands interact very weakly with each other. As seen in Fig. 5c, quiescent (i.e. diffusive) patterns dominate the dynamics in regime I whereas propagating (i.e. ballistic) events proliferate and eventually solely account for the dynamics in regime II. The latter proliferation of propagating events is directly illustrated through ultrasonic imaging at σ = 80 Pa in Fig. 15 in Appendix A 5, see also Supplemental Video 2. This proliferation also goes along with a decorrelation of the global shear rate relative to local velocities, see Fig. 16. The fact that the statistics ofγ become Gaussian in regime II also supports a picture where propagating bands are independent random events.

IV. DISCUSSION AND INTERPRETATION IN TERMS OF UNSTABLE VORTICITY BANDS
Our results demonstate that the unsteady dynamics observed in the DST of a dense cornstarch suspension originates from localized, intermittent propagating events that proliferate as the shear stress is increased. While vertical (or S-shaped) sections of the flow curve are usually attributed to steady-state banding [21], we do not observe any shear localization involving solidlike arrested regions, unlike in experiments of Refs. [20,26,28]. This is most probably because these previous works explored the shear-jammed region of the flow phase diagram [15] either due to higher volume fractions [20,28] or due to imposing the shear rate in a wide-gap Taylor-Couette geometry which led to particle migration and flow separation [26]. Here, we observe that the flow always remains homogeneously sheared in the gradient direction. Moreover we do not observe any signature of strong local variations of the volume fraction. The curvature of the velocity profiles seen in Fig. 4e, indicative of apparent shear-thinning, could be attributed to a mild migration towards the fixed outer cylinder, consistently with recent concentration measurements through X-ray imaging [40]. Still, migration-induced local jamming cannot account for our experimental findings. (c) Coefficients α (triangles) and β (squares) of the best fit of δγ 2 by a linear combination α δγ 2 q +β δγ 2 p of the variances of the two patterns shown in Fig. 4c. In the fitting procedure, the variance of the increments is normalized by the variance of the shear rate so that the coefficients can be compared from one shear stress to the other.
Based on the above observations, we may thus interpret the present results in terms of a homogeneous DST system that essentially undergoes an unsteady yet uniform shear rate. In such a case, the only physical parameter that can be invoked to explain instability is the local stress through, e.g., the presence of vorticity bands or stress-bearing structures. Of course, our experimental setup does not provide access to local stresses in the bulk suspension but we shall argue below that the observed propagating events are fully consistent with traveling vorticity bands that were predicted and observed numerically very recently [22].
Theoretical arguments rule out the possibility of steady-state vorticity bands for hard non-Brownian particles [20]. Indeed, were vorticity bands stable, both the particle pressure and the solvent pressure would need to balance separately in order to stabilize the band interface. This would lead to identical shear stresses and thus to the absence of bands. Rather, DST flows were predicted to be always unsteady and a two-dimensional extension of the Wyart and Cates model was shown to yield oscillatory and chaotic-like dynamics in frictional granular particles [29] calling for three-dimensional approaches. Such three-dimensional simulations were performed recently by Chacko et al. [22] who also devised a  one-dimensional instability model for vorticity bands. Both the model and the simulations in Ref. [22] predict shear-thickened bands bearing a shear stress larger than average that travel along the vorticity direction at a given speed c. Simulations show that the band moves by about 800 particle radii per strain unit. With our cornstarch grains of mean diameter 15 µm (Appendix A 1) and a typical strain rate of 2 s −1 , this corresponds to c 12 mm s −1 in striking quantitative agreement with Fig. 5a. Although simulations still remain to be made fully comparable to experiments, in particular concerning boundary conditions, this agreement strongly suggests that our propagating events correspond to the traveling bands of Ref. [22].
Another prediction by Chacko et al. [22] is that band propagation goes along with small local variations of the vertical velocity v z and of the volume fraction φ. Although we cannot demonstrate the presence of non-zero vertical velocities due to the low resolution of our ultrasonic imaging setup in the z-direction, we may get some hints for local variations of the density by focusing on the local ultrasonic intensity I. Indeed, using the same setup on various non-Brownian suspensions in the dilute and semidilute regimes, we have shown that I directly relates to φ, yet in a complex manner for suspensions above a few volume percent [44]. While intensity maps remain featureless in regime I at σ = 12 Pa (see Fig. 13d), a convincing correlation between the tangential velocity and the local intensity is reported in Fig. 6 for σ = 80 Pa in regime II. There, propagating events are more pronounced and we observe that the intensity varies in the band concomitantly with the local velocity, see also Fig. 15d for more events. Therefore, our experiments are consistent with the simulations that display an increased volume fraction at the front of the vorticity bands followed by a small depletion within the band. More quantitative local measurements of the volume fraction, e.g. through time-resolved X-ray imaging, are in line to strengthen the present qualitative observations.
To push the discussion further, we note that the interpretation of propagating events in terms of unstable vorticity bands is also supported by the behaviors of the local velocity and wall slip shown in Fig. 4 during an event. If the volume fraction is slightly larger at the outer cylinder, as indicated by the curvature of the velocity profiles, then one expects the network contacts to exert a larger stress there in the case of sliding friction at the walls. Equivalently, in the case of lubricated friction, one expects the lubricating layer to be compressed and thinner at the fixed wall. This would result in a lower slip velocity at the stator during the propagation of a band bearing a higher local stress. Thus we may speculate that slippage of the suspension indirectly probes the local stress. We also note that the events unveiled in the present study clearly differ from large-amplitude stickslip oscillations that may couple to the rheometer inertia as reported in other works focusing on shear-jammed DST suspensions [25,45]. However, wall slip during DST of dense suspensions definitely requires more attention in future work especially to elucidate its microscopic origin and its possible interplay with the unstable bulk dynamics.
Finally, our results raise a number of interesting open issues. First, we observed that propagating events go along with sudden acceleration and deceleration of the rotating cylinder. These global peaks are likely to correspond to fast elastic waves generated during the nucleation of the bands or from their interactions with the bottom wall and/or the top surface. Conversely, the fact that these peaks are not systematically associated with a propagating event detected at the probe location suggests that the bands do not span the whole circumference of the cell. This points to the need to simulate extended yet bounded three-dimensional systems in order to investigate how the bands precisely nucleate and what sets their spatial extent. 3D simulations of dense suspensions were recently performed based on a fluid dynamics model for dilatant fluids where the inertia of the fluid plays a central role [27,46]. This approach accounts for the large and fast "shear thickening oscillation" observed in the wide-gap Taylor-Couette flow of potato-starch suspensions [47] and predicts localized shear-thickened bands that bear negative pressures but that do not seem to organize nor propagate along the vorticity direction [27]. The case of a small-gap Taylor-Couette flow with much smaller inertia and more homogeneous shear rate and volume fraction conditions still remains to be explored numerically as well as experimentally through local wall or particle pressure measurements [48][49][50].
Second, we successfully rationalized the complex DST dynamics of cornstarch based on two different spatiotemporal patterns extracted from a statistical analysis inspired by turbulence studies. This approach was only possible thanks to long experiments and to local insights into the flow. We do not however report any clear oscillatory dynamical regime such as the periodic and chaoticlike states in the experiments of Ref. [20] or the "locally oscillating bands" in the model and simulations of Ref. [22]. This could be due to differences in the location of the system within the shear-thickening phase diagram. Flow curves measured at different cornstarch weight fractions are shown in Fig. 17 in Appendix A 7 in order to facilitate the comparison with other results. We also note that the systems investigated in Refs. [20,22] are more confined than in the present work. The possibility of truly chaotic dynamical regions in the DST phase diagram should thus be ascertained together with their dependence on geometry. To stimulate such a future study, we show in Fig. 18 in Appendix A 8 that fluctuations similar to those reported above are also observed in two other geometries. Yet, it is highly probable that the nature of the instability and its spatiotemporal signature depend on a number of physical degrees of freedom including shearing geometry, sample size and boundary conditions. Third, our statistical approach has revealed a specific range of timescales where the dynamics is diffusive-like in regime I. Such a separation of timescales is reminiscent of multiscale processes, e.g. hydrodynamic turbulence [41] or imbibition processes [51]. In the present case of shear-thickening flows, however, it is still unclear whether some kind of self-similarity is involved. In particular, the traveling bands unveiled here could represent an analogue of large-scale coherent structures in turbulent flows but the available statistics does not allow us to tell whether they present some hierarchy of scales. Moreover, although the short timescale was related to the duration of global peaks in the mechanical response, the nucleation mechanism of the bands and what exactly sets this timescale are unknown. Along the same lines, the timescale above which the statistics become fully Gaussian is probably linked to the lifetime of the bands, to their propagation speed and/or to their occurrence frequency but the processes that select such characteristics are still to be discovered.

V. CONCLUSION
By carefully investigating a dense cornstarch suspension with a combination of rheometry and ultrasound imaging over long time scales, we have pinned the origin of unsteady dynamics in DST to the existence of transient localized bands that travel along the vorticity direction. The global dynamics can be decomposed into ballistic phases that correspond to such traveling bands and diffusive phases that correspond to a background of random kicks applied to the moving tool. As the shear stress is increased in the DST regimes, propagating bands progressively dominate over diffusive phases, leading to Gaussian fluctuations. These key observations urge to introduce stress spatiotemporal correlations in models and to simulate dense three-dimensional systems together with the shearing boundaries. In particular, the striking similarities between our experiments and recent model and simulations strongly support an explanation of shear-thickening dynamics based on localized vorticity bands bearing high shear stresses. Further experimental work on DST suspensions should now focus on the measurements of bulk stresses to directly confirm this picture as well as more microscopic observations including sliding of the material at the walls.
Our results concretize recent scientific effort to unravel shear-thickening in dense suspensions and to understand its control through the addition of plasticizers [25,52,53], nanofibers [54] or large particles [55] and the application of magnetic fields [56] or transverse shear [57]. Our work may also shed new light on undesirable effects in industrial processes such as stress surges [23,58], granulation and slurry fracture [59,60] and help to further de-velop promising applications such as fabrics impregnated with shear-thickening fluids for improved resistance to impact [61][62][63][64]. Dense cornstarch suspensions are obtained by dispersing 41% wt. of cornstarch (Sigma-Aldrich, CAS 9005-25-8, S4126-2KG) at room temperature in a density-matched solvent composed of 46% wt. water and 54% wt. cesium chloride (CsCl) [38]. Figures 7 and 8 respectively show density measurements of water-cesium chloride mixtures and the cornstarch size distribution inferred from optical microscopy. Cornstarch particles have a mean diameter a = 15 µm and the standard deviation of the diameter distribution is 7 µm. Cornstarch is progressively added to the solvent and lumps of starch are broken using a spatula. The suspension is then centrifuged at 200 g for two minutes before being poured into the Taylor-Couette cell. Following Ref. [38], the weight fraction used in the present work corresponds to a volume concentration of about 47%. To remove trapped air bubbles that may affect the rheology of the cornstarch suspension, preshear is applied for at least 5 h under a low shear stress σ = 5 Pa. As seen in Fig. 9, the suspension initially scatters ultrasound very strongly due to air bubbles. Shearing the suspension for several hours in the shear-thinning regime below shear-thickening allows us to progressively remove air bubbles, as confirmed by the slow decrease in the backscattered ultrasound intensity. The speckle signal obtained after ∼ 5 h is homogeneous and corresponds to backscattering by the cornstarch grains with a negligible contribution from air bubbles (if any).    Figure 10 shows a picture of our experimental setup. Rheological measurements are carried out in a concentriccylinder (or Taylor-Couette) geometry driven by a stress-controlled rheometer (TA Instruments ARG2). The fixed outer cylinder (stator) has an inner radius R o = 25 mm and the rotating inner cylinder (rotor) an outer radius R i = 23 mm, leaving a gap e = R o − R i = 2 mm between the two cylinders. In such a small-gap Taylor-Couette geometry, the curvature induces a relative decrease of the shear stress by δσ/σ = (R o /R i ) 2 − 1 = 0.18 from the inner cylinder to the outer cylinder. Let us emphasize that the stress heterogeneity inherent to our cell is twice smaller than in the Taylor-Couette geometry used by Hermes et al. [20], for which δσ/σ = 0.36, and ten times smaller than in the wide-gap geometry of Fall et al. [26] where δσ/σ = 1.8, i.e. where the stress is almost three times smaller at the stator than at the rotor. Shear is even more heterogeneous in a parallel-plate geometry where the shear rate goes from zero on the axis of rotation to its maximum value at the periphery [20]. Therefore, our small-gap Taylor-Couette geometry ensures a good homogeneity of the shear field, which helps a lot in mitigating particle migration due to shear gradients.

Experimental setup
The gap width e 130 a is still large compared to the mean particle size such that no significant effect of confinement is expected. Moreover, our cell has a large aspect ratio thanks to its height H = 63 mm 30e 4,000 a . This large aspect ratio corresponds to a small ratio of the free surface of the suspension to the rotating cylinder surface over which the total torque measured by the rheometer is integrated, S free /S tool e/H 0.03. For typical cone-andplate or parallel plate devices with radius 20 mm and gap width 1 mm at the periphery, the ratio is S free /S tool 0.1, which makes those geometries more sensitive to instabilities of the free surface.
Both cylinders are made of smooth Delrin (polyoxymethylene). This material was chosen because it leads to limited slippage compared to, e.g., smooth PMMA surfaces. However, as discussed in the main text, wall slip has a deep connection to the DST bulk dynamics and certainly deserves more attention in itself. In particular, rather than a mere artifact that has to be eliminated, it should be treated as a complex yet interesting physical phenomenon that may carry key information on the system under study. Moreover the inertia due to the geometry and to the rheometer can be neglected in all experiments. Indeed our setup has a total moment of inertia I tot = 54 µN ms 2 and the maximum strain acceleration was measured to beγ max 8 s −2 so that the corresponding stress is at most σ i,max = eIγ max /(2πHR 3 i ) 0.18 Pa, always below 2% of the imposed shear stress. The fluid has a moment of inertia I fluid = πρH(R 4 o − R 4 i )/2 18 µN ms 2 and brings an even smaller contribution to inertial stresses. The Reynolds number, Re = ρe 2γ /η, never exceeds 10 −2 . We may therefore neglect inertia in our stress measurements.
Finally, the stator is closed by a lid as sketched in Fig. 1. A groove machined in the upper surface of the rotor is filled with water to act as a solvent trap and minimize evaporation. This allows us to perform reproducible experiments on a time span of more than 20 hours with the same loading of the cell. We used the ARG2 Auxiliary Sample utility to retrieve the applied stress and the shear rate response as a function of time with a sampling frequency of 500 Hz.   Fig. 2).

Signal analysis
Long experiments under a constant stress in the DST regime yield statistically stationary shear rate responsesγ(t). Figure 11 shows a selection of such responses on more than 2,000 s. To characterize the temporal fluctuations of the shear rate, we use the following tools borrowed from studies of turbulent signals [41,42].
• The probability distribution function (PDF) of the increments δγ(t, ∆t) =γ(t + ∆t) −γ(t) is the probability P(δγ, ∆t) to find an event of amplitude δγ in the increment time series for a given ∆t. By definition, for any given value of ∆t, the mean of δγ(t, ∆t) is zero.
• The second moment of the PDF is the variance δγ 2 (∆t) = δγ(t, ∆t) 2 t , where the brackets · · · t denote the average over time t. As explained in the main text, whatever the applied stress in DST, the variance δγ 2 can be fitted by linear combinations of the variances of two elementary processes δγ 2 q and δγ 2 p extracted just above DST onset. Such fits are shown in Fig. 12.
• Related to the third moment of the PDF, we introduce the skewness defined as S = 0 indicates that the PDF is symmetric. S > 0 (S < 0 resp.) means that the distribution is concentrated on the right (left resp.) side of δγ .
• Finally, the fourth moment of the PDF is used to define the logarithm of the normalized kurtosis: The normalization is such that K = 0 for a Gaussian distribution. K > 0 (K < 0 resp.) indicates that outlier events are more (less resp.) probable than in a Gaussian distribution. FIG. 12. Variance δγ 2 (∆t) of the shear rate for various imposed stresses (colored symbols) together with the fits by linear combinations of δγ 2 q and δγ 2 p as detailed in the main text (colored lines). The coefficients α and β used for the fits are shown in Fig. 5c. For clarity, successive data sets were shifted vertically by 2 starting from σ = 13 Pa.

Ultrasound imaging
Rheological measurements are synchronized with ultrafast ultrasonic imaging [43]. Our echography technique relies on a custom-made high-frequency scanner driving an array of 128 piezoelectric transducers that sends short plane ultrasonic pulses with 15 MHz center frequency across the gap of the Taylor-Couette cell. As displayed in Fig. 10, the ultrasonic probe is immersed in a large water tank for acoustic impedance matching and is set vertically at about 25 mm from the stator. The water tank is connected to a circulation bath (Huber Ministat 125) that regulates the temperature to T = 25 ± 0.1 • C. The full specifications of this rheo-ultrasound setup can be found in Ref. [43]. During its propagation through the suspension, a plane ultrasonic pulse gets scattered by the cornstarch particles. The backscattered pressure signals constitute an ultrasonic "speckle" that is recorded by the transducer array, sampled at 160 MHz and further post-processed into images of the echogeneous structure.
Cross-correlating successive images yields the tangential velocity of the sample v(r, z, t) at time t and position (r, z) in cylindrical coordinates, r denoting the distance to the stator and z the position along the vertical direction pointing downwards with the origin taken at transducer #1 (see Fig. 10b). The ultrasonic images cover 32 mm in height, i.e. about half of the Taylor-Couette cell in the vorticity direction. Moreover, the characteristics of the ultrasonic beam corresponds to an azimuthal span of 300 µm. The spatial resolution along the z-direction is given by the spacing of 250 µm between two adjacent transducers and the resolution in the radial direction r is 75 µm. The frame rate for ultrasonic images is fixed to 50 fps in all the experiments presented here.
Velocity measurements from ultrasound speckle images rely on the phase of the backscattered pressure signal. We can also measure the amplitude P (r, z, t) of the backscattered signal through a Hilbert transform. As shown in Ref. [44] for a number of suspensions of non-Brownian particles with diameters ranging from 20 to 80 µm, P (r, z, t) is a monotonic function of the local volume fraction φ(r, z, t). Here, we shall only use P (r, z, t) as an indication of local variations of φ(r, z, t) since an absolute determination of the local volume fraction requires careful calibration and strong additional assumptions on scattering and attenuation by the suspension [44]. 5. Spatiotemporal diagrams from v(r, z, t) and P (r, z, t) In order to put emphasis on the dynamics of the local tangential velocity either along the vorticity direction z or along the gradient direction r, spatiotemporal diagrams are computed by averaging v(r, z, t) along the r-direction or the z-direction, which yields respectively the v(z, t) and v(r, t) maps shown in Figs. 4a,b, 6a, 13b,e and 15b,e. Similarly, we present spatiotemporal diagrams I(z, t) of the local ultrasonic intensity in Figs. 6b, 13d and 15d. There, in order to remove small yet systematic discrepancies in the intensity response of the different transducers, the intensity I(z, t) is computed from P (z, t), the r-average of P (r, z, t), as I(z, t) = (P (z, t) − P (z, t) t )/σ P (z) , where σ P (z) denotes the standard deviation of P (z, t) taken over t.
Typical velocity and intensity spatiotemporal diagrams are displayed over the course of about 1 min in Fig. 13 for σ = 12 Pa (regime I) and in Fig. 15 for σ = 80 Pa (regime II). For the lower shear stress, we observe that the propagating event detected at t 575 s is preceded and followed by global peaks in the shear rateγ, see Fig. 13a. These peaks also show as vertical lines along the z and r axes in the spatiotemporal maps of Fig. 13b,e. Based on the presence of various other peaks inγ, we infer that during this specific time sequence, a few other propagating events take place yet outside the ultrasound region of interest. Global peaks disappear from the spatiotemporal representations when the local velocity is normalized by the rotor velocity v Rheo (t), see Fig. 13c,f. This confirms that the propagating events are local while the sudden acceleration phases before and after the events are global.
For the larger shear stress, the local dynamics are decorrelated from the global ones. Yet, locally, propagating events are still delineated in time by maxima of the velocity that occur simultaneously for all z and r. The fact that normalizing the local velocity by the rotor velocity does not significantly affect the spatiotemporal diagrams in this case (compare Fig. 15b,e with Fig. 15c,f) suggests that many different propagating events take place independently at the same time throughout the sample and average out in the global signals v Rheo (t) andγ(t).

Analysis of local velocity profiles
Local velocity profiles show significant wall slip both at the rotor and at the stator. In order to compare the dynamics of the local velocity with the global dynamics inferred from the shear rateγ(t) measured by the rheometer, we consider the z-averaged velocity profiles v(r, t) and estimate the sample velocities v S (t) and v R (t) respectively at the stator and at the rotor through linear extrapolations of v(r, t) (see red lines in Figs. 14a and 16a). Slip velocities are then simply given by v SS (t) = v S (t) at the stator and v SR (t) = v Rheo (t) − v R (t) at the rotor, where v Rheo (t) is the velocity of the rotor at time t. v Rheo is linked to the global shear rate through v where the last approximation results from the small-gap limit e R i . The "bulk" velocity profile is then obtained by subtracting the slip velocity at the stator to v(r, t): v B (r, t) = v(r, t)−v SS (t). Taking the average over r, we get v B (t) from which we define the effective shear rate felt by the bulk material asγ Figures 14 and 16 show the results of the analysis of the velocity profiles corresponding to the time sequences of Figs. 13 and 15. In regime I, we observe a very strong positive correlation between v Rheo , v SS and v B and a fairly positive correlation between v Rheo and v SR , see Fig. 14c. In regime II, however, the various velocities show very different dynamics, see Fig. 16b,c. The fact that v SS + v SR + 2 v B does not exactly match v Rheo can be ascribed to the downward curvature of the velocity profiles that leads to an average value v B smaller than (v R − v S )/2.

Flow curves for various concentrations
It is well established that the rheological properties of cornstarch not only depend on its concentration but also on its polydispersity, on its porosity, on the ambient humidity and on the supplier [24,38,39]. Moreover different authors have used different suspending fluids, e.g. pure water [26], water-glycerol mixtures [20] or water-CsCl mixtures as in the present work [33,38,39]. Therefore comparisons with previous or future studies require to locate our sample within the shear-thickening phase diagram of cornstarch. Although a systematic investigation of the influence of the concentration on the effects studied here is out of the scope of the present paper, Fig. 17 shows flow curves measured on our cornstarch system at different weight fractions in a parallel-plate geometry. The transition from continuous to discontinuous shear thickening occurs at about 38% wt., above which the flow curves clearly show a sigmoidal part. The corresponding volume fraction has been noted φ c or φ DST in previous works [15,20,22,26]. Further tests at 43% wt. reveal strong surface instability at low shear stress so that such a large weight fraction most probably corresponds to a volume fraction above the jamming point for frictional particles, noted φ m or φ µ J in Refs. [15,20,22]. We conclude that the weight fraction of 41% wt. studied in details in the present paper lies close to full jamming but still belongs to the DST region of the phase diagram.

DST in other geometries
In order to test for the robustness of the dynamics observed in our smooth Taylor-Couette geometry, we conduct additional experiments in a rough Taylor-Couette geometry and in a rough parallel-plate geometry, see Fig. 18. The rough Taylor-Couette geometry uses a rotor of radius R i = 23.9 mm and height H = 61 mm covered with sandpaper (P-320 grade, grain size 46 µm) and a stator of radius R o = 25.0 mm made of sandblasted PMMA (typical roughness of a few micrometers). The gap width of this Taylor-Couette cell is thus e = 1.1 mm. Strong scattering of ultrasound due to the roughness of the stator makes it impossible to perform ultrasonic imaging simultaneously to rheometry with this Taylor-Couette cell. Consequently we could not check for wall slip or propagating events in this case. The parallel-plate device consists of two plates of diameter 29 mm covered with the same sandpaper as the rough rotor and separated by a gap of 1 mm.
The data in Fig. 18 reveal the same phenomenology as that found with the smooth Taylor-Couette geometry of gap 2 mm. Both flow curves show the distinctive vertical portion characteristic of DST. Within the DST regime, the shear rate measured under an imposed stress becomes unsteady. The critical shear rate and shear stress at DST onset inferred from stress sweeps are slightly less than those in Fig. 2a, which could indicate some degree of sensitivity to the geometry and/or sample preparation. The fact that, in the parallel plates, long experiments under constant stress record shear rates that are increasingly larger than in the preceding, rather fast stress ramp is most probably due to humidity absorption by the suspension. As suggested in Ref. [20], particle migration is also more likely to occur in the parallel-plate geometry. In any case, the shear rateγ(t) always slowly drifts to larger values such that the shear rate response in parallel plates is not strictly statistically stationary. This precludes a detailed statistical analysis of long signals as described in the main text. Still, the fluctuations observed in Fig. 18 are qualitatively similar to those of Fig. 2b. Therefore, they are likely to be due to the same local dynamical phenomena, although this remain to be fully investigated and demonstrated from local measurements.