Metastability at the Yield-Stress Transition in Soft Glasses

We study the solid-to-liquid transition in a two-dimensional fully periodic soft-glassy model with an imposed spatially heterogeneous stress. The model we consider consists of droplets of a dispersed phase jammed together in a continuous phase. When the peak value of the stress gets close to the yield stress of the material, we find that the whole system intermittently tunnels to a metastable"fluidized"state, which relaxes back to a metastable"solid"state by means of an elastic-wave dissipation. This macroscopic scenario is studied through the microscopic displacement field of the droplets, whose time statistics displays a remarkable bimodality. Metastability is rooted in the existence, in a given stress range, of two distinct stable rheological branches as well as long-range correlations (e.g., large dynamic heterogeneity) developed in the system. Finally, we show that a similar behavior holds for a pressure-driven flow, thus suggesting possible experimental tests.


I. INTRODUCTION
Soft amorphous materials-such as emulsions, microgels, foams and colloidal suspensions-display a solidto-liquid transition for sufficiently large values of an external forcing: they are solid at rest and able to store energy via elastic deformations, whereas they flow whenever the stress is above a critical threshold known as the yield stress [1]. The complex spatio-temporal behavior shown by soft-glasses at the yield-stress transition has been the subject of intense scrutiny in the recent years [2][3][4][5]. Some materials, often denoted as "simple" yield-stress fluids [5] (e.g. microgels [6], nonadhesive emulsions [7,8]) exhibit yielding properties which are rather homogeneous in space: For any imposed shear rate, even a small one, there is always a stress at which these materials can fluidify homogeneously; the steady flow dynamics is also typically preceded by a nontrivial transient behavior [9,10]. In other materials with thixotropic properties [5], like adhesive emulsions [11], a specific kind of heterogeneous flow can be steadily established: If an imposed shear rate is smaller than a given threshold, the system may decompose in two distinct spatial regions, showing a solid and fluidized behavior respectively. By changing the shear rate value, the widths of the two regions are changed, whereas the shear stress remains constant. This phenomenon is known as shear banding [12][13][14][15][16][17][18]. Here, the term "shear banding" as a form of heterogeneous flow characterized by shear localization independently of any stress heterogeneity [5]. This differs from the shear lo-calization induced by stress heterogeneity, where part of the material is above yield and part below; it also differs from the shear localization emerging in presence of slippage at the walls.
From the theoretical point of view, different phenomenological models have been proposed to capture the fundamental physics underlying soft-glasses behaviors. In some cases [such as the soft-glassy-rheology (SGR) model [19][20][21] or shear-transformation-zone (STZ) theory [22]] the notion of "effective temperature" provides a useful way to describe the onset of the plastic flow in soft glasses. Such "temperature" is actually thought of as a quantification of the mechanical noise induced by the flow itself [19][20][21] and triggers activated hopping through the energy landscape of the system. Moreover, it has been clearly demonstrated both experimentally [23][24][25] and numerically [26] that soft glasses exhibit a nontrivial size dependence. This may give rise to "nonlocal" rheological effects [7,8] parameterized by a cooperativity length [6][7][8]27] estimating the typical size of the region involved in plastic rearrangements of the constituents following local elastic deformations. A recent proposal [18] has also linked cooperativity effects and nonlocal rheology to the emergence of shear-banding configurations. From a more general perspective, the shear-banding phenomenon has often been interpreted as the signature of a dynamic transition with a "phase coexistence" of two distinct states in space [27][28][29]: a jammed solid state and a fluidized state. A common explanation is to assume an underlying nonmonotonous rheological curve relating the stress to the shear rate [12,15,16], with two stable branches separated by an unstable branch. This nonmonotonicity has also been linked to the competition between different timescales related to different physical processes [29][30][31][32] (e.g. aging vs. flow-induced rejuvenation in Ref. [29] or restructuring time vs. stress-release time in Ref. [31]). When the minimum of the rheological curve occurs at very small shear rates, one can draw a "simple" picture of coexisting branches [5]: a solid branch described by zero shear (S = 0) and stress σ in the interval [0, σ st ], where σ st is referred to as the static yield stress; and a fluidized branch characterized by an Herschel-Bulkley (HB) relation of the type σ = σ Y + AS n [33], with σ Y < σ st denoted as the dynamic yield stress. For stress values σ ∈ [σ Y , σ st ] the shear rate is multivalued, hence the phase coexistence in space. For shear rate S greater than the critical shear S c = ((σ st − σ Y )/A) 1/n , the rheology of the system is described uniquely by the HB relation and no shear-banding is observed. This scenario has been explored and discussed in glassy models and numerical simulations [28,[34][35][36][37].
In this paper we want to look at the statistical properties of the yield-stress transition when σ Y < σ st from a different point of view. Permanent shear bands are often observed by applying an external velocity difference, say ∆U on a system of size L [18]. For ∆U/L < S c the system shows a homogeneous stress in space and splits into two shearing regions (a solid and a fluidized band) which permanently persist in time. Now, let us consider the same system under an imposed space-dependent stress ranging, say, from 0 to some value σ p close to σ st . In this case, we have two solutions linked to the two possible branches. If the rate of plastic rearrangements is large enough, the system can perform activated processes and transitions in time between the two solutions may be observed. In other words, for a relatively narrow range of values of the imposed shear stress peak σ p , one should be able to observe a clear bimodality in the probability distribution of a global rheological variable, like the space-averaged velocity, or some other convenient observable. Hence, we expect a time bimodality because of the repeated (back-and-forth) transitions between two different states which are unimodal in space. Such transitions are expected to be enhanced by the choice of a heterogeneous stress field which reduces the extent of the spatial region in which transitions take place. Based on numerical simulations of a soft-glassy model [38][39][40][41][42][43] (see Sec. II) we aim at providing a clear evidence that the above scenario holds.
In Sec. III we will analyze the rheological response at "large scales" and analyze the signatures of bimodality in the time evolution of the flow; then, in Sec. IV we enrich these observations with a comprehensive analysis of the rheological response at "small scales", i.e., by studying the statistical properties of the displacement field of the microstructural constituents. When bimodality is observed, we also observe that the overlap-overlap correlation length (see Sec. V) becomes of the same order of the system size. We argue that a long-range correlation function among plastic events is necessary in order to observe transitions in time from one state to the other. Preliminary investigations for a pressure-driven flow (see Sec. VI) will also support the same scenario, thus suggesting an experimental setup that could be used to test the predictions of numerical simulations. Some conclud-ing remarks will be given in Sec. VII. We believe that our results open a new perspective in the phenomenology of shear-banding in soft glasses.

II. MODEL
We simulated a soft-glassy model by means of a lattice Boltzmann (LB) equation which allows the simulations of droplets of one component dispersed in another component [38][39][40][41][42][43][44]. Droplets are stabilized against coalescence (see Fig. 1) by the combined effect of attractive and repulsive interactions [44]. In previous publications we showed that the model displays many of the wellknown properties observed for soft glasses. Importantly, for shear-controlled experiments (i.e., in Couette geometry) it behaves as a non-Newtonian fluid displaying a dynamic yield stress σ Y [38], a nonlinear HB rheology with HB exponent n 0.5 [41], elastic shear waves and plastic rearrangements [45]. For values of the stress σ larger than σ Y , our model shows quantitative agreement with nonlocal rheology theories [39,42] which have been used to rationalize the flow of concentrated emulsions [7,8,27] and other yield-stress fluids [6,42,46] in confined geometries. The values of the cooperativity scale ξ extracted from the model [41,47,48] are in agreement with experimental observations [7,49]. Recently, the model has been used in synergy with experiments on real emulsions in order to quantify the impact of the fluidization induced by the roughness of microchannels on the flow behavior of the emulsion [47,48]. As shown in the reminder of the paper, for this "model emulsion" the static and dynamic yield-stress values are found to differ. Looking at the literature on real emulsions, we know that "pure" emulsions do not show this behavior, whereas loaded "attractive" emulsions actually do [37,50,51]. Hence, in terms of the yielding properties, our model bears similarities with the behavior of an "attractive" emulsion with σ Y < σ st . Hereafter, we present all our numerical results by rescaling the LB units in such a way that the flowing rheological branch for a Couette forcing is given by σ/σ Y = 1+S 1/2 , where S is the shear. The system we consider is two dimensional, with x and y being the streamwise and spanwise coordinates respectively. We will study the rheology of our model by imposing a space-dependent stress. For this purpose, we consider fully periodic boundary conditions with a space-dependent forcing imposing the xy component of the stress (Kolmogorov flow): where L is the system size which has the same value in both directions and σ p is the peak value for the stress (see Fig. 1). A very similar setting has been used in previous experimental [26] and numerical [52,53] works. The choice of a fully periodic setup is initially taken in order to avoid possible wall effects and dependence on boundary conditions, which may alter the rheological response of the system [49,[54][55][56][57]. Later, in section VI we will discuss some preliminary simulations for a pressure-driven flow. In the fully periodic setup, for a Newtonian fluid with constant viscosity η, the streamwise component of the stationary velocity field induced by the stress would read where the peak value for the Newtonian velocity profile u N 0 = L 2πη σ p is a constant. In the model, an important control parameter is the quantity R = 2δ √ N /L where δ is the average thickness of the continuous phase, N is the number of droplets and L is the system size. Such a quantity is a measure of the ratio between the interface area and the area occupied by the droplets. Note that 1 − R should be considered proportional to the packing fraction in our system. The numerical simulations for the Kolmogorov flow have been performed with L = 1024 grid points, N = 512 droplets and R = 0.09 which implies a packing fraction well above the jamming point.

III. RHEOLOGICAL RESPONSE AT "LARGE SCALES"
The simplest way to measure the rheology in our system is to compute the characteristic shear S as a function of σ p . The value of S is computed using the average streamwise velocity profile u x (y, t) = L −1 x u x (x, y, t) at time t and performing its projection onto the viscous profile in Eq. (2) From u s (t) we compute s(t) = 2πu s (t)/L, whose time average provides the value of the shear S. In the top panel of Fig. 2 we show the rheological curve obtained in our system: Starting from σ p /σ Y 1 we perform a series of numerical simulations (red bullets) by increasing stepwise the peak stress σ p (i.e. "ramp-up" protocol). At relatively large values of the forcing (namely for σ p /σ Y 2.2) the system is completely fluidized. Next we reduce the forcing (i.e. "ramp-down" protocol) using exactly the same values σ p /σ Y of the "ramping-up" simulations: As observed in the top panel of Fig. 2 a clear (although small) hysteresis loop is observed. In the top panel of Fig. 2 the black continuous line refers to the same quantities for a simple HB fluid whose parameters are the same as those observed for our model in a Couette geometry [41] while the blue connected crosses refer to the same HB fluid supplemented with cooperativity effects, obtained using the steady nonlocal fluidity model [7,8,27]. Finally, in the lower panel of Fig. 2 we show the average velocity profiles u x (y) observed at σ p /σ Y 1.5 for the "ramp-up" simulation (purple squares) and "ramp-down" simulation (green triangles). Velocity profiles u x (y) are obtained from an average in time of u x (y, t). The results shown in Fig. 2 clearly demonstrate the existence in our system of two rheological branches with a dynamical yield stress smaller than the static one. Moreover, looking at the top panel of Fig. 2 we can immediately observe that the yielding point is above the yielding threshold evaluated in homogeneous conditions, i.e. σ p /σ Y 1.4. Qualitatively, we can argue that this is a consequence of the nonlocality in the  . Rheological data extracted from simulations are compared to the results obtained from the fluidity model [7,8,27] (blue crosses) and the Herschel-Bulkley (HB) [33] fit in a Couette geometry (black continuous line). Neither models can describe the solid branch and the transition to the plastic one; however, the fluidity model better describes the flowing regime. The purple empty square and the green empty triangle signal the position on the rheological curve of the corresponding profiles shown in the bottom panel. Bottom panel: Velocity profiles for the Kolmogorov flow at fixed peakstress value σp/σY 1.5 obtained from two different protocols: ramping up from lower peak-stress values (purple empty squares) and ramping down from larger peak-stress values (green empty triangles). This is clear evidence of the existence of two stable rheological branches signaling that the static yielding threshold is above the dynamic one [28]. flow coupled to the heterogeneity of the stress. Indeed, for the flow to occur, the peak stress needs to be above σ Y in a spatial region of the order of the cooperativity length [40,53]. The net effect of this is to increase the yielding threshold. However, a closer quantitative inspection reveals that the nonlocal model works very well only when the peak stress is well above the yield stress, while it fails to describe the transition point for σ p /σ Y 1.4. We indeed observe an abrupt transition in the rheological response that neither the simple HB model nor the stationary nonlocal fluidity model are able to capture. This contrasts with previous observations in yield-stress flu-ids subject to heterogeneous stress distribution [53]. We are therefore interested in investigating the nature and properties of this transition.
To get an intuitive picture on the system behavior at the transition, we show in Fig. 3 the time behavior of u s (t) for three different values of σ p . All of the following simulations have been performed using the ramp-up protocol unless explicitly stated otherwise (see Section VI). For relatively small σ p (top panel) the system intermittently tries to flow with an average value of u s close to zero; at large σ p (lower panel) the system is fluidized, and the signal corresponds to a plastic flow, as expected. The interesting point is the behavior of the system at σ p /σ Y 1.4 (middle panel): The system persists for relatively long time in a fluidized state and then goes back in a "solid" state. We also notice in the upper and middle panels of Fig. 3 strong periodic oscillations of u s (t). These oscillations are due to elastic waves generated in the system. The signal shown in the upper panel recalls the "stick-slip" behavior observed near the yield-stress transition in shear controlled systems [28]: since we impose the stress, the shear (or the velocity) shows intermittent bursts of activity. It is much less immediate, however, to understand the physics behind the behavior  . For small forcing, σp/σY 0.8, the system responds elastically and dissipates mainly through elastic waves visible from the periodic oscillations around 0. At σp/σY 1.4, the systems is at the middle point between the two branches (see Fig. 2) and it intermittently switches between an elastic response and a plastic flowing regime for which us(t) > 0. At σp/σY 1.6 the system is plastically flowing [7,8,27]. of u s (t) shown at σ p /σ Y 1.4. Since the intermittency in u s (t) is due to plastic rearrangements occurring in the system, it is important to inspect the system behavior at the scales of the microstructural constituents in order to get a deeper insight about the nature of the observed transition.
IV. RHEOLOGICAL RESPONSE AT "SMALL SCALES" Plastic rearrangements are localized topological changes in the droplets configurations. In our system, we can identify plastic rearrangements, corresponding to topological changes in the Voronoi tessellation of the centers of mass, by using its dual Delaunay triangulation (see Fig. 1): A plastic event happens whenever a link in the triangulation flips [59]. Next, we need to measure the droplet displacement during plastic rearrangements and try to understand whether this measure can be correlated to the observations discussed in Fig. 3. For this purpose, we start by looking at the displacement ∆ i (t) of the droplets defined as where x i (t) is the position of the center of mass of the ith droplet at time t and δt is a given time interval which in our simulations is set to be δt = 100 simulation time steps. This choice corresponds roughly to δt = t drop /10, where t drop = η R /γ is the droplet time, with R the average radius and γ the surface tension. As expected, | ∆ i (t)| is a highly intermittent quantity both in i (space) and time: It fluctuates around a small value when there are no plastic rearrangements, while it becomes large and strongly localized in space when a plastic rearrangement occurs somewhere in the system. For this reason, we consider as a quantitative measure of plastic activity in the system. The behavior in time of ∆ s (t) is shown in Fig. 4 for the same values of the peak stress discussed in Fig. 3. Quite remarkably (but not surprisingly), the behavior of ∆ s (t) is qualitatively similar to the one shown by u s (t). However, an important difference must be stressed: ∆ s is not affected by the presence of elastic waves. This difference can be understood in a simple way: the displacement due to elastic waves is relatively small and it is coherent in space (all droplets oscillate); in contrast the displacement due to plastic rearrangements is rather large and not coherent in space. Therefore, our quantity ∆ s is not sensitive to elastic waves.  and σ p /σ Y 1.6 for which the system is plastically flowing [7,8,27]. Inspection of Fig. 5 suggests that we should consider the probability distribution P (log(∆ s )), in agreement with the approach to intermittent fluctuations in dynamical systems theory [60]. We remark that, upon writing Z = log(∆ s ), it is easily shown that P (Z) = ∆ s P (∆ s ), i.e. the peak in the probability distribution of P (log(∆ s )) corresponds to the relevant value of ∆ s contributing to the average ∆ s [61].
The probability distributions P (log(∆ s )) are shown in Fig. 6 for the three different peak stresses already considered before: at small σ p , P (log(∆ s )) is peaked at small values and shows a rather long tail; at large σ p , P (log(∆ s )) is peaked at large values corresponding to the plastic flow previously discussed. Remarkably, at the transition point σ p /σ Y 1.4, the probability distribution P (log(∆ s )) is bimodal, i.e. the system shows transitions in time between two states with small (solid) and large (fluidized) values. Hence, we observe bimodality in time of two states that are unimodal in space. Now, we go back to the results shown in Fig. 3. The results discussed in terms of ∆ s suggest that transitions from the solid branch to the fluidized branch should be observed for u s as well. As already remarked, however, u s (t) is strongly perturbed by elastic waves which makes it impossible to observe the same bimodality unless the damped oscillations near u s = 0. Knowing the period and the dissipation time of the elastic wave [41], it is possible to fit the damped oscillations rather well as shown from the blue dashed line in the lower panel. We then obtain the filtered signalû s (t) by computing a running averageû where 2T el is the oscillation period of the elastic waves. In Fig. 8 we report the probability distribution for log |û s | for different peak stresses. We consider the log |û s | for the same reasons previously discussed for ∆ s . Comparing Figs. 6 and 8, once the elastic waves are filtered from the original signals, the probability distributions of log |û s | display the same features as P (log(∆ s )) at the same forcing.

V. DISCUSSION
Bimodal distributions and/or metastability have already been reported in the literature of amorphous systems [52,[61][62][63][64][65][66]. Regarding bimodal distributions, a recent theoretical work on amorphous solids by Jaiswal et al. [65] showed bimodality for an order parameter ad hoc constructed to see how much the system is correlated to the initial condition after an athermal, quasi-static (AQS) shearing protocol is applied. An experimental study on colloidal glasses by Chikkadi et al. [64] reported bimodality for the spatial distribution of an order parameter constructed with the time-integrated mean-square displacement of particles. It is also worth to recall some other studies on glasses under shear [61][62][63], in which a nontrivial statistics has been observed in the nonaffine displacements of particles, whose probability distribution exhibits peaks in different displacement ranges dependently on the observation time. The present investigation differs from previous studies in an important way: The results displayed in Fig. 6 and Fig. 8 show the succession in time of two metastable states at σ p /σ Y 1.4 corresponding to different rheological branches. In other words, the whole system spends roughly the same amount of time in both the elastic and fluidized phases, constantly tunneling back and forth from one state to the other. Transitions are due to plastic events which eventually drive the system from the solid to the fluidized branch. Once the system reaches the fluidized branch, it flows plastically with a large number of plastic rearrangements (see Fig. 5). Plastic flow dissipates energy quite efficiently, and eventually, the power input due to the forcing is not able to sustain the energy dissipation due to plastic flow and the system goes back to the solid branch. Last but not least, we argue that the choice of heterogeneous stress enhances the probability to perform transitions between the two branches because this choice reduces the region (in physical space) where the system may switch form a flowing regime to a solid/elastic state (and vice versa). This phenomenology differs from the bimodality discussed by Chikkadi et al. [64], since that is related to bimodality "in space" of the underlying shear. Our transitions in time, between elastic and fluidized states, also differ qualitatively from the observations of Refs. [65,66] on the intermittent periods of elastic loadings displayed in the failure of amorphous solids. Indeed, the loading and the failure take place on remarkably different timescales, which leads to a power-law distribution of the displacement field rather than a bimodal distribution (see Ref. [67] for a study of our model under shear flow). It must be also emphasized that in Ref. [65] bimodality is reported for the overlap variable describing how well the system remembers its initial configuration as a function of the applied quasistatic deformation: Such a choice would not allow to probe whether or not a given system repeatedly tunnels from a jammed to a flowing state and back, since the overlap is measured with respect to the starting configuration; thus, attaining a high overlap after a low value is reached is highly improbable. On the other hand, we observe bimodality for the time evolution of a rheological observable, signaling repeated transitions in time from a jammed to a flowing state and back, both states being unimodal in space. The presence of bimodality in time, for both log[∆ S (t)] and log[û s (t)], should be related to long-range space correlations of plastic events, of the order of the domain size. In fact, for systems with a short-range space correlation, the effect of a single plastic rearrangement is unable to develop a cascade (in space and time) of other plastic events and trigger the transition of the whole system from the metastable solid branch to the metastable fluidized branch. A similar reasoning applies for the reversed transition: Once plastic rearrangements stop occurring in some part of the system, the flow ceases locally and the transition to the solid branch for the whole system necessitates a correlation length that allows to cover the entire system size. This picture is actually borne out by a direct calculation of the correlation. A simple and intuitive way to look at space correlations is to compute the overlap-overlap correlation G(r) that was already used in Ref. [41]: We follow the analysis presented in Ref. [68], based on the idea of Ref. [69]. The physical meaning of G(r) is rather clear. In a nutshell we can say that small values of G(r) indicate that a part of the system moves somewhere while some other parts do not; large values of G(r) mean the opposite, implying that different parts of the system move or not move at the same time. In other words, for large values of dr G(r) (also known as dynamic heterogeneity [70]) the system either moves everywhere or does not move almost everywhere. We compute G(r) as follows: we consider two times t and t + T q and at each time we define the field . . x stands for space average and ρ A , ρ B are the densities of the continuous and dispersed phases. Then, we define the overlap q(x, y, t, t + T q ) as: Using Eq. (7) we define the overlap-overlap correlation function, centered in the middle of the channel at y = L/2: where . . . t,x stands for time and x averages and T q is chosen to be of the order of the time needed to perform a plastic rearrangement [41]. In Fig. 9 we show G c (r) (the connected part of G(r)) for σ p /σ Y 1.2, 1.4 and 1.6. Clearly, at the transition point σ p /σ Y 1.4, G c (r) is very large everywhere in the system. It is crucial to remark that the correlation length observed for σ p /σ Y 1.4 differs from the cooperative scale ξ of the system, the latter being equal to few droplet diameters [41]. When the system is in the fluidized branch for σ p /σ Y 1.6, the function G(r) decays to zero with a correlation scale of the order of ξ [7,8]. These features in our model have already been observed in conditions of imposed shear in Ref. [41].
Here they are confirmed in a setup with imposed heterogeneous stress. Moreover, stress-controlled experiments (like the one we propose) somehow offer valid alternatives to the shear controlled ones [28,37] in order to investigate the presence of multiple rheological branches.
Indeed, if the system shows long-range correlations [27] among plastic events, it may well be that in a shearcontrolled experiment, the shear bands (if they form) are strongly fluctuating both in time and space. Eventually, these fluctuations would simply disappear in the average flow profile and one should rather observe a complex dynamics in time of the shear stress characterized by a strong intermittency of the time derivative of the stress, a phenomenology well reminiscent of the stick-slip behavior [28,29,71]. In favor of this argument, we can mention the study by Varnik et al. on a model glass [28], where the authors find that long-lived shear bands are replaced by the emergence of stick-slip phenomena with intermittent bursts; this happens at very low shear rates, i.e., at the point of discontinuity between the solid branch at S = 0 and the fluid branch. We also mention the study by Pignon et al. [71] and by Picard et al. [29] where the two regimes of stick-slip and shear bands are observed for different apparent shear rates; however, one has to notice that the stick-slip observed here is rather an oscillatory flow with undetectable intermittency. In the specific case of the theoretical model by Picard et al. [29] this may be possibly related to the minimalistic nature of the model; i.e., no noise is added [18]. In order to stress the combined role of multiple rheological branches and space correlations, it is worthwhile to further connect our observations with some other results presented in the literature [53,73]. A recent work by Chaudhuri et al. [53] studied the interplay between the system size and the cooperative length in the flow arrest. Specifically, the model is that of soft-jammed repulsive disks (the "Durian" model [74]) in a periodic flow setup with heterogeneous stress, very similar to our stress profile. Upon decreasing the driving force, the authors determine the yielding threshold at which the flow ceases: Interestingly, under the conditions of periodic flow [53], when the cooperative length becomes of the order of the system size, the authors find that the yielding threshold is increased with respect to the yield stress σ Y , somehow in line with our findings (see Fig. 2). However, although an increased intermittency is reported at the onset of flow, the authors in Ref. [53] do not report any signature of metastable states like the one we observe, whereas simulation results are well predicted by the stationary fluidity model [27]. This contrasts with our observations. The interplay between system size and cooperative scale was also highlighted in another work by Chaudhuri & Horbach [73], studying the transition to the flowing regime in a pressure-driven flow for a Yukawa binary fluid [75,76]. When the cooperative length is of the order of the system size, it is shown that (in the long time limit) the system fluidizes nearly homogeneously. This behavior bears similarities with the transition from the solid-to-fluidized branch that we observe (see Fig. 3), with an important difference: The study by Chaudhuri & Horbach [73] does not report the existence of metastable states; i.e., once the fluidized state is reached it is shown to persist for the whole simulation time. However, the time spent by the system in the solid phase is remarkably long, much longer than the time that would be observed for an unstable state. In other words, one may argue that in Ref. [73] two metastable branches coexist although the possibility of transition between the branches has not been inves-   Fig. 6). The integral of G(r) is usually known as dynamic heterogeneity [72]. tigated in detail. According to the results shown in the previous section and to the overlap-overlap correlation function shown in Fig. 9, we identify two conditions that should be satisfied for a clear signature of metastable states: there must exist a difference between the static and the dynamic yield-stress values (i.e., there must exist two rheological branches) and there must be long-range correlations among plastic events. In Refs. [53,73] it is unknown whether one or both requirements are not met. We may argue that the model used by the authors in Ref. [53] is rather a model for a nonadhesive emulsion [77], and the difference between the static and dynamic yield stress is so small [37] that metastability between two different rheological branches cannot be observed. The Durian [74] model has also been used recently by Kawasaki & Berthier [52] to study the yielding transition under oscillatory flow. By analyzing the displacement fields of the particles the authors report a rather discontinuous transition at the yield stress: While above yield the fluctuations in the displacement fields are persistent in time (fluidized state), below the yield stress they are metastable and cease after some time. The possibility of transitions back to the fluidized state has not been studied in detail when changing the stress protocol and/or for longer simulation times, but again we argue that it would not be observed because of the model used [53]. All these considerations suggest that studies regarding the presence of shear bands and "stick-slip" should be consistently accompanied with measurement of the correlation functions. Correlations of the microscopic strain field were actually measured by Chikkaddi et al. [23] in colloidal glasses showing the formation of shear bands; however, such results were only obtained for the two bands separately.
Further analysis in our numerical simulations is also stimulated by a direct comparison of the phenomenology that we observe to that of glassy models [28,34,35] and, in particular, finite size p-spin models [34]. The nontrivial and interesting point is the observation that the system spontaneously develops two stable branches in its phase-space dynamics, similarly to the two rheological branches needed to describe the formation of shear bands. Such systems are also known to display a dynamic transition at some temperature T d . For T < T d the system is trapped in a large number of states, which grows as the exponential of its size. Upon applying an external force, the system shows a dynamic transition similar to a yield-stress transition. For a finite number of spins, the system exhibits bursts of activity, i.e., the activated process, which show self-similarity in size and time [34,78]. The probability distribution of the trapping time τ , namely the time between two successive bursts, shows a scaling behavior P (τ ) ∼ τ −a with a = 1 + T /T d . This behavior is qualitatively similar to the one described by SGR theories [19][20][21] based on the trap model [79]. Going back to our results, for the case where P (log(∆ s )) is bimodal (σ p /σ Y 1.4), we can define the trapping time τ spent by the system in the solid branch: We use the value of log(∆ s ) at the local minimum (see Fig. 6) as a threshold to condition the data. We expect τ to be a random variable and we look at the probability distribution P (τ ) shown in Fig. 10. The probability distribution P (τ ) behaves as a scaling function of τ , i.e. P (τ ) ∼ τ −α with α ∼ 1 showing the existence of nontrivial time correlations. In the bottom panel of Fig. 10 we show the running average S R (t) of the shear S(t) = 2πu s (t)/L [see Eq. 3 and below], normalized to its maximum S M R for the bimodal forcing σ p /σ Y 1.4. The running average S R (t) is computed when the system is in the flowing phase, i.e., when u s (t) belongs to the larger peak shown in the middle panel of Fig. 8; thus the value of the minimum of log |û s | is used as a cutoff. From the bottom panel of Fig. 10 it is possible to see that the time evolution of S R (t) is consistent with a logarithmic decay. Indeed, this is a further characterization of our results that can be verified in non-homogeneous stress experiments such as the one that we outline in the following section. The distribution can be well fit by a power law P (τ ) ∼ τ −α with α ∼ 1 indicating a self-similar structure in the intermittent transitions from the solid to the fluidized state. Such self-similar distribution has also been measured in spin glasses [34]. In the lower panel we show the running average SR(t) of the shear S(t) = 2πus(t)/L, normalized to its maximum S M R (see text for details). The behavior of SR(t) is consistent with a logarithmic decay.

VI. PRESSURE-DRIVEN FLOWS
The results discussed in the previous sections refer to fully periodic boundary conditions. In this section we want to comment about the possibility to obtain the same results in the case of realistic boundary conditions. In particular we consider the case of a pressure-driven flow in a two-dimensional channel and streamwise periodic boundary conditions. Since the system is driven with a constant force (pressure gradient) in the streamwise direction, the stress is a linear function of the coordinate y (see Fig. 1) and its absolute value reaches the maximum σ p at the boundaries. The system shows a rather clear apparent slip [49,57,80] at the smooth boundaries and this goes together with a nonzero mean flow; hence, the analysis in terms of ∆ s is no longer suitable. Furthermore, because of strong localization of plastic events at the boundary, there is less energy available to switch from one rheological branch to the other for the whole system. This implies that the characteristic "trapping" time becomes much longer with respect to the one observed in periodic boundary conditions. To perform long-time numerical simulations, we choose a square system with side L = 512 lattice points. In Fig. 11 we show the most interesting information obtained from our simulations. We choose σ p /σ Y 1.4 and we run simulations imposing a pressure gradient on a configuration picked from a lower forcing steady state (i.e., ramp-up protocol). The interesting variable to look at is the velocity flux u(t) defined as the space average at time t of the stream-wise velocity. In the upper panel of Fig. 11 we show u(t) (thick red line) for about 9 × 10 3 shear times. The system shows a nonzero average velocity (due to the slip at the boundaries) with superimposed bursts of larger values, similar to a stick-slip behavior. The probability distribution of u is shown in the middle panel, while the average velocity profile is shown in the bottom panel. Next, we increase σ p so that the system reaches a fluidized state (not shown). Once the statistical properties in the fluidized state could be considered stationary, we reduced the pressure gradient (i.e., ramp-down protocol) and perform a new numerical simulation at the same value of the peak stress σ p 1.4 already discussed. For this new simulation the results are reported with the thin blue line in Fig. 11. It is quite clear that the system shows transitions in the rheological behavior, characterized by small and large values of u (see the probability distribution). The qualitative picture is similar to the one discussed in the previous section although the time scale is much longer.
The results shown in Fig. 11 can be considered a preliminary investigation for systems with realistic boundary conditions. The point we want to highlight here is that the existence of two metastable states, discussed in the previous section, can be observed numerically and (most importantly) experimentally with long-time statistics (order 10 4 shear times of the system) and with a fine scanning of the forcing parameters. Moreover, further 1.4. Different line thicknesses (indicated also with different colors) correspond to different system preparations: Data reported using the thick red line refer to a system previously driven from a lower forcing (i.e., ramp-up protocol), whereas data represented by a blue thin line refer to a system previously driven at a larger forcing (i.e., ramp-down protocol). Top panel: Velocity flux u(t) as a function of time normalized by the shear time τshear. The thick red line displays some intermittent spikes while the blue one shows a sequence of transitions. Middle panel: Probability distributions for the velocity flux u displaying a single peak for the thick red line and a bimodal character for the thin blue line. Bottom panel: Velocity profile for ux(y) averaged over time and along the streamwise direction x, the thick red curve indicates a plug-flow dynamics dominated by an elastic bulk, while the thin blue one shows a developed velocity gradient near the boundaries.
analysis is required to investigate hysteresis effects.

VII. CONCLUDING REMARKS
Based on numerical simulations of a soft-glassy model we have studied its rheological response with an imposed space-dependent stress in an ideal fully periodic setup. The rheological properties of the model under study show the existence of multiple rheological branches with a difference between static and dynamic yield stress. The peak value σ p of the imposed stress is set close to the static yield stress of the material. We observe that the time dynamics of the system is remarkably non-steady, as it tunnels intermittently between two different states, a "solid" state and a "fluidized" one. Numerical simulations [38][39][40][41][42][43] allow to bridge the rheological response at large scales to the behavior displayed deeper down at small scales, where we observe a bimodal probability distribution of the largest value of the displacement field describing recurrent transitions in time between two unimodal states in space. Our results highlight the role of plastic rearrangements as the mechanical trigger for the hopping between the two states as well as the role of long-range correlations for the hopping to occur. Preliminary investigations have shown that such scenario holds for the more realistic case of a flow driven by a constant pressure gradient. Hence, such nonsteady yielding dynamics with recurrent transitions can be put to test by laboratory experiments.
From a general perspective, we point out that the existence of multiple rheological branches has been often introduced to explain the formation of permanent shear bands in soft-glasses [28,34]. From this point of view, the formation of shear bands can be considered as a "phase separation" in space, allowing the space coexistence of solid and fluidized regions [27]. Our observations somehow take a broader perspective and generalize the idea of coexistence in the time domain. For this coexistence in time, both long-range space correlations and the stress protocol are crucial. We indeed argue that a spatial correlation length of the order of the system size is crucial to trigger transitions between states and establish the time coexistence. Moreover, we argue that the choice of heterogeneous stress enhances the probability to perform transitions between the two branches because this choice reduces the region (in physical space) where the system may switch form a flowing regime to a solid/elastic state (and vice versa).
Given the role of space correlations in our system, it is then natural to comment on their expected role in a "classical" shear-banding scenario, i.e., when heteroge-neous flow is observed in presence of a shear-controlled experiments with homogeneous stress [5]. We indeed argue that the "phase coexistence" in space can be observed only if short-ranged correlations are present, whereas in presence of long-ranged correlations one would rather expect a stick-slip behavior. Noteworthy, preliminary simulations of our system under the conditions of imposed shear flow do not show permanent shear bands [67]. The above discussion may suggest that some complex dynamic and rheological properties observed in some softglasses, namely stick-slip behavior [10,29,71,81] and formation of permanent shear bands [13,15], can somehow be unified within the same theoretical framework, dependently on the range of space correlations. Given this view, it could be interesting to revisit our recent proposal [18] where cooperativity effects have been linked to the formation of permanent shear bands. One could add to the model a tunable correlation between plastic rearrangements and explore the consequences on the formation of the bands.
We remark that our findings share many features with the analysis performed on p-spin glasses near the dynamic transition at the temperature T = T d . The analysis performed in Ref. [34] shows that for T < T d the system develops two stable rheological branches. Moreover, the trapping time in the solid branch shows a power-law distribution which is also observed in our system. Finally, the theoretical analysis in Ref. [72] shows that, near the critical temperature, the system displays bimodality in the order parameter and long-range correlations in space (i.e., diverging dynamic heterogeneity), because of the spinodal character of the transition. All the above features are observed in our simulations.